diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py index 7051075..df47c03 100755 --- a/Raphael/DWBA_ZR.py +++ b/Raphael/DWBA_ZR.py @@ -1,474 +1,53 @@ #!/usr/bin/env python3 -from boundState import BoundState -from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot -import matplotlib.pyplot as plt - -# boundState = BoundState(16, 8, 1, 0, 1, 0, 0.5, -3.273) -# boundState = BoundState(16, 8, 1, 0, 0, 2, 2.5, -4.14) -# boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.30) -# boundState.FindPotentialDepth(-75, -50, 0.1) -# # boundState.PrintWF() -# boundState.PlotBoundState() - -# exit() - -from distortedWave import DistortedWave - -# dw = DistortedWave("60Ni", "p", 30) -# dw.ClearPotential() -# dw.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.20, 0.669), False) -# dw.AddPotential(WS_SurfacePot(-6.878j, 1.28, 0.550), False) -# dw.AddPotential(SpinOrbit_Pot(-5.250 + 0.162j, 1.02, 0.590), False) -# dw.AddPotential(CoulombPotential(1.258), False) - -# dw = DistortedWave("60Ni", "d", 60) -# dw.PrintInput() -# dw.ClearPotential() -# dw.AddPotential(WoodsSaxonPot(-81.919, 1.15, 0.768), False) -# dw.AddPotential(WoodsSaxonPot(-4.836j, 1.33, 0.464), False) -# dw.AddPotential(WS_SurfacePot(-8.994j, 1.373, 0.774), False) -# dw.AddPotential(SpinOrbit_Pot(-3.557, 0.972, 1.011), False) -# dw.AddPotential(CoulombPotential(1.303), False) - - -# dw.CalScatteringMatrix() -# dw.PrintScatteringMatrix() - -# dw.PlotDCSUnpolarized(180, 1) - -# exit() - -# for i in range(1, 19): -# theta = 10*i -# # ruth = dw.RutherFord(theta) -# # coulAmp = dw.CoulombScatterintAmp(theta) -# dw.CalLegendre(theta) -# nuAmp1 = dw.NuclearScatteringAmp(-0.5, 0.5, 14) -# nuAmp2 = dw.NuclearScatteringAmp(0.5, -0.5, 14) -# # dsc = dw.DCSUnpolarized(theta, 14) -# # print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}, {dsc:10.6f}, {ruth:10.6f}") -# print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}") - - -import sys, os -import re -import numpy as np -from scipy.integrate import simpson -import matplotlib.pyplot as plt - -sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) -from IAEANuclearData import IsotopeClass - -from clebschGordan import clebsch_gordan, quantum_factorial, obeys_triangle_rule -from sympy.physics.quantum.cg import wigner_9j - -# Woods-Saxon -v = 0 -r0 = 0 -a = 0 -vi = 0 -ri0 = 0 -ai = 0 -# Woods-Saxon Surface -vsi = 0 -rsi0 = 0 -asi = 0 -# Spin-orbit -vso = 0 -rso0 = 0 -aso = 0 -vsoi = 0 -rsoi0 = 0 -asoi = 0 -# Coulomb -rc0 = 0 - -def AnCai(A : int, Z : int, E : float): - global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0 - - A3 = A**(1./3.) - v = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3 - r0 = 1.152 - 0.00776 / A3 - a = 0.719 + 0.0126 * A3 - - vi = 1.104 + 0.0622 * E - ri0 = 1.305 + 0.0997 / A3 - ai = 0.855 - 0.1 * A3 - - vsi = 10.83 - 0.0306 * E - rsi0 = 1.334 + 0.152 / A3 - asi = 0.531 + 0.062 * A3 - - vso = 3.557 - rso0 = 0.972 - aso = 1.011 - - vsoi = 0.0 - rsoi0 = 0.0 - asoi = 0.0 - - rc0 = 1.303 - -def Koning(A : int, Z : int, E : float, Zproj : float): - global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0 - - N = A-Z - A3 = A**(1./3.) - - vp1 = 59.3 + 21.*(N-Z)/A - 0.024*A - vn1 = 59.3 - 21.*(N-Z)/A - 0.024*A - - vp2 = 0.007067 + 0.00000423*A - vn2 = 0.007228 - 0.00000148*A - - vp3 = 0.00001729 + 0.00000001136 * A - vn3 = 0.00001994 - 0.00000002 * A - - vp4 = 7e-9 # = vn4 - vn4 = vp4 - - wp1 = 14.667 + 0.009629*A - wn1 = 12.195 + 0.0167*A - - wp2 = 73.55 + 0.0795*A # = wn2 - wn2 = wp2 - - dp1 = 16 + 16.*(N-Z)/A - dn1 = 16 - 16.*(N-Z)/A - - dp2 = 0.018 + 0.003802/(1 + np.exp((A-156.)/8)) # = dn2 - dn2 = dp2 - - dp3 = 11.5 # = dn3 - dn3 = dp3 - - vso1 = 5.922 + 0.003 * A - vso2 = 0.004 - - wso1 = -3.1 - wso2 = 160 - - epf = -8.4075 + 0.01378 *A - enf = -11.2814 + 0.02646 *A - - rc = 1.198 + 0.697/pow(A3,2) + 12.995/pow(A3,5) - vc = 1.73/rc * Z / A3 - - v = vp1*(1 - vp2*(E-epf) + vp3*pow(E-epf,2) - vp4*pow(E-epf,3)) + vc * vp1 * (vp2 - 2*vp3*(E-epf) + 3*vp4*pow(E-epf,2)) - #neutron - if Zproj == 0 : - v = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3)) - - r0 = 1.3039 - 0.4054 / A3 - a = 0.6778 - 0.000148 * A - - vi = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2)) - if Zproj == 0 : - vi = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2)) - - ri0 = 1.3039 - 0.4054 / A3 - ai = 0.6778 - 0.000148 * A - - vsi = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf)) - if Zproj == 0 : - vsi = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf)) - - rsi0 = 1.3424 - 0.01585 * A3 - asi = 0.5187 + 0.0005205 * A - if Zproj == 0: - asi = 0.5446 - 0.0001656 * A - - vso = vso1 * np.exp(-vso2 * (E-epf)) - if Zproj == 0: - vso = vso1 * np.exp(-vso2 * (E-enf)) - - rso0 = 1.1854 - 0.647/A3 - aso = 0.59 - - vsoi = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2)) - if Zproj == 0 : - vsoi = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2)) - - rsoi0 = 1.1854 - 0.647/A3 - asoi = 0.59 - - rc0 = rc - -def ConvertLSym(LSym :str) -> int: - if LSym == "s" : - return 0 - elif LSym == "p" : - return 1 - elif LSym == "d" : - return 2 - elif LSym == "f" : - return 3 - elif LSym == "g" : - return 4 - elif LSym == "h" : - return 5 - elif LSym == "i" : - return 6 - elif LSym == "j" : - return 7 - elif LSym == "k" : - return 8 - else : - return -1 - -#========================================== - -nu_A = "16O" -nu_a = "d" -nu_b = "p" -nu_B = "17O" -ELabPreU = 10 # MeV/u -Ex = 0.87 -J_B = "1/2+" -orbital = "1s1/2" - import time -start_time = time.time() # Start the timer +from dwba_zr import DWBA_ZR -iso = IsotopeClass() +haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10) -A_A, Z_A = iso.GetAZ(nu_A) -A_a, Z_a = iso.GetAZ(nu_a) -A_b, Z_b = iso.GetAZ(nu_b) -A_B, Z_B = iso.GetAZ(nu_B) +# haha.boundState.PlotBoundState() -A_x = abs(A_a - A_b) -Z_x = abs(Z_a - Z_b) +haha.PrintRadialIntegral() +# haha.PlotRadialIntegral() -mass_A = iso.GetMassFromSym(nu_A) -mass_a = iso.GetMassFromSym(nu_a) -mass_b = iso.GetMassFromSym(nu_b) -mass_B = iso.GetMassFromSym(nu_B) +# haha.PlotDistortedWave(True, 2, 1, 20) +# haha.PlotDistortedWave(False, 2, 1.5, 20) -mass_x = iso.GetMassFromAZ( A_x, Z_x) +j = 1/2 +sa = 1 +sb = 1/2 -if A_A < A_B : # (d,p) - A_c = A_A - Z_c = Z_A - BindingEnergy = mass_B - mass_A - mass_x + Ex -else: #(p,d) - A_c = A_B - Z_c = Z_B - BindingEnergy = mass_A - mass_B - mass_x +# hehe = haha.Gamma(0, 1, 0, 0.5, 0, 1, 0.5) +# print(hehe) -sym_A = iso.GetSymbol(A_A, Z_A) -sym_B = iso.GetSymbol(A_B, Z_B) +# haha.PreCalLegendreP(10) +# print(haha.legendrePArray) -spin_A_str = iso.GetJpi(A_A, Z_A) -spin_B_str = J_B +# jaja = haha.Beta(-1, 1, -0.5) * haha.ffactor +# print(jaja) -spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str))) -spin_B = float(eval(re.sub(r'[+-]', '', J_B))) +# lala = haha.AngDist(10) +# print(lala) -if A_a == 2 and Z_a == 1: - spin_a = 1.0 - spin_b = 0.5 -else: - spin_a = 0.5 - spin_b = 1.0 +angList = [] +xsec = [] -s = 0.5 # spin of x, neutron or proton +start_time = time.time() +for i in range(0, 181, 5): + angList.append(i) + kaka = haha.AngDist(i) + xsec.append(kaka) + print(i, kaka) -#=================== digest orbital -match = re.search(r'[a-zA-Z]', orbital) # Find first letter -if match: - index = match.start() # Get position of the first letter - -node = int(orbital[:index]) -l_sym = orbital[index:index+1] -j_sym = orbital[index+1:] -j = eval(j_sym) -l = ConvertLSym(l_sym) +stop_time = time.time() +print(f"Total time {(stop_time - start_time) :.2f} sec") -#==== check the angular conservasion -passJ = False -if obeys_triangle_rule(spin_A, spin_B, j): - passJ = True -else: - print(f"the orbital spin-J ({j}) does not consver J({nu_A}) + J({nu_B}) = {spin_A} + {spin_B}.") +import matplotlib.pyplot as plt -passS = False -if obeys_triangle_rule(spin_a, spin_b, s): - passS = True -else: - print(f"the orbital spin-s ({s})does not consver S({nu_a}) + S({nu_b}) = {spin_a} + {spin_b}.") - -passl = False -if obeys_triangle_rule(j, s, l): - passl = True -else: - print(f"the orbital spin-l ({l})does not consver J({j}) + J({s}).") - -if passJ == False or passS == False or passl == False : - print("Fail angular momentum conservation.") - exit() - -reactionStr = f"{nu_A}({spin_A_str})({nu_a},{nu_b}){nu_B}({Ex:.3f}|{spin_B_str}, {orbital}) @ {ELabPreU:.1f} MeV/u" - -print("==================================================") -print(reactionStr) - -Q_value = mass_A + mass_a - mass_b - mass_B - Ex -print(f"Transfer Orbtial : {orbital}") -print(f"Q-value : {Q_value:10.6f} MeV") -print(f"Binding : {BindingEnergy:10.6f} MeV") - -#=================== find the maximum L for partial wave -mass_I = mass_A * mass_a / (mass_A + mass_a) # reduced mass of incoming channel -hbarc = 197.3269788 # MeV.fm -k_I = np.sqrt(2*mass_I * A_a * ELabPreU)/hbarc # wave number of incoming channel -touching_Radius = 1.25*(A_A**(1./3) + A_a**(1./3)) + 10 # add 10 fm -maxL = int(touching_Radius * k_I) # maximum partial wave -print(f"max L : {maxL}") - -#================== Bound state -print("====================== Bound state ") -boundState = BoundState(A_c, Z_c, A_x, Z_x, node, l, j, BindingEnergy) -boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.30) -boundState.FindPotentialDepth(-70, -55, 0.1) -# # boundState.PrintWF() -# boundState.PlotBoundState() - -# exit() - -#================== incoming wave function -print("====================== Incoming wave function ") -AnCai(A_A, Z_A, A_a * ELabPreU) - -dwI = DistortedWave(nu_A, nu_a, ELabPreU * A_a) -dwI.maxL = maxL -dwI.PrintInput() -dwI.ClearPotential() -dwI.AddPotential(WoodsSaxonPot(-v, r0, a), False) -dwI.AddPotential(WoodsSaxonPot(-1j*vi, ri0, ai), False) -dwI.AddPotential(WS_SurfacePot(-1j*vsi, rsi0, asi), False) -dwI.AddPotential(SpinOrbit_Pot(-vso , rso0, aso), False) -dwI.AddPotential(SpinOrbit_Pot(- 1j* vsoi, rsoi0, asoi), False) -dwI.AddPotential(CoulombPotential(rc0), False) -dwI.PrintPotentials() - -sm_I, wfu_I = dwI.CalScatteringMatrix() - -dwI.PrintScatteringMatrix() - -# dwI.PlotDistortedWave(1, 1, 20) -# dwI.PlotScatteringMatrix() - -#================= outgoing wave function -print("====================== Outgoing wave function ") -Koning(A_B, Z_B, A_a*ELabPreU + Q_value - Ex, Z_b) - -dwO = DistortedWave(nu_B, nu_b, ELabPreU * A_a + Q_value - Ex) -dwO.maxL = maxL -dwO.ClearPotential() -dwO.AddPotential(WoodsSaxonPot(-v, r0, a), False) -dwO.AddPotential(WoodsSaxonPot(-1j*vi, ri0, ai), False) -dwO.AddPotential(WS_SurfacePot(-1j*vsi, rsi0, asi), False) -dwO.AddPotential(SpinOrbit_Pot(-vso , rso0, aso), False) -dwO.AddPotential(SpinOrbit_Pot(- 1j* vsoi, rsoi0, asoi), False) -dwO.AddPotential(CoulombPotential(rc0), False) -dwO.PrintPotentials() - -sm_O, wfu_O = dwO.CalScatteringMatrix() - -dwO.PrintScatteringMatrix() - -# dwO.PlotDistortedWave(1, 1.5, 20) - -end_time = time.time() # End the timer -print(f"Time used {(end_time - start_time) * 1000:.2f} milliseconds") - -#=================== Calculate radial integral -print("====================== Calculating Radial integrals") - -def FormatSpin(spin : float) -> str: - if int(2*spin) % 2 == 0 : - return f"{int(spin):+d}" - else: - return f"{int(2*spin):+d}/2" - -spin_a = dwI.spin_a -spin_b = dwO.spin_a - -radialInt = np.zeros((maxL+1, int(2*spin_a+1), int(2*spin_b+1)), dtype=complex) - -bs = boundState.GetBoundStateWF() - -for L in range(0, maxL+1): - for index1 in range(0, len(wfu_I[L])): - wf1 = wfu_I[L][index1] - for index2 in range(0, len(wfu_O[L])): - wf2 = wfu_O[L][index2] - # if L == 0 and index1 == 2 and index2 == 1 : - # for i in range(0, len(bs)): - # if i%50 == 0 : - # print(bs[i], wf1[i], wf2[i], bs[i]* wf1[i]* wf2[i]) - pf1 = np.exp(1j*dwI.CoulombPhaseShift(L)) - pf2 = np.exp(1j*dwI.CoulombPhaseShift(L)) - integral = simpson (bs*wf1*wf2, dx=boundState.dr) - radialInt[L][index1][index2] = integral * pf1 * pf2 - -#print radial integral -for index1 in range(0, int(2*spin_a) + 1): - for index2 in range(0, int(2*spin_b) + 1): - print(f"======================= J1 = L{FormatSpin(index1-spin_a)}, J2 = L{FormatSpin(index2-spin_b)}") - for L in range(0, maxL+1): - J1 = L + index1 - spin_a - J2 = L + index2 - spin_b - print(f"{L:2d}, {J1:4.1f}, {J2:4.1f}, {np.real(radialInt[L][index1][index2]):12.4e} + {np.imag(radialInt[L][index1][index2]):12.4e}I") - -# Plot radial integral -fig, axes = plt.subplots(int(2*spin_b+1), int(2*spin_a+1), figsize=(6*int(2*spin_a+1), 4*int(2*spin_b+1))) - -for index2 in range(0, int(2*spin_b) + 1): - for index1 in range(0, int(2*spin_a) + 1): - haha = [] - l_list = [] - for L in range(0, maxL+1): - J1 = L + index1 - spin_a - J2 = L + index2 - spin_b - if J1 < 0 or J2 < 0 : - continue - l_list.append(L) - haha.append(radialInt[L][index1][index2]) - axes[index2, index1].plot(l_list, np.real(haha), label="Real", marker='o') - axes[index2, index1].plot(l_list, np.imag(haha), label="Imag", marker='x') - axes[index2, index1].legend() - axes[index2, index1].set_xlabel('L') - axes[index2, index1].set_ylabel('Value') - axes[index2, index1].set_title(f'Radial Int. vs L for Spin J1 = L{FormatSpin(index1-spin_a)}, J2 = L{FormatSpin(index2-spin_b)}.') - axes[index2, index1].set_xlim(-1, maxL+1) - axes[index2, index1].grid() - -plt.tight_layout() +plt.plot(angList, xsec) +plt.xlim(-1, 181) +plt.yscale("log") +plt.grid() plt.show(block=False) + input("Press Enter to continue...") - - -def Gamma(L1, J1, L2, J2, m, ma, mb): - if int(L1 + L2 + l)%2 != 0: - return 0 - else: - fact0 = wigner_9j(j, l, s, J1, L1, spin_a, J2, L2, spin_b) - if fact0 == 0: - return 0 - else: - fact1 = pow(-1, m) * np.power(1j, L1-L2-l) * (2*L2+1) * np.sqrt((2*l+1)*(2*s+1)*(2*L1+1)*(2*J2+1)) - fact2 = np.sqrt( quantum_factorial(L2-m) / quantum_factorial(L2+m) ) - fact3 = clebsch_gordan(J2, mb-m, j, m-mb+ma, J1, ma) - fact4 = clebsch_gordan(L1, 0, spin_a, ma, J1, ma) - fact5 = clebsch_gordan(L2, -m, spin_b, mb, J2, mb-m) - fact6 = clebsch_gordan(L1, 0, l, 0, L2, 0) - return fact0 * fact1 * fact2 * fact3 * fact4 * fact5 * fact6 - - -def Beta(m, ma, mb, theta_deg): - return 0 - - diff --git a/Raphael/boundState.py b/Raphael/boundState.py index 9ff1ff8..22db364 100755 --- a/Raphael/boundState.py +++ b/Raphael/boundState.py @@ -39,6 +39,7 @@ class BoundState(SolvingSE): def FindPotentialDepth(self, Vmin, Vmax, Vstep=1, isPathWhittaker = True): start_time = time.time() # Start the timer + print(f"Potential Depth Search from {Vmin:.2f} to {Vmax:.2f} MeV, step {Vstep:.2f}") V0List = np.arange(Vmin, Vmax, Vstep) lastSolU = [] minLastSolU = 0 diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py index c03e623..4a9ff30 100755 --- a/Raphael/distortedWave.py +++ b/Raphael/distortedWave.py @@ -53,16 +53,18 @@ class DistortedWave(SolvingSE): self.ScatMatrix = [] self.distortedWaveU = [] + tempZeroList = np.zeros(len(self.rpos), dtype=np.complex128) + for L in range(0, maxL+1): - sigma = self.CoulombPhaseShift() + # sigma = self.CoulombPhaseShift() temp_ScatMatrix = [] temp_distortedWaveU = [] for J in np.arange(L-self.S, L + self.S+1, 1): - if J < 0: + if J < 0 or (L==0 and J != self.S): temp_ScatMatrix.append(0) - temp_distortedWaveU.append(0) + temp_distortedWaveU.append(tempZeroList) continue self.SetLJ(L, J) diff --git a/Raphael/dwba_zr.py b/Raphael/dwba_zr.py new file mode 100755 index 0000000..673859d --- /dev/null +++ b/Raphael/dwba_zr.py @@ -0,0 +1,372 @@ +#!/usr/bin/env python3 +import sys, os +import re +import numpy as np +from scipy.integrate import simpson +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt +import time +from sympy import S +from sympy.physics.quantum.cg import wigner_9j + +sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) +from IAEANuclearData import IsotopeClass + +from assLegendreP import associated_legendre_array +from clebschGordan import clebsch_gordan, quantum_factorial, obeys_triangle_rule +from boundState import BoundState +from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot +from distortedWave import DistortedWave +import opticalPotentials as op + +class DWBA_ZR: + def __init__(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float): + start_time = time.time() + iso = IsotopeClass() + + A_A, Z_A = iso.GetAZ(nu_A) + A_a, Z_a = iso.GetAZ(nu_a) + A_b, Z_b = iso.GetAZ(nu_b) + A_B, Z_B = iso.GetAZ(nu_B) + + self.ELab = A_a * ELabPerU + + mass_A = iso.GetMassFromSym(nu_A) + mass_a = iso.GetMassFromSym(nu_a) + mass_b = iso.GetMassFromSym(nu_b) + mass_B = iso.GetMassFromSym(nu_B) + self.ExB = ExB + + # sym_A = iso.GetSymbol(A_A, Z_A) + # sym_B = iso.GetSymbol(A_B, Z_B) + + spin_A_str = iso.GetJpi(A_A, Z_A) + self.spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str))) + self.spin_B = float(eval(re.sub(r'[+-]', '', JB))) + + if A_a == 2 and Z_a == 1: + self.spin_a = 1.0 + self.spin_b = 0.5 + else: + self.spin_a = 0.5 + self.spin_b = 1.0 + + #====== transfering nucleon + self.s = 1/2 # spin of x, neutron or proton + A_x = abs(A_a - A_b) + Z_x = abs(Z_a - Z_b) + + mass_x = iso.GetMassFromAZ( A_x, Z_x) + + #======== core + if A_A < A_B : # (d,p) + A_c = A_A + Z_c = Z_A + BindingEnergy = mass_B - mass_A - mass_x + self.ExB + else: #(p,d) + A_c = A_B + Z_c = Z_B + BindingEnergy = mass_A - mass_B - self.ExB - mass_x + + #=================== digest orbital + match = re.search(r'[a-zA-Z]', orbital) # Find first letter + if match: + index = match.start() # Get position of the first letter + + node = int(orbital[:index]) + l_sym = orbital[index:index+1] + j_sym = orbital[index+1:] + self.j = eval(j_sym) + self.l = op.ConvertLSym(l_sym) + + passJ = False + if obeys_triangle_rule(self.spin_A, self.spin_B, self.j): + passJ = True + else: + print(f"the orbital spin-J ({self.j}) does not consver J({nu_A}) + J({nu_B}) = {self.spin_A} + {self.spin_B}.") + + passS = False + if obeys_triangle_rule(self.spin_a, self.spin_b, self.s): + passS = True + else: + print(f"the orbital spin-s ({self.s}) does not consver S({nu_a}) + S({nu_b}) = {self.spin_a} + {self.spin_b}.") + + passL = False + if obeys_triangle_rule(self.j, self.s, self.l): + passL = True + else: + print(f"the orbital spin-L ({self.l}) does not consver J({self.j}) + S({self.s}).") + + self.isSpinBalanced = passJ * passS * passL + if self.isSpinBalanced == False : + print("Fail angular momentum conservation.") + return + else: + print("All Spin are balance.") + + self.reactionStr = f"{nu_A}({spin_A_str})({nu_a},{nu_b}){nu_B}({ExB:.3f}|{JB}, {orbital}) @ {ELabPerU:.1f} MeV/u" + print("==================================================") + print(self.reactionStr) + + self.Q_value = mass_A + mass_a - mass_b - mass_B - ExB + print(f"Transfer Orbtial : {orbital}") + print(f"Q-value : {self.Q_value:10.6f} MeV") + print(f"Binding : {BindingEnergy:10.6f} MeV") + + print("====================== Bound state ") + self.boundState = BoundState(A_c, Z_c, A_x, Z_x, node, self.l, self.j, BindingEnergy) + self.boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.30) + self.boundState.FindPotentialDepth(-70, -45, 0.5) + + 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) + self.dwI.AddPotential(WoodsSaxonPot(-1j*op.vi, op.ri0, op.ai), False) + self.dwI.AddPotential(WS_SurfacePot(-1j*op.vsi, op.rsi0, op.asi), False) + self.dwI.AddPotential(SpinOrbit_Pot( -op.vso, op.rso0, op.aso), False) + self.dwI.AddPotential(SpinOrbit_Pot(-1j*op.vsoi, op.rsoi0, op.asoi), False) + self.dwI.AddPotential(CoulombPotential( op.rc0), False) + self.dwI.PrintPotentials() + + 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 + + sm_I, wfu_I = self.dwI.CalScatteringMatrix() + + self.dwI.PrintScatteringMatrix() + + Ecm_I = self.dwI.Ecm + Ecm_O = Ecm_I + self.Q_value + Eout = ((Ecm_O + mass_b + mass_B + self.ExB)**2 - (mass_b + mass_B + ExB)**2)/2/mass_B + + print("====================== Outgoing wave function ") + op.Koning(A_B, Z_B, self.ELab + self.Q_value - ExB, Z_b) + + self.dwO = DistortedWave(nu_B, nu_b, Eout) + self.dwO.spin_A = self.spin_B + self.dwO.maxL = self.maxL + self.dwO.PrintInput() + self.dwO.ClearPotential() + self.dwO.AddPotential(WoodsSaxonPot( -op.v, op.r0, op.a), False) + self.dwO.AddPotential(WoodsSaxonPot(-1j*op.vi, op.ri0, op.ai), False) + self.dwO.AddPotential(WS_SurfacePot(-1j*op.vsi, op.rsi0, op.asi), False) + self.dwO.AddPotential(SpinOrbit_Pot( -op.vso, op.rso0, op.aso), False) + self.dwO.AddPotential(SpinOrbit_Pot(-1j*op.vsoi, op.rsoi0, op.asoi), False) + self.dwO.AddPotential(CoulombPotential( op.rc0), False) + + self.dwO.PrintPotentials() + + sm_O, wfu_O_temp = self.dwO.CalScatteringMatrix() + + #============ rescale the outgoing wave + rpos_O_temp = self.dwO.rpos * A_B/A_A + self.rpos_O = [] + rpos_O_filled = False + self.wfu_O = [] + for L in range(0, self.maxL+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') + temp_wfu = [] + for r in self.dwI.rpos: + if r > 20 : + break + if rpos_O_filled == False: + self.rpos_O.append(r) + temp_wfu.append(wfu_O_inter_real(r) + 1j * wfu_O_inter_imag(r)) + 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) + + bs = self.boundState.GetBoundStateWF() + + for L1 in range(0, self.maxL+1): + for index1 in range(0, len(wfu_I[L1])): + 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)])): + 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 + + mass_I = self.dwI.mu + k_I = self.dwI.k + mass_O = self.dwO.mu # reduced mass of outgoing channel + k_O = self.dwO.k # wave number of outgoing channel + + D0 = 1.55e+4 # for (d,p) + + self.massFactor = A_B/A_A + self.ffactor = np.sqrt(4*np.pi)/k_I /k_O * A_B/A_A + + 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) + + stop_time = time.time() + print(f"Total time {(stop_time - start_time) * 1000:.2f} msec") + + #========== end of contructor + + def FormatSpin(self, spin : float) -> str: + if int(2*spin) % 2 == 0 : + return f"{int(spin):+d}" + 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 PlotRadialIntegral(self): + spin_b = self.spin_b + spin_a = self.spin_a + l = int(self.l) + maxL = self.maxL + + 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.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') + axes[int(2*l+1)*index2 + indexL2, index1].plot(l_list, np.imag(haha), label="Imag", marker='x') + 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].grid() + + plt.tight_layout() + 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 Gamma(self, L1:int, J1, L2:int, J2, m, 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, S(2*self.l)/2, S(2*self.s)/2, S(2*J1)/2, S(2*L1)/2, S(2*self.spin_a)/2, S(2*J2)/2, S(2*L2)/2, S(2*self.spin_b)/2).evalf() + 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) ) + 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, int(self.l), 0, L2, 0) + return fact0 * fact1 * fact2 * fact3 * fact4 * fact5 * fact6 + + def Beta(self, m:int, ma, mb): + 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): + + if not obeys_triangle_rule(J1, self.j, J2): + continue + if not(abs(m) <= L2): + continue + 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.legendrePArray[int(L2)][int(abs(m))] + if m < 0 : + lp *= (-1)**m * quantum_factorial(int(L2)+m)/ quantum_factorial(int(L2)-m) + 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}") + + result += gg * lp * ri + + return result + + def PreCalLegendreP(self, theta_deg:float, maxL:int = None, maxM:int = None): + if maxL is None: + maxL = self.maxL + if maxM is None: + maxM = int(self.j + self.spin_b + self.spin_a) + self.legendrePArray = associated_legendre_array(maxL, maxM, theta_deg) + + def AngDist(self, theta_deg): + xsec = 0 + self.PreCalLegendreP(theta_deg) + for ma in np.arange(-self.spin_a, self.spin_a + 1, 1): + for mb in np.arange(-self.spin_b, self.spin_b + 1, 1): + for m in np.arange(-self.j + mb - ma, self.j + mb -ma + 1, 1): + haha = self.Beta(m, ma, mb) + xsec += np.abs(haha)**2 + + return xsec * self.xsecScalingfactor * 10 # factor 10 for fm^2 = 10 mb + + + + + + + diff --git a/Raphael/opticalPotentials.py b/Raphael/opticalPotentials.py new file mode 100755 index 0000000..d159c9c --- /dev/null +++ b/Raphael/opticalPotentials.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python3 + +#need to use import + +import numpy as np + +# Woods-Saxon +v = 0 +r0 = 0 +a = 0 +vi = 0 +ri0 = 0 +ai = 0 +# Woods-Saxon Surface +vsi = 0 +rsi0 = 0 +asi = 0 +# Spin-orbit +vso = 0 +rso0 = 0 +aso = 0 +vsoi = 0 +rsoi0 = 0 +asoi = 0 +# Coulomb +rc0 = 0 + +def AnCai(A : int, Z : int, E : float): + global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0 + + A3 = A**(1./3.) + v = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3 + r0 = 1.152 - 0.00776 / A3 + a = 0.719 + 0.0126 * A3 + + vi = 1.104 + 0.0622 * E + ri0 = 1.305 + 0.0997 / A3 + ai = 0.855 - 0.1 * A3 + + vsi = 10.83 - 0.0306 * E + rsi0 = 1.334 + 0.152 / A3 + asi = 0.531 + 0.062 * A3 + + vso = 3.557 + rso0 = 0.972 + aso = 1.011 + + vsoi = 0.0 + rsoi0 = 0.0 + asoi = 0.0 + + rc0 = 1.303 + +def Koning(A : int, Z : int, E : float, Zproj : float): + global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0 + + N = A-Z + A3 = A**(1./3.) + + vp1 = 59.3 + 21.*(N-Z)/A - 0.024*A + vn1 = 59.3 - 21.*(N-Z)/A - 0.024*A + + vp2 = 0.007067 + 0.00000423*A + vn2 = 0.007228 - 0.00000148*A + + vp3 = 0.00001729 + 0.00000001136 * A + vn3 = 0.00001994 - 0.00000002 * A + + vp4 = 7e-9 # = vn4 + vn4 = vp4 + + wp1 = 14.667 + 0.009629*A + wn1 = 12.195 + 0.0167*A + + wp2 = 73.55 + 0.0795*A # = wn2 + wn2 = wp2 + + dp1 = 16 + 16.*(N-Z)/A + dn1 = 16 - 16.*(N-Z)/A + + dp2 = 0.018 + 0.003802/(1 + np.exp((A-156.)/8)) # = dn2 + dn2 = dp2 + + dp3 = 11.5 # = dn3 + dn3 = dp3 + + vso1 = 5.922 + 0.003 * A + vso2 = 0.004 + + wso1 = -3.1 + wso2 = 160 + + epf = -8.4075 + 0.01378 *A + enf = -11.2814 + 0.02646 *A + + rc = 1.198 + 0.697/pow(A3,2) + 12.995/pow(A3,5) + vc = 1.73/rc * Z / A3 + + v = vp1*(1 - vp2*(E-epf) + vp3*pow(E-epf,2) - vp4*pow(E-epf,3)) + vc * vp1 * (vp2 - 2*vp3*(E-epf) + 3*vp4*pow(E-epf,2)) + #neutron + if Zproj == 0 : + v = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3)) + + r0 = 1.3039 - 0.4054 / A3 + a = 0.6778 - 0.000148 * A + + vi = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2)) + if Zproj == 0 : + vi = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2)) + + ri0 = 1.3039 - 0.4054 / A3 + ai = 0.6778 - 0.000148 * A + + vsi = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf)) + if Zproj == 0 : + vsi = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf)) + + rsi0 = 1.3424 - 0.01585 * A3 + asi = 0.5187 + 0.0005205 * A + if Zproj == 0: + asi = 0.5446 - 0.0001656 * A + + vso = vso1 * np.exp(-vso2 * (E-epf)) + if Zproj == 0: + vso = vso1 * np.exp(-vso2 * (E-enf)) + + rso0 = 1.1854 - 0.647/A3 + aso = 0.59 + + vsoi = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2)) + if Zproj == 0 : + vsoi = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2)) + + rsoi0 = 1.1854 - 0.647/A3 + asoi = 0.59 + + rc0 = rc + +def ConvertLSym(LSym :str) -> int: + if LSym == "s" : + return 0 + elif LSym == "p" : + return 1 + elif LSym == "d" : + return 2 + elif LSym == "f" : + return 3 + elif LSym == "g" : + return 4 + elif LSym == "h" : + return 5 + elif LSym == "i" : + return 6 + elif LSym == "j" : + return 7 + elif LSym == "k" : + return 8 + else : + return -1 diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index 1a782af..b39bc93 100755 --- a/Raphael/solveSE.py +++ b/Raphael/solveSE.py @@ -197,7 +197,7 @@ class SolvingSE: self.eta = 0 else: self.eta = self.Z * self.ee * self.k /2 /self.Ecm - self.maxL = int(self.k * (1.4 * (self.A_A**(1/3) + self.A_a**(1/3)) + 5)) + self.maxL = int(self.k * (1.4 * (self.A_A**(1/3) + self.A_a**(1/3)) + 10)) # def SetA_ExSpin(self, ExA, sA ): # self.ExA = ExA