subroutine TNUCTH (kp, k, k0, cthnuc, MN, ApB, AmB, CpD, CmD, E) ************************************************************************ c *** See R.H. Landau's program LPOTT, 1981 c *** TNUCTH calculates the theta dependent amplitudes in p-nucleus frame c with angle dependent potential c *** tofth(j) code, t = ..., for j=... c 1 retp-spinindep 2 imtp-spin indep c 3 retn- spin indep 4 imtn-spin indep c 5 retp-spin dependt 6 imtp-spin dep c 7 retn-soin dept 8 imtn spon-depend c code for mApB(i),AmB(i),mAmB(i),etc..i=1 real part prot c i=2, real part neutrons. i=3 img part prot,i=4 im neut. c the subroutine changes also mAmB to AmB etc. implicit real*8 (a-h, i, k, m, o-z) dimension nifty(20) dimension tofth(8), pl(10), plp(10) dimension km(4), qim(4), em(4), qm(2), pm(4) dimension kapm(4), alph(4), beta(4) dimension ApB(4), AmB(4), CpD(4), CmD(4), E(4) dimension MApB(4), MAmB(4), MCpD(4), MCmD(4), ME(4) common /params/ hbarc, pi, Mp, xMN, nz, na, nes, nwaves common /switch/ nifty common /Tmats/ MApB, MAmB, MCpD, MCmD, ME c *** Statement FUNCTIONS now follow Ep(p) = SQRT( Mp2 + p*p ) EN(p) = SQRT( MN2 + p*p) c *** the p-nucleon c.o.m.. momemtum**2 for a given s kcm2(s) = (s - (Mp + MN)**2) * (s - (Mp - MN)**2)/4./s c *** internal nucleon energy, new optimal definition Eint(p,pp) = SQRT( MN2 + (rhl1**2) * p * p * 0.25 + $ (rhl2**2) * pp * pp * 0.25 - $ 0.5 * rhl1 * rhl2 * p * pp * cthnuc ) c *** softh is the new, optimal, theta-dependent s softh(p,pp) = Mp2 + MN2 + rhl1 * p * p - rhl2 * p * pp * cthnuc + $ 2. * Ep(p) * Eint(p,pp) c *** End of FUNCTION definitions *** c >>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< nit = 1 amass = na * MN/xMN Ma = na * xMN rhl1 = (amass + 1.)/amass rhl2 = (amass - 1.)/amass MN2 = MN * MN Mp2 = Mp * Mp m12 = (Ma - MN)**2 rhl3 = Mp2 - MN2 Ma2 = Ma * Ma if (nifty(9) .ne. 0) then do ii = 1, 4 ApB(ii) = 0. AmB(ii) = 0. CpD(ii) = 0. CmD(ii) = 0. E(ii) = 0. enddo return endif sin = softh( k, kp) sout = softh( kp, k) ENin = Eint(k,kp) ENout = Eint(kp,k) IF ((nifty(5) .gt. 5) .AND. (nifty(5) .ne. 9)) go to 20 c *** Calculate on-shell kappa with 2 body energy son = softh( k0, k0) kapon2 = kcm2( son ) kapon = SQRT( kapon2 ) kap2 = kcm2( sin ) kap = SQRT( kap2 ) kapp2 = kcm2( sout ) kapp = SQRT( kapp2 ) IF (kapp .le. 0.) kapp = 0.0001 IF (kap .le. 0.) kap = 0.0001 cthn = ( ( Ep(kap) * Ep(kapp)- Ep(k) * Ep(kp) )/kapp/kap ) $ + k * kp * cthnuc/kap/kapp IF (nifty(5) .eq. 2) cthn = cthnuc IF (abs(cthn) .le. 2. .OR. nifty(5) .eq. 9) go to 10 10 continue gamth = SQRT( Ep(kap) * Ep(kapp) * EN(kap) * EN(kapp) $ /(Ep(k) * Ep(kp) * ENin * ENout) ) costh = cthn factor = gamth c *** Choose defintion of energy c *** nifty(5) = 0 % use softh(on shell k) c *** 1 use sin of theata (incoming momentum) c *** 9 use softh(on shell k)plus aay angles c *** 2 softh (on shell), no angle t.f. kape = kapon IF (nifty(5) .eq. 1) kape = kap ekape = Ep( kape ) + EN( kape ) IF (nifty(5) .lt. 5) go to 70 IF (nifty(5) .eq. 9) go to 40 c *** 3 body definition of energy c *** nifty(5)= 5: e3b with nonfolded or w0-folded t + old,optimal angle c *** = 7: same as above yet with aay magic angle t.f.(standard) c *** c *** N.B. n(5) = 7 is the standard case, whereas 6 is for trial c *** don't use 6 if have w0-folded t s 20 continue spinuc = Ep( k0 ) + SQRT( Ma2 + k0*k0 ) q2 = kp*kp + k*k - 2.*kp*k*cthnuc qk = k*kp*cthnuc - k*k be = 5. * nifty(2) pkap2 = rhl2 * rhl2 * ( k*k + q2/4. + 174. * 174. + qk) ekape = spinuc - be - SQRT( pkap2 + m12) ekape = ekape * ekape - pkap2 c *** E3b-folded t s being used c *** non relativistic definition for folding if (nifty(5) .ne. 6) go to 30 twomu = 2. * (Mp + MN) * (Ma - MN)/(Mp + Ma) ekape = spinuc - Ma + MN + 5 - be - (k*k + (q2 - k0*k0)/4. + qk) $ * rhl2 * rhl2/twomu 30 rhl = 1. IF (ekape .lt. 0.) rhl = -1. ekape = rhl * SQRT( rhl * ekape ) IF (nifty(5) .eq. 5) go to 70 c *** Nifty(5) = 6, 7, 9 use magic vector prescription(aay) c *** 1st two components x, y of kin, 3rd and 4th of x, y of kout 40 rhl = SQRT( 1. - cthnuc**2 ) ki0 = Ep( k ) + Eint( k, kp ) kf0 = Ep( kp ) + Eint( kp, k ) km(1) = k km(2) = 0. km(3) = kp * cthnuc km(4) = kp * rhl do 50 j = 1,2 j2 = j + 2 qm(j) = km(j2) - km(j) pm(j) = -km(j)/amass + rhl2 * qm(j)/2. pm(j2) = pm(j) - qm(j) kapm(j) = km(j) + pm(j) kapm(j2) = kapm(j) alph(j) = 0. alph(j2) = 0. qim(j) = 0.5 * ( km(j) - pm(j) - alph(j) * kapm(j) ) qim(j2) = 0.5 * ( km(j2) - pm(j2) - alph(j2) * kapm(j2) ) beta(j) = qim(j) * kapm(j)/ki0/(ki0 + SQRT( sin )) beta(j2) = qim(j2) * kapm(j2)/kf0/(kf0 + SQRT( sout )) IF (j .ne. 2) go to 50 beta(1) = beta(1) + beta(2) beta(2) = beta(1) beta(3) = beta(3) + beta(4) beta(4) = beta(3) 50 continue do 60 j = 1,4 em(j) = qim(j) - beta(j) * kapm(j) 60 continue kap = SQRT( em(1)**2 + em(2)**2 ) kapp = SQRT( em(3)**2 + em(4)**2 ) costh = ( em(1) * em(3) + em(2) * em(4) )/kap/kapp factor = SQRT( Ep(kap) * Ep(kapp) * EN(kap) * EN(kapp) $ /(Ep(k) * Ep(kp) * ENin * ENout) ) c *** c *** Proton part, see eqtns 4.30-4.40 of Paez thesis c *** pmj = (3. * amass + 1.)/(4. * amass) pmb = pmj - beta(1) * rhl2/2. pmc = -rhl2/4. - beta(1) * rhl2/2. pmbp = pmj - beta(3) * rhl2/2. pmcp = -rhl2/4. - beta(3) * rhl2/2. pmj1 = pmb * pmbp - pmc * pmcp pmj1sq = pmj1 * pmj1 c *** Magnitude of vector kap+kapp 70 continue pmj3 = 2. * kap * kapp * costh kaplkp = SQRT( kap**2 + kapp**2 + pmj3 ) c *** Magnitude of vector kapp-kap kapmsk = SQRT( kap**2 + kapp**2 - pmj3) c *** Magnitude of vector k+kp pmj4 = 2. * k * kp * cthnuc kplkp = SQRT( k**2 + kp**2 + pmj4 ) c *** Magnitude of vector kp-k kpmsk = SQRT( k**2 + kp**2 - pmj4 ) c *** Run with shifted p-N subenergy( for cex) 700 eshift = ekape IF (MN .lt. 1000.) then call TNOFF2 (kapp, kap, eshift, kp, k, costh) endif c *** c *** call tnoff2.....results tfed via sectr c *** lap = 1, 2, 3, 4, .... re/im p/n with angle dependence c *** a + b a - b amplitudes c *** do 51 lap = 1,4 ApB(lap) = MApB(lap) * factor AmB(lap) = MAmB(lap) * factor E(lap) = 0.0 51 continue pmj1rt = SQRT( pmj1sq ) IF (pmj1rt .eq. 0.) go to 52 do 102 lys = 1,4 E(lys) = ME(lys) * factor * pmj1/pmj1rt 102 continue 52 deth = 0. eth = 0. geth = 0.5 * (pmb - pmbp + pmcp - pmc) * kpmsk/kaplkp efeth = 0.5 * (pmb + pmbp + pmc + pmcp) * kplkp/kaplkp if (kapmsk .eq. 0.) go to 62 deth = 0.5 * (pmb + pmbp - pmc - pmcp) * kpmsk/kapmsk eth = 0.5 * (pmb - pmbp + pmc - pmcp) * kpmsk/kapmsk 62 do 63 li = 1,4 CpD(li) = (MCpD(li) * (deth**2) + MCmD(li) * (geth**2)) *factor CmD(li) = (MCmD(li) * (efeth**2) + MCpD(li) * (eth**2)) *factor 63 continue if (nifty(6) .ne. 3) return c *** Nifty(6) = 3 protons with spin 0 * 0.5 the old way x = costh sinth = SQRT( abs( 1. - cthnuc**2 ) ) if (sinth .eq. 0.) sinth = 0.0001 do 81 j = 1,4 AmB(j) = ApB(j) CpD(j) = E(j)/sinth 81 continue return c *** FORMATS 100 format(1h ,21h cos thea pi-n gt 2= , e15.3) 250 format(1h ,4h em ,4e13.5) end