c subroutine FfInt(npts,P,Zk,Uborn,Urop,Uct,Ucs, & zUborn,zUrop,zUct1,zUct2,zUct0,zUcs1,zUcs2,zUcs0) c c*********************************************************************** c FfInt calc form factor for each grid k by num integration. c ff at onshell point is calc'ed in 2nd part. c c npts: the number of grid points c rmev: radius of the bag in unit of 1/Mev c Uborn: vertex form factor array in Born diagrams c Uct: 3D array of integral in time component c Ucs: 3D array of integral in space component c clmax: the maximal channel number c lmax: the maximal orbital angular momentum c ns2: the normalization constant for ground state c z*: corresponding onshell point matrix c CDM and CBM share same matrix var, their actual values diff c because of diff inte upper limit c*********************************************************************** c LOCAL variables integer first,ich,jch,nch,ipt,jpt,l,nmax,lmax integer nx,NB,i,nbag,nfor,iq parameter (nmax=101,nx=100,nfor=200,lmax=1,NB=1000) double precision rmev, vs, ns2, q, tmp, Utmp double precision m1, m2, c1r, c2r, hbarc, Urop0(nmax),Vrop0(nmax) double precision P0(nmax), Uborn0(nmax), Uct0(nmax), Ucs0(nmax) double precision Q0(nmax), Vborn0(nmax), Vct0(nmax), Vcs0(nmax) double precision Y2b(nmax), Y2t(nmax), Y2s(nmax), Y2r(nmax) double precision j0s(NB), j1s(NB), jk(NB,64,0:2) double precision Xs(NB),Ws(NB), Xx(nx),Wx(nx), Xq(nfor),Wq(nfor) double precision fb(NB), ft(NB), fs(NB), fr(NB), ak, ak2 double precision ic0, ic1, is0, is2, twoopi, mid, upper, FFgau double complex CFfgau, cjk(NB,2,0:2), Csbj, zq2 double complex ztmp, zc0, zc1, zs0, zs2 c SHARED variables integer Nifty, npts, Fidx double precision Uborn(64), Uct(64,64,0:1), Ucs(64,64,0:1) double precision Urop(64), P(64), zero, SBj double precision f0, mgb, alpha, c1, c2, Mbare, rcbm double complex Zk(2), zUborn(2), zUrop(2) double complex zUct1(2,64,0:1),zUcs1(2,64,0:1) double complex zUct2(64,2,0:1),zUcs2(64,2,0:1) double complex zUct0(2,2,0:1),zUcs0(2,2,0:1) c parameter (zero=0.0d0,vs=2.0427869d0,hbarc=197.3d0) parameter (twoopi=0.63661977d0) c common /control/ Nifty(15) c bag.h Version 1.4 contains the CCDM parameters. first 3 independent common /bag/ rcbm,f0,mgb,alpha,c1,c2,Mbare(3),Fidx(3) c data first /1/ c*********************************************************************** c nch = Nifty(2) if (first .eq. 1) then open (1, file='cdmff', status='unknown') if (Nifty(1).ne.0) then c include Roper read(1,*) alpha, mgb, c1, c2, c1r, c2r else c donot include Roper read(1,*) alpha, mgb, c1, c2 endif c if (nifty(5).eq.1) then c Compute hardron bare mass from Phatak's parameterization c Mbare(1) = (c1 - c2*alpha)*mgb c Mbare(2) = (c1 + c2*alpha)*mgb c else c use input file to scale m(nucleon) and m(delta) directly c for CBM nif14=3, mgb and alpha doesnot mean anything mgb = (Mbare(2) + Mbare(1))/c1/2.0 alpha = (Mbare(2) - Mbare(1))/c2/2.0/mgb c endif Mbare(3) = (c1r - c2r*alpha)*mgb c c if (Nifty(14).eq.1) then c***** Read in CDM form factor computed by Dr. Phatak's code. ***** c Interpolate for gauss grids from 101 equally distributed points. c nmax is actual total num of momentum points. do ipt = 1, nmax if (Nifty(1).ne.0) then read (1,*) Q0(ipt),Vct0(ipt),Vborn0(ipt),Vcs0(ipt),Vrop0(ipt) else read (1,*) Q0(ipt),Vct0(ipt),Vborn0(ipt),Vcs0(ipt) endif c back to the original ff ak = (ipt-1.)/25. c ak2 = ak**2/4. c ak2 = -ak**2/4. ak2 = 0. Vct0(ipt) = Vct0(ipt) *exp(ak2) Vborn0(ipt) = Vborn0(ipt) *exp(ak2) Vcs0(ipt) = Vcs0(ipt) *exp(ak2) Vrop0(ipt) = Vrop0(ipt) *exp(ak2) c P0(ipt) = mgb * P0(ipt) c Ucs0(ipt) = Ucs0(ipt)/mgb c divided by hbarc?? enddo close(1) endif endif c if ((Nifty(14).eq.1) .and. & (first.eq.1 .or. m1.ne.Mbare(1) .or. m2.ne.Mbare(2))) then first = 0 m1 = Mbare(1) m2 = Mbare(2) mgb = (Mbare(2) + Mbare(1))/c1/2.0 alpha = (Mbare(2) - Mbare(1))/c2/2.0/mgb Mbare(3) = (c1r - c2r*alpha)*mgb c c convert to actual unit, then multiplied by pi ff to each leg c here rpi = 0.4fm do ipt = 1, nmax P0(ipt) = mgb * Q0(ipt) Ucs0(ipt) = Vcs0(ipt)/mgb ak2 = exp(-P0(ipt)**2/(40.*hbarc**2)) Uct0(ipt) = Vct0(ipt) * ak2*ak2 Uborn0(ipt) = Vborn0(ipt) * ak2 Ucs0(ipt) = Ucs0(ipt) * ak2*ak2 Urop0(ipt) = Vrop0(ipt) * ak2 c write to the file 'debug' write(22,221) P0(ipt), Uborn0(ipt) enddo 221 format(1x,'dg-ff1',2e12.3) call Spline(P0, Uct0, nmax, 0.d0, 0.0d0, Y2t) call Spline(P0, Uborn0, nmax, 0.d0, 0.0d0, Y2b) call Spline(P0, Ucs0, nmax, 0.d0, 0.0d0, Y2s) call Spline(P0, Urop0, nmax, 0.d0, 0.0d0, Y2r) c c gauss grids for coordinate integration 0-4fm nbag = 400 mid = 1.5/hbarc upper = 3.0/hbarc call Gauss(nbag, 1, mid, upper, Xs, Ws) c call Gauss(nbag, 0, zero, upper, Xs, Ws) c gauss grid for fourier transf from p to r space c call Gauss(nfor, 2, zero, 252.d0, Xq, Wq) c note that maximum Xq(nfor) must be less than maximum P0(nmax) c phatak' input range:0 -> 0.451e+4 MeV. c diff mgb may result in a smaller P0(nnmax). call Gauss(nfor, 0, zero, 4500.d0, Xq, Wq) c calc fourier transf of ff, convert p to r space write(8,*) 'Fourier transform: i : x(i) : fbi : ft : fs:fr' do i = 1, nbag fb(i) = zero ft(i) = zero fs(i) = zero fr(i) = zero do iq = 1, nfor tmp = Xq(iq)**2*Wq(iq)*SBj(0,Xq(iq)*Xs(i)) call Splint(P0, Uborn0, Y2b, nmax, Xq(iq), Utmp) fb(i)=fb(i) + tmp*Utmp call Splint(P0, Uct0, Y2t, nmax, Xq(iq), Utmp) ft(i)=ft(i) + tmp*Utmp call Splint(P0, Ucs0, Y2s, nmax, Xq(iq), Utmp) fs(i)=fs(i) + tmp*Utmp call Splint(P0, Urop0, Y2r, nmax, Xq(iq), Utmp) fr(i)=fr(i) + tmp*Utmp enddo fb(i) = fb(i)*twoopi ft(i) = ft(i)*twoopi fs(i) = fs(i)*twoopi fr(i) = fr(i)*twoopi if (Nifty(4).ne.3 .and. mod(i,10).eq.0) & write(8,901) i, Xs(i),fb(i),ft(i),fs(i),fr(i) enddo 901 format(i4, 5e15.5) c do i = 1, nbag do l = 0, lmax+1 do ipt = 1, npts jk(i,ipt,l) = SBj(l, P(ipt)*Xs(i)) enddo enddo enddo c c Born(N and Delta) and Roper ff do ipt = 1, npts Uborn(ipt) = zero Urop(ipt) = zero do i = 1, nbag Uborn(ipt) = Uborn(ipt) + Xs(i)**2*Ws(i) & *fb(i)*jk(i,ipt,0) Urop(ipt) = Urop(ipt) + Xs(i)**2*Ws(i) & *fr(i)*jk(i,ipt,0) enddo enddo c Contact ff do jpt = 1, npts do ipt = 1, jpt ic0 = zero ic1 = zero is0 = zero is2 = zero do i = 1, nbag ic0=ic0 + Xs(i)**2*Ws(i)*ft(i)*jk(i,ipt,0)*jk(i,jpt,0) ic1=ic1 + Xs(i)**2*Ws(i)*ft(i)*jk(i,ipt,1)*jk(i,jpt,1) c is0=is0 + Xs(i)**2*Ws(i)*fs(i)*jk(i,ipt,0)*jk(i,jpt,0) is2=is2 + Xs(i)**2*Ws(i)*fs(i)*jk(i,ipt,2)*jk(i,jpt,2) enddo Uct(ipt,jpt,0) = ic0 Uct(ipt,jpt,1) = ic1 Ucs(ipt,jpt,0) = zero Ucs(ipt,jpt,1) = P(ipt)*P(jpt)*(is0-is2)/3.0d0 enddo enddo c do i = 1,nbag c write(22,223) i, xs(i),ws(i),fb(i),jk(i,1,0) c enddo c do i = 1, nmax c write(22,224) i, p0(i), uborn0(i) c enddo c do i = 1, nfor c write(22,225) i, xq(i), wq(i) c enddo 222 format(1x,'dg-ff2',i5, 2e12.3) 223 format(1x,'dg-ff3',i5, 4e12.3) 224 format(1x,'dg-ff4',i5, 2e12.3) 225 format(1x,'dg-ff4',i5, 2e12.3) c c else if (Nifty(14).eq.2 .and. first.eq.1) then c***** GAUSS form factor ***** first = 0 do ipt = 1, npts UBorn(ipt) = Ffgau(Nifty(12), 2, (P(ipt)/mgb)**2) enddo c c contact term ff by integrating over Xx=cos(angle) call Gauss(nx, 0, -1.0d0, 1.0d0, Xx, Wx) c do jpt = 1, npts do ipt = 1, jpt Uct(ipt, jpt, 0) = zero Uct(ipt, jpt, 1) = zero Ucs(ipt, jpt, 0) = zero Ucs(ipt, jpt, 1) = zero do i = 1, nx q = sqrt(P(ipt)**2+P(jpt)**2-2.*P(ipt)*P(jpt)*Xx(i)) Utmp = Ffgau(Nifty(12), 1, (q/mgb)**2) Uct(ipt,jpt,0) = Uct(ipt,jpt,0) + 0.5*Wx(i)*Utmp Uct(ipt,jpt,1) = Uct(ipt,jpt,1) + 0.5*Xx(i)*Wx(i)*Utmp c Utmp = Ffgau(Nifty(12), 3, (q/mgb)**2)/mgb Ucs(ipt,jpt,1) = Ucs(ipt,jpt,1) + P(ipt)*P(jpt)* & 0.25*(1.-Xx(i)**2)*Wx(i)*Utmp enddo enddo enddo c c else if (Nifty(14).eq.3 .and. first.eq.1) then c***** CBM Born ff (Eq 2.28 of Veit ea,PRD33,1859(1986)) ***** first = 0 c calc bessel func at grid for later use. nbag = 200 rmev = rcbm/hbarc call Gauss(nbag, 0, 0.0d0, rmev, Xs, Ws) c do i = 1, nbag j0s(i) = SBj(0, vs*Xs(i)/rmev) j1s(i) = SBj(1, vs*Xs(i)/rmev) do l = 0, lmax do ipt = 1, npts jk(i, ipt, l) = SBj(l, P(ipt)*Xs(i)) enddo enddo enddo c c Born vertex function do ipt = 1, npts tmp = P(ipt)*rmev Uborn(ipt) = SBj(1, tmp)/tmp * vs/(vs-1.d0) enddo c c Normalization constant for ground state ns2 = 1.d0/(2.d0*SBj(0,vs)**2*rmev**3) * vs/(vs-1.d0) c c contact interaction: time component Uct; space component Ucs do l = 0, lmax do jpt = 1, npts do ipt = 1, jpt Uct(ipt, jpt, l) = zero Ucs(ipt, jpt, l) = zero do i = 1, nbag Uct(ipt, jpt, l) = Uct(ipt, jpt, l) + & Xs(i)*Xs(i)*(j0s(i)*j0s(i)+j1s(i)*j1s(i)) & *jk(i,ipt,l)*jk(i,jpt,l)*Ws(i)*ns2 c Ucs(ipt,jpt,l) = Ucs(ipt,jpt,l)+Ws(i)*Xs(i) & *2.d0*j0s(i)*j1s(i)*jk(i,ipt,l)*jk(i,jpt,l)*ns2 enddo enddo enddo enddo endif c c first = 0 c endif c endif first, E-independent part, except Mbare changes c c************************************************ c**** On-Shell point ff matrix **** c************************************************ cwrite(*,1001) Zk(1),Zk(2) 1001 format('ffint: onshl momentum',4e12.3) Mbare(3) = (c1r - c2r*alpha)*mgb if (Nifty(14) .eq. 1) then c **** CDM **** c Set up onshell complex bessel func do i = 1, nbag do l = 0, lmax+1 do ich = 1, nch cjk(i,ich,l) = Csbj(l,Zk(ich)*Xs(i)) enddo enddo enddo c c BORN(N and Delta) and Roper ff do ich = 1, nch zUborn(ich) = zero zUrop(ich) = zero do i = 1, nbag zUborn(ich) = zUborn(ich) + Xs(i)**2*Ws(i) & *fb(i)*cjk(i,ich,0) zUrop(ich) = zUrop(ich) + Xs(i)**2*Ws(i) & *fr(i)*cjk(i,ich,0) enddo enddo c c CONTACT ff: k'=0, zUct1; k=0, zUct2 do ich = 1, nch do ipt = 1, npts zc0 = zero zc1 = zero zs0 = zero zs2 = zero do i = 1, nbag zc0=zc0 + Xs(i)**2*Ws(i)*ft(i)*cjk(i,ich,0)*jk(i,ipt,0) zc1=zc1 + Xs(i)**2*Ws(i)*ft(i)*cjk(i,ich,1)*jk(i,ipt,1) c zs0=zs0 + Xs(i)**2*Ws(i)*fs(i)*cjk(i,ich,0)*jk(i,ipt,0) zs2=zs2 + Xs(i)**2*Ws(i)*fs(i)*cjk(i,ich,2)*jk(i,ipt,2) enddo zUct1(ich,ipt,0) = zc0 zUct1(ich,ipt,1) = zc1 zUct2(ipt,ich,0) = zUct1(ich,ipt,0) zUct2(ipt,ich,1) = zUct1(ich,ipt,1) zUcs1(ich,ipt,0) = zero zUcs2(ipt,ich,0) = zero zUcs1(ich,ipt,1) = Zk(ich)*P(ipt)*(zs0-zs2)/3.0d0 zUcs2(ipt,ich,1) = zUcs1(ich,ipt,1) enddo enddo c c onshell k'=k=k0, zUct0 do ich = 1, nch do jch = 1, nch zc0 = zero zc1 = zero zs0 = zero zs2 = zero do i = 1, nbag zc0 = zc0 + Xs(i)**2*Ws(i)*ft(i)*cjk(i,ich,0)*cjk(i,jch,0) zc1 = zc1 + Xs(i)**2*Ws(i)*ft(i)*cjk(i,ich,1)*cjk(i,jch,1) zs0 = zs0 + Xs(i)**2*Ws(i)*fs(i)*cjk(i,ich,0)*cjk(i,jch,0) zs2 = zs2 + Xs(i)**2*Ws(i)*fs(i)*cjk(i,ich,2)*cjk(i,jch,2) enddo zUct0(ich,jch,0) = zc0 zUct0(ich,jch,1) = zc1 zUcs0(ich,jch,0) = zero zUcs0(ich,jch,1) = Zk(ich)*Zk(jch)*(zs0-zs2)/3.0d0 enddo enddo c c***** GAUSS form factor ***** else if (Nifty(14) .eq. 2) then do ich = 1, nch zUborn(ich) = CFfgau(Nifty(12), 2, (Zk(ich)/mgb)**2) enddo c c k'=0, zUct1; k=0, zUct2 do ich = 1, nch do ipt = 1, npts zUct1(ich, ipt, 0) = zero zUct1(ich, ipt, 1) = zero zUcs1(ich, ipt, 0) = zero zUcs1(ich, ipt, 1) = zero do i = 1, nx zq2 = Zk(ich)**2+P(ipt)**2-2.*Zk(ich)*P(ipt)*Xx(i) ztmp = CFfgau(Nifty(12), 1, zq2/mgb**2) zUct1(ich,ipt,0) = zUct1(ich,ipt,0) + 0.5*Wx(i)*ztmp zUct1(ich,ipt,1) = zUct1(ich,ipt,1) + 0.5*Xx(i)*Wx(i)*ztmp c ztmp = CFfgau(Nifty(12), 3, zq2/mgb**2)/mgb zUcs1(ich,ipt,1) = zUcs1(ich,ipt,1) + Zk(ich)*P(ipt)* & 0.25*(1.-Xx(i)**2)*Wx(i)*ztmp enddo zUct2(ipt,ich,0) = zUct1(ich,ipt,0) zUct2(ipt,ich,1) = zUct1(ich,ipt,1) zUcs2(ipt,ich,0) = zUcs1(ich,ipt,0) zUcs2(ipt,ich,1) = zUcs1(ich,ipt,1) enddo enddo c c k'=k=k0, zUct0 do ich = 1, nch do jch = 1, nch zUct0(ich, jch, 0) = zero zUct0(ich, jch, 1) = zero zUcs0(ich, jch, 0) = zero zUcs0(ich, jch, 1) = zero do i = 1, nx zq2 = Zk(ich)**2+Zk(jch)**2-2.*Zk(ich)*Zk(jch)*Xx(i) ztmp = Cffgau(Nifty(12), 1, zq2/mgb**2) zUct0(ich,jch,0) = zUct0(ich,jch,0) + 0.5*Wx(i)*ztmp zUct0(ich,jch,1) = zUct0(ich,jch,1) + 0.5*Xx(i)*Wx(i)*ztmp c ztmp = Cffgau(Nifty(12), 3, zq2/mgb**2)/mgb zUcs0(ich,jch,1) = zUcs0(ich,jch,1) + Zk(ich)*Zk(jch)* & 0.25*(1.-Xx(i)**2)*Wx(i)*ztmp enddo enddo enddo c else if (Nifty(14) .eq. 3) then c **** CBM **** do ich = 1, nch ztmp = Zk(ich)*rmev zUborn(ich) = Csbj(1, ztmp)/ztmp * vs/(vs-1.d0) enddo c c Set up onshell complex bessel func do l = 0, lmax do ich = 1, nch do i = 1, nbag ztmp = Zk(ich)*Xs(i) cjk(i,ich,l) = Csbj(l,ztmp) enddo enddo enddo c c k'=0, zUct1; k=0, zUct2 do l = 0, lmax do ich = 1, nch do ipt = 1, npts zUct1(ich, ipt, l) = zero zUct2(ipt, ich, l) = zero zUcs1(ich, ipt, l) = zero zUcs2(ipt, ich, l) = zero do i = 1, nbag zUct1(ich,ipt,l) = zUct1(ich,ipt,l) & +Ws(i)*Xs(i)*Xs(i)*(j0s(i)**2+j1s(i)**2) & *cjk(i,ich,l)*jk(i,ipt,l) * ns2 zUcs1(ich,ipt,l) = zUcs1(ich,ipt,l) & +Ws(i)*Xs(i)*(2.0d0*j0s(i)*j1s(i)) & *cjk(i,ich,l)*jk(i,ipt,l) * ns2 enddo zUct2(ipt,ich,l) = zUct1(ich,ipt,l) zUcs2(ipt,ich,l) = zUcs1(ich,ipt,l) enddo enddo enddo c c k'=k=k0, zUct0 do l = 0, lmax do ich = 1, nch do jch = 1, nch zUct0(ich, jch, l) = zero zUcs0(ich, jch, l) = zero do i = 1, nbag zUct0(ich,jch,l) = zUct0(ich,jch,l) & +Ws(i)*Xs(i)*Xs(i)*(j0s(i)**2+j1s(i)**2) & *cjk(i,ich,l)*cjk(i,jch,l) * ns2 zUcs0(ich,jch,l) = zUcs0(ich,jch,l) & +Ws(i)*Xs(i)*(2.0d0*j0s(i)*j1s(i)) & *cjk(i,ich,l)*cjk(i,jch,l) * ns2 enddo enddo enddo enddo c endif return end c c Append gaussian form factor, real and complex case2 double precision function Ffgau(nif12, num, p2) integer nif12, num, first double precision ut0, ub0, us0, wt, wb, ws, t0, b0, s0, p2 parameter (ut0=.9999969,ub0=.8580252,us0=1.471414) data first/1/ c if (first .eq. 1) then if (nif12 .eq. 1) then wt = 0.0 wb = 0.0 ws = 0.0 t0=5.9528714 b0=6.3627147 s0=4.5191895 else wt= 1.0/.7**2 wb=1.0/.66**2 ws=1.0/.86**2 t0= 3.4103130 b0= 3.4369686 s0= 2.9414983 endif first = 0 endif c if (num .eq. 1) then Ffgau = ut0*(1.0-wt*p2)*exp(-t0*p2) else if (num .eq. 2) then Ffgau = ub0*(1.0-wb*p2)*exp(-b0*p2) else Ffgau = us0*(1.0-ws*p2)*exp(-s0*p2) endif return end c double complex function Cffgau(nif12, num, p2) integer nif12, num, first double precision ut0, ub0, us0, wt, wb, ws, t0, b0, s0 double complex p2 parameter (ut0=.9999969,ub0=.8580252,us0=1.471414) data first/1/ c if (first .eq. 1) then if (nif12 .eq. 1) then wt = 0.0 wb = 0.0 ws = 0.0 t0=5.9528714 b0=6.3627147 s0=4.5191895 else wt= 1.0/.7**2 wb=1.0/.66**2 ws=1.0/.86**2 t0= 3.4103130 b0= 3.4369686 s0= 2.9414983 endif first = 0 endif c if (num .eq. 1) then Cffgau = ut0*(1.0-wt*p2)*exp(-t0*p2) else if (num .eq. 2) then Cffgau = ub0*(1.0-wb*p2)*exp(-b0*p2) else Cffgau = us0*(1.0-ws*p2)*exp(-s0*p2) endif return end