program lpotp ******************************************************************************* c *** Paez modification for proton scattering - revised by RHL c Further revisions by M. Sagen c Minor changes to output format by L.Berge 26 Nov 88 c Now will only use Nucleons as the projectiles. All c references to the Pion code have been removed. c *** Converted to double precision by M. Sagen (August, 1986) c *** Calculation of R-matrix and T-matrix by matrix inversion with c potential in momentum space for spin 0 and 1/2, at any energy. c *** Energies, momenta, and masses in MeV units throughout program c *** Input switches - Nifty(1) thru Nifty(20) follow c c nifty(1) = -2, (p,n) cex after p,n c = 5, protons c = 9, neutron c nifty(2)*5MeV = binding energy (-5. * nifty(2) used) for 3-body energy c nifty(3) = 0, default case--takes no action. c 1, single scatt (when n(6)used else) c 2, set t=r (nonunitarize) but with ms c 3, double scattering for 1/2*1/2 c nifty(4) = 0, Regular Triplet-Singlet Mixing c 1, Triplet Singlet Mixing turned off c nifty(5) = 0, use softh(k0) c 1, sin(theat) c 2, so(on shell),no angle tf c 5, e3b energy c 6, e3b energy with folded tmatrice c 7, e3b with AAY magic vector prescription for theta c 9, soth(ko), AAY c nifty(6) = 0, nospin flip and no calc of ds and ts (0 x 0) c 1, calculate R matrix up to and include ds c 2, calculate R matrix up to and include ts c 3, include spin flip terms (0 x 1/2) c 4, no scatt (u=0) as check of wf c 5, R matrix for only ss c 6, Ustrong = 0 c 7, Imaginary part of U = 0 c 8, spin 1/2 x 1/2 with 6 spin states c nifty(7) = 0, Set H, I, J, and K NN-partial waves to zero c = 1, Retain all NN-partial waves. c nifty(8) = 0, don't read in phases, use stored t c = 1, read in phases c nifty(9) = 0, Not used by the Proton code. c nifty(10) = 0, no Coulomb amplitude c 1, Coulomb amplitude and phase change included c 2, Coulomb amplitude, but not phase change c 3, Exact coulomb, spherical app. c 4, Exact coulomb matching, realistic shape c *** Nifty(11) may be used as temporary switch c nifty(11) = 0, Schroedinger theory c = 1, Relativistic (Dirac equation) theory c (12) = 0, for off-shell GRAZ structure c = 1, for off-shell determined via f-ratios c = 2, for PEST interaction c (13) = 0, not used for proton code c (14) = 0, not used for proton code c nifty(15) = 0, divide out proton ff from He ff c 1, dont divide out proton ff c nifty(16) = 0, modified aussian rho and rho**2 & for p* z=n c nucleus will also use t(spin)*rho(matter) c 1, Wood Saxon rho,..but Gaussian rho**2 c 2, Wood Saxon rho and rho**2 c 3, special form for He or Deuteron c 4, He3..fnsp=fnmat,fpsp=0 c 5, He3..fnsp=fnmat=fpmt, fpsp=0 c 6, use fmatter*tspin for 0*.5 (non-Gaussian nucleus) c 7, uspin=(tsp*fsp+tspn*fmt) for 0*.5 c 8, use Ray's form factor for C13 c 9, use form factor of 1 (constant) c nifty(17) = 0, KMT potential c 1, Watson potential, no A/(A-1) c nifty(18) = 0, dont turn off any 1/2*1/2 amps c 1, no E terms c 2, no A+B terms c 3, no E and no A+B terms c 4, no F term in pA amplitude c nifty(19) = 0, nifty(16) controls ff's c 1, Hadjimichael impulse approximation form factors (no MEC) c 2, full Hadjimichael (i.e. with exchange currents) c nifty(20) = 0, reduced output, NO plot file to tape10 c = 1, full output, NO plot file to tape10 c = 2, reduced output, write plot file to tape10 c = 3, full output, write plot file to tape10 c nsex =0-normal/ =1-1st time thru/ =2-2nd time thru/ =3 diff cex c *** plot up lab cross sections if nang was read in as a negative number c *** U(nk,1) = Uplus for spin flip, U(nk, 2)= Uminus c *** set = for nonflip (same code for tr,ti,trc,tic,tbr,tbi) c ---------------------------------------------- include '/usr/local/include/fpvm3.h' implicit real*8 (a-h, o-z) integer mytid, consti, tids, bufid real*8 k0 , ko , ko2 , kk real*8 Mnuc , Mnuc2 , MN , Mp , Mp2 complex*16 zf c *** N.B. need to match dimensions and # of entries in data statement c *** t's enlarged 3/84 for cex (m.s.) c *** dimension of u/f depends on ngp c %%% dimension of ur,ui (ngp+1,ngp+1,6) - LB 10 Dec 88 c %%% dimension of fr,fi (2*ngp+2,2*ngp+2) - LB 10 Dec 88 c %%% dimension of den (2*ngp+2) - LB 10 Dec 88 dimension tr(100,12),ti(100,12) dimension ttr(12), tti(12), tids(97) dimension kk(97) dimension ur(43,43,8) , ui(43,43,8) dimension fff(43,88,7,43) c dimension zf(194,194) dimension zf(86,86) dimension den(194) , gp(96) , wt(96) , sigbrn(361) dimension theta(361) , dflip(361) , dnof(361) dimension x1(361), x2(361), x11(361), x2p(361), x3(361), $ x5(361), x4(361) dimension tva(361), tvb(361), tvc(361), tvd(361), $ tve(361), tvf(361), tma(361), tmb(361), $ tmc(361), tmd(361), tme(361), tmf(361) dimension dsig(361), chart(61,61), polar(361) dimension scoul(100) dimension nifty(20), constr(10), consti(10) c common /nlspfl/ ur, ui common /parallel/ fff common /foptp/ zf common /kinemt/ El , Ecm , k0, pl , s , Mp2 , Mnuc , Mnuc2 common /params/ hbarc, pi, Mp, MN, nz, na, nes, nwaves common /switch/ nifty common /sizes/ achp, acmp, wsp, achn, acmn, wsn common /ranges/ rcoul, rcut common /inputs/ Tlab , b , ymin1 , ymin2 , xnang, $ kode , lxmax , ngp , nr common /Rcomn/ rr , ri , rrb , rib common /spins/ nspina , nspinb , nspinc , nspind , nspine common /Tcoul/ Trc , Tic , Tr3c , Ti3c , Tr5c , Ti5c, xgam common /Tcomn/ tr , ti ******************************************************************************** hbarc = 197.3289e+0 pi = 3.141593e+0 icex = 1 call GetIns() c nplot = 0 if (xnang .lt. 0.) nplot = 1 xnang = abs(xnang) amass = na zcm = acmp zch = achp zws = wsp nsex = 0 nifty1 = nifty(1) c *** charge exchange, do both p-He3 and n-He3, then subtract if ( nifty1 .eq. -2 ) then nsex = 1 nifty(1) = 5 endif 60 continue c *** cex(proton) mod (m.s.) if (nifty(1) .eq. 5) Mp = 938.2796 if (nifty(1) .eq. 9) Mp = 939.5731 c *** Mnuc here refers to nucleus mass, MN is nucleon MN = (nz * 938.2796 + (na - nz) * 939.5731)/na Mnuc = MN * amass c write(6,780) amass, Mp, MN Mnuc2 = Mnuc * Mnuc Mp2 = Mp * Mp 80 continue if (nr .le. 0) write(6,1000) c *** Now k0 for non-relativistic case is the same as for relativistic, c *** only the wave equation differs El = Tlab + Mp Ecm = SQRT( Mnuc2 + Mp2 + 2. * Mnuc * El ) k0 = SQRT( (El * El - Mp2) * Mnuc2/Ecm/Ecm ) pl = SQRT( El * El - Mp2 ) c write(6,960) Tlab, El, Ecm, k0, pl c write(6,980) if (nr .gt. 4) write(6,990) c write(6,940) c write(6,950) Tlab, k0 n1 = ngp + 1 n2 = n1 + n1 a = k0 c *** special check if b=0, use different distrb of points if (b .eq. 0.) a = 200. call GAUSS2 (ngp, kode, 0.0d0, a, gp, wt) do 90 i1 = 1 , ngp kk(i1) = gp(i1) 90 continue kk(n1) = k0 ko = k0 ko2 = ko * ko emnos = sqrt(ko2 + Mnuc2) epios = sqrt(ko2 + Mp2) if (nr.eq.0) then emnos = Mnuc epios = Mp endif s = (emnos + epios)**2 c *** Now set up the denominator for either relativistic or non-relativistic c *** case call DENOM( k0, Mnuc, emnos, epios, psfac, kk, wt, den) c *** c *** set lBorn and lmax if needed c *** lclass = ko * 1.33 * (na**(.333333))/hbarc ldmax = lxmax if (lxmax .eq. 0) ldmax = 2.5 * lclass + 2.5 lBorn = 2 * lclass if (lBorn .gt. ldmax) lBorn = ldmax xgam = 0. lBorn = ldmax icex = 2 write(6,730) ldmax, lBorn, lxmax c write(6,4020) c should set lborn to lmax so born approx not made,rhl 11/88 write(6,730) ldmax, lBorn, lxmax c *** Call to PRAMPS at present necessary for He-3 in order to get proper c *** initialization though it is not clear why this should be so. c *** include coulomb, setup - EVEN FOR NEUTRON *** zp = 1.0 c set charge to zero for neutrons or for no coulomb case if ( nifty(1) .eq. 9 .or. nifty(10).eq.0 ) zp = 0.0 xgam = zp * nz/137.036/(ko * (emnos + epios)/emnos/epios) alph = (nz - 2)/3. alph = alph/2./(2. + 3. * alph) if (na .lt. 4) alph = 0. c *** set up sigl *** rhoc = ko * rcut/hbarc lcmax = ldmax c if changed from eq to le 9/6/90 if (lcmax .le. 1) lcmax = 2 eta = xgam call SIGCL (eta, scoul, lcmax, 100) sig0 = scoul(1) c *** c *** do loop over l for this e and state c *** c *** *** Calculate Optical potential with subroutine call *** *** *** c *** *** The Phase I parallelization is in this subroutine *** *** *** call optp1 (kk, n1, ldum, ldmax, lBorn, aovera) c **** \/\/\/\/ Begin PARALLEL Phase II\/\/\/\/ ************ c ******** Ready to send slaves to solve LS eq. and Coul. match. *** newdim = n1 * n1 * 8 constr(1) = psfac consti(2) = n2 constr(2) = aovera consti(3) = n1 constr(3) = xgam consti(4) = ldmax constr(4) = k0 constr(5) = rcut constr(6) = hbarc consti(5) = newdim call pvmfmytid( mytid ) call pvmfcatchout(1) call pvmfspawn('Conquistador',PVMDEFAULT,'*',ldmax,tids,numt) do 420 ldum = 1,ldmax ld = ldum l = ldum - 1 consti(1) = ldum c*************************************************************************** c declare dimensions here to pass along to optp c*************************************************************************** c %%% changed to allow 48 integration points - LB 11 Dec 88 %%% ************** tvm 10/26/92 add aovera to call optp ****************** call optp2 (n1, ldum, ldmax, lBorn, aovera) call pvmfinitsend(PVMRAW,bufid) call pvmfpack( INTEGER4, consti(1), 5, 1, info) call pvmfpack( REAL8, constr(1), 6, 1, info) call pvmfpack( INTEGER4, nifty(1), 20, 1, info) call pvmfpack( REAL8, den(1), 194, 1, info) call pvmfpack( REAL8, scoul(1), 100, 1, info) call pvmfpack( REAL8, ur(1,1,1), newdim, 1, info) call pvmfpack( REAL8, ui(1,1,1), newdim, 1, info) call pvmfsend( tids(ldum), 1, info) 420 continue c ****** ^^End send loop, \/\/ Begin recall loop do 425 ld = 1,ldmax call pvmfrecv(tids(ld), ld, bufid) call pvmfunpack( REAL8, ttr(1), 8, 1, info) call pvmfunpack( REAL8, tti(1), 8, 1, info) do 422 nsp = 1, 8 tr(ld,nsp) = ttr(nsp) ti(ld,nsp) = tti(nsp) 422 continue testy2 = tr(ld,5) call coulomb2(scoul, sigl, aovera, ld, njump ) 425 continue call pvmfexit(info) c **** ^^^^ End Parallel, LS eq. solved. T-matrix received. *** c *** Coulomb matching executed and phase adjusted. ********* c *** End of do loop over partial waves(after remove coulomb) *** 430 continue do 434 ll = 1,ldmax c ************* tvm 9/2/92 TS mixing Why fight? 2 - 7 interchange in nspin *************** convention from this point on. ********************** tmdum = tr(ll,2) tvdum = ti(ll,2) tr(ll,2) = tr(ll, 7) ti(ll,2) = ti(ll, 7) tr(ll,7) = tmdum ti(ll,7) = tvdum ********* Fix iteration convention in line below, -1 tvm 7/28/92 **** c write(6,1031) ll-1, (tr(ll,np2), np2=1,8), (ti(ll,np2), np2=1,8) 434 continue 435 continue if (nifty(6) .ne. 8) write(6,1090) if (nifty(6) .eq. 3) write (6,1110) c write(7, 7213) c7213 format(1h1,'th_c.m.', 4x, 'Aoolm', 7x, 'Aooml', c $ 7x, 'Aoomm', 7x, 'Aooll') c-----calculate angular disrtbn --------------------------------------- call xsects( theta, dsig, sigbrn, dnof, dflip, polar, $ scoul, alph, xgam, sig0, ldmax, noangs, $ nspin, nsex ) if (nplot .ne. 1) go to 581 c *** plot up cross sections on semilog 5 cycle paper 561 nspin = 0 if (nifty(6) .eq. 3 .OR. nifty(6) .eq. 8) nspin = 1 c *** variables for plotting titles chart(1,1) = Tlab chart(2,1) = achp chart(3,1) = acmp chart(4,1) = wsp chart(5,1) = achn chart(6,1) = acmn chart(7,1) = wsn chart(8,1) = nz chart(9,1) = amass do 570 ii = 1,20 j = ii + 9 chart(j,1) = nifty(ii) 570 continue if (nifty(6) .ne. 8) then call logplt (theta, dsig, dnof, noangs, nspin, chart, ymin1) endif if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) then call logplt (theta, dnof, dflip, noangs, 1, chart, ymin2) endif if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) go to 581 c dsig/domega c call logplt (theta, dsig, dnof, noangs, 0, chart, ymin1) c call logplt (theta, dsig, dnof, noangs, 0, chart, ymin2) if(nifty(20) .lt.2) go to 576 c read in plot data from tape10 rewind 10 do 575 i = 1, noangs if (na.eq.13) then read(10,*) theta(i), x1(i), x2(i), x11(i), x2p(i), x3(i), $ x4(i), x5(i) endif 575 continue 576 continue c 2,2' 11 if(na.eq.13) go to 577 c call logplt (theta, polar, dnof, noangs, -1, chart, ymin2) c call logplt (theta, dflip, dnof, noangs, -1, chart, ymin2) c call logplt (theta, dnof, dnof, noangs, -1, chart, ymin2) 577 continue c rhl add for c13 if(na.ne.13) go to 579 do 578 i=1,noangs 578 theta(i)= 3.*theta(i) call logplt (theta, dsig, dnof, noangs, 0, chart, ymin1) if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) then call logplt (theta, dnof, dflip, noangs, 0, chart, ymin2) endif if (nifty(6) .ne. 3 .AND. nifty(6) .ne. 8) go to 581 call logplt (theta, dsig, dflip, noangs, 1, chart, ymin2) call logplt (theta, polar, dflip, noangs, -1, chart, ymin2) if (nifty(20).lt. 2) go to 579 c plot up plotdata stored on tape10 c call logplt (theta, x2, dflip, noangs, -1, chart, ymin2) call logplt (theta, x11, dflip, noangs, -1, chart, ymin2) call logplt (theta, x2p, dflip, noangs, -1, chart, ymin2) call logplt (theta, x3, dflip, noangs, -1, chart, ymin2) call logplt (theta, x4, dflip, noangs, -1, chart, ymin2) call logplt (theta, x5, dflip, noangs, -1, chart, ymin2) 579 continue 581 if (nsex .eq. 0) go to 620 if (nsex .eq. 3) go to 620 if (nsex .eq. 2) go to 610 c *** nsex = 1 pass done, setup 2nd pass c *** proton done, now do neutron nsex = 2 nifty(1) = 9 c______________________________________________________________________________ c Now reinitialize t's by calling TNCM for charge exchange c______________________________________________________________________________ energy=-200.e+00 call tncm (energy) c *** go to 60 c *** nsex = 2, two passes for cex done,now do last 610 nsex = 3 ymin1 = ymin1 - 1 go to 430 620 continue write(6,841) c stop c _______________ formats c *** c227 format(i3, 8e15.6 / 200(3x, 8e15.6 / )) 700 format (15h p, n, then cex) 730 format (1h0,'lmax(calculated) = ',i4,5x,'lBorn = ',i4,5x, $ 'lmax(read) = ',i4) 740 format (1h1,5i4) 780 format (1h0,' A=',f10.4,' Mp=',f10.4,' m(nucleon)=',f10.4//7x, $ ' Tlab', 9x, 'Eplab ', 8x, 'Ecm(sqrts)', 4x, $ 'kcm', 11x, 'plab') c820 format (47h R matrix calculated with ms series up to order,i3) 841 format(12h eof in main ,i3) c928 format(1h0,33h re/im f matrix before inversion ) c930 format (1h ,26h normalized r(k,k..k)= ,24x,2e14.7) 940 format (1h0,10x,8h Energy,14x,15h Momentum ,/) 950 format (1h ,7x,f14.4,10x,f14.4) 960 format (1h ,5f14.4) 980 format (1h0,31h relativistic calculation) 990 format (1h ,10x,33happroximate klein gordon equation) 1000 format (1h0,47h nonrelativstc calc--but with same k0 as reltv,40x 1,32h i.e. only the wave eqtndiffers) c1010 format (1h+,85x,4heta=,f10.4,8h delta=,f10.4) c1020 format (41h0********** eta greater than 1 **********,//) c1030 format (1h ,15h t(k,k,k)= ,35x,4e14.7) 1031 format (1h ,16h t(k,k,k),for l=,i3,/,10(8e10.3,/,8e10.3,/)) c1040 format (38h now with exact coulomb after matching) c1050 format (20x,40hcoulomb phase modified nuclear amplitude/16h t(k, c 1k.k) = ,4e14.7) 1070 format (1h0,30h----- angular momentum = ,i4,5h-----) c1080 format (1h ,8x,23h born approximation ) 1090 format (17h1th-c.m. cos(cm),3x,10hdsig/dw-cm,2x,11hds/dw(b.a.), 13x,1ht,4x,6hth-lab,2x,8hcos(lab),2x,5hds/dw,4x,11hds/dw(b.a.),2x, 210href(th)-cm,1x,3himf,4x,6hsignuc,4x,7hsigcoul//) c1091 format(1h1,' th-c.m. cos(cm),3x,4hq2fm,7x,10hdsig/dw-cm , c 1 1x,10hds/dw(nuc) ,1x,11hds/dw(b.a.) ,2x,5hpolar ,4x, c 2 7hnoflip ,4x,4hflip ,5x,4hpnoo ,4x,5hkolmo ,4x, c 35hclmoo ,2x,7hsigcoul /1h ,18x, c 420hza ,zb ,zc ,zd ,ze /1h ,18x, c 520hzab,zbb,zcb,zdb,zeb //) 1091 format(1h1,'th-c.m. q ',4x,'ds/dw',6x,'P=Aoono',6x, $ 'Dnono',7x,'Dloso',7x,'Dsoso',7x,'Dlolo',7x,'Aoonn',7x,'Aooon') c1093 format(1h ,f7.0,2f8.3,3e11.3,f8.3,2e11.3,3f9.3,e11.3) c1100 format(1h ,f7.0,f8.3,2e11.3,2f8.1,f10.3,2e11.3,2f9.3,e12.4 c 1,e10.2) 1110 format (1h+,21x,6hnoflip,7x,8hnoflip-b,5x,5hpolar,9x,4hflip,9x, 18hflip-b,r,//) c1120 format (1h ,18x,10e11.3) c1498 format(' ???main ur(lid,lid,1)= ',e10.4,' ui(same)= 'e10.4) c1921 format(1h ,7h sigot=,e13.5,7h sigit=,e13.5, c 27h sig2t=,e13.5,2he=,f8.2,5hxgam=,e13.5) 1979 format(1h ,17h ur ui(k0,k0,k0)= ,2e13.5) 2211 format(1h ,' .......... nspin= ',i4,' ...... ') c3335 format(1h ,19h tr ti (nspin=3-6)=,8e13.5) c4000 format(' ', f7.2, 3e11.3) c4010 format(1h ,f7.2,2e11.3) 4020 format(1h2) end