diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py index 92a7417..edcfb1e 100755 --- a/Raphael/DWBA_ZR.py +++ b/Raphael/DWBA_ZR.py @@ -1,15 +1,9 @@ #!/usr/bin/env python3 -from solveSE import WS, Coulomb, SO, WSSurface, SolvingSE - - -boundState = SolvingSE(16, 8, 1, 0, -4.14) -boundState.SetRange(0, 0.1, 300) -boundState.PrintInput() - -boundState.ClearPotential() -boundState.AddPotential(WS(-40, 1.10, 0.65)) - -boundState.SetLJ(0, 0.5) - -print(boundState.GetPotentialValue(0)) \ No newline at end of file +from boundState import BoundState + +boundState = BoundState(16, 8, 1, 0, 0, 2, 2.5, -4.14) +boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65) +boundState.FindPotentialDepth(-70, -60, 0.5) +# boundState.PrintWF() +boundState.PlotBoundState() \ No newline at end of file diff --git a/Raphael/boundState.py b/Raphael/boundState.py new file mode 100755 index 0000000..df20038 --- /dev/null +++ b/Raphael/boundState.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 + +from solveSE import WSPotential, CoulombPotential, SpinOrbitPotential, SolvingSE + +from scipy.optimize import curve_fit +from scipy.interpolate import interp1d +from scipy.integrate import simpson +import numpy as np +import matplotlib.pyplot as plt +from mpmath import whitw +import mpmath + +mpmath.mp.dps = 15 # Decimal places of precision + +def SevenPointsSlope(data, n): + return (-data[n + 3] + 9 * data[n + 2] - 45 * data[n + 1] + 45 * data[n - 1] - 9 * data[n - 2] + data[n - 3]) / 60 + + +class BoundState(SolvingSE): + def __init__(self, A, ZA, a, Za, node, L, J, BE): + super().__init__(A, ZA, a, Za, BE) + self.Ecm = BE + self.BE = BE + + self.SetRange(0, 0.1, 300) # default range + self.PrintInput() + + self.node = node # number of nodes of the wave function r > 0 + self.SetLJ(L, J) + + self.FoundBounfState = False + + def SetPotential(self, r0, a0, Vso, rso, aso, rc = 0.0): + self.r0 = r0 + self.a0 = a0 + self.Vso = Vso + self.rso = rso + self.aso = aso + self.rc = rc + self.ClearPotential() + self.AddPotential(WSPotential(-60, r0, a0), False) # not use mass number of a + self.AddPotential(SpinOrbitPotential(Vso, rso, aso), False) # not use mass number of a + if rc > 0 and self.Z_a > 0: + self.AddPotential(CoulombPotential(rc), False) # not use mass number of a + + def FindPotentialDepth(self, Vmin, Vmax, Vstep=1, isPathWhittaker = True): + V0List = np.arange(Vmin, Vmax, Vstep) + lastSolU = [] + minLastSolU = 0 + maxLastSolU = 0 + for V0 in V0List: + self.ClearPotential() + self.AddPotential(WSPotential(V0, self.r0, self.a0), False) + self.AddPotential(SpinOrbitPotential(self.Vso, self.rso, self.aso), False) # not use mass number of a + if self.rc > 0 and self.Z_a > 0: + self.AddPotential(CoulombPotential(self.rc), False) + wf = self.SolveByRK4() + lastSolU.append(np.real(wf[-1])) + if lastSolU[-1] < minLastSolU: + minLastSolU = lastSolU[-1] + if lastSolU[-1] > maxLastSolU: + maxLastSolU = lastSolU[-1] + + if minLastSolU >= 0 or maxLastSolU <= 0: + print("No bound state found in the range") + print("min Last SolU = ", minLastSolU) + print("max Last SolU = ", maxLastSolU) + plt.plot(V0List, lastSolU, marker="x") + plt.grid() + plt.xlabel('V0 (MeV)') + plt.show(block=False) + input("Press Enter to exit.") + + return + + f = interp1d(lastSolU, V0List, kind='cubic') + self.V_BS = f(0) + print(f"Potential Depth = {self.V_BS:15.11f} MeV") + + # recaluclate the wave function for V_BS + self.ClearPotential() + self.AddPotential(WSPotential(self.V_BS, self.r0, self.a0), False) + self.AddPotential(SpinOrbitPotential(self.Vso, self.rso, self.aso), False) # not use mass number of a + if self.rc > 0 and self.Z_a > 0: + self.AddPotential(CoulombPotential(self.rc), False) + + self.SolveByRK4() + self.SolU = np.real(self.SolU) + + # find number of node in self.SolU with in to potenital range + # find how many time self.SolU change sign + detNode = 0 + for i, r in enumerate(self.rpos): + if r > self.r0 * (pow(self.A_A, 1/3) + pow(self.A_a, 1/3)) + 5 * self.a0: + break + if self.SolU[i] * self.SolU[i-1] < 0: + detNode += 1 + + if detNode != self.node: + print(f"\033[91mNumber of node is not matched, expected : {self.node}, found : {detNode}\033[0m") + return + + self.FoundBounfState = True + + # normalize the wave function + norm = simpson(self.SolU**2) * self.dr + self.SolU /= np.sqrt(norm) + self.wf = np.zeros_like(self.SolU) + self.wf[1:] = self.SolU[1:] / self.rpos[1:] + self.wf[0] = self.SolU[0] # Handle the first element separately if needed + + #extrapolate wf from 0.2, 0.1 to 0.0 + def func(x, a, b, c): + return a * x**2 + b * x + c + popt, pcov = curve_fit(func, self.rpos[1:6], self.wf[1:6]) + self.wf[0] = func(0, *popt) + + if isPathWhittaker: + # patch the long range function to Whittaker function + kappa = np.sqrt(2*self.mu*np.abs(self.BE))/self.hbarc + charge = self.Z + self.eta = kappa * charge * self.ee / 2 / self.BE + + R = min(self.r0 * (pow(self.A_A, 1/3) + pow(self.A_a, 1/3)) + 5 * self.a0, max(self.rpos)) + rIndex = self.NearestPosIndex(R) + print(f"replacing wave function with Whittaker after R : {R:5.1f}, index : {rIndex}") + if rIndex < len(self.rpos): + R = self.rpos[rIndex] + W_values = [float(whitw(-1j*self.eta, self.L + 0.5, 2*kappa * r))/(kappa * r) for r in self.rpos] + + self.ANC = self.wf[rIndex] / W_values[rIndex] + print(f"ANC : {self.ANC:10.6e}") + self.wf[rIndex:] = self.ANC * np.array(W_values[rIndex:]) + + def PlotBoundState(self, maker=None): + if not self.FoundBounfState: + plt.plot(self.rpos[1:], self.SolU[1:]/self.rpos[1:]) + else: + if maker is None: + plt.plot(self.rpos, self.wf) + else: + plt.plot(self.rpos, self.wf, marker=maker) + plt.grid() + plt.xlabel('r (fm)') + plt.ylabel('u(r)/r') + plt.show(block=False) + input("Press Enter to exit.") + + def PrintWF(self, func = None): + if func is None: + for r, wf in zip(self.rpos, self.wf): + print(f"{r:5.1f}, {wf:10.6e}") + else: + for r, wf in zip(self.rpos, func): + print(f"{r:5.1f}, {wf:10.6e})") + diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index 2efd270..1b8aca0 100755 --- a/Raphael/solveSE.py +++ b/Raphael/solveSE.py @@ -8,7 +8,7 @@ import sys, os sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) from IAEANuclearData import IsotopeClass -class Coulomb: +class CoulombPotential: def __init__(self, rc): self.rc = rc self.id = 0 @@ -29,7 +29,7 @@ class Coulomb: else: return (self.Charge * self.ee) / (2 * self.Rc) * (3 - (x / self.Rc)**2) -class WS: +class WSPotential: def __init__(self, V0, r0, a0) : self.V0 = V0 self.r0 = r0 @@ -45,7 +45,7 @@ class WS: def output(self, x): return self.V0/(1 + math.exp((x-self.R0)/self.a0)) -class SO: +class SpinOrbitPotential: def __init__(self, VSO, rSO, aSO) : # the LS factor is put in the SolvingSE Class self.VSO = VSO @@ -89,10 +89,8 @@ class SolvingSE: dr = 0.05 nStep = 600*5 rpos = np.arange(rStart, rStart+nStep*dr, dr) - SolU = [] + SolU = [] # raidal wave function maxSolU = 0.0 - WF = np.empty([nStep], dtype=float) # radial wave function - maxWF = 0.0 #constant mn = 939.56539 #MeV/c2 @@ -111,8 +109,8 @@ class SolvingSE: potential_List = [] def PrintInput(self): - print(f" A : ({self.A:3d}, {self.ZA:3d})") - print(f" a : ({self.a:3d}, {self.Za:3d})") + print(f" A : ({self.A_A:3d}, {self.Z_A:3d})") + print(f" a : ({self.A_a:3d}, {self.Z_a:3d})") print(f" Elab : {self.Energy : 10.3f} MeV") print(f" mu : {self.mu: 10.3f} MeV/c2") # print(f" Ecm : {self.Ecm: 10.3f} MeV") @@ -125,10 +123,10 @@ class SolvingSE: def __init__(self, A, ZA, a, Za, Energy): - self.A = A - self.a = a - self.ZA = ZA - self.Za = Za + self.A_A = A + self.A_a = a + self.Z_A = ZA + self.Z_a = Za self.Z = ZA * Za self.sA = 0 @@ -139,14 +137,13 @@ class SolvingSE: self.J = 0 haha = IsotopeClass() - self.mass_A = haha.GetMassFromAZ(self.A, self.ZA) - self.mass_a = haha.GetMassFromAZ(self.a, self.Za) + self.mass_A = haha.GetMassFromAZ(self.A_A, self.Z_A) + self.mass_a = haha.GetMassFromAZ(self.A_a, self.Z_a) #self.mu = (A * a)/(A + a) * self.amu self.mu = (self.mass_A * self.mass_a)/(self.mass_A + self.mass_a) self.Energy = Energy - self.Ecm = self.Energy # self.E_tot = math.sqrt(math.pow((a+A)*self.amu,2) + 2 * A * self.amu * eng_Lab) # self.Ecm = self.E_tot - (a + A) * self.amu # self.k = math.sqrt(self.mu * 2 * abs(self.Ecm)) / self.hbarc @@ -175,10 +172,7 @@ class SolvingSE: self.rStart = rStart self.dr = dr self.nStep = nStep - self.WF=np.empty([nStep], dtype=float) self.rpos = np.arange(self.rStart, self.rStart+self.nStep*dr, self.dr) - self.WF = np.empty([self.nStep], dtype=float) - self.maxWF = 0.0 self.SolU = [] self.maxSolU = 0.0 @@ -189,15 +183,15 @@ class SolvingSE: if pot.id == 0: pot.setCharge(self.Z) if useBothMass: - pot.setAa(self.A, self.a) + pot.setAa(self.A_A, self.A_a) else: - pot.setA(self.A) + 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: + if pot.id == 2 and self.L > 0: value = value + self.LS() * pot.output(x) else: value = value + pot.output(x) @@ -240,7 +234,11 @@ class SolvingSE: self.SolU.append(y + dy) dSolU.append(z + dz) - return self.SolU + if abs(self.SolU[-1]) > self.maxSolU: + self.maxSolU = abs(self.SolU[-1]) + + return self.SolU + + def NearestPosIndex(self, r): + return min(len(self.rpos)-1, int((r - self.rStart) / self.dr)) - def normalize_boundState(self): - pass