From e47bd2dc02b05dfa718f778f537ad549bd464fb6 Mon Sep 17 00:00:00 2001 From: "Ryan@Home" Date: Sun, 23 Feb 2025 18:46:37 -0500 Subject: [PATCH] fix many bug for angular momentum transfer not zero --- Raphael/DWBA_ZR.py | 32 +++++-- Raphael/dwba_zr.py | 216 ++++++++++++++++++++++++++------------------- 2 files changed, 149 insertions(+), 99 deletions(-) diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py index 0319557..d73bd20 100755 --- a/Raphael/DWBA_ZR.py +++ b/Raphael/DWBA_ZR.py @@ -4,20 +4,38 @@ import time import matplotlib.pyplot as plt from dwba_zr import DWBA_ZR -haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10) -haha.FindBoundState() -haha.CalRadialIntegral() +# haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10) -# haha.PrintRadialIntegral() +haha = DWBA_ZR("16O", "d", "p", "17O", "5/2+", "0d5/2", 0.0, 10) +haha.FindBoundState() # haha.boundState.PlotBoundState() +haha.CalRadialIntegral() + +# haha.PlotScatteringMatrix(False) + +# haha.PlotIncomingDistortedWave(2, 1, 20) +# haha.PlotOutgoingDistortedWave(15, 15.5, 20) + +# haha.PrintRadialIntegral() # haha.PlotRadialIntegral() +# haha.PlotRadialIntegralSigle(1, 1, -0.5) -# haha.PlotDistortedWave(True, 2, 1, 20) -# haha.PlotDistortedWave(False, 2, 1.5, 20) -haha.CalAngDistribution() + +haha.CalAngDistribution(0, 180, 1) +haha.PrintAngDist() haha.PlotAngDist() +###################################### + +# print(haha.GetPreCalNineJ(1, 1, 3, 3.5)) + +# print(haha.Gamma(1, 1, 3, 3.5, 0, 1, 0.5)) + + + +# haha.PreCalLegendreP(5) +# haha.Beta(0, 1, 0.5) diff --git a/Raphael/dwba_zr.py b/Raphael/dwba_zr.py index 0d74c1c..bd7a303 100755 --- a/Raphael/dwba_zr.py +++ b/Raphael/dwba_zr.py @@ -125,10 +125,8 @@ class DWBA_ZR: print("====================== Incoming wave function ") op.AnCai(A_A, Z_A, self.ELab) - self.maxL = 15 self.dwI = DistortedWave(nu_A, nu_a, self.ELab) - self.dwI.maxL = self.maxL self.dwI.PrintInput() self.dwI.ClearPotential() self.dwI.AddPotential(WoodsSaxonPot( -op.v, op.r0, op.a), False) @@ -141,7 +139,9 @@ class DWBA_ZR: self.mass_I = self.dwI.mu # reduced mass of incoming channel self.k_I = self.dwI.k # wave number of incoming channel - self.maxL = self.dwI.maxL + self.maxL1 = self.dwI.maxL + + self.maxL2 = self.maxL1 + self.l Ecm_I = self.dwI.Ecm @@ -153,7 +153,7 @@ class DWBA_ZR: self.dwO = DistortedWave(nu_B, nu_b, Eout) self.dwO.spin_A = self.spin_B - self.dwO.maxL = self.maxL + self.dwO.maxL = self.maxL2 self.dwO.PrintInput() self.dwO.ClearPotential() self.dwO.AddPotential(WoodsSaxonPot( -op.v, op.r0, op.a), False) @@ -175,9 +175,9 @@ class DWBA_ZR: D0 = 1.55e+4 # for (d,p) self.massBoverMassA = A_B/A_A - self.ffactor = np.sqrt(4*np.pi)/k_I /k_O * A_B/A_A + self.ffactor = np.sqrt(4*np.pi)/k_I /k_O - self.xsecScalingfactor = A_B / A_A *D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1) + self.xsecScalingfactor = D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1) self.PreCalNineJ() @@ -189,31 +189,30 @@ class DWBA_ZR: else: return f"{int(2*spin):+d}/2" - - def PrintRadialIntegral(self): - for index1 in range(0, int(2*self.spin_a) + 1): - for index2 in range(0, int(2*self.spin_b) + 1): - print(f"======================= J1 = L{self.FormatSpin(index1-self.spin_a)}, J2 = L{self.FormatSpin(index2-self.spin_b)}") - for L1 in range(0, self.maxL+1): - print("{", end="") - for L2 in np.arange(abs(L1 - self.l), L1 + self.l + 1): - J1 = L1 + index1 - self.spin_a - J2 = int(L2) + index2 - self.spin_b - indexL2 = int(L2 - abs(L1-self.l)) - print(f"{{{L1:2d}, {J1:4.1f}, {int(L2):2d}, {J2:4.1f}, {np.real(self.radialInt[L1][index1][indexL2][index2]):12.4e} + {np.imag(self.radialInt[L1][index1][indexL2][index2]):12.4e}I}},", end="") - print("},") - print("=========================== end of Radial Integrals.") - - def FindBoundState(self): self.boundState.FindPotentialDepth(-70, -45, 0.5) + def ConvertLJ2RadialIndex(self, L1:int, J1, L2:int, J2): + index1 = int(J1 - L1 + self.spin_a) + indexL2 = int(L2 - L1 + self.l) + index2 = int(J2 - L2 + self.spin_b) + return [L1, index1, indexL2, index2] + + def ConvertRadialIndex2LJ(self, in1:int, in2:int, in3:int, in4:int): + L1 = in1 + J1 = L1 + in2 - self.spin_a + L2 = in3 + L1 - self.l + J2 = L2 + in4 - self.spin_b + return [L1, J1, L2, J2] + +########################################################### def CalRadialIntegral(self): start_time = time.time() sm_I, wfu_I = self.dwI.CalScatteringMatrix() self.dwI.PrintScatteringMatrix() sm_O, wfu_O_temp = self.dwO.CalScatteringMatrix() + self.dwO.PrintScatteringMatrix() #============ rescale the outgoing wave print("====================== Scaling the outgoing wave") @@ -221,14 +220,14 @@ class DWBA_ZR: self.rpos_O = [] rpos_O_filled = False self.wfu_O = [] - for L in range(0, self.maxL+1): + for L2 in range(0, self.maxL2 + 1): temp_wfu_L = [] for k in range(0, int(2*self.spin_b)+1): - wfu_O_inter_real = interp1d(rpos_O_temp, np.real(wfu_O_temp[L][k]), kind='cubic') - wfu_O_inter_imag = interp1d(rpos_O_temp, np.imag(wfu_O_temp[L][k]), kind='cubic') + wfu_O_inter_real = interp1d(rpos_O_temp, np.real(wfu_O_temp[L2][k]), kind='cubic') + wfu_O_inter_imag = interp1d(rpos_O_temp, np.imag(wfu_O_temp[L2][k]), kind='cubic') temp_wfu = [] for r in self.dwI.rpos: - if r > 20 : + if r > max(rpos_O_temp) : break if rpos_O_filled == False: self.rpos_O.append(r) @@ -236,29 +235,51 @@ class DWBA_ZR: rpos_O_filled = True temp_wfu_L.append(temp_wfu) self.wfu_O.append(temp_wfu_L) - - self.dwO.PrintScatteringMatrix() print("====================== Calculating Radial integrals") - - self.radialInt = np.zeros((self.maxL+1, int(2*self.spin_a+1), int(2*self.l+1), int(2*self.spin_b+1)), dtype=complex) - + self.radialInt = np.zeros((self.maxL1+1, int(2*self.spin_a+1), int(2*self.l+1), int(2*self.spin_b+1)), dtype=complex) bs = self.boundState.GetBoundStateWF() - for L1 in range(0, self.maxL+1): - for index1 in range(0, len(wfu_I[L1])): + for L1 in range(0, self.maxL1+1): + for J1 in np.arange(L1-self.spin_a, L1 + self.spin_a + 1, 1): + if J1 < 0 : + continue + index1 = int(J1 - L1 + self.spin_a) wf1 = wfu_I[L1][index1] - for L2 in np.arange(abs(L1 - self.l), L1 + self.l + 1, 1): - for index2 in range(0, len(self.wfu_O[int(L2)])): + for L2 in np.arange(L1 - self.l, L1 + self.l + 1, 1): + if L2 < 0 : + continue + for J2 in np.arange(L2 - self.spin_b, L2 + self.spin_b + 1, 1): + if J2 < 0: + continue + index2 = int(J2 - L2 + self.spin_b) wf2 = self.wfu_O[int(L2)][index2] pf1 = np.exp(1j*self.dwI.CoulombPhaseShift(L1)) pf2 = np.exp(1j*self.dwO.CoulombPhaseShift(L2)) - integral = simpson (bs[:200]*wf1[:200]*wf2[:200], dx=self.boundState.dr) - indexL2 = int(L2 - abs(L1-self.l)) - self.radialInt[L1][index1][indexL2][index2] = integral * pf1 * pf2 + integral = simpson (bs*wf1*wf2, dx=self.boundState.dr) + indexL2 = int(L2 - L1 + self.l) + product = integral * pf1 * pf2 * self.massBoverMassA + self.radialInt[L1][index1][indexL2][index2] = product + # if J1 == L1 + self.spin_a and L2 == L1 + 1 and J2 == L2 - self.spin_b: + # print(f"{L1:2d}, {J1:4.1f}({index1}), {L2:2d}({indexL2}), {J2:4.1f}({index2}), {integral * pf1 * pf2 * self.massBoverMassA:.6f}") stop_time = time.time() print(f"Total time for distorted wave and radial intergal {(stop_time - start_time) * 1000:.2f} msec") +########################################################### + + def PrintRadialIntegral(self): + for index1 in range(0, int(2*self.spin_a) + 1): + for index2 in range(0, int(2*self.spin_b) + 1): + print(f"======================= J1 = L{self.FormatSpin(index1-self.spin_a)}, J2 = L{self.FormatSpin(index2-self.spin_b)}") + for L1 in range(0, self.maxL1+1): + print("{", end="") + for L2 in np.arange(L1 - self.l, L1 + self.l + 1, 1): + J1 = L1 + index1 - self.spin_a + J2 = int(L2) + index2 - self.spin_b + indexL2 = int(L2 - L1 + self.l) + print(f"{{{L1:2d}, {J1:4.1f}, {int(L2):2d}, {J2:4.1f}, {np.real(self.radialInt[L1][index1][indexL2][index2]):11.4e}+{np.imag(self.radialInt[L1][index1][indexL2][index2]):11.4e}I}}, ", end="") + print("},") + print("=========================== end of Radial Integrals.") def PlotRadialIntegral(self): if self.radialInt is None: @@ -267,21 +288,20 @@ class DWBA_ZR: spin_b = self.spin_b spin_a = self.spin_a l = self.l - maxL = self.maxL + maxL1 = self.maxL1 fig, axes = plt.subplots(int(2*spin_b+1)*int(2*l+1), int(2*spin_a+1), figsize=(6*int(2*spin_a+1), 4*int(2*spin_b+1)*int(2*l+1))) for index2 in range(0, int(2*spin_b) + 1): for index1 in range(0, int(2*spin_a) + 1): - l_list = [] for indexL2 in range(0, int(2*l) + 1): haha = [] - for L1 in range(0, maxL+1): - J1 = L1 + index1 - spin_a - L2 = int(abs(L1-l) + indexL2) - J2 = L2 + index2 - spin_b - if J1 < 0 or J2 < 0 : - continue + l_list = [] + for L1 in range(0, maxL1+1): + # J1 = L1 + index1 - spin_a + # L2 = int(L1 - l + indexL2) + # J2 = L2 + index2 - spin_b + # [dummy, J1, L2, J2] = self.ConvertRadialIndex2LJ(L1, index2, indexL2, index2) l_list.append(L1) haha.append(self.radialInt[L1][index1][indexL2][index2]) axes[int(2*l+1)*index2 + indexL2, index1].plot(l_list, np.real(haha), label="Real", marker='o') @@ -289,85 +309,104 @@ class DWBA_ZR: axes[int(2*l+1)*index2 + indexL2, index1].legend() axes[int(2*l+1)*index2 + indexL2, index1].set_xlabel('L1') axes[int(2*l+1)*index2 + indexL2, index1].set_ylabel('Value') - axes[int(2*l+1)*index2 + indexL2, index1].set_title(f'Radial Int. vs L for Spin J1 = L{self.FormatSpin(index1-spin_a)}, L2 = L1{indexL2-l:+d}, J2 = L{self.FormatSpin(index2-spin_b)}.') - axes[int(2*l+1)*index2 + indexL2, index1].set_xlim(-1, maxL+1) + axes[int(2*l+1)*index2 + indexL2, index1].set_title(f'Radial Int. vs L for Spin J1 = L1{self.FormatSpin(index1-spin_a)}, L2 = L1{indexL2-l:+d}, J2 = L2{self.FormatSpin(index2-spin_b)}.') + axes[int(2*l+1)*index2 + indexL2, index1].set_xlim(-1, maxL1+1) + axes[int(2*l+1)*index2 + indexL2, index1].set_xticks(np.arange(0, maxL1+1, 2)) axes[int(2*l+1)*index2 + indexL2, index1].grid() plt.tight_layout() plt.show(block=False) input("Press Enter to continue...") + def PlotRadialIntegralSigle(self, dJ1, dL, dJ2): + if self.radialInt is None: + print("Radial integral not computed.") + return + haha = [] + l_list = [] + for L1 in range(0, self.maxL1+1): + l_list.append(L1) + [dummy, index1, indexL2, index2] = self.ConvertLJ2RadialIndex(L1, L1 + dJ1, L1 + dL, L1 + dL + dJ2) + haha.append(self.radialInt[L1][index1][indexL2][index2]) + print(f"{L1:2d}, {L1 + dJ1:4.1f}({index1}), {L1 + dL:.0f}({indexL2}), {L1 + dL + dJ2:4.1f}({index2}), {haha[-1]:.6f}") + plt.plot(l_list, np.real(haha), label="Real", marker='o') + plt.plot(l_list, np.imag(haha), label="Imag", marker='x') + plt.xlabel("L1") + plt.title(f'Radial Int. vs L for Spin J1 = L1{self.FormatSpin(dJ1)}, L2 = L1{dL:+d}, J2 = L2{self.FormatSpin(dJ2)}.') + plt.grid() + plt.show(block=False) + input("Press Enter to continue...") + def PlotScatteringMatrix(self, isIncoming): if isIncoming : self.dwI.PlotScatteringMatrix() else: self.dwO.PlotScatteringMatrix() - def PlotDistortedWave(self, isIncoming, L, J, maxR = None): - if isIncoming: - self.dwI.PlotDistortedWave(L, J, maxR) - else: - plt.plot(self.rpos_O, np.real(self.wfu_O[L][int(J-L + self.spin_b)]), label="Real") - plt.plot(self.rpos_O, np.imag(self.wfu_O[L][int(J-L + self.spin_b)]), label="Imag") - plt.title(f"Radial wave function for L={L} and J={J}") - if maxR != None : - plt.xlim(-1, maxR) - plt.legend() - plt.grid() - plt.show(block=False) - input("Press Enter to continue...") + def PlotIncomingDistortedWave(self, L, J, maxR = None): + self.dwI.PlotDistortedWave(L, J, maxR) + + + def PlotOutgoingDistortedWave(self, L, J, maxR = None): + plt.plot(self.rpos_O, np.real(self.wfu_O[L][int(J-L + self.spin_b)]), label="Real") + plt.plot(self.rpos_O, np.imag(self.wfu_O[L][int(J-L + self.spin_b)]), label="Imag") + plt.title(f"Radial wave function for L={L} and J={J}") + if maxR != None : + plt.xlim(-1, maxR) + plt.legend() + plt.grid() + plt.show(block=False) + input("Press Enter to continue...") def approximate_to_half_integer(self, value): return round(value * 2) / 2 def PreCalNineJ(self): - self.NineJ = np.zeros((self.maxL+1, int(2*self.spin_a+1), (2*self.l+1), int(2*self.spin_b+1)), dtype=float) - for L1 in range(0, self.maxL+1): + self.NineJ = np.zeros((self.maxL1+1, int(2*self.spin_a+1), (2*self.l+1), int(2*self.spin_b+1)), dtype=float) + for L1 in range(0, self.maxL1+1): for ind1 in range(0, int(2*self.spin_a+1)): for indL2 in range(0, 2*self.l+1): - for ind2 in range(0, int(2*self.spin_B+1)): + for ind2 in range(0, int(2*self.spin_b+1)): J1 = self.approximate_to_half_integer(L1 + ind1 - self.spin_a) L2 = int(L1 + indL2 - self.l) J2 = self.approximate_to_half_integer(L2 + ind2 - self.spin_b) self.NineJ[L1, ind1, indL2, ind2] = wigner_9j(self.j, self.l, self.s, J1, L1, self.spin_a, J2, L2, self.spin_b) def GetPreCalNineJ(self, L1:int, J1, L2:int, J2): - ind1 = int(J1 - L1 + self.spin_a) - indL2 = int(L1 - L2 + self.l) - ind2 = int(J2 - L2 + self.spin_b) + [dummy, ind1, indL2, ind2] = self.ConvertLJ2RadialIndex(L1, J1, L2, J2) return self.NineJ[L1, ind1, indL2, ind2] def Gamma(self, L1:int, J1, L2:int, J2, m:int, ma, mb): - if int(L1 + L2 + self.l)%2 != 0: #check if the sum of L1 + L2 + l is even return 0 else: - # fact0 = wigner_9j(S(2*self.j)/2, self.l, S(2*self.s)/2, S(2*J1)/2, L1, S(2*self.spin_a)/2, J2, S(2*L2)/2, S(2*self.spin_b)/2).evalf() - # fact0 = wigner_9j(self.j, self.l, S(2*self.s)/2, S(2*J1)/2, L1, S(2*self.spin_a)/2, J2, S(2*L2)/2, S(2*self.spin_b)/2).evalf() - # fact0 = wigner_9j(self.j, self.l, self.s, S(2*J1)/2, L1, S(2*self.spin_a)/2, J2, S(2*L2)/2, S(2*self.spin_b)/2).evalf() - # fact0 = wigner_9j(self.j, self.l, self.s, S(2*J1)/2, L1, self.spin_a, J2, S(2*L2)/2, S(2*self.spin_b)/2).evalf() - # fact0 = wigner_9j(self.j, self.l, self.s, S(2*J1)/2, L1, self.spin_a, J2, S(2*L2)/2, self.spin_b).evalf() fact0 = self.GetPreCalNineJ(L1, J1, L2, J2) if fact0 == 0: return 0 else: fact1 = pow(-1, m) * np.power(1j, L1-L2-self.l) * (2*L2+1) * np.sqrt((2*self.l+1)*(2*self.s+1)*(2*L1+1)*(2*J2+1)) - fact2 = np.sqrt( quantum_factorial(L2-m) / quantum_factorial(L2+m) ) + fact2 = np.sqrt( quantum_factorial(L2-abs(m)) / quantum_factorial(L2 + abs(m)) ) fact3 = clebsch_gordan(J2, mb-m, self.j, m-mb+ma, J1, ma) fact4 = clebsch_gordan(L1, 0, self.spin_a, ma, J1, ma) fact5 = clebsch_gordan(L2, -m, self.spin_b, mb, J2, mb-m) - fact6 = clebsch_gordan(L1, 0, self.l, 0, L2, 0) + fact6 = clebsch_gordan(L2, 0, self.l, 0, L1, 0) + # print(f"{fact1:.5f}, {fact2:.5f}, {fact3:.5f}, {fact4:.5f}, {fact5:.5f}, {fact6:.5f}") return fact0 * fact1 * fact2 * fact3 * fact4 * fact5 * fact6 def Beta(self, m:int, ma, mb): if self.radialInt is None : return result = 0 - for L1 in np.arange(0, self.maxL+1): - for J1 in np.arange(abs(L1 - self.spin_a), L1 + self.spin_a + 1, 1): - for L2 in np.arange(abs(L1 - self.l), L1 + self.l + 1, 1): - for J2 in np.arange(abs(L2 - self.spin_b), L2 + self.spin_b + 1, 1): - + for L1 in np.arange(0, self.maxL1+1): + for J1 in np.arange(L1 - self.spin_a, L1 + self.spin_a + 1, 1): + if J1 < 0 : + continue + for L2 in np.arange(L1 - self.l, L1 + self.l + 1, 1): + if L2 < 0: + continue + for J2 in np.arange(L2 - self.spin_b, L2 + self.spin_b + 1, 1): + if J2 < 0 : + continue if not obeys_triangle_rule(J1, self.j, J2): continue if not(abs(m) <= L2): @@ -375,14 +414,13 @@ class DWBA_ZR: if int(L1 + L2 + self.l) % 2 != 0: continue - index1 = int(J1 - L1 + self.spin_a) - index2 = int(J2 - L2 + self.spin_b) - indexL2 = int(L1 - abs(L1 - self.l)) gg = self.Gamma(L1, J1, L2, J2, m, ma, mb) if gg == 0: continue - lp = self.GetPreCalLegendreP(L2, m) + lp = self.legendrePArray[L2][int(abs(m))] + + [dummy, index1, indexL2, index2] = self.ConvertLJ2RadialIndex(L1, J1, L2, J2) ri = self.radialInt[int(L1)][index1][indexL2][index2] # print(f"{L1:2d}, {J1:4.1f}({index1:d}), {L2:2d}({indexL2:d}), {J2:4.1f}({index2:d}), {gg:10.6f}, {ri *self.ffactor :.10f}, {lp:10.6f}") @@ -392,16 +430,10 @@ class DWBA_ZR: def PreCalLegendreP(self, theta_deg:float, maxL:int = None, maxM:int = None): if maxL is None: - maxL = self.maxL + maxL = max(self.maxL1, self.maxL2) if maxM is None: maxM = int(self.j + self.spin_b + self.spin_a) self.legendrePArray = associated_legendre_array(maxL, maxM, theta_deg) - - def GetPreCalLegendreP(self, L2:int, m:int): - lp = self.legendrePArray[int(L2)][int(abs(m))] - if m < 0 : - lp *= (-1)**m * quantum_factorial(int(L2)+m)/ quantum_factorial(int(L2)-m) - return lp def AngDist(self, theta_deg:float) -> float: xsec = 0 @@ -425,7 +457,7 @@ class DWBA_ZR: for i in np.arange(angMin, angMax + angStep, angStep): self.angList.append(i) self.angDist.append(self.AngDist(i)) - if time.time() - progress_time > 2: + if time.time() - progress_time > 1: elapsed_time = time.time() - start_time print(f"\r Time elapsed: {elapsed_time:.2f} sec, Progress: {100 * (i - angMin) / (angMax - angMin):.1f}%", end="") progress_time = time.time()