program lpotp c *** Calculation of R-matrix and T-matrix by matrix inversion with c potential in momentum space for spin 1/2 x 1/2, at any energy. c-------modifications for PVM shown, other parts edited 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 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 approx 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 *** 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 *** 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 c not good for exact Coulomb 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 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 %%% c ************** 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 435 continue c-----calculate angular disrtbn --------------------------------------- ......lines left out................ 620 continue end