program bsmain c 1=output,tape4=64,tape8=64,tape7,debug=output) c bound states via det(1-gv)=0 c ridge32 version july 85 c nes =no 2 body energies for t s, or or swtch for k-p potents c nwaves lt use non reltv reduced mass in sig pi channel c zero, then reset c nifty(19)= 0 real*8 k2 in v2......dont redefine k2 in v2 as e change c 1 .. do ... c 2 comple .... do ... implicit real*8 (a-h,o-z) external detcal real*8 xeri(2),w(31),drr(10),dri(10),detri(2) common/det/ar11,alfanz c the open stmts are c'ed out so files spefied on run card open(5,file='runbs',status='unknown') open(6,file='outbs',status='unknown') open(7,status='scratch',access='sequential', 1 form='unformatted') i=signal(119,dumsub,1) write(6,879) write(6,879) write(6,879) 879 format(' BOUND STATE calculation via det(1-GV)=0, ridge32, ', 1'8/86, program bsdetr, detcal') 10 read(5,*,end=233)drr(1),eim,(drr(i),i=2,7) read(5,*)tol,aitmax 880 format(8f10.6) write(6,882)drr(1),eim,(drr(i),i=2,7) write(6,882)tol,aitmax 882 format(8f15.6) maxit = aitmax nrei=2 if(eim.ge.0.d+00)nrei=1 l=0 c calculate point coulomb b.e. for sch eqtn or kge c no dirac yet, i.e. no spin c fj=l +- 1/2....i.e.fj=l for these cases fj=l ekge=0. 226 format(24h nbohr,e(se),e(kge1,2)= ,i3,3(1pe20.10),6h eta= , 1 e20.8) c--------do loop over energies do 213 ix=1,7 if(drr(ix).eq.0.d+00) go to 213 xeri(1)=drr(ix) xeri(2)=eim/ix write(6,880)xeri(1),xeri(2) call snsqe(detcal,jac,2,2,xeri,detri,tol,-1,info,w,31,maxit) if (info .eq. 0)then write(6,100) 100 format(1h0,'improper input parameters in search') else if (info .eq. 1) then write(6,101) tol 101 format(1h0,'normal termination, relative error is at most', 1 f13.10) else if (info .eq. 2)then write(6,102) maxit 102 format(1h0,'iterations exceed ',i4) else if (info .eq. 3) then write(6,103) tol 103 format(1h0,f13.10,'tol too small.no further 1improvement possible in xeri') else if (info .eq. 4) then write(6,104) 104 format(1h0,'iteration not making good pr 1ogress') endif nmin=l+1 if(ix.ne.1) go to 211 do 223 nbohr=nmin,4 esch=(-ar11/2.)*(alfanz/nbohr)**2 dri(nbohr)=esch eta=nbohr-(fj+.5d+00)+sqrt((fj+.5d+00)**2-(alfanz)**2) wron1=(alfanz/eta)**2 eta=sqrt(1.+wron1) wron2=ar11*(-.5*wron1+wron1**2/8.-5.*wron1**3/16.) ekge=ar11*(1./eta-1.) write(6,226) nbohr,esch,ekge,wron2,eta 223 continue 211 epsev=(dri(ix)-xeri(1)) gamev=-xeri(2)*2. esch=-epsev write(6,881) xeri(1),xeri(2),epsev,gamev,esch,xeri(2) 881 format('0E(MeV)=',2f16.14,' eps,gamma(eV)=',2(6pf10.1)/ 1 ' Del E(ev)=',2(6pf10.1)) 213 continue go to 10 233 write(6,841) 841 format(' **************eof in main***********') stop end C================================================================= subroutine detcal(nrei,xeri,detri,iflag) c version modified sept85 for antiproton bs, 8/86 for reltv implicit real*8 (a-h,o-y) integer mypvt(384) complex*16 detin(2),work(384) complex*16 h, x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12, 1 x13,x14,x15,x16,x17 complex*16 zi,zw(10), za(384,384),zdet 1 ,den(384,3),sum1,sum2,sum3,zk22,zk33,zk11,zk2,zk1,zk3 2 ,ze,zs,ar11,ar22,ar33 real*8 kapg,mn,mn2,k0,ko,ko2,k,k2,kk,mpi,mpi2,mr,k02 1, imtij,bnuc,bn,retij,xbarc,xi,xpi,xn,xgam,zcm,zch,theta, 1dsig,bnucf,bnf,dflip,dnof,wsp,wsn,achp,achn,acmp,acmn,ymin1, 2ymin2 c NB NEED MATCH DIMEN AND NO ENTRIES IN DATA STM FOR IBM dimension m1(100),m2(100),y(50),gin(14,100),h(15,20),kapg(100 1), x1(15),x2(1),x3(4), 2 x4(3),x5(10),x6(8),x7(6),x8(6),x9(2),x10(9),x11(1), 3 x12(1),x13(1),x14(1),x15(4),x16(4),x17(2),intger(384),mm1(384), 4 mm2(384),xeri(2),detri(2) dimension kk(130),u(36864,2),f(36864),gp(66),wt(66) 1,ptptw(66),wa(130) dimension theta(95), yy(50), dflip(95), dnof(95) dimension nifty(20),bnuc(20),bn(16),retij(14),imtij(14),bnucf 1(20), dsig(95),chart(61,61),polar(61),denom(14,25),ecmden(25), 2 bnf(16), 3 reul(25,2), aimul(25,2) C NEW DIMENSIONS AND EQUIVL FOR BS,BE CAREFUL WITH U NEED C F CAN BE USED.nb order u(11) u(21)...u(n1) u(12) u(22) ..u(n2) c common /nlspfl/ gin,kapg,denom,ecmden,fmn,rew(14,100),aimw(14,100) 1,akapd(100),nkapgs,necms,dummy(1)/foptp/f common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,xbarc,xi,xpi,xn,nz,nes 1,nwaves,nifty,na common /sec4/ achp,acmp,wsp,achn,acmn,wsn/sec6/rcoul,rcut common /absp/ b0r,b0i,c0r,c0i common/det/mr,alfanz equivalence(u(1,1),y(1)),(u(51,1),yy(1)),(m2(1),dsig(1)),(chart(1, 11),u(1,1)), (gin(1,1),u(1,1)), (aimw(12,34),dnof(1)), (aimw(14,41) 2,dflip(1)),(aimw(2,49),polar(1)),(gin(1,1),reul(1,1)),(gin(1,5), 3aimul(1,1)),(x1(1),h(1,1)),(x2(1),h(1,2)),(x3(1),h(1,3)), 4(x4(1),h(1,4)),(x5(1),h(1,5)),(x6(1),h(1,6)),(x7(1),h(1,7)), 5(x8(1),h(1,8)),(x9(1),h(1,9)),(x10(1),h(1,10)),(x11(1),h(1,11)), 6(x12(1),h(1,12)),(x13(1),h(1,13)),(x14(1),h(1,14)),(x15(1),h(1,15) 7),(x16(1),h(1,16)),(x17(1),h(1,17)),(ko,k0),(ko2,k02), 9(m1(1),f(800)),(m2(1),f(1000)),(intger(1),f(1200)),(wa(1),mm1(1)) a,(m1(1),mm1(1)),(m2(1),mm2(1)),(za(1,1),u(1,1)), b(h(1,1),u(1,1)),(zw(1),f(1)) c-1 ibm version with (u,h) equivalence removed to avoid block data c c revised constants and masses nov 81 for bound s problem data alfa /7.29735d-3/,an,ap,akmn,ak0,asig,ndim/ 1 939.5731d+00,938.2796d+00,493.699d+00,497.67d+00,1193.35d+00, 2 384/,ix/1/ data hbarc,api,pi/197.3289d+0,139.5669d+0,3.1415926535d+0/,x1/ 1 10h(p,n) ,10h(k+,ko) ,10hpi+ ,10hpi- , 2 10hpi+,pi- ,10hpi cex ,10hpio ,10hproton , 3 10hcex(+pio) ,10hk+ ,10hk0 ,10hneutron , 3 10hk- ,10hk0bar ,10hp bar /, 4x2 /10h*5=be /,x3 /10hbs not res ,10hres not bs, 5 10hna ,10hna /,x4 /10htpin=const, 6 10hnlsp-g(p) ,10hlocal lapp /,x5 /10hso,kapo , 7 10hsin,kap0 ,10hso,no ang ,2*10hundef ,10he3b,no aay 8,10hundef ,10he3b, aay ,10hundefined ,10hso,aay /,x6/ 910hno spin ms,10hss+ds ,10hss+ds+ts ,10hspin,ms , a 10hus=uc=0 ,10hss+unitary,10hus=0,ucne0, a 10huim = 0 / data x7/ b 10hno abs/s-f,10hno abs/s-s,10hno abs/s-p ,10hvabs-old , c10hvabs-yool ,10habs-old,sp/,x8 /10hno del/fld,10hfld,no del , d10hdel,no fld,10hdel in,fld,10hk+ martin ,10hmartin,fld /,x9/ e 10hno pauli ,10hpauli tpin /,x10/10hno coul ,10h+fc+phase , f10h+fc ,10hexact coul ,10hexact+f(q) ,10hbs,pt coul, f 10hbs,sph col,10hbs,real cl,10hbs,no coul/,x11/ g10hpi0,pi-,k+/,x12 /10hpi+,ko,k0b/,x13 / h 10hpi-channel /,x14 /10hpi+channel/,x15 /10h1 chanel=p, h 10h2 chan=pip,10h3 chn=ppin, i10h2 chan=pn /,x16 /10hrhog,rho2g,10hrhows,rh2g,10hrhows,r2ws, j10hhe/d /,x17 /10hkmt ,10hwatson / if(ix.ne.1) go to 150 10 continue C MN HERE REFERS TO NUCLEUS MASS,XN (IN COMMON) TO NUCLEOn C READ MOMENTA IN CENTRE _OF_MASS C NR.LE.0 IS NONRELATIVISTIC NR .GT. 0 > RELATIVISTIC CASe read(5,*,end=650) nr,lxmax write (6,740) nr,lxmax read (5,*) ngp,kode,b,nang,ymin1,ymin2 if (nang.eq.0) nang = 3 write (6,760) ngp,kode,b,nang,ymin1,ymin2 if (ngp.le.0) go to 650 write (6,790) ngp C -----READ IN DATA read (5,*) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut write (6,880) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut read (5,*) nz,na,(nifty(n),n=1,20) write (6,830) nz,na,(nifty(n),n=1,20) nifty(15)=abs(nifty(15)) C PRINT OUT NIFTY do 25 nif=1,17 jj=0 j = nifty(nif)+1 if (nif.eq.1) j = nifty(nif)+3 if (nif.eq.2.or.nif.eq.11.or.nif.eq.12.or.nif.eq.13.or.nif.eq. 1 14 )jj=1 if (jj.eq.1) j=1 write (6,710) nif,nifty(nif),h(j,nif) if (jj.eq.1) write(6,711) 25 continue read (5,*) nes,nwaves,b0r,b0i,c0r,c0i write (6,910) nes,nwaves,b0r,b0i,c0r,c0i if (b0i.gt.1.d-03.or.c0i.gt.1.d-03) go to 30 c use old values if nothing is read in b0r = -0.04d+00 b0i = 0.04d+00 c0r = 0.d+00 c0i = 0.08d+00 write (6,910) nes,nwaves,b0r,b0i,c0r,c0i if (nwaves.lt.0)write(6,884)nwaves 884 format(' Nwaves lt 0 for % partial absorption, Nwaves=',i4) 30 if (na.ne.3) go to 40 c square input size params for he3(prevoius read in squaed) achp = achp*achp acmp = acmp*acmp wsp = wsp*wsp write(6,880)achp,acmp,wsp 40 continue mpi=api if ((nifty(1).eq.7).or.(nifty(1).eq.8).or.(nifty(1).eq.-1)) mpi = 1 akmn if(nifty(1).eq.10.or.nifty(1).eq.11)mpi=akmn if (nifty(1).eq.12)mpi=ap zcm = acmp amass=na zch = achp zws = wsp nifty1 = nifty(1) c for pi 0(ko) , use no coul if not b.s. if(nifty(10).gt.4) go to 45 if (nifty1.eq.4.or.nifty1.eq.8.or.nifty1.eq.11) nifty(10) = 0 C FOR B S ALWAYS USE WATSON POTEMNTIAl 45 if(nifty(10).ge.5)nifty(17)=1 c c DEFINE MASSES WHICH DONT CHANGE WITH ENERGY c c mn here refers to nucleus mass,xn is nucleon mpi2 = mpi*mpi zi=(0.d+00,1.d+00) xn = (nz*ap+(na-nz)*an)/amass mn=xn*amass mn2=mn*mn xbarc = hbarc xi = pi xpi = mpi mr = mn*mpi/(mn+mpi) ar11=mr ar112=ar11*ar11 alfanz=nz*alfa ar2=mr*mr c need readjust mass for other nuclei C SPECIAL COUPLED CHANNELS SECTION e12=akmn+ap-api-asig e13 = akmn + ap -ak0 - an if (nifty(1).eq.12)then e12=0. e13=ap +ap -an -an ak0=an akmn=ap end if c for n.r. case use haw's reltivis sig pi channel reduced mass am1=akmn+ap am2=asig+api am3=ak0+an am12=am1**2 am22=am2**2 am32= am3**2 am1m=(akmn-ap)**2 am2m=(asig-api)**2 am3m =(ak0 -an )**2 a21=asig**2 a22=api**2 a31=ak0**2 a32=an**2 ar22=(am1+am2)*(am12-am2m)/8.d+00/am12 c set booby trap on ar22 for pbar if(nifty(1).eq.12)ar22=0. ar222=ar22*2.d+00 api2=api*api ar33=ak0*an/(ak0+an) ar332=ar33*2.d+00 ak02=ak0**2.d+00 an2=an**2.d+00 asig2=asig**2.d+00 65 continue if(nr.le.0.) write(6,1000) if(nr.gt.0.) write(6,980) write (6,830) nz,na,(nifty(n),n=1,20) write (6,780) mpi,amass,mn write (6,781)ar11,ar22,ar33 rhl2 =4.d+00*mr/pi rhl22=4.d+00*ar22/pi rhl23=4.d+00*ar33/pi rhl3 = 2.d00/pi rhl4 = 2.*rhl3 wron1=2.d+00*mn wron2=2.d+00*mpi wron3=mn/10.d+00 wron4=mpi/10.d+00 ldmax=lxmax if(lxmax.eq.0)ldmax=1 lborn=ldmax xgam=0.d+00 c include coulomb, set uop zp = 1 nzp = nifty(1)+1 go to (150,130,140,140,140,150,140,150,140,140,130,140,130), nzp 130 zp = -1 go to 150 140 zp = 0. c-----------come here directly on all but 1st call----------------- 150 continue C bs ENERGY IS ke IN com FOR BOTH rEL AND nR c need redefine for scatt where KE is lab KE e=xeri(1) eimin =xeri(2) zdet=zi*eimin ze=e+zdet e2=e+e12 e3=e+e13 if(nr.gt.0)goto 156 c calc complex*16 e s and momenta for chan 2 zk11= (e+zdet)*2.*mr zk22=(e2+zdet)*ar222 zk33=(e3+zdet)*ar332 go to 157 c RELATIVISTIC ENERGY c take sqrt(s) = e + m1 +m2 as def of E, same in all channels c zk11 is now k2 in channel 1, NOT energy, ar is complex mr 156 ze = e + zdet +mn +mpi zs = ze**2 zk11 = (zs-am12)*(zs-am1m)/4./zs zk22 = (zs-am22)*(zs-am2m)/4./zs zk33 = (zs-am32)*(zs-am3m)/4./zs ar11 = sqrt((zk11 + mn2)*(zk11+mpi2))/ze ar22 = sqrt((zk22 + a21)*(zk22+a22))/ze ar33 = sqrt((zk33 + a31)*(zk33+a32))/ze if(ix.eq.1)write(6,781) ar11,ar22,ar33 157 continue c place complex momentum in quadrant II for BS analytic continuation c place complex momentum in quadrant IV for Resonance (nif(3)=1) x = zk11 yi = -zi*zk11 th = atan2(yi,x) + 2.*pi if(nifty(3).eq.1) th = atan2(yi,x) zk1 = sqrt(sqrt(x**2 + yi**2))*(cos(th/2.) + zi*sin(th/2.)) x = zk22 yi = -zi*zk22 th = atan2(yi,x) +2.*pi if(nifty(3).eq.1) th = atan2(yi,x) zk2 = sqrt(sqrt(x**2 + yi**2))*(cos(th/2.) + zi*sin(th/2.)) x = zk33 yi = -zi*zk33 th = atan2(yi,x) +2.*pi if(nifty(3).eq.1) th = atan2(yi,x) zk3 = sqrt(sqrt(x**2 + yi**2))*(cos(th/2.) + zi*sin(th/2.)) if(ix.eq.1) 1 write(6,108) e,eimin,e2,eimin,e3,eimin,zk1,zk2,zk3 c theta(1),(2) used in call to optpbs if whant complex momentum theta(1)=zk2 theta(2)=-zi*zk2 c C-----------SET UP GRID POINTS--------------------------------------- c n1 = ngp+1 n2 = n1+n1 a = zk2 a=abs(a) c special check if b=0,use different distrb of points if(ix.gt.1)go to 89 if (b.eq.0.) a = 200. if(kode .lt. 0) go to 84 if(ngp.gt.48) go to 82 call gauss (ngp,kode,a,b,gp,wt) go to 83 c ngp gt 48,redistrb pts, 0-ko ,k0-inf(half up to 10*ko) 82 npt1=24 npt2=ngp-npt1 b=a call gauss(npt1,011,0.d+00,b,gp,wt) if(npt2.gt.48)write(6,107) b=8.*a call gauss(npt2,041,b,a,dflip,dnof) ngp=npt1+npt2 write(6,730)ngp,npt1,npt2 do 81 i=1,npt2 ii=i + npt1 gp(ii)=dflip(i) 81 wt(ii)=dnof(i) go to 83 107 format(20h NPT2 GT 48 IN MAIN ) c special bs point distb for kode lt 0- 84 cpmax=10.d+00**(-kode) catom=b write(6,883) kode,catom 883 format(' kode < 0, use bs points,cpmax=10**-kode, ' 1 ,'catom <,>0 :gausbs-rhl mod of kt', 2 'gausbs-kt version ,kode,catom=',i4,e12.4) nnuc=nang if(nang.eq.3)nnuc=0 cnucl=ymin2 csize=ymin1 write(6,104) ngp,nnuc,catom,csize,cnucl,cpmax 104 format(40h Ngp,Nnuc,Catom,Csize,Cnucl,Cpmax(MeV)= ,2i4,4e12.4/) if(catom.ge.0.) 1call gausb2(ngp,nnuc,catom,csize,cnucl,cpmax,gp,wt,ptptw) if(catom.lt.0.) 1call gausbs(ngp,nnuc,-catom,csize,cnucl,cpmax,gp,wt,ptptw) 83 continue write(6,105)(gp(i),i=1,ngp) write(6,106)(wt(i),i=1,ngp) 105 format(8h points= ,10e12.4) 106 format(8h weight= ,10e12.4) 89 continue sum1=0. sum2=0. sum3=0. if (nr.gt.0) go to 100 c ----------------------------------------------------------------- c NONRELATIVISTIC DENOMINATOR c ----------------------------------------------------------------- do 90 i1=1,ngp k = gp(i1) kk(i1) = k k2 = k*k ak2w=k2*wt(i1) den(i1,1)=ak2w/(k2-zk11)*rhl2 if(nifty(15).eq.0) go to 90 den(i1,2)=ak2w/(k2-zk22)*rhl22 den(i1,3)=ak2w/(k2-zk33)*rhl23 c 2 closed channels, replace den2 with den3 and dont use channel 3 c closed channels: (K-p,Kobar-n) or (p-pbar,n-nbar) if(nifty(15).eq.3)den(i1,2)=den(i1,3) if(nifty(15).eq.3)den(i1,3)=0. c include Prin Value sum in denom for open channel if(e .gt.0.)sum1 = sum1 +wt(i1)/(k2-zk11) if(e2.gt.0.)sum2 = sum2 +wt(i1)/(k2-zk22) if(e3.gt.0.)sum3 = sum3 +wt(i1)/(k2-zk33) 90 continue c SET DEN(N1) WITH DEN2 CONTAINING BOTH I EPS AND p SUBTRN c can use this procedure for any real/complex energy, here use reE>0 switch if(e .lt.0.)den(n1,1)=0. if(e2.lt.0.)den(n1,2)=0. if(e3.lt.0.)den(n1,3)=0. kk(n1)=zk1 if(nifty(15).eq.0.and.e.lt.0.) go to 120 c channel 2 (usually the only open one) if(e2.gt.0.)then c channel 2 open kk(n1)=zk2 zw(1)= - sum2*rhl22*zk22 zw(2)= + zi*ar222*zk2 den(n1,2)= zw(1) +zw(2) endif if(e3.gt.0.)then c channel 3 open zw(1)= - sum3*rhl23*zk33 zw(2)= + zi*ar332*zk3 den(n1,3)= zw(1) +zw(2) endif if(e.gt.0.)then c channel 1 open zw(1)= - sum1*rhl2*zk11 zw(2)= + zi*ar2*zk1 den(n1,1)= zw(1) +zw(2) endif c for closed channels set den(n1,channel)=0 c for 2 channel case (nbar-n or kobar-n) switch channel 3 to 2 if(nifty(15).eq.3)then den(n1,2)=den(n1,3) den(n1,3)=0. endif go to 120 c --------------------------------------------------------------- c RELATIVISTIC DENOMINATOR c --------------------------------------------------------------- 100 continue do 110 i1=1,ngp k = gp(i1) kk(i1) = k k2 = k*k ak2w=k2*wt(i1) emn = sqrt(k2+mn2) epi = sqrt(k2+mpi2) c special low momentum expansion, maybe do for other channels ? if(k.lt.wron3)emn=(k2/wron1)*(1.-k2/wron1/wron1)+mn if(k.lt.wron4)epi=(k2/wron2)*(1.-k2/wron2/wron2)+mpi den(i1,1)=ak2w* rhl3/(emn+epi-ze) if(nifty(15).eq.0)goto 110 den(i1,2)=ak2w*rhl3/(sqrt(k2+a21)+sqrt(k2+a22)-ze) den(i1,3)=ak2w*rhl3/(sqrt(k2+a31)+sqrt(k2+a32)-ze) c 2 closed channels, replace den2 with den3 and dont use channel 3 c closed channels: (K-p,Kobar-n) or (p-pbar,n-nbar) if(nifty(15).eq.3)den(i1,2)=den(i1,3) if(nifty(15).eq.3)den(i1,3)=0. c include Prin Value sum in denom for open channel if(e .gt.0.)sum1 = sum1 +wt(i1)/(k2-zk11) if(e2.gt.0.)sum2 = sum2 +wt(i1)/(k2-zk22) if(e3.gt.0.)sum3 = sum3 +wt(i1)/(k2-zk33) 110 continue c set den(n1) with den2 containing both i eps and P subtrn if(e .lt.0.)den(n1,1)=0. if(e2.lt.0.)den(n1,2)=0. if(e3.lt.0.)den(n1,3)=0. kk(n1)=zk1 if(nifty(15).eq.0.and.e.lt.0.) go to 120 if(e2.gt.0.)then c channel 2 open (the usual case for k-) kk(n1)=zk2 den(n1,2)= (-sum2*rhl4*zk22 +zi*2.*zk2)*ar22 endif if(e3.gt.0.)then c channel 3 open den(n1,3)= (-sum3*rhl4*zk33 +zi*2.*zk3)*ar33 endif if(e.gt.0.)then c channel 1 open den(n1,1)= (-sum1*rhl4*zk11 +zi*2.*zk1)*ar11 endif c for closed channels set den(n1,channel)=0 c for 2 channel case (nbar-n or kobar-n) switch channel 3 to 2 if(nifty(15).eq.3)then den(n1,2)=den(n1,3) den(n1,3)=0. endif c---------------end relativistic denom------------------------------- 120 if(ix.eq.1.or.e.gt.0..or.e2.gt.0..or.e3.gt.0.) 1 write(6,91) (den(n1,i),i=1,3) if(ix.ne.1)go to 203 ii=3 if(nifty(15).eq.0)ii=1 do 202 i=1,ii 202 write(6,102) (den(i1,i),i1=1,n1) 102 format(' Zden(3/1,i)=',10e12.4) c C----- DO LOOP OVER L FOR THIS E AND STATE c 203 l = ldmax-1 if(ix.eq.1)write (6,1070) l c C BOUND STATE OR EIGENENERGY CALCULATION ( A* RESONANCES ) c calc bound state e's only for ldmax c theta(3)=10. call optpbs(kk,n1,wt,ldmax,ldmax,lborn,theta,e,eimin) c for three coupled channels possible, need 3*maxngp as dimen c nunit=ngp or n1 is size of block matrix nunit=ngp if(nifty(15).eq.1.or.nifty(15).eq.2)nunit=ngp+1 nn=nunit ng2=nunit*2 if(nifty(15).eq.2)nn=3*nunit if(nifty(15).eq.1.or.nifty(15).eq.3)nn=ng2 c WRITE OUT MATRICs if(ix.ne.1)go to 208 write(6,229) i=nn-1 if (nn .lt.6) then do 9995 ii = 1,nn 9995 write(6,227) ii,(za(ii,ij),ij=1,nn) else do 217 ii=1,nn,i 217 write(6,227) ii,(za(ii,ij),ij=1,nn) endif ii=nn 228 format(1h0,' 1-G(complex*16)*V matrix before determinant' 1 ,' found') 229 format(1h0,' complex*16 potential matrix before 1-GV formed') 227 format(1h0,i3,8e15.6/200(3x,8e15.6/)) 208 continue C DETERMINE CHANNEL FOR DEN + ROW OF SUBMATRX irow=1 do 219 ii=1,nn if(ii.gt.nunit)irow=2 if(ii.gt.ng2)irow=3 do 216 ij=1,nn jjj=ij if(ij.gt.nunit)jjj=ij-nunit if(ij.gt.ng2)jjj=ij-ng2 icol=1 c revised version with g1 and g2 switched if(ij.gt.nunit)icol=2 if(ij.gt.ng2)icol=3 za(ii,ij)=za(ii,ij)*den(jjj,icol) 2140 continue if(ii.eq.ij)za(ii,ii)=za(ii,ii)+1. 216 continue 219 continue C----------DO LOOP OVER II,IJ ENDs c write out matrics if(ix.gt.1)go to 298 write(6,228) i=nn-1 if (nn .lt. 6)then do 9996 ii=1,nn 9996 write(6,227) ii,(za(ii,ij),ij=1,nn) else do 297 ii=1,nn,i 297 write(6,227) ii,(za(ii,ij),ij=1,nn) endif if(nifty(6).eq.6.or.nifty(6).eq.7)go to 298 298 continue 212 continue C SEARCH FOR DET=0 WITH LANDE SUBRTN IN VCOUL call cgefa(za,384,nn,mypvt,info) call cgedi(za,384,nn,mypvt,detin,work,10) c calculate deter zdet=(1.,0.) do 223 i=1,nn ipvt=wa(i) if(ipvt.ne.i)zdet=-zdet zdet=zdet*za(i,i) 223 continue if(ix.eq.1)ix=2 rpow = detin(2) zdet = detin(1)*(10.d00**rpow) detri(1)=zdet detri(2)=-zi*zdet write(6,224)zk1,e,eimin,zdet 224 format(' zk1=',2(1pE11.3),' Er,Ei=', 1 2(1pe18.8),' zdet=',2(e15.3)/) return C _______________ FORMATS 650 write(6,841) stop c 108 format(' zE1, zE2, zE3 =',3(2e12.4,4x)/ 1 ' zk1, zk2, zk3 =',3(2e12.4,4x)) 91 format(' den(n1,1-3) (P, Del parts) =',3(2e12.4,4x)) 710 format(7h Nifty(,i2,2h)=,i3,10x,a8,a2) 711 format(1h+,36x,5hshift ) 730 format (4x,4i4) 740 format (' Nr Lmax =',5i4) 760 format (2i5,f10.2,i5,2f10.0) 780 format(1h ,3x,' Mproj, At, MA =', 1 5x,3(f8.3,11x)) 781 format(1h ,3x,' ar11, ar22, ar33=', 1 5x,3(2f8.3,3x)/) 790 format (1h ,18hNo of grid points=,i4) 830 format (1h ,2i5,1x,i2,i2,8i1,10i3) 841 format(12h EOF IN MAIN ,i3) 850 format (i3,f15.5) 851 format(29h g(p)-unit4/d(e)-unit2 print ,/i3,f15.5) 860 format (f10.4/4(6f13.6,/),4f13.6) 861 format(1x ,150(f10.4/6f13.6/6f13.6/2f13.6/)) 870 format (e10.4/8e10.3/6e10.3) 871 format(1x ,150(f10.4/8e10.3/6e10.3/)) 880 format (8f10.5) 910 format (28h Nes,Nwaves,b0r,b0i,c0r,c0i=,2i5,4f10.5) 940 format (1h0,10x,8h Energy,10x,15h Momentum ,/) 950 format (1h ,7x,f14.4,10x,f14.4/) 960 format (1h ,5f14.4) 980 format (1h0,31h RELATIVISTIC CALCULATION) 1000 format (1h0,' NONRELATIVSTC CALC') 1070 format (1h0,30h----- Angular Momentum = ,i4,5h-----) end C================================================================= subroutine gausb2(npt,nnuc,catom,csize,cnucl,cpmax,p,wp,ppw) c rhl mod of /tabakin + kwon trig scaled grid c bound state problem implicit real*8 (a-h,o-z) dimension p(132),wp(132),ppw(132),x(128),wx(128) if(nnuc.ne.0) rata=catom/csize if(nnuc.eq.0) rata=0. c rata = catom/csize natom=npt-nnuc c atomic region if (natom.eq.0) go to 12 call gauss(natom,010,-1.d+00,1.d+00,x,wx) do 10 i=1,natom z1=x(i) p(i)=catom*(1.+z1)/(1.-(1.-2.*rata)*z1) dpdx=2.*catom*(1.-rata)/((1.-(1.-2.*rata) 1 *z1)**2) wp(i)=dpdx*wx(i) 10 ppw(i)=p(i)*p(i)*wp(i) if(nnuc.eq.0) go to 50 12 continue c nuclearr region grid nuc=nnuc/2 if(nuc.ne.natom) call gauss(nuc,010,-1.d+00,1.d+00,x,wx) n1=natom+nuc n2=natom+1 do 20 j=n2,n1 p(j)=(csize+cnucl+(cnucl-csize)*x(j-natom))/2. dpdx=(cnucl-csize)/2. wp(j)=dpdx*wx(j-natom) ppw(j)=p(j)*p(j)*wp(j) 20 continue n3=natom+2*nuc if(n3.ne.npt) write(6,900)npt,nnuc,nuc 900 format(' nnuc not div by 2 in gausb2,npt,nnuc,nuc=',3i3) do 30 j=n1+1,n3 p(j)=(cnucl+cpmax+(cpmax-cnucl)*x(j-n1))/2. dpdp=(cpmax-cnucl)/2. wp(j)=dpdp*wx(j-n1) ppw(j)=p(j)*p(j)*wp(j) 30 continue 50 continue return end C================================================================= subroutine gausbs(npt,nnuc,catom,csize,cnucl,cpmax,p,wp,ppw) c tabakin + kwon trig scaled grid c bound state problem implicit real*8 (a-h,o-z) dimension p(132),wp(132),ppw(132),x(128),wx(128) pid4=3.14159265/4. if(nnuc.ne.0) rata=catom/csize if(nnuc.eq.0) rata=0. c rata = catom/csize ratb=cnucl/cpmax natom=npt-nnuc c atomic region if (natom.eq.0) go to 12 call gauss(natom,010,-1.d+00,1.d+00,x,wx) do 10 i=1,natom z1=pid4*(1.+x(i)) p(i)=catom*sin(z1)/(cos(z1) + rata*sin(z1)) dpdx=catom*pid4/(cos(z1)+rata*sin(z1))**2 wp(i)= dpdx*wx(i) 10 ppw(i)=p(i)*p(i)*wp(i) if(nnuc.eq.0) go to 50 12 continue c nuclearr region grid if(nnuc.ne.natom) call gauss(nnuc,010,-1.d+00,1.d+00,x,wx) do 20 j=natom+1,npt z2=pid4*(1.+x(j-natom)) p(j)=csize+cnucl*sin(z2)/(cos(z2)+ratb*sin(z2)) dpdx=cnucl*pid4/(cos(z2)+ratb*sin(z2))**2 wp(j)=dpdx*wx(j-natom) ppw(j)=p(j)*p(j)*wp(j) 20 continue 50 continue return end C================================================================= * subroutine gauss(npts, kode, a, b, xs, wts) c head.h Version 1.1 c*********************************************************************** c c Oregon State University Nuclear Theory Group c Multi-Channel Anti Kaon -- Nucleon Code c Guangliang He & Rubin H. Landau c c*********************************************************************** c fortran.h Version 1.4 c*********************************************************************** c gauss.F Version 1.3 c Latest Revision Date 10:35:54 8/16/91 c*********************************************************************** c subroutine gauss(npts,job,a,b,x,w) c c*********************************************************************** c rescale rescals the gauss-legendre grid points and weights c c npts number of points c kode = 0 rescalling uniformly between (a,b) c 1 for integral (0,b) with 50% points inside (0, ab/(a+b)) c 2 for integral (a,inf) with 50% inside (a,b+2a), with c rational scaling c kode = 3 for integral (a, inf) with 50% inside (a, b+2a), with c atan() scaling. (a+b) has to great than zero. c c x, w output grid points and weights. c*********************************************************************** c integer npts,m,i,j,kode double precision x(npts),w(npts),a,b,xi double precision t,t1,pp,p1,p2,p3,aj double precision eps,pi,zero,two,one,half,quarter,pio4 parameter (pi = 3.14159265358979323846264338328, eps = 3.0d-15) parameter (zero=0.0d0,one=1.0d0,two=2.0d0) parameter (half=0.5d0,quarter=0.25d0) c c *** FIRST EXECTUABLE ************************************************* c kode=0 m=(npts+1)/2 do 10020 i=1,m t=cos(pi*(i-quarter)/(npts+half)) 10000 continue p1=one p2=zero aj=zero do 10010 j=1,npts p3=p2 p2=p1 aj=aj+one p1=((two*aj-one)*t*p2-(aj-one)*p3)/aj 10010 continue pp=npts*(t*p1-p2)/(t*t-one) t1=t t=t1-p1/pp c if(abs(t-t1).gt.eps) goto 10000 c x(i)=-t x(npts+1-i)=t w(i)=two/((one-t*t)*pp*pp) w(npts+1-i)=w(i) 10020 continue c c*********************************************************************** c rescale the grid points c*********************************************************************** c if (kode.eq.0) then c scale to (a,b) uniformly do 10030 i=1,npts x(i)=x(i)*(b-a)/two+(b+a)/two w(i)=w(i)*(b-a)/two 10030 continue c elseif (kode.eq.1) then c scale to (0,b) with 50% points inside (0,ab/(a+b)) do 10040 i=1,npts xi=x(i) x(i)=a*b*(one+xi)/(b+a-(b-a)*xi) w(i)=w(i)*two*a*b*b/((b+a-(b-a)*xi)*(b+a-(b-a)*xi)) 10040 continue c elseif (kode.eq.2) then c scale to (a,inf) with 50% points inside (a,b+2a) do 10050 i=1,npts xi=x(i) x(i)=(b*xi+b+a+a)/(one-xi) w(i)=w(i)*two*(a+b)/((one-xi)*(one-xi)) 10050 continue else if (kode.eq.3) then c scale to (a, inf) with 50% points inside (a, b+2a) pio4 = atan(one) if (a+b .le. zero) stop 'gauss: a+b less than zero' do 10060 i = 1, npts xi = x(i) x(i) = a+(a+b)*tan(pio4*(one+xi)) w(i) = w(i)*(a+b)*pio4 & /(cos(pio4*(one+xi))*cos(pio4*(one+xi))) 10060 continue else stop 'gauss: Wrong value of kode' endif c return end C================================================================= subroutine optpbs (kk,n1,wt,lpia,lmaxin,lborn,fx,er,eim) c coupled channels version...see dimens + equiv c modified sept85 for anti protons, n.b. ze now passed c *** see r h landau, program lpott, 1981 implicit real*8 (a-h,k,o-z) real*8 imtij,mpi,mn,k,kp,kk,imvabs,k1 complex*16 vkp(3,3),zi,zk,zkp,zk2,zu(384,384),ze dimension rho(4,50), bi(100), vcoull(30), revabs(30,2), imvabs(30, 12) ,wt(128) dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij 1(14),imtij(14),kk(128),u(36864,2),f(18432,2),fx(30), 2 ui(384,384),ur(384,384),vr(6),vi(6) dimension ff(715,6) dimension np(5) common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,mn,nz,nes 1,nwaves,nifty,na common/nlspfl/u/foptp/f common /sec4/ achp,acmp,wsp,achn,acmn,wsn common /sec5/ drho(4,50),rho,xkp,xk,imax/sec6/rcoul,rcut equivalence (vkp(1,1),bnuc(1)),(iset,nes) ,(ur(1,1),u(1,1)), 1(ui(1,1),u(1,2)),(f(1,1),vr(1)),(f(7,1),vi(1)),(f(1,1),ff(1,1)) equivalence(zu(1,1),u(1,1)) data vcoull,revabs,imvabs/150*0.0e+00/,ikp,zi/0,(0.e+00,1.e+00)/ 1 ,nint,eri,eimi/0,0.,0./ c c need reset dimensns and equivs as max no grid increase c should use n=ngp +1 for scattering c-----1st equiv---------------------------------------------------------- c ur+ui(3n,3n) , u( (3n)**2,2) c-----2nd equivalence---------------------------------------------------- c f(n(n+1)*2)-in main f(n(n+1),2)-in optpbs c ff(ndimf/6, 6) also needed ,and equivalenced to f(1) c for equiv need length+1 to overlay c also see intialize statements which now change with dimensions 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 c comple k2 needed for kp channel 2-open,fx(1),(2)=re,im k c fx(3) lt 1 define ur,ui ge 1 define complex*16 zu c zk2=fx(1)+zi*fx(2) nx=fx(3) nint=nint+1 if(er.eq.eri.and.eim.eq.eimi)goto 9 c energy has changed, recalculate u for e-dependent potential c eri=er eimi=eim go to 8 9 if (lpia.ne.1.or.nint.ne.1) go to 250 c 1st call to optpbs, calculate all u s c intialize (cant use data for this due to commoms) 8 do 7 i=1,2 do 7 j=1,18432 7 f(j,i)=0. do 6 j=1,36864 u(j,1)=0. 6 u(j,2)=0. if(kk(n1).gt.0.)k1=kk(n1) if(nifty(10).gt.4)k1=kk(1) if (nint.eq.1) write(6,410) k1 k = -200. c ikp = k- p switch for special sub for opt pot if(nifty(1).ge.10.and.na.eq.1)ikp=1 if(ikp.eq.1) go to 5 c if(nifty(6).ne.4.and.nifty(6).ne.6)call tncm (k) 5 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. if(nz.ne.0)ap = (nz-2.)/6./nz*(acmp/achp)**2 if (nz.le.2) ap = 0. if(nn.ne.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.le.2) aovera = 1. c set u(strong)=0 if ( nifty(6).eq.4) aovera = 0. xhl = aovera*2.*pi*pi aparam = 1. i = nifty(1)+1 go to (30,10,20,20,20,30,20,30,20,20,10,20,10), i 10 aparam = -1. go to 30 20 aparam = 0. 30 if(rcoul.lt.1.e-10 )rcoul=1.e-6 aparam = (1.5*nz*aparam/rcoul)*hbarc*7.29735e-3 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) lmax = ldmax+5 if(nifty(6).eq.6) go to 110 if (nifty(16).ne.0) go to 90 if(ikp.gt.0)go to 110 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. c 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) go to 60 if ((achp.ne.achn).or.(acmp.ne.acmn)) go to 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 go to 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. c 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) go to 110 c spin for c13 nff = 4 zs = k*kp*rhl5 xs = (k*k+kp*kp)*rhl5/2. c 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).and.nint.eq.1) 1 write (6,350) k,(rho(4,i),i=1,lmax) go to 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 c call rhokff (nff) if (nff.ne.1) go to 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).lt.3.or.nifty(10).eq.8) go to 120 c include cutoff coulomb potential in momentum space c nb vl contains scale factor ,i.e. it is v(j=l+1/2) ak = k akp = kp c dec83 add,no coul for open chan if(i1.eq.n1.or.i2.eq.n1)goto120 if(ak.le.0..or.akp.le.0.)go to 120 if(i1.ne.i2) go to 118 c k = kp, need special vcoul for bs so remove singularity if(nifty(10).lt.5) go to 118 call vcdiag(ak,i2,aparam,xhl1,xhl2,vcoull,ldmax, 1 nif,n1,kk,wt) go to 120 118 call vcoulb (ak,akp,aparam,xhl1,xhl2,vcoull,ldmax,nif) 120 if (nifty(7) .lt. 3 .and. nifty(7).ne.5) go to 130 c include absorption on 2 nucleons c 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 c if(nifty(6).ne.6.and.ikp.ne.1)call tnuc (kp,k,kk(n1),xmn) c do loop over l opt c do 230 ldum=1,ldmax if ((ldum.gt.lborn).and.((i1+i2).ne.(2*n1))) go to 230 j1 = ldum urc = 0. vc=0 uic = 0. urs = 0. uis = 0. if (nifty(6).eq.3) rhl1 = j1*(j1-1.) if(nifty(6).eq.6.or.ikp.eq.1) go to 205 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 c******************* rhl = cgc2(j3,j2,j1) c spin independt potential if (nifty(1).eq.1.or.nifty(1).eq.8) go to 140 c pi plus scatterinf jpr = 4*j2-3 jpi = jpr+1 jnr = jpr+2 jni = jpr+3 go to 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) go to 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) go to 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)) 180 if (nifty(6).ne.3) go to 190 c ---- ---spin dependnet potential if (j1.eq.1) go to 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 c rhl=(2pi**2/2l+1)(a-1/a) 200 continue 205 continue c form potentials for state of definite j rhl = xhl/(2.*j1-1.) vc = 0. if (nifty(10).lt.3.or.nifty(10).eq.8) go to 210 vc = vcoull(j1)*aovera if (nifty(6).eq.6) vc = vcoull(j1) 210 if (nifty(7).ne.3.and.nifty(7).ne.4) go to 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 (nifty(7).ne.4.or.nifty(6).ne.3) go to 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 c nb the f is - imv as stored, latter reversed npot1 = 2*ldum f(npot1-1,1) = rhl*(urc+(j1-1)*urs)+vc f(npot1,1) = -rhl*(uic+(j1-1)*uis) ff1=f(npot1-1,1) ff2=f(npot1,1) c pure coulomb or non k-p if(nifty(6).eq.6 .or. ikp. ne .1) go to 226 c c special k-p/or pbar p/section coupled or uncoupled c nes=iset=switch for model (k-p) or channel(pbar) c only l=0 for k- 225 zkp=kp zk= k 227 continue call vkbarp(zkp,zk) f(npot1-1,1)=vkp(1,1)+f(npot1-1,1) ff1=f(npot1-1,1) c nb reverse sign of v here for storage as f f(npot1 ,1)=+zi*vkp(1,1)+f(npot1 ,1) ff2=f(npot1,1) if(nifty(15).eq.0) go to 226 ff(npot1-1,1)=ff1 ff(npot1,1)= ff2 c coupled channels kp c store matrices with new, nonlinear forms c not -imv any longer so need add 2 imv to correct the above f c nifty(15)=0 ..uncoupled 1...k-p+sigma-pi c 2...k-p,sig pi,ko/ n 3...k-p, k0/ n c nspin for coupled kp sigma-pi kobarn c 1---v11 2---v22 3---v33 c 4---v12 5---v13 6---v23 c ff2=ff2-2.*zi*vkp(1,1) ff(npot1,1)=ff2 i=nifty(15) go to (231,231,232),i c k-p plus sigma pi (coul only for k-p) 231 ff(npot1-1,2)=vkp(2,2) ff(npot1,2)=-zi*vkp(2,2) ff(npot1-1,4)=vkp(1,2) ff(npot1,4)=-zi*vkp(1,2) if(nifty(15).eq.1) go to 226 c kip, sigma-pi, kobar n 232 ff(npot1-1,3)=vkp(3,3) ff(npot1,3)=-zi*vkp(3,3) ff(npot1-1,5)=vkp(1,3) ff(npot1,5)=-zi*vkp(1,3) if (nifty(15).eq.3) go to 226 ff(npot1-1,6)=vkp(2,3) ff(npot1,6)=-zi*vkp(2,3) 226 nnn = npot1-1 if(kp.eq.k1.and.kp.eq.k.and.ldum.eq.1.and.nint.eq.1) 1 write(6,410)k1 if(k.eq.k1.and.ikp.eq.1.and.lpia.eq.2.and.nint.eq.1 1 .and.ldum.eq.2)goto228 if(k.eq.k1.and.ikp.eq.1.and.lpia.eq.2.and.ldum.eq.1 1 .and.ldum.eq.1)goto402 if(k.ne.k1.or.ldum.ne.1.or.nint.ne.1) go to 402 c write out potential matrices on printer (half offshell) 228 write (6,400)kp,ff1 1 ,vc,revabs(j1,1),revabs(j1,2),ff2,imvabs(j1,1), 2 imvabs(j1,2) n=nifty(15) if(i1.ne.1.and.i1.ne.8) go to 402 if(n.eq.1.or.n.eq.2)write(6,401)(i,ff(nnn,i),ff(npot1,i),i=2,4,2) if(n.eq.2.or.n.eq.3)write(6,401)(i,ff(nnn,i),ff(npot1,i),i=3,5,2) if(n.eq.2 )write(6,401)(i,ff(nnn,i),ff(npot1,i),i=6,6,1) 401 format(1x,i2,19x,e19.4,41x,e19.4) 402 if (nifty(6).ne.3) go to 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.and.nifty(15).eq.0) write (7) 1 (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) if(nifty(15).ne.0)write(7) ((ff(npot2-1,i),ff(npot2,i),i=1,6 1 ),npot2=2,npot1m,2) c do loop over kp and k ends 240 continue c c----------------read u from unit7 on nonintial calls c 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 c mod 9/85 so dont need lpia=1 to set up 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(15).ne.0)nspin=6 if (nifty(6).ne.3) go to 260 c spin-dependent case, now 4 v s adjacent for each l vs 2 if nf nspin = 2 npot = 2*npot c coupled channels kn, 12 v s adjacent,non imv sign reversal 260 ldum = 2*(lpia-1)*nspin if(nifty(15).ne.0) go to 345 if (lpia.gt.lborn) go to 290 c do loop over npts (i.e. records) c nb ldum contains both no of l's and nspin(pot types) to skip do 280 npt=1,npts if (lpia.ne.1) go to 270 read (7) (f(npt,i),i=1,npot) go to 280 c space to l value needed 270 read(7)(rhl,i=1,ldum),(f(npt,i),i=1,npot) 280 continue go to 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) (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 if(nifty(10).ge.5)go to 331 npot1 = (i2-1)*n2+i1 320 npot2 = npot1+n1*(n2+1) c imv:npot4 -imv:npot3 (the one stored for uncoupled channels) npot3 = npot2-n1 npot4 = npot1+n1 u(npot1,1) = f(npt,1) u(npot2,1) = u(npot1,1) c rhl change feb82 for bs if(nifty(6).eq.7)f(npt,2)=0. u(npot3,1) = f(npt,2) u(npot4,1) = -f(npt,2) c spin-dependnt potentials if (nifty(6).ne.3) go to 330 c special turn off of spin for bs to avoid mstrix def write(6,999) 999 format(' ********** spin section in opt no good at 999 ') c u(npot1,2) = f(npt,3) c u(npot2,2) = u(npot1,2) c u(npot3,2) = f(npt,4) if(nifty(6).eq.7) u(npot3,2)=0. 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)) go to 340 npot1 = (i1-1)*n2+i2 iflag = 1 go to 320 c bound state case, set up seperate ur,ui 331 if(nx.lt.1)ur(i1,i2)=f(npt,1) if(nx.lt.1)ui(i1,i2)=-f(npt,2) if(nx.gt.1)zu(i1,i2)=f(npt,1)-zi*f(npt,2) if(i1.eq.i2) go to 340 c off diagnol symmetrty if(nx.lt.1)ur(i2,i1)=ur(i1,i2) if(nx.lt.1)ui(i2,i1)=ui(i1,i2) if(nx.gt.1)zu(i2,i1)=zu(i1,i2) 340 continue return c c special coulpled channels calc....now for k-p,kobar-n, sigma-pi c 1st fill lower half of full matrix,then latter upper c then set ur,ui to super matrix (ur,ui equiv to the above u) c v1 v4 v5 c v4 v2 v6 c v5 v6 v3 c npt is order written out/read in from storage c here only for b.s. with on shell point now included as n1 345 ngp=n1 if(nifty(15).eq.3)ngp=n1-1 ng2=2*ngp ng3=3*ngp c npt is locatn in triangle matrix stored linearly c rearrang with (i,j) index and above symmetry if(lpia.gt.lborn) write(6,420) lpia,lborn 420 format(' lpia gt lborn in optpbs bs********',2i3) do 347 i1=1,ngp c skip on shell point if ko k- if(i1.eq.n1.and.nifty(15).eq.3) go to 347 npot1=i1+ngp npot3=i1+ng2 do 346 i2=1,i1 npot2=i2+ngp npot4=i2+ng2 if(lpia.eq.1)read(7)(vr(i),vi(i),i=1,6) if(lpia.ne.1)read(7)(rhl,i=1,ldum),(vr(i),vi(i),i=1,6) if(nifty(6).eq.7)then c set imv = 0 do 425 i=1,6 425 vi(i)=0. endif c off diagnol block matrices,build in internal symmetry if(nifty(15).eq.3) go to 421 if(nx.lt.1)ur(npot1,i2)=vr(4) if(nx.lt.1)ui(npot1,i2)=vi(4) if(nx.lt.1)ur(npot3,i2)=vr(5) if(nx.lt.1)ui(npot3,i2)=vi(5) if(nx.lt.1)ur(npot3,npot2)=vr(6) if(nx.lt.1)ui(npot3,npot2)=vi(6) if(nx.gt.1)zu(npot1,i2)=vr(4)+zi*vi(4) if(nx.gt.1)zu(npot3,i2)=vr(5)+zi*vi(5) if(nx.gt.1)zu(npot3,npot2)=vr(6)+zi*vi(6) go to 422 c n(15)=3 st up 2 block matrix 421 zu(npot1,i2)=vr(5)+zi*vi(5) 422 continue if(i1.eq.i2) go to 348 c off diagnol block submatrix in lower lh side,fill its c non diagnol elements via symmetry if(nifty(15).eq.3) go to 423 if(nx.lt.1)ur(npot2,i1)=vr(4) if(nx.lt.1)ui(npot2,i1)=vi(4) if(nx.lt.1)ur(npot4,i1)=vr(5) if(nx.lt.1)ui(npot4,i1)=vi(5) if(nx.lt.1)ur(npot4,npot1)=vr(6) if(nx.lt.1)ui(npot4,npot1)=vi(6) if(nx.gt.1)zu(npot2,i1)=vr(4)+zi*vi(4) if(nx.gt.1)zu(npot4,i1)=vr(5)+zi*vi(5) if(nx.gt.1)zu(npot4,npot1)=vr(6)+zi*vi(6) go to 348 423 if(nx.gt.1)zu(npot2,i1)=vr(5)+zi*vi(5) c block diagnol matrices (only lower part as rest to follow) 348 if(nx.lt.1)ur(i1,i2)=vr(1) if(nx.lt.1)ui(i1,i2)=vi(1) if(nx.lt.1)ur(npot1,npot2)=vr(2) if(nx.lt.1)ui(npot1,npot2)=vi(2) if(nx.lt.1)ur(npot3,npot4)=vr(3) if(nx.lt.1)ui(npot3,npot4)=vi(3) if(nx.gt.1)zu(i1,i2)=vr(1)+zi*vi(1) if(nifty(15).eq.3) go to 424 if(nx.gt.1)zu(npot1,npot2)=vr(2)+zi*vi(2) if(nx.gt.1)zu(npot3,npot4)=vr(3)+zi*vi(3) go to 346 424 if(nx.gt.1)zu(npot1,npot2)=vr(3)+zi*vi(3) 346 continue 347 continue istart=ng3 if(nifty(15).eq.1.or.nifty(15).eq.3) istart=ng2 c build up upper rh part of super matrix via symmetry do 349 i2=2,istart i22=i2-1 do 349 i1=1,i22 if(nx.gt.1)zu(i1,i2)=zu(i2,i1) if(nx.lt.1)ur(i1,i2)=ur(i2,i1) 349 if(nx.lt.1)ui(i1,i2)=ui(i2,i1) c redefine on shell v for channel 2 (others mult by 0 in main) c only second row of block matrices changed if(nint.eq.1.or.nifty(19).eq.0)go to 245 if(nifty(15).eq.3) go to 245 write(6,244) 244 format(' define on shell v for channnel2 ala e') istart=n1+1 iend=ng2 iend2=1 zk=kk(n1) if(nifty(19).eq.2)zk=zk2 np(1)=n1 np(2)=ng2 np(3)=ng3 do 242 ii=istart,iend idum=ii-n1 zkp=kk(idum) if(ii.eq.iend) iend2=n1 do 242 ij=1,iend2 if(ii.ne.iend) go to 243 c last lint, fill in entire n1 th row for all block 2 zk=kk(ij) if(ij.eq.n1.and.nifty(19).eq.2)zk=zk2 np(1)=ij np(2)=ij+n1 np(3)=ij+ng2 243 call vkbarp(zkp,zk) do 242 i=1,3 jjj=np(i) if(nx.lt.1)ur(ii,jjj)=vkp(2,i) if(nx.lt.1)ui(ii,jjj)=-zi*vkp(2,i) if(nx.gt.1)zu(ii,jjj)=vkp(2,i) if(nint.eq.2.and.nx.lt.1) 1write(6,246)ii,jjj,zkp,zk,ur(ii,jjj),ui(ii,jjj) 246 format(' i,j,zkp,zk,u=',2i5,6e15.6) 242 continue 245 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/6e16.7) 400 format (3x,3e19.4,2e11.3,e19.4,2e11.3) 410 format(28h0opt poten for k =k0(k1-bs)=,e18.5,13hand l=0(1-bs)/ 1 20x,2hkp, 2 13x,3hrev,16x,2hvc,15h rabs(n0 flip) ,11h rabs(flip),17x, 3 4h-imv,27h imvabs(noflip),imabs(flip)/) end c======================================================================= subroutine vcdiag(ak,i2,aparam,xhl1,xhl2,vl,ldmax,nif,n1,kk,wt) c vcdiag calculates the diagnol term of point coulomb pot c including the subtractions implicit real*8 (a-h,k,m,o-z) dimension vl(30),vcoull(30),ail(30),pl(30),kk(128),wt(128) data nit/0/,ail/0.,1.,1.2247449,1.3228549,1.3777022,1.4127050, 1 1.4369752,1.4547898,1.4684210,1.4791869,1.4879047,19*1.5/ if(nit.ne.0) go to 10 c 1st call intialize nit=1 pi2=3.14159265/2. ngp=n1-1 if(ldmax .gt. 11) write(6,900)ldmax 900 format(' ldmax too large for il s in vcdiag ,ldmax=',i3//) c c set up sl(p)...diff norm than tabakin slp=+(pi2)*aparam*xhl1/1.5 lmax=ldmax-1 10 ak2=ak*ak twoak=2.*ak do 5 l=1,ldmax 5 vl(l)=0. c do loop over grid points do 100 npt=1,ngp if(npt.eq.i2)go to 100 kn=kk(npt) zmn=(ak2+kn*kn)/twoak/kn call vcoulb(ak,kn,aparam,xhl1,xhl2,vcoull,ldmax,nif) c legendre poly and vl for all l s up to lmax pl(1)=1. pl(2)=zmn do 20 l=1,ldmax if (l.eq.1)go to 15 g=zmn*pl(l) pl(l+1)=2.*g-pl(l-1)-(g-pl(l-1))/l 15 vl(l)=vl(l) + vcoull(l)*wt(npt)/pl(l) 20 continue 100 continue do 25 l=1,ldmax 25 vl(l)=(slp*(pi2-ail(l))/ak- vl(l))/wt(i2) return end C================================================================= subroutine vcoulb(k,kp,a,rc,r,vl,nls,nif) c point coulomb potential....to add to reg vcoul implicit real*8 (a-h,k,m,o-z) dimension vl(30) zee=+(a/1.5)*rc/2. if (nif.ge.5)go to 100 c put reg sub here continue 100 continue c c point coulomb potenetial/non diagnol 110 fac=zee/k/kp z=(k*k+kp*kp)/2./k/kp call qlfun(nls,z,vl) do 120 i=1,nls 120 vl(i)=fac*vl(i) return end C================================================================= subroutine vkbarp(kp,k ) c sub to calculate the henley alberg willets (haw) kn and sigma pi c two body potentials using their params yet with the lpott conventions c also the schick-gibson potentials c so that the potential is added directly to the optpbs potent c nnb potentials are complex*16, have 3 channels, and are in /sec2/ c nb only s waves for the present form, no spin dependeece c nes in /sec2/ is iset, switch for the 3 potent sets(i=1) c c channel labels: 1 k-p 2 sigma pi (only i=0 part for haw) c 3 kobar-n c iset lt 0 is a switch for only k-p in this sub(used as prog check) c cc iset switch : for i=1, 1,2,3, refer to sets a, b, c (nes=switch) c iset= 30...use bag model i=0 via sub vkbag implicit real*8 (a,b,c-h,m,k, o-u,w-z) complex*16 vkp,l,bsq,l0,l1,b1,b0,bo,lo,g ,l0haw,l1haw,zi 1 ,kp,k,k2,kp2,rhl c 2 ,b1i.l1i dimension l0(3),lo(3),l1(5),b1(5),vkp(3,3),fl(3) dimension nifty(20),bnuc(20),bnucf(20),bn(16),bnf(16),retij(14), 1 aimtij(14) common/sec2/bnuc,bn,bnucf,bnf,retij,aimtij,hbarc,pi,mpi,mn,nz,nes 1 ,nwaves,nifty,na equivalence(l0,lo),(b0,bo),(bnuc(1),vkp(1,1)),(iset,nes), 1 (fl(1),l0(1)) data l0/(-326.,0.),(-112.,0.),(-412.,0.)/,b0,zi/(5.56,0.),(0.,1.)/ 1 , l1/(-232.,-64.9),(-388.,-175.),(-30.8,-22.4),2*(0.,0.)/, 2b1/(0.824,2.29),(4.,0.),(2.,0.),2*(0.,0.)/,nit /0 / g(bsq,l)= l/(bsq+k2)/(bsq+kp2) c intialize...convert to mev + do mult so no need repeat ea time if(nit.ne.0) go to 10 fourpi=4.*pi icoul=0 rt2=sqrt(2.) i1chan=1 acoup=1. if (iset.eq.6) then read(5,800) b0 read(5,800) l0(1) read(5,800) l0(2) read(5,800) l0(3) read(5,800) b1(1) read(5,800) l1(1) c ****** convert units to those of HAW 1980 do 15 i=1,3 l0(i)=l0(i)/b0 15 continue l1(1)=l1(1)/hbarc b0=b0/hbarc b1(1)=b1(1)/hbarc c ****** treat as HAW iset=1 endif c switch for pure elastic kp only if(iset .lt.0)i1chan=0 iset=abs(iset) ibag=0 iset=3 20 if(iset.eq.5)go to 7 if(iset.eq.4)go to 6 if(iset.ge.100)icoul=1 if(iset.eq.99)icoul=2 if(iset.gt.100)iset=iset-100 if(iset.le.3)go to 4 if(nwaves.le.0)acoup=-nwaves/100. c set pot real*8 for iset=25 if(iset.eq.25)a=l1(3) if(iset.eq.25)l1(3)=a write(6,901)ibag, iset,acoup 901 format(' iset gt 3 in vkbarp,set=3 ibag,iset,acoup=',2i3,f10.3) iset=3 go to 4 c schick-gibson complex*16 single channel potentials(kim or bea) c c kim 6 b0=0.180 lo(1)=(-6.45,-0.5928) b1(4)=0.180 l1(4)=(-3.855,-3.101) go to 8 c bea 7 b0=0.07643 l0(1)=(-77.89,-3.289) b1(5)=0.07643 l1(5)=(-54.94,-16.55) 8 l0(2)=0. l0(3)=0. 4 continue c special run,pure i=0 with 2x strength for b.s. if(i1chan.eq.0)l0(2)=l0(2)*rt2 if(i1chan.eq.0)l0(1)=2.*l0(1) if(i1chan.eq.0)l1(iset)=0. write(6,900)ibag, 1 iset,b0,acoup,i1chan,(l0(i),i=1,3),b1(iset),l1(iset) 900 format(' ibag,iset,b0=',2i3,2d15.5,'acoup,i1chanl=',f10.3,i3/ 1 ' l0/b0(1,2,3 )=',6d15.5/' b1,l1(iset)=',4d15.5/) nit=1 if(iset.eq.4.or.iset.eq.5)go to 9 b0=b0*hbarc do 5 i=1,3 l1(i)=hbarc*l1(i)*fourpi c reduce im part of potential by amount reduce coupling if(acoup.eq.1) go to 12 b=l1(i) a=-acoup*zi*l1(i) l1(i)=b+zi*a 12 l0(i)=l0(i)*fourpi*b0 5 b1(i)=(b1(i)*hbarc)**2 b0=b0**2 go to 11 c intialize and convert s-g pot into lpt convention 9 b0=(hbarc/b0)**2 b1(iset)=(hbarc/b1(iset))**2 l0haw=l0(1)*pi/2. l1haw=l1(iset)*pi/2. l0(1)=l0haw*(hbarc**2)*fourpi l1(iset)=l1haw*(hbarc**2)*fourpi 11 write(6,899) 899 format(' now in mev units,b=b*b,l=4pi*l--for lpt pot') write(6,900)ibag, 1 iset,b0,acoup,i1chan,(l0(i),i=1,3),b1(iset),l1(iset) c compute the l=0 optical potetial, lpt conventins 10 kp2=kp*kp 21 k2=k*k if(ibag.ne.1)rhl=g(b0,l0(1)) vkp(1,1)=0.5*(rhl + g(b1(iset),l1(iset))) bsq=5.569836 l=-0.1626465 if(icoul.eq.1)vkp(1,1)=vkp(1,1)+g(bsq,l) if(icoul.eq.2)vkp(1,1)=g(bsq,l) if(ibag.ne.1)vkp(2,2)= g(b0,l0(3)) if(ibag.eq.1)vkp(2,2)=fl(2) vkp(3,3)=vkp(1,1) if(ibag.ne.1)vkp(1,2)= acoup*g(b0,l0(2))/rt2 c the equivalent separable coulomb pot of rcb vkp(2,1)=vkp(1,2) vkp(3,2)=-vkp(1,2) vkp(2,3)=vkp(3,2) vkp(1,3)=acoup*0.5*(g(b1(iset),l1(iset)) - rhl) vkp(3,1)=vkp(1,3) c 800 format(2(2x,e15.8)) return end C=================================================================