subroutine TNCM(kappa) ******** This subroutine inputs the NN phase shifts as input data **** *********************************************************************** c *** Subroutine TNCM calcs averaged on-shell pN c.m. amplitudes c *** c *** kappa(or vector kap) = p-nucleon c.m. momentum c *** c *** Initialize if phase shifts or nucleus (Z, A) changes implicit real*8 (a-h, i, k, m, o-z) c ********* tvm 2/24/92 change all dimensions of 100 to 160 ******** dimension nifty(20), Ecm(160), Tlab(160), kap(160), plab(160) dimension d(40,160), eta(40,160), eps(8,160) dimension t(40), tx(40), rewave(40,160), imwave(40,160) dimension bn(40), bnuc(40), bnf(40), bnucf(40) CCC increase dimensionality to include other off-diagonal CCC ra(7)--> ra(14); r(7)-->r(14) : pjd 4/4/93 c dimension ta(7), tax(7), ra(7), rax(7) dimension ta(7), tax(7), ra(14), rax(14) CCC dimension tawave(7,160), taxwve(7,160), retij(40), imtij(40) common /tjmats/ t, tx, ra, rax common /bandt/ bnuc, bn, bnucf, bnf, retij, imtij common /params/ hbarc, pi, Mp, Mn, nZ, nA, nes, nwaves common /switch/ nifty common /tape/ rewave, imwave, tawave, taxwve, kap equivalence (rewave(1,1), d(1,1)), (imwave(1,1), eta(1,1)) c *** Flip is more complicated, see Alexander's notes c *** Nwave code (order for read in changed from Alex) c *** c *** I = 1 c *** c *** 1) 1S0 2) 1D2 3) 1G4 4) 3P0 5) 3P1 6) 3P2 c *** 7) 3F2 8) 3F3 9) 3F4 10) 3H4 11) 3H5 12) 3H6 c *** c *** c *** I = 0 c *** c *** 13) 1S0 14) 1P1 15) 1F3 16) 1H5 17) 3S1 18) 3D1 c *** 19) 3D2 20) 3D3 21) 3G3 22) 3G4 23) 3G5 c *** c *** c *** I = 1 c *** c *** 24) 1I6 25) 3J6 26) 3J7 c *** c *** c *** I = 0 c *** c *** 27) 3I5 28) 3I6 29) 3I7 30) 1J7 31) 3K7 32) 3K8 c *** c *** After read in phases, the mixing parameters follow c *** epsilon (eps) 1,2,3,..... c *** N.B. This is modified form from Alexander c>>>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Mln = Mp nucln = 1 ndata = 0 IF (kappa .le. -100.0) ndata = 1 kappa = abs(kappa) IF (ndata .eq. 0) go to 2001 IF (nifty(8) .ne. 1) go to 2000 c *** Nifty(8) = 1, then read in phase shifts. c *** Tmatrix (unfolded) for this nucleus. c *** Amplitudes from input phase shifts and energy. c *** c *** Tlab in MeV, delta (d) read in in degrees, f calculated in MeV-1, c *** eta = exp(-2 * Im(delta)), if read in as 0 then set to 1. c *** open(1, file='tape1', status='unknown') nei = 1 rewind 1 read(1,919) ne, nwaves c *** Proton or neutron do 8011 ne = nei, nes read(1,9001) Tlab(ne), $ (eta(nwave,ne), d(nwave,ne), nwave=1,nwaves) read(1,9001) (eps(j,ne), j=1,8) c *** Write out phases c if( nifty(20) .eq. 1 .OR. nifty(20) .eq. 3 ) then c write(6,9003) Tlab(ne), c $ (eta(nwave,ne), d(nwave,ne), nwave=1,nwaves) c write(6,9003) Tlab(ne), (eps(j,ne), j=1,8) c else if( ne .eq. 1 .OR. ne .eq. nes ) then c write(6,9003) Tlab(ne), c $ (eta(nwave,ne), d(nwave,ne), nwave=1,nwaves) c write(6,9003) Tlab(ne), (eps(j,ne), j=1,8) c endif c *** Lab KE read in, convert to Ecm Ecm(ne) = sqrt(2. * Mln * (Tlab(ne) + 2. * Mln)) 8011 continue rewind 1 do 548 nu = 1,7 ta(nu) = 0. tax(nu) = 0. 548 continue c------- Calculate kinematics do 1000 ne = 1,nes s = Ecm(ne) * Ecm(ne) kap(ne) = SQRT((s - (Mln + Mp) * (Mln + Mp)) $ * (s - (Mln - Mp) * (Mln - Mp))/4./s) plab(ne)= kap(ne) * SQRT(s)/Mln k = kap(ne) Epofk = SQRT( Mp*Mp + k*k ) ENofk = SQRT( Mln*Mln + k*k ) ftot = -(Epofk + ENofk)/4./pi/pi/Epofk/ENofk factor = ftot/k c----------- Compute alphas and convert deltas c *** c *** See LPT p. 311 c *** See also Stapp, Ypsilantis, Metropolis Phys Rev 105 #1,305/57 c *** c *** Code for ta, tax p. 305 Stapp et al.. c *** ta(1) = Re(alpha1/2i), ta(2) = Re(alpha2/2i) etc. c *** tax for imaginary parts of the same alphas. c *** A factor of 1/2, not in the formula, included to agree c *** with rewave imwave in TNCM then one factor of 1/2 not included c *** in the formulas for t11, t1m1, t00, and tss. ds = (d(18, ne) + d(17, ne)) * pi/180. fp = (d( 7, ne) + d( 6, ne)) * pi/180. gd = (d(21, ne) + d(20, ne)) * pi/180. hf = (d(10, ne) + d( 9, ne)) * pi/180. ig = (d(27, ne) + d(23, ne)) * pi/180. ajh = (d(25, ne) + d(12, ne)) * pi/180. aki = (d(31, ne) + d(29, ne)) * pi/180. do 544 nwave = 1,nwaves IF (eta(nwave, ne) .le. 0.) eta(nwave, ne) = 1. 544 continue eta21 = SQRT( eta(18, ne) ) eta01 = SQRT( eta(17, ne) ) eta32 = SQRT( eta( 7, ne) ) eta12 = SQRT( eta( 6, ne) ) eta43 = SQRT( eta(21, ne) ) eta23 = SQRT( eta(20, ne) ) eta54 = SQRT( eta(10, ne) ) eta34 = SQRT( eta( 9, ne) ) eta65 = SQRT( eta(27, ne) ) eta45 = SQRT( eta(23, ne) ) eta76 = SQRT( eta(25, ne) ) eta56 = SQRT( eta(12, ne) ) eta87 = SQRT( eta(31, ne) ) eta67 = SQRT( eta(29, ne) ) et1s = eta21 * eta01 et2s = eta32 * eta12 et3s = eta43 * eta23 et4s = eta54 * eta34 et5s = eta65 * eta45 et6s = eta76 * eta56 et7s = eta87 * eta67 c *** To make s = sin(2 * epsilon * pi/180) s1 = SIN( eps(1, ne) * pi/90. ) s2 = SIN( eps(2, ne) * pi/90. ) s3 = SIN( eps(3, ne) * pi/90. ) s4 = SIN( eps(4, ne) * pi/90. ) s5 = SIN( eps(5, ne) * pi/90. ) s6 = SIN( eps(6, ne) * pi/90. ) s7 = SIN( eps(7, ne) * pi/90. ) x1 = s1 * et1s x2 = s2 * et2s x3 = s3 * et3s x4 = s4 * et4s x5 = s5 * et5s x6 = s6 * et6s x7 = s7 * et7s tax(1) = x1 * SIN(ds) tax(2) = x2 * SIN(fp) tax(3) = x3 * SIN(gd) tax(4) = x4 * SIN(hf) tax(5) = x5 * SIN(ig) tax(6) = x6 * SIN(ajh) tax(7) = x7 * SIN(aki) ta(1) = x1 * COS(ds) ta(2) = x2 * COS(fp) ta(3) = x3 * COS(gd) ta(4) = x4 * COS(hf) ta(5) = x5 * COS(ig) ta(6) = x6 * COS(ajh) ta(7) = x7 * COS(aki) c *** To make compatible with LPT convention for rewave, imwave do 545 ni = 1,7 tawave(ni, ne) = 0. taxwve(ni, ne) = 0. tawave(ni, ne) = 0.5 * factor * ta(ni) taxwve(ni, ne) = 0.5 * factor * tax(ni) 545 continue do 500 nwave = 1,nwaves etadum = eta(nwave, ne) ddum = d(nwave, ne) * pi/180. IF (etadum .le. 0.) etadum = 1. rewave(nwave, ne) = 0. imwave(nwave, ne) = 0. IF (nwave .gt. nwaves) go to 500 c *** Nucleon case, mix partial waves, at least for MAW convention IF ((nwave .eq. 6) .OR. (nwave .eq. 7)) then etadum = etadum * COS( 2.0 * eps(2, ne) * pi/180. ) endif IF ((nwave .eq. 9) .OR. (nwave .eq. 10)) then etadum = etadum * COS( 2.0 * eps(4, ne) * pi/180. ) endif IF ((nwave .eq. 17) .OR. (nwave .eq. 18)) then etadum = etadum * COS( 2.0 * eps(1, ne) * pi/180. ) endif IF ((nwave .eq. 20) .OR. (nwave .eq. 21)) then etadum = etadum * COS( 2.0 * eps(3, ne) * pi/180. ) endif IF ((nwave .eq. 23) .OR. (nwave .eq. 27)) then etadum = etadum * COS( 2.0 * eps(5, ne) * pi/180. ) endif IF ((nwave .eq. 12) .OR. (nwave .eq. 25)) then etadum = etadum * COS( 2.0 * eps(6, ne) * pi/180. ) endif IF ((nwave .eq. 29) .OR. (nwave .eq. 31)) then etadum = etadum * COS( 2.0 * eps(7, ne) * pi/180. ) endif rewave(nwave, ne) = 0.5 * etadum * SIN(2. * ddum) * factor imwave(nwave, ne) = 0.5 * (1. - etadum * COS(2. * ddum)) $ * factor 500 continue 1000 continue ndata = 0 c-----Amplitudes all stored, interpolate on kappa, using 2 point lagrange c *** Amplitudes all stored on tape 3, read back 2000 continue ndata = 0 c *** Change do loop so UBC wont mess up the read 2001 continue x = kappa **************************************************************** IF (kappa .lt. kap(1)) x = kap(1) c *** Don't want scattering length for nucleon projectile. c *** Here t and tx refer to interpolated values of rewave and imwave c *** tvm 2/24/92 replace all 100 's with 160's in next 4 lines **** 2100 call LAGRNG (x, kap, t, rewave, nes, nwaves, 2, 160, 40) call LAGRNG (x, kap, tx, imwave, nes, nwaves, 2, 160, 40) call LAGRNG (x, kap, ra, tawave, nes, 7, 2, 160, 7) call LAGRNG (x, kap, rax, taxwve, nes, 7, 2, 160, 7) t(13) = 0. tx(13) = 0. CCC For on-shell, the off-diagonal components are equal, therefore CCC set ra(i) = ra(i+7) etc. (pjd 4/4/93) do 2002 ip=1,7 ra(ip+7) = ra(ip) 2002 rax(ip+7) = rax(ip) CCC c *** Special case H, I, J, and K NN-partial waves set = 0. IF (nifty(7) .eq. 1) go to 8006 do 2424 nwave = 1,nwaves IF ( nwave .eq. 10 $ .OR. $ nwave .eq. 11 $ .OR. $ nwave .eq. 12 $ .OR. $ nwave .eq. 16 $ .OR. $ nwave .ge. 24 ) $ then t(nwave) = 0.0 tx(nwave) = 0.0 endif 2424 continue do 2525 nica = 4,7 ra(nica) = 0.0 rax(nica) = 0.0 CCC As above (+7) ra(nica+7) = 0.0 rax(nica+7) = 0.0 2525 continue 8006 call MIXUPN (bn, bnf, t, tx) do 410 n = 1,32 retij(n) = t(n) imtij(n) = tx(n) 410 continue return c-----formats 902 format('0','# of energies, and # of deltas for each Ecm =', $ 2i10,//) 919 format(2i5) 921 format(' ******nes > nstore=', i5,'nes set = store *****,',/) 923 format('0', ' Tmatrices (TNCM)', $ ' calculated from Arndt SM83', $ ' phases nstore(nes) nwaves =', 2i5/) 940 format('0','kap Ecm plab Tlab first line after rewave', $ ' imwave',/ $ ' ', 13x, '1S0', 23x, '1D2', 23x, '1G4', 23x, '3P0', $ 23x, '3P1',/ $ ' ', 13x, '3P2', 23x, '3F2', 23x, '3F3', 23x, '3F4', $ 23x, '3H4',/ $ ' ', 13x, '3H5', 23x, '3H6', 23x, '1S0', 23x, '1P1', $ 23x, '1F3',/ $ ' ', 13x, '1H5', 23x, '3S1', 23x, '3D1', 23x, '3D2', $ 23x, '3D3',/ $ ' ', 13x, '3G3', 23x, '3G4', 23x, '3G5', 23x, '1I6', $ 23x, '3J6',/ $ ' ', 13x, '3J7', 23x, '3I5', 23x, '3I6', 23x, '3I7', $ 23x, '1J7',/ $ ' ', 13x, '3K7', 23x, '3K8') 941 format('0', 4(f10.3),/, 100( 10e13.5,/)) 9001 format(8f10.6) 9003 format(1x,f10.3/120(8f13.5/)) end