fix many bug for angular momentum transfer not zero
This commit is contained in:
Ryan@Home
parent
9d2a251106
commit
e47bd2dc02
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
@ -237,28 +236,50 @@ class DWBA_ZR:
|
|||
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,17 +430,11 @@ 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
|
||||
self.PreCalLegendreP(theta_deg)
|
||||
|
|
@ -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()
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user