From 3b67db12e247fa8e7a735bb81fd075c35235a830 Mon Sep 17 00:00:00 2001 From: "Ryan@Home" Date: Sat, 22 Feb 2025 20:26:56 -0500 Subject: [PATCH] snapshot, finial step for DWBA_ZR --- .gitignore | 3 + Cleopatra/IAEANuclearData.py | 2 +- PyGUIQt6/DWBA | 4 +- Raphael/DWBA_ZR.py | 187 +++++++++++++++++++++++++++++++---- Raphael/boundState.py | 9 ++ Raphael/clebschGordan.py | 125 +++++++++++++++++++++++ Raphael/distortedWave.py | 39 ++++---- Raphael/solveSE.py | 146 ++++++++++++++++----------- 8 files changed, 418 insertions(+), 97 deletions(-) diff --git a/.gitignore b/.gitignore index 708c0b2..9bd0176 100644 --- a/.gitignore +++ b/.gitignore @@ -14,5 +14,8 @@ IAEA_NuclearData.csv *.in *.out +*.prof + +test.py DWUCK4 diff --git a/Cleopatra/IAEANuclearData.py b/Cleopatra/IAEANuclearData.py index 59ae5e7..a09100f 100644 --- a/Cleopatra/IAEANuclearData.py +++ b/Cleopatra/IAEANuclearData.py @@ -136,7 +136,7 @@ class IsotopeClass: except: return "unknown" - def GetJpi(self, A : int, Z : int): + def GetJpi(self, A : int, Z : int) -> str: try: dudu = self.data[(self.data['z']==Z) & (self.data['A']==A)] return dudu['jp'].iloc[0] diff --git a/PyGUIQt6/DWBA b/PyGUIQt6/DWBA index c573061..9bc925f 100644 --- a/PyGUIQt6/DWBA +++ b/PyGUIQt6/DWBA @@ -43,4 +43,6 @@ #10Be(t,p)12Be 0 1L=0 0+ 0.000 5MeV/u lA #two-nucleon_transfer #32Si(t,p)34Si 0 0L=0 0+ 0.000 8MeV/u lA #two-nucleon_transfer #133Sb(t,3He)133Sn 7/2 0g7/2 0+ 0.000 8.5MeV/u Ax .... cannot cal -#12C(6Li,d)16O 0 nL=2 2+ 0.0 10MeV/u 6K \ No newline at end of file +#12C(6Li,d)16O 0 nL=2 2+ 0.0 10MeV/u 6K + +16O(d,p)17O 0 1s1/2 1/2+ 0.87 10MeV/u AK \ No newline at end of file diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py index 2e11ec3..7051075 100755 --- a/Raphael/DWBA_ZR.py +++ b/Raphael/DWBA_ZR.py @@ -1,13 +1,13 @@ #!/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.SetPotential(1.25, 0.65, -6, 1.10, 0.65, 1.30) -# boundState.FindPotentialDepth(-75, -40, 0.1) +# 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() @@ -15,12 +15,12 @@ import matplotlib.pyplot as plt 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", "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() @@ -32,12 +32,12 @@ dw.AddPotential(CoulombPotential(1.258), False) # dw.AddPotential(CoulombPotential(1.303), False) -dw.CalScatteringMatrix() +# dw.CalScatteringMatrix() # dw.PrintScatteringMatrix() -dw.PlotDCSUnpolarized(180, 1) +# dw.PlotDCSUnpolarized(180, 1) -exit() +# exit() # for i in range(1, 19): # theta = 10*i @@ -54,10 +54,15 @@ exit() 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 @@ -220,9 +225,12 @@ nu_b = "p" nu_B = "17O" ELabPreU = 10 # MeV/u Ex = 0.87 -J_B = 0.5 +J_B = "1/2+" orbital = "1s1/2" +import time +start_time = time.time() # Start the timer + iso = IsotopeClass() A_A, Z_A = iso.GetAZ(nu_A) @@ -252,6 +260,12 @@ else: #(p,d) 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) +spin_B_str = J_B + +spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str))) +spin_B = float(eval(re.sub(r'[+-]', '', J_B))) + if A_a == 2 and Z_a == 1: spin_a = 1.0 spin_b = 0.5 @@ -259,10 +273,7 @@ else: spin_a = 0.5 spin_b = 1.0 -Q_value = mass_A + mass_a - mass_b - mass_B - Ex - -print(f"Q-value : {Q_value:10.6f} MeV") -print(f"Binding : {BindingEnergy:10.6f} MeV") +s = 0.5 # spin of x, neutron or proton #=================== digest orbital match = re.search(r'[a-zA-Z]', orbital) # Find first letter @@ -275,6 +286,39 @@ j_sym = orbital[index+1:] j = eval(j_sym) l = ConvertLSym(l_sym) +#==== 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}.") + +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 @@ -284,17 +328,22 @@ 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.25, 0.65, -6, 1.10, 0.65, 1.30) +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) @@ -302,12 +351,17 @@ 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) @@ -319,7 +373,102 @@ 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() \ No newline at end of file +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.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 3069309..9ff1ff8 100755 --- a/Raphael/boundState.py +++ b/Raphael/boundState.py @@ -9,6 +9,7 @@ import numpy as np import matplotlib.pyplot as plt from mpmath import whitw import mpmath +import time mpmath.mp.dps = 15 # Decimal places of precision @@ -21,6 +22,7 @@ class BoundState(SolvingSE): self.PrintInput() self.node = node # number of nodes of the wave function r > 0 self.FoundBounfState = False + self.wf = None def SetPotential(self, r0, a0, Vso, rso, aso, rc = 0.0): self.r0 = r0 @@ -36,6 +38,7 @@ class BoundState(SolvingSE): self.AddPotential(CoulombPotential(rc), False) # not use mass number of a def FindPotentialDepth(self, Vmin, Vmax, Vstep=1, isPathWhittaker = True): + start_time = time.time() # Start the timer V0List = np.arange(Vmin, Vmax, Vstep) lastSolU = [] minLastSolU = 0 @@ -121,6 +124,12 @@ class BoundState(SolvingSE): print(f"ANC : {self.ANC:10.6e}") self.wf[rIndex:] = self.ANC * np.array(W_values[rIndex:]) + end_time = time.time() # End the timer + print(f"Finding Potential Depth and Bound state took {(end_time - start_time) * 1000:.2f} milliseconds") + + def GetBoundStateWF(self): + return self.wf + def PlotBoundState(self, maker=None): if not self.FoundBounfState: plt.plot(self.rpos[1:], self.solU[1:]/self.rpos[1:]) diff --git a/Raphael/clebschGordan.py b/Raphael/clebschGordan.py index 8b000e2..ec5927b 100755 --- a/Raphael/clebschGordan.py +++ b/Raphael/clebschGordan.py @@ -13,6 +13,25 @@ from scipy.special import gamma import numpy as np from math import sqrt +def KroneckerDelta(i, j): + if i == j: + return 1 + else: + return 0 + +def obeys_triangle_rule(j1, j2, j3): + """Check if j1, j2, j3 obey the vector summation rules.""" + # Ensure non-negativity (optional if inputs are guaranteed positive) + if j1 < 0 or j2 < 0 or j3 < 0: + return False + # Triangle inequalities + if (j3 < abs(j1 - j2) or j3 > j1 + j2): + return False + # Check if j1 + j2 + j3 is an integer (for half-integer j, this is automatic) + if (j1 + j2 + j3) % 1 != 0: + return False + return True + def quantum_factorial(n): """ Calculate factorial for integer or half-integer numbers using gamma function. @@ -87,3 +106,109 @@ def clebsch_gordan(j1, m1,j2, m2, j, m): return prefactor * sum_result + +#============ don;t use, very slow, use the sympy package +def threej(j1, m1, j2, m2, j3, m3): + if m1 + m2 + m3 != 0: + return 0 + if obeys_triangle_rule(j1, j2, j3) == False: + return 0 + + cg = clebsch_gordan(j1, m1, j2, m2, j3, -m3) + norm = pow(-1, j1-j2-m3)/(2*j3+1)**0.5 + return norm * cg + +def sixj(j1, j2, j3, j4, j5, j6): + """Compute the 6j symbol using Clebsch-Gordan coefficients.""" + # Check triangle conditions + if not (obeys_triangle_rule(j1, j2, j3) and + obeys_triangle_rule(j1, j5, j6) and + obeys_triangle_rule(j4, j2, j6) and + obeys_triangle_rule(j4, j5, j3)): + return 0.0 + + sixj_value = 0.0 + + # Ranges for m values + m1_range = range(-j1, j1 + 1) + m2_range = range(-j2, j2 + 1) + m4_range = range(-j4, j4 + 1) + m5_range = range(-j5, j5 + 1) + + # Sum over m values + for m1 in m1_range: + for m2 in m2_range: + m3 = - m1 - m2 + for m4 in m4_range: + for m5 in m5_range: + m6 = m2 + m4 + + if m3 + m5 not in m4_range or m1 + m6 not in m5_range: + continue + + # cg1 = threej(j1, -m1, j2, -m2, j3, -m3) + + cg1 = (-1)**(j1-j2+m3) * clebsch_gordan(j1, -m1, j2, -m2, j3, m3) / (2*j3+1)**0.5 + + cg2 = threej(j1, m1, j5, -m5, j6, m6) + cg3 = threej(j4, m4, j2, m2, j6, -m6) + cg4 = threej(j4, -m4, j5, m5, j3, m3) + + norm = pow(-1, j1-m1 + j2-m2 + j3-m3 + j4-m4 + j5-m5 + j6-m6) + + sixj_value += cg1 * cg2 * cg3 * cg4 * norm + + return sixj_value + +def ninej(j1, j2, j3, j4, j5, j6, j7, j8, j9): + """Compute the 9j symbol using 6j symbols.""" + # Check triangle conditions for rows + if not (obeys_triangle_rule(j1, j2, j3) and + obeys_triangle_rule(j4, j5, j6) and + obeys_triangle_rule(j7, j8, j9)): + return 0.0 + + # Check triangle conditions for columns + if not (obeys_triangle_rule(j1, j4, j7) and + obeys_triangle_rule(j2, j5, j8) and + obeys_triangle_rule(j3, j6, j9)): + return 0.0 + + ninej_value = 0.0 + + # Determine the range of intermediate angular momentum x + x_min = max(abs(j1 - j9), abs(j4 - j8), abs(j2 - j6)) + x_max = min(j1 + j9, j4 + j8, j2 + j6) + + # Sum over x (must be integer or half-integer depending on inputs) + step = 1 if all(j % 1 == 0 for j in [j1, j2, j3, j4, j5, j6, j7, j8, j9]) else 0.5 + for x in [x_min + i * step for i in range(int((x_max - x_min) / step) + 1)]: + # if not (obeys_triangle_rule(j1, j4, j7) and + # obeys_triangle_rule(j1, j9, x) and # j1 j9 + # obeys_triangle_rule(j8, j9, j7) and + # obeys_triangle_rule(j8, j4, x) and # j8 j4 + # obeys_triangle_rule(j2, j5, j8) and + # obeys_triangle_rule(j2, x, j6) and # j2 j6 + # obeys_triangle_rule(j4, j5, j6) and + # obeys_triangle_rule(j4, x, j8) and # j4 j8 + # obeys_triangle_rule(j3, j6, j9) and + # obeys_triangle_rule(j3, j1, j2) and + # obeys_triangle_rule( x, j6, j2) and # j2 j6 + # obeys_triangle_rule( x, j1, j9)): # j1 j9 + # continue + + if not (obeys_triangle_rule(j1, j9, x) and # j1 j9 + obeys_triangle_rule(j8, j4, x) and # j8 j4 + obeys_triangle_rule(j2, x, j6)): # j1 j9 + continue + + sixj1 = sixj(j1, j4, j7, j8, j9, x) + sixj2 = sixj(j2, j5, j8, j4, x, j6) + sixj3 = sixj(j3, j6, j9, x, j1, j2) + + phase = (-1) ** int(2 * x) # Phase factor + weight = 2 * x + 1 # Degeneracy factor + + ninej_value += phase * weight * sixj1 * sixj2 * sixj3 + + return ninej_value \ No newline at end of file diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py index 44a74ed..c03e623 100755 --- a/Raphael/distortedWave.py +++ b/Raphael/distortedWave.py @@ -17,20 +17,8 @@ def SevenPointsSlope(data, n): def FivePointsSlope(data, n): return ( data[n + 2] - 8 * data[n + 1] + 8 * data[n - 1] - data[n - 2] ) / 12 -# from sympy.physics.quantum.cg import CG -# from sympy import S -# def clebsch_gordan(j1, m1, j2, m2, j, m): -# cg = CG(S(j1), S(m1), S(j2), S(m2), S(j), S(m)) -# result = cg.doit() -# return np.complex128(result) - -def KroneckerDelta(i, j): - if i == j: - return 1 - else: - return 0 - -from clebschGordan import clebsch_gordan +from clebschGordan import clebsch_gordan, KroneckerDelta +import time ############################################################ class DistortedWave(SolvingSE): @@ -57,6 +45,8 @@ class DistortedWave(SolvingSE): return np.angle(gamma(L+1+1j*eta)) def CalScatteringMatrix(self, maxL = None, verbose = False): + start_time = time.time() # Start the timer + if maxL is None: maxL = self.maxL @@ -101,16 +91,28 @@ class DistortedWave(SolvingSE): temp_ScatMatrix.append(ScatMatrix) dwU = np.array(self.solU, dtype=np.complex128) - dwU *= np.exp(-1j*sigma)/(B-A*1j) + #dwU *= np.exp(-1j*sigma)/(B-A*1j) + dwU *= 1./(B-A*1j) temp_distortedWaveU.append(dwU) self.ScatMatrix.append(temp_ScatMatrix) self.distortedWaveU.append(temp_distortedWaveU) + end_time = time.time() # End the timer + print(f"Calculate Scattering Matrixes took {(end_time - start_time) * 1000:.2f} milliseconds") + return [self.ScatMatrix, self.distortedWaveU] def PrintScatteringMatrix(self): + print("======================= Scattering Matrix") for L in range(0, len(self.ScatMatrix)): + + if L == 0 : + print(" ", end="") + for i in range(0, len(self.ScatMatrix[L])): + print(f"{{{'L':>2s},{'J':>4s}, {'Real':>10s} + {'Imaginary':>10s}}}, ", end="") + print("") + print("{", end="") for i in range(0, len(self.ScatMatrix[L])): print("{", end="") @@ -125,11 +127,14 @@ class DistortedWave(SolvingSE): return self.ScatMatrix[L][J-L+self.S] def GetDistortedWave(self, L, J): - return self.distortedWaveU[L][J-L+self.S] + return self.distortedWaveU[L][int(J-L+self.S)] - def PlotDistortedWave(self, L, J): + def PlotDistortedWave(self, L, J, maxR = None): plt.plot(self.rpos, np.real(self.GetDistortedWave(L, J)), label="Real") plt.plot(self.rpos, np.imag(self.GetDistortedWave(L, J)), label="Imaginary") + 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) diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index f9c86c1..1a782af 100755 --- a/Raphael/solveSE.py +++ b/Raphael/solveSE.py @@ -8,79 +8,100 @@ import sys, os sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) from IAEANuclearData import IsotopeClass -class CoulombPotential: - def __init__(self, rc): - self.rc = rc - self.id = 0 - self.ee = 1.43996 # MeV.fm - - def setA(self, A): - self.Rc = self.rc * math.pow(A, 1/3) - - def setAa(self, A, a): - self.Rc = self.rc * (math.pow(A, 1/3) + math.pow(a, 1/3)) +class PotentialForm: + def __init__(self): + self.V0 = -10 + self.r0 = 1.3 + self.a0 = 0.65 + self.id = -1 def setCharge(self, Z): self.Charge = Z - def output(self, x): - if x >self.Rc: - return (self.Charge * self.ee) / (x + 1e-20) # Add a small value to avoid division by zero - else: - return (self.Charge * self.ee) / (2 * self.Rc) * (3 - (x / self.Rc)**2) + # def setA(self, A): + # self.R0 = self.r0 * math.pow(A, 1/3) -class WoodsSaxonPot: + def setAa(self, A, a): + self.R0= self.r0 * (math.pow(A, 1/3) + math.pow(a, 1/3)) + + def output(self, x): + return 0 + + def printPot(self, msg:str): + print(f"{msg:20s} : V0 {np.real(self.V0):7.3f} + {np.imag(self.V0):7.3f} I, R0 {self.R0:6.4f}({self.r0:6.4f}), a0 {self.a0:6.4f}, ") + +class CoulombPotential(PotentialForm): + def __init__(self, rc): + self.V0 = 0 + self.r0 = rc + self.a0 = 0 + self.id = 0 + self.ee = 1.43996 # MeV.fm + + def output(self, x): + if self.Charge == 0 : + return 0 + else: + if x >self.R0: + return (self.Charge * self.ee) / (x + 1e-20) # Add a small value to avoid division by zero + else: + return (self.Charge * self.ee) / (2 * self.R0) * (3 - (x / self.R0)**2) + + def printPot(self): + return super().printPot("Coulomb") + +class WoodsSaxonPot(PotentialForm): def __init__(self, V0, r0, a0) : self.V0 = V0 self.r0 = r0 self.a0 = a0 self.id = 1 - def setA(self, A): - self.R0 = self.r0 * math.pow(A, 1/3) - - def setAa(self, A, a): - self.R0 = self.r0 * (math.pow(A, 1/3) + math.pow(a, 1/3)) - def output(self, x): - return self.V0/(1 + math.exp((x-self.R0)/self.a0)) + if self.V0 == 0.0: + return 0 + else: + return self.V0/(1 + math.exp((x-self.R0)/self.a0)) -class SpinOrbit_Pot: + def printPot(self): + return super().printPot("Woods-Saxon") + +class SpinOrbit_Pot(PotentialForm): def __init__(self, VSO, rSO, aSO) : # the LS factor is put in the SolvingSE Class - self.VSO = VSO - self.rSO = rSO - self.aSO = aSO + self.V0 = VSO + self.r0 = rSO + self.a0 = aSO self.id = 2 - def setA(self, A): - self.RSO = self.rSO * math.pow(A, 1/3) - - def setAa(self, A, a): - self.RSO = self.rSO * (math.pow(A, 1/3) + math.pow(a, 1/3)) - def output(self, x): - if x > 0 : - return 4*(self.VSO * math.exp((x-self.RSO)/self.aSO))/(self.aSO*math.pow(1+math.exp((x-self.RSO)/self.aSO),2))/x - else : - return 4*1e+19 + if self.V0 == 0.0 : + return 0 + else: + if x > 0 : + return 4*(self.V0 * math.exp((x-self.R0)/self.a0))/(self.a0*math.pow(1+math.exp((x-self.R0)/self.a0),2))/x + else : + return 4*1e+19 -class WS_SurfacePot: + def printPot(self): + return super().printPot("Spin-Orbit") + +class WS_SurfacePot(PotentialForm): def __init__(self, V0, r0, a0): self.V0 = V0 self.r0 = r0 self.a0 = a0 self.id = 3 - def setA(self, A): - self.R0 = self.r0 * math.pow(A, 1/3) - - def setAa(self, A, a): - self.R0 = self.r0 * (math.pow(A, 1/3) + math.pow(a, 1/3)) - def output(self, x): - exponent = (x - self.R0) / self.a0 - return 4* self.V0 * math.exp(exponent) / (1 + math.exp(exponent))**2 + if self.V0 == 0 : + return 0 + else: + exponent = (x - self.R0) / self.a0 + return 4* self.V0 * math.exp(exponent) / (1 + math.exp(exponent))**2 + + def printPot(self): + return super().printPot("Woods-Saxon Surface") #======================================== class SolvingSE: @@ -205,23 +226,30 @@ class SolvingSE: def ClearPotential(self): self.potential_List = [] - def AddPotential(self, pot, useBothMass : bool = False): - if pot.id == 0: - pot.setCharge(self.Z) - if useBothMass: - pot.setAa(self.A_A, self.A_a) - else: - pot.setAa(self.A_A, 0) - self.potential_List.append(pot) + def AddPotential(self, pot : PotentialForm, useBothMass : bool = False): + if isinstance(pot, PotentialForm): + if pot.id == 0: + pot.setCharge(self.Z) + if useBothMass: + pot.setAa(self.A_A, self.A_a) + else: + pot.setAa(self.A_A, 0) + self.potential_List.append(pot) def __PotentialValue(self, x): value = 0 for pot in self.potential_List: - if pot.id == 2 and self.L > 0: - value = value + self.LS() * pot.output(x) - else: - value = value + pot.output(x) + if isinstance(pot, PotentialForm): + if pot.id == 2 and self.L > 0: + value = value + self.LS() * pot.output(x) + else: + value = value + pot.output(x) return value + + def PrintPotentials(self): + for pot in self.potential_List: + if isinstance(pot, PotentialForm): + pot.printPot() def GetPotentialValue(self, x): return self.__PotentialValue(x)