Subroutine optp (kk,n1,lpia,lmaxin,lborn,u,f,fx) c *** see RH Landau, Program lpott, 1981 Implicit Real*8 (a-h,o-z) Real*8 imtij,mpi,mn,k,kp,kk,imvabs Dimension voptp1(17,17,4,2),voptp2(17,17,4,2) Common /twob/ t2body(17,17,2,14) Dimension rho(4,50), bi(100), vcoull(30), revabs(30,2), imvabs(30, 12) Dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij 1(14), imtij(14), kk(50), u(4356,2), f(1089,4), fx(1) Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,mn,nz,nes 1,nwaves,nifty,na Common /sec4/ achp,acmp,wsp,achn,acmn,wsn Common /sec5/ drho(4,50),rho,xkp,xk,imax/sec6/rcoul,rcut Common /nocheck/ nocheck data nocheckfirst/1/ data vcoull,revabs,imvabs/150*0.0e+00/ c c kk= grid point vector, n1= size of matrix,kk and u size set by ca c calculate partial wave decomp of rho- the form factor c rho(nff,i) is the i-1 th patial wave ampltd of the form factor c for a value of k a nd kp. nff determines which f.f. is used c nff=1 proton matter 2 neutron matter c 3 protone spin 4 neutron spin c Do all lvalues on 1stcall,then store on file7, ft is for temp c If (lpia.ne.1) GoTo 250 c 1st call to optp, calculate all u s Write(6,410) kk(n1) k = -200. Call tncm (k) n2 = 2*n1 nn = na-nz c make sure lmax=1 case Does not run into indef on intial runs ldmax = lmaxin If (ldmax.eq.1) ldmax = 2 rhl3 = (achp/hbarc)**2/2. rhl4 = (achn/hbarc)**2/2. c addeded terms for c13 spin...may want to change size here rhl5 = (wsn/hbarc)**2/2. as = 1./6. ap = (nz-2.)/6./nz*(acmp/achp)**2 If (nz.le.2) ap = 0. an = (nn-2.)/6./nn*(acmn/achn)**2 If (nn.le.2) an = 0. c special c13 , add in p1/2 n and rescale If (na.eq.13) an = (nn-3.)/6./(nn-1.)*(acmn/achn)**2 c Do loop over all grid points c size of u must be 4*(ngp+1)**2,2 Rewind 7 If (nifty(6).eq.6) nifty(17) = 1 aovera = (na-1.)/na If (nifty(17).eq.1) aovera = 1. If (na.eq.2) aovera = 1. If ((nifty(6).eq.4).or.(nifty(6).eq.6)) aovera = 0. xhl = aovera*2.*pi*pi aparam = 1. i = nifty(1)+1 GoTo (30,10,20,20,20,30,20,30,20,20), i 10 aparam = -1. GoTo 30 20 aparam = 0. 30 aparam = (1.5*nz*aparam/rcoul)*hbarc/137.036 xhl1 = rcoul/hbarc xhl2 = rcut/hbarc nif = nifty(10) Do 240 i1=1,n1 Do 240 i2=1,i1 k = kk(i2) kp = kk(i1) jj1 = i1 jj2 = i2 lmax = ldmax+5 If (nifty(16).ne.0) GoTo 90 c fo shell model nuclesu(c,0) use analytic expression with il*s c different neutron and protone densiteies possible zp = k*kp*rhl3 xp = (k*k+kp*kp)*rhl3/2. Call beslia (zp,bi) b = 0. If ((zp-xp).gt.-150.) b = exp(zp-xp) Do 40 i=1,lmax rho(3,i) = 0. rho(1,i) = b*(2.*i-1)*((1.-4.*ap*(xp-i+1))*bi(i)+4.*zp*ap 1 *bi(i+1)) 40 Continue If (na.eq.13) GoTo 60 If ((achp.ne.achn).or.(acmp.ne.acmn)) GoTo 60 c equal n and p densities Do 50 i=1,lmax rho(4,i) = 0. rho(2,i) = rho(1,i) 50 Continue nff = 1 GoTo 110 c unequal n and p densites in h.o. model 60 nff = 2 zn = k*kp*rhl4 xn = (k*k+kp*kp)*rhl4/2. Call beslia (zn,bi) b = 0. If ((zn-xn).gt.-150.) b = exp(zn-xn) Do 70 i=1,lmax rho(2,i) = b*(2.*i-1)*((1.-4.*an*(xn-i+1))*bi(i)+4.*zn*an 1 *bi(i+1)) rho(4,i) = 0. 70 Continue If (na.ne.13) GoTo 110 c spin for c13 nff = 4 zs = k*kp*rhl5 xs = (k*k+kp*kp)*rhl5/2. Call beslia (zs,bi) b = 0. If ((zs-xs).gt.-150.) b = exp(zs-xs) c normalize to 1/nn so can mult by nn and have 1 neutron b = b/nn Do 80 i=1,lmax rho(4,i) = b*(2.*i-1)*((1.-4.*as*(xs-i+1))*bi(i)+4.*zs*as 1 *bi(i+1)) rho(2,i) = ((nn-1.)/nn)*rho(2,i)+rho(4,i) 80 Continue If ((i1.eq.i2).and.(i1.eq.n1)) Write (6,350) k,(rho(4,i),i=1 1 ,lmax) GoTo 110 c for wood saxon typr(or arbitrary) nuclera n density,calc rho nume 90 xkp = kp xk = k imax = ldmax+5 lmax = imax c deteermine number of f.f.*s needed,if nifty(6)=3, need spin in nff = 2 If ((achp.eq.achn).and.(acmp.eq.acmn).and.(wsp.eq.wsn)) nff 1 = 1 If (nifty(6).eq.3) nff = 4 c for rho**2 in ca If (nifty(16).eq.2) nff = 4 Call rhokff (nff) If (nff.ne.1) GoTo 110 c nff=1, seet n and p matter densites=, all spin dens=0 Do 100 i=1,lmax rho(2,i) = rho(1,i) rho(3,i) = 0. rho(4,i) = 0. 100 Continue 110 Continue If ((nifty(10).ne.3).and.(nifty(10).ne.4)) GoTo 120 c include cutoff Coulomb potential in momentum space ak = k akp = kp Call vcoul (ak,akp,aparam,xhl1,xhl2,vcoull,ldmax,nif) 120 If (nifty(7) .lt. 3 .and. nifty(7).ne.5) GoTo 130 c include absorption on 2 nucleons Call vabs (k,kp,kk(n1),revabs,imvabs,ldmax) 130 Continue c new version....outter Do loopn over kp,k,,innerover all l s xmn = mn Call tnuc (jj1,jj2,kp,k,kk(n1),xmn) c---------------------------------------------------------------------- c Do loop over l opt c Do 230 ldum=1,ldmax If ((ldum.gt.lborn).and.((i1+i2).ne.(2*n1))) GoTo 230 j1 = ldum urc = 0. uic = 0. urs = 0. uis = 0. If (nifty(6).eq.3) rhl1 = j1*(j1-1.) c Do over pi-n angular momentum;s,p,d,f,g, Do 200 j2=1,5 j4 = iabs(j1-j2)+1 j5 = j1+j2-1 If (nifty(6).eq.3) rhl2 = j2*(j2-1.) Do 190 j3=j4,j5 rhl = cgc2(j3,j2,j1) c spin independt potential If (nifty(1).eq.1.or.nifty(1).eq.8) GoTo 140 c pi plus scatterinf jpr = 4*j2-3 jpi = jpr+1 jnr = jpr+2 jni = jpr+3 GoTo 150 140 Continue c pi minus scatterin, rearrange code jnr = 4*j2-3 jni = jnr+1 jpr = jnr+2 jpi = jnr+3 150 Continue If (nifty(1).ne.4) GoTo 160 c pi zero scattering, take average of p and n amplitudes bnuc(jpr) = (bnuc(jpr)+bnuc(jnr))/2. bnuc(jnr) = bnuc(jpr) bnuc(jpi) = (bnuc(jpi)+bnuc(jni))/2. bnuc(jni) = bnuc(jpi) If (nifty(6).ne.3) GoTo 160 bnucf(jpr) = (bnucf(jpr)+bnucf(jnr))/2. bnucf(jnr) = bnucf(jpr) bnucf(jpi) = (bnucf(jpi)+bnucf(jni))/2. bnucf(jni) = bnucf(jpi) 160 Continue urc = urc+rhl*(nz*rho(1,j3)*bnuc(jpr)+nn*rho(2,j3)* 1 bnuc(jnr)) uic = uic+rhl*(nz*rho(1,j3)*bnuc(jpi)+nn*rho(2,j3)* 1 bnuc(jni)) If ((nifty(3).ne.2).and.(nifty(3).ne.3)) GoTo 180 c print out off-shell tmatrix onto tape 3 If (ldum.ne.1) GoTo 180 If ((i1.ne.i2).or.(i1.ne.1)) GoTo 170 c first time thru, print out identifiers Rewind 3 Write (3,360) nz,na,(nifty(n),n=1,20) Write (3,370) achp,acmp,wsp,achn,acmn,wsn,rcoul, 1 rcut ngp = n1-1 Write (3,360) ngp Write (3,380) (kk(n),n=1,n1) c print out pi n t matrix in pi-a com for k' gtr k 170 If (i1.lt.i2) GoTo 180 Write (3,390) i1,i2,bnuc(jpr),bnuc(jpi),bnuc(jnr), 1 bnuc(jni) If (nifty(6).eq.3) Write (3,380) bnucf(jpr),bnucf 1 (jpi),bnucf(jnr),bnucf(jni) 180 If (nifty(6).ne.3) GoTo 190 c ---- ---spin dependnet potential If (j1.eq.1) GoTo 190 rhl = rhl*(rhl2+rhl1-j3*(j3-1.))/(2.*rhl1) urs = urs+rhl*(nz*rho(3,j3)*bnucf(jpr)+nn*rho(4,j3) 1 *bnucf(jnr)) uis = uis+rhl*(nz*rho(3,j3)*bnucf(jpi)+nn*rho(4,j3) 1 *bnucf(jni)) 190 Continue 200 Continue c form potentials for state of definite j c********* mr lu stores optp here and prints it later ********* c pure unclear part If (ldum .le. 4) then voptp1(i1,i2,ldum,1) = urc voptp1(i1,i2,ldum,2) = uic Endif rhl = xhl/(2.*j1-1.) vc = 0. If ((nifty(10).ne.3).and.(nifty(10).ne.4)) GoTo 210 vc = vcoull(j1)*aovera If (nifty(6).eq.6) vc = vcoull(j1) 210 If (nifty(7).ne.3.and.nifty(7).ne.4) GoTo 220 c add in vabs to central potential c aovera scaling alReady Done for vabs in sub vabs urc = urc+revabs(j1,1)/rhl uic = uic+imvabs(j1,1)/rhl if (ldum .le. 4) then voptp2(i1,i2,ldum,1) = revabs(j1,1)/rhl voptp2(i1,i2,ldum,2) = imvabs(j1,1)/rhl Endif If (nifty(7).ne.4.or.nifty(6).ne.3) GoTo 220 c include spin-dependent absorption urs = urs+revabs(j1,2)/rhl uis = uis+imvabs(j1,2)/rhl 220 Continue c set up potential needed for this l in integral eqtn npot1 = 2*ldum f(npot1-1,1) = rhl*(urc+(j1-1)*urs)+vc f(npot1,1) = -rhl*(uic+(j1-1)*uis) nnn = npot1-1 If(kp.eq.kk(n1).and.kp.eq.k.and.ldum.eq.1) 1 Write(6,410)kk(n1) If ((kp.eq.kk(n1)).and.(ldum.eq.2)) Write (6,400) k,f(nnn 1 ,1),vc,revabs(j1,1),revabs(j1,2),f(npot1,1),imvabs(j1,1), 2 imvabs(j1,2) If (nifty(6).ne.3) GoTo 230 f(npot1-1,2) = rhl*(urc-j1*urs)+vc f(npot1,2) = -rhl*(uic-j1*uis) c Do loop over ldum ends 230 Continue npot1m = 2*lborn If (ldmax.lt.lborn) npot1m = 2*ldmax If ((i1.eq.i2).and.(i1.eq.n1)) npot1m = 2*ldmax c new unformatted Read for time saving If (nifty(6).ne.3) Write (7) (f(npot1,1),npot1=1,npot1m) If (nifty(6).eq.3) Write (7) ((f(npot2-1,i),f(npot2,i),i=1,2 1 ),npot2=2,npot1m,2) c Do loop over kp and k ends 240 Continue 250 Rewind 7 c pots all stored,Read in the values for cuurrent lpia c npts =? of i1,i2 combos(followed by u s for all l s) npts = (n1*(1+n1))/2 If (lpia.ne.1) GoTo 260 c first call for Read, setup params c n.b....the values of vl(i,j) for all l but fixed c i,j are on one record via unformatted Write c need skip to new record to get new ij, nspin = 1 npot = 2 If (nifty(6).ne.3) GoTo 260 c spin-dependent case, now 4 v s adjacent for each l vs 2 If nf nspin = 2 npot = 2*npot 260 ldum = 2*(lpia-1)*nspin If (lpia.gt.lborn) GoTo 290 c Do loop over npts (i.e. records) Do 280 npt=1,npts If (lpia.ne.1) GoTo 270 Read (7,end=280) (f(npt,i),i=1,npot) GoTo 280 c space to l value needed 270 Read (7,end=280)(rhl,i=1,ldum),(f(npt,i),i=1,npot) 280 Continue GoTo 310 c lpia .gt. lborn, Read only last on-shell point 290 nskip = npts-1 Do 300 i=1,nskip Read (7) 300 Continue Read (7,end=310) (rhl,i=1,ldum),(f(npts,i),i=1,npot) c set up potential matrix 310 istart = 1 If (lpia.gt.lborn) istart = n1 Do 340 i1=istart,n1 Do 340 i2=istart,i1 Iflag = 0 npt = ((i1-1)*i1)/2+i2 npot1 = (i2-1)*n2+i1 320 npot2 = npot1+n1*(n2+1) npot3 = npot2-n1 npot4 = npot1+n1 u(npot1,1) = f(npt,1) u(npot2,1) = u(npot1,1) u(npot3,1) = f(npt,2) u(npot4,1) = -u(npot3,1) c spin-dependnt potentials If (nifty(6).ne.3) GoTo 330 u(npot1,2) = f(npt,3) u(npot2,2) = u(npot1,2) u(npot3,2) = f(npt,4) u(npot4,2) = -u(npot3,2) c set up potntial with u(kp,k)=u(k,kp)..symmetric u 330 If ((i1.eq.i2).or.(iflag.eq.1)) GoTo 340 npot1 = (i1-1)*n2+i2 Iflag = 1 GoTo 320 340 Continue Return c 350 format (30h c13 neutron spin ff for k=kp=,e15.5/10(10x,6e15.6/)) 360 format (2i5,2x,i1,i2,8i1,10i3,f10.3,e13.6) 370 format (8f10.4) 380 format (5e16.7) 390 format (2i5/5e16.7) 400 format (3x,3e19.4,2e11.3,e19.4,2e11.3) 410 format(28h0opt poten values for kp=k0=,f10.3,7hand l=1/20x,1hk, 1 13x,3hrev,16x,2hvc,15h rabs(n0 flip) ,11h rabs(flip),3himv, 217x,27h imvabs(noflip),imabs(flip)/) End