c/* %W% latest revision %G% %U% */ subroutine TPN (x, ApB, AmB, CpD, CmD, E) ******************************************************************************* c *** Computes the coefficients A, B, C, D, and E in Bistricky, Lehar c *** and Winternitz Journal de Physique 39, page 1(1978) eq 2.14 c *** The convention for the associated Legendre polynomials c *** is the same as in Jackson's E&M which differs in c *** convention used by Dr Landau in Pl1 (not Pl prime) in c *** a negative sign. Not included here in sine = -sqrt(1.-x*x) c *** because the antisymmetrization was taken by multi- c *** plying by 2 rather than taking f(theta) + f(pi-theta). c *** Code for ApB(i), AmB(i), CpD(i), CmD(i), E(i) c *** i = 1 for real part protons, i = 2 real part neutrons c *** i = 3 for im part protons, i = 4 im part neutrons c *** REVISION HISTORY *** c 2/18/88 (M.S.) calculate E amplitude directly form Tspin c output rather than Mixupn c 3/1/88 (M.S.) Remove 1/2 from calculation of E term so will c agree with old calculation from Mixupn implicit real*8 (a-h, i, k, o-z) complex*16 zi complex*16 Ap, An, Bp, Bn, Cp, Cn, Dp, Dn, Ep, En CCC Added for off-diagonal problem dimension t(40), tx(40), ra(14), rax(14) c dimension t(40), tx(40), ra(7), rax(7) CCC dimension t1m1(40), t11(40), tss(40), t00(40), t01(40), t10(40) dimension br(40), brf(40) dimension brfa(36) dimension tssa(36), t11a(36), t1m1a(36), t00a(36) dimension t01a(36), t10a(36) dimension ApB(4), AmB(4), CpD(4), CmD(4), E(4) dimension Pl(0:9), Plp(0:9), Pldp(0:9) common /tjmats/ t, tx, ra, rax common /amplit/ Ap, An, Bp, Bn, Cp, Cn, Dp, Dn, $ Ep, En c******************************************************************************* call LEGPOL (x, Pl, 10) call PLPRME (x, Plp, 10) call PLDBLP (x, Pldp, 10) call TSPIN (t1m1, t11, t00, t01, t10, tss, t, tx, ra, rax) call MIXUPN (br, brf, t, tx) sinsq = 1.0 - x*x mm = 0 sine = SQRT( 1.0 - x*x ) zi = (0.0, 1.0) Rebnfp = 0. Imbnfp = 0. Rebnfn = 0. Imbnfn = 0. ReAp = 0. ImAp = 0. ReAn = 0. ImAn = 0. ReBp = 0. ImBp = 0. ReBn = 0. ImBn = 0. ReDp = 0. ImDp = 0. ReDn = 0. ImDn = 0. ReCp = 0. ImCp = 0. ReCn = 0. ImCn = 0. do 10 l = 0,8 do 20 m = 1,4 mjp = mm + m brfa(mjp) = brf(mjp) * Plp(l) * sine t01a(mjp) = t01(mjp) * Plp(l) * sine t10a(mjp) = t10(mjp) * Plp(l) * sine t00a(mjp) = t00(mjp) * Pl(l) tssa(mjp) = tss(mjp) * Pl(l) t11a(mjp) = t11(mjp) * Pl(l) t1m1a(mjp) = t1m1(mjp) * Pldp(l) * sinsq 20 continue mm = mm + 4 10 continue c *** Form separate sums for Re and Im bn for protons and c *** neutrons and Re and Im bnf for protons and neutrons Rebnfp = brfa(1) Imbnfp = brfa(2) Rebnfn = brfa(3) Imbnfn = brfa(4) c *** For A = 0.5 * (M11 + M00 - M1m1) ReAp = 0.5 * (t11a(1) + t00a(1) - t1m1a(1)) ImAp = 0.5 * (t11a(2) + t00a(2) - t1m1a(2)) ReAn = 0.5 * (t11a(3) + t00a(3) - t1m1a(3)) ImAn = 0.5 * (t11a(4) + t00a(4) - t1m1a(4)) c *** For B = 0.5 * (M11 + Mss + M1m1) ReBp = 0.5 * (t11a(1) + t1m1a(1) + tssa(1)) ImBp = 0.5 * (t11a(2) + t1m1a(2) + tssa(2)) ReBn = 0.5 * (t11a(3) + t1m1a(3) + tssa(3)) ImBn = 0.5 * (t11a(4) + t1m1a(4) + tssa(4)) c *** First terms for D = 1/2 * (-M11 + M00 + M1m1)/cos(theta) ReDp = 0.5 * (-t11a(1) + t1m1a(1) + t00a(1)) ImDp = 0.5 * (-t11a(2) + t1m1a(2) + t00a(2)) ReDn = 0.5 * (-t11a(3) + t1m1a(3) + t00a(3)) ImDn = 0.5 * (-t11a(4) + t1m1a(4) + t00a(4)) c *** For C = 0.5 * (M11 - Mss + M1m1) ReCp = 0.5 * (t11a(1) - tssa(1) + t1m1a(1)) ImCp = 0.5 * (t11a(2) - tssa(2) + t1m1a(2)) ReCn = 0.5 * (t11a(3) - tssa(3) + t1m1a(3)) ImCn = 0.5 * (t11a(4) - tssa(4) + t1m1a(4)) c *** For E = 0.5 * 1/sqrt(2) * ( M10 - M01 ) *** ReEp = 1./sqrt(2.) * (t10a(1) - t01a(1)) ImEp = 1./sqrt(2.) * (t10a(2) - t01a(2)) ReEn = 1./sqrt(2.) * (t10a(3) - t01a(3)) ImEn = 1./sqrt(2.) * (t10a(4) - t01a(4)) do 30 ml = 5,33,4 Rebnfp = Rebnfp + brfa( ml ) Imbnfp = Imbnfp + brfa( ml+1 ) Rebnfn = Rebnfn + brfa( ml+2 ) Imbnfn = Imbnfn + brfa( ml+3 ) c *** For A ReAp = ReAp + 0.5 * (-t1m1a(ml) + t00a(ml) + t11a(ml)) ImAp = ImAp + 0.5 * (-t1m1a(ml+1) + t00a(ml+1) + t11a(ml+1)) ReAn = ReAn + 0.5 * (-t1m1a(ml+2) + t00a(ml+2) + t11a(ml+2)) ImAn = ImAn + 0.5 * (-t1m1a(ml+3) + t00a(ml+3) + t11a(ml+3)) c *** For B ReBp = ReBp + 0.5 * (t11a(ml) + t1m1a(ml) + tssa(ml)) ImBp = ImBp + 0.5 * (t11a(ml+1) + t1m1a(ml+1) + tssa(ml+1)) ReBn = ReBn + 0.5 * (t11a(ml+2) + t1m1a(ml+2) + tssa(ml+2)) ImBn = ImBn + 0.5 * (t11a(ml+3) + t1m1a(ml+3) + tssa(ml+3)) c *** For D ReDp = ReDp + 0.5 * (-t11a(ml) + t1m1a(ml) + t00a(ml)) ImDp = ImDp + 0.5 * (-t11a(ml+1) + t1m1a(ml+1) + t00a(ml+1)) ReDn = ReDn + 0.5 * (-t11a(ml+2) + t1m1a(ml+2) + t00a(ml+2)) ImDn = ImDn + 0.5 * (-t11a(ml+3) + t1m1a(ml+3) + t00a(ml+3)) c *** For C ReCp = ReCp + 0.5 * (t11a(ml) - tssa(ml) + t1m1a(ml)) ImCp = ImCp + 0.5 * (t11a(ml+1) - tssa(ml+1) + t1m1a(ml+1)) ReCn = ReCn + 0.5 * (t11a(ml+2) - tssa(ml+2) + t1m1a(ml+2)) ImCn = ImCn + 0.5 * (t11a(ml+3) - tssa(ml+3) + t1m1a(ml+3)) c *** for E ReEp = ReEp + 1./sqrt(2.) * (t10a(ml) - t01a(ml)) ImEp = ImEp + 1./sqrt(2.) * (t10a(ml+1) - t01a(ml+1)) ReEn = ReEn + 1./sqrt(2.) * (t10a(ml+2) - t01a(ml+2)) ImEn = ImEn + 1./sqrt(2.) * (t10a(ml+3) - t01a(ml+3)) 30 continue c *** a factor of 1/sqrt(2) not written because a factor c *** of sqrt(2) was not included in t, tx (they would cancel) Ap = ReAp + zi * ImAp An = ReAn + zi * ImAn Bp = ReBp + zi * ImBp Bn = ReBn + zi * ImBn Cp = ReCp + zi * ImCp Cn = ReCn + zi * ImCn Ep = ReEp + zi * ImEp En = ReEn + zi * ImEn c *** To avoid division by zero take into account that c *** D is equal to 0 at x = 0, theta = 90 degrees IF (x .eq. 0.) then ReDp = 0.0 ImDp = 0.0 ReDn = 0.0 ImDn = 0.0 Dp = (0.0, 0.0) Dn = (0.0, 0.0) else ReDp = ReDp/x ImDp = ImDp/x ReDn = ReDn/x ImDn = ImDn/x Dp = ReDp + zi * ImDp Dn = ReDn + zi * ImDn endif c *** Code for ApB(i), AmB(i) etc c *** i = 1 for real part protons c *** i = 2 for real part neutrons c *** i = 3 for imag part protons c *** i = 4 for imag part neutrons c *** This is to form ApB = 0.5(A+B), AmB = 0.5(A-B) c *** CpD = 0.5(C+D), CmD = 0.5(C-D) and E in Bistricky and Lehar ApB(1) = 0.5 * (ReAp + ReBp) ApB(2) = 0.5 * (ReAn + ReBn) ApB(3) = 0.5 * (ImAp + ImBp) ApB(4) = 0.5 * (ImAn + ImBn) AmB(1) = 0.5 * (ReAp - ReBp) AmB(2) = 0.5 * (ReAn - ReBn) AmB(3) = 0.5 * (ImAp - ImBp) AmB(4) = 0.5 * (ImAn - ImBn) CpD(1) = 0.5 * (ReCp + ReDp) CpD(2) = 0.5 * (ReCn + ReDn) CpD(3) = 0.5 * (ImCp + ImDp) CpD(4) = 0.5 * (ImCn + ImDn) CmD(1) = 0.5 * (ReCp - ReDp) CmD(2) = 0.5 * (ReCn - ReDn) CmD(3) = 0.5 * (ImCp - ImDp) CmD(4) = 0.5 * (ImCn - ImDn) c--------RHL change 4/83 multiply by i to get E------ c *** note E is the same as bnf E(3) = ReEp E(4) = ReEn E(1) = -ImEp E(2) = -ImEn return end