Program TestPiNdel c Reads pi-nucleon phase shifts from file delal c Produces disk file: TestPiN.dat with two columns: c COM energy and eta or d for selected nwave c 1<=nwave <= 14 , eta :real part and d: imaginary part c nwave code c 1-----s31 2-----p31 3-----p33 4-----s11 c 5-----p11 6-----p13 7-----d33 8-----d35 c 9-----f35 10----f37 11----d13 12----d15 c 13----f15 14----f17 c Implicit Real*8 (a-h,o-z) Dimension ecm(100), eta(14,100), d(14,100), retf(14,100) Dimension aimtf(14,100) pi = 3.141592 nwaves = 14 nes=96 !25, RSL phaseshifts, 71 A-L nei=26 Open(6, File='TestPiNdel.dat', Status='Unknown') Write(6,*) nes, nwaves Open(1,file='delal',status='unknown') Read (1,170) (ecm(ne),(eta(nwave,ne),d(nwave,ne),nwave=1,nwaves), 1 ne=nei,nes) c Get first 25 phase shifts from RSL formula come=ecm(26) Call rosalan(come,eta,d,ecm) c d tabulated in degrees, now convert to radians Do nwave=1,14 Do ne=1,nes d(nwave,ne)=d(nwave,ne)*pi/180.0 End Do End Do c Call Argand Diagram calc Call argand (ecm, eta,d,retf,aimtf,nes) Call crossect (ecm,retf,aimtf,nes) 170 Format (13f6.1/8f6.1/8f6.1) Close(1) Close(6) Stop End Subroutine argand(ecm, eta,d,retf,aimtf,nes) c c Output is file "TestPiNdel.dat" c with Ecm, Re f. Im F c f=-0.5i(eta*exp(2id -1)) c select below for desired nwave Implicit Real*8 (a-h,o-z) Real*8 mpi Dimension eta(14,100),d(14,100),retf(14,100) Dimension aimtf(14,100), ecm(100) pi=3.141592 mpi=139.578 c Open(6,file='TestPiNdel.dat', status='unknown') Do nwave=1,14 Do ne=1,nes twodel=2.0*d(nwave,ne) retf(nwave,ne)=0.5*eta(nwave,ne)*sin(twodel) aimtf(nwave,ne)=0.5*(1.0-eta(nwave,ne)*cos(twodel)) c retf(nwave,ne)=0.5*sin(twodel) ! eta=1 c aimtf(nwave,ne)=0.5*(1.0-cos(twodel)) !eta=1 End Do End Do nwave=3 Do ne=1,nes-1 Write (6,100) retf(nwave,ne), aimtf(nwave,ne), ecm(ne) !argand dia End Do 100 Format(1h , e15.5, e15.4, f8.2) c Close(6) Return End Subroutine kinemat(come,kappa,pp) c input: center of mass energy c output: kappa, pp :com pion momentum and lab pion momentum Implicit Real*8 (a-h,o-z) Real*8 mn,mpi,kappa mpi=139.578 !pion mass mn=938.272 !proton mass epi=(come**2-mpi**2-mn**2)/(2*mn) pp=sqrt(epi**2-mpi**2) s = come**2 kappa = Sqrt((s-(mn+mpi)*(mn+mpi))*(s-(mn-mpi)*(mn-mpi))/4./s) Return End Subroutine crossect(ecm,retf,aimtf,nener) c input: ecm aimtf and c nener: number of com energies and aimmtf's needed c output file "TestPiNdel.dat" c that contains either the total cross sect. versus pi-lab momemt. c or the contribution of any partial wave to the cross section c or tpilab and reft for a wave c or tplilab and aimtf for a wase c depending of your choice of the i value c write(6,100) statement Implicit Real*8 (a-h,o-z) Real*8 kappa, mpi Dimension ecm(100), retf(14,100), aimtf(14,100) pi=3.141592 hbarc=197.3289 mpi=139.578 ! pion mass c fact : to give total cross section in mb c Open(6,file='TestPiNdel.dat', status='unknown') c call kinemat gives com enery (come) and c gets kappa: com pion momentum c pp: lab. pion momentum c tpilab: kinetic energy of the pion in lab system c i=0 tpilab and total cross section c i=1 tpilab and cntrbtn of any partial wave to crossection c i=2 tpilab and retf for a wave c i=3 tpilab and imtf for a wave i=0 Do ne=1,nener come = ecm(ne) Call kinemat(come,kappa,pp) tpilab=sqrt(pp*pp+mpi*mpi)-mpi fact=10.0*4.0*pi*(hbarc/kappa)**2 swve=aimtf(1,ne) pwve=2.0*aimtf(3,ne)+aimtf(2,ne) dwve=3.0*aimtf(8,ne)+2.0*aimtf(7,ne) fwve=4.0*aimtf(10,ne)+3.0*aimtf(9,ne) sigtot=fact*(swve+pwve+dwve+fwve) c first write for sigma total c second write for contribution of a wave to sigma total c uncomment one and comment the other according to your c choice if(i.eq.0)Write(6,100) tpilab,sigtot c if(i.eq.1)Write(6,100)tpilab,fact*2.0*aimtf(3,ne) c if(i.eq.2)Write(6,100)tpilab,retf(3,ne) c if(i.eq.3)Write(6,100)tpilab,aimtf(3,ne) End Do c Close(6) 100 Format(1h ,e15.5,e15.4) Return End Subroutine rosalan(come,eta,d,ecm) c Rowe, Salomon, Landau phase shifts c input: come = 1st energy, then does 25 below this c output: eta(14,100), d(14,100) and ecm(1-25) c taken from routine tpicm from lpott c the output is a file "TestPiNdel.dat" c select wave (1 to 14), 3 in this case Implicit Real*8 (a-h,o-z) Real*8 kappa,kaph,dumkap Dimension eta(14,100),d(14,100),ecm(100) c Open(6,file='TestPiNdel.dat', status='unknown') pi=3.141592 Call kinemat(come,kappa,pp) If (kappa.le.5.) kappa = 6. kaph = (kappa-5.)/20. kappa = 0. c Do nkap=1,25 ! as in tpicm c comment last line and uncomment next to plot phashifts Do nkap=1,51 !for com pi momentum, up to 51 for P33graph kappa = kappa+kaph If (nkap.le.6) kappa = (nkap*1.5) Do ni=1,14 !for partial waves Call tpirsl (kappa,ni,ddum,w,tr,ti) c w is given by last call is com energy eta(ni,nkap) = 1. d(ni,nkap) = ddum*180./pi c tangent jumps at 90 prom + to -, correct to sum phases c for ni=3 if((ni.eq.3).and.(d(3,nkap).lt.0.0)) d(3,nkap)= 1 180.0+d(3,nkap) End Do ecm(nkap) = w c input w in kinemat and use pp: lab pion kin. ener. Call kinemat(w,dumkap,pp) c To plot S31 P31,P33, P11 phaseshits in degrees c Write(6,*)pp,d(3,nkap) !here End Do c Close(6) Return End Subroutine tpirsl (kappa,ni,ddum,w,tr,ti) c *** see RH Landau, Program lpott, 1981 c calc of pin phases and amplitudes with rowe,salomon,Landau c formula, or Landau s for higher l Implicit Real*8 (a-h,i,k,m,o-z) Dimension a(14), b(14), c(14), xx(14), w0(14), q0(14), gam(14) c use rowe,salomon,Landau for low e s,p,d13,d15;Landau for others c data mpi,mn,pi,kold/139.578e+0,938.903,3.141593e+0,0./, 1 a/-11.2,-2.91,11.4,16.8,-5.7,-1.31,11.2e-13,-3.15e-12, 2 8.31e-18,5.46e-17,0.109,0.112,-38.26e-18,7.43e-18/, 3 b/-30.7,3.45,-15.4,-35.4,25.4,1.22,-0.369e-14,0.118e-13, 4 -0.363e-19,-0.236e-18,-0.031,-0.270,0.205e-18,-0.222e-19/, 5 c/21.,-1.5,7.2,27.,-29.,-0.4,0.,0.,0.,0.,0.003,0.19,0.,0./, 6 xx/.31,.22,.99,.44,.61,.23,0.,0.,0.,0.,.54,.43,0.,0./, 7 w0/1655.,1850.,1233.,1550.,1435.,1815.,0.,0.,0.,0., 8 1525.,1670.,0.,0./, 9 q0/550.,678.,228.,477.,393.,656.,0.,0.,0.,0.,459.,560.,0.,0./, a gam/170.,200.,116.,105.,230.,255.,0.,0.,0.,0.,125.,155.,0.,0./ If (abs(kappa-kold).lt.1.e-03) GoTo 10 kold = kappa rhl = kappa*kappa If (kappa.eq.0.) kappa = 1.e-10 epikap = Sqrt(mpi*mpi+rhl) enkap = Sqrt(mn*mn+rhl) w = epikap+enkap s = w*w factor = (-1./kappa)*w/4./pi/pi/epikap/enkap q = kappa/mpi q2 = q*q q4 = q2*q2 c determine l value for formula 10 If ((ni.eq.1).or.(ni.eq.4)) l = 0 If ((ni.eq.2).or.(ni.eq.3).or.(ni.eq.5).or.(ni.eq.6)) l = 1 If ((ni.eq.7).or.(ni.eq.8).or.(ni.eq.11).or.(ni.eq.12)) l = 2 If ((ni.eq.9).or.(ni.eq.10).or.(ni.eq.13).or.(ni.eq.14)) l = 3 jj = 2*l+1 If ((ni.eq.7).or.(ni.eq.8).or.(ni.eq.9).or.(ni.eq.10).or.(ni.eq.13 1).or.(ni.eq.14)) GoTo 20 rhl = q0(ni)/mpi del = (q**jj)*(a(ni)*1.e-2+b(ni)*q2*1.e-3+c(ni)*q4*1.e-4+xx(ni)* 1gam(ni)*w0(ni)/(rhl**jj)/(w0(ni)*w0(ni)-s)) ddum = atan(del) GoTo 30 c no salomon fit, use old rhl fit 20 ddum = (a(ni)+b(ni)*kappa)*(kappa**jj)*pi/180. 30 tr = 0.5*factor*sin(2.*ddum) ti = 0.5*factor*(1.-cos(2.*ddum)) Return End