added boundState.py class

Raphael/DWBA_ZR.py

@@ -1,15 +1,9 @@
 #!/usr/bin/env python3
 
-from solveSE import WS, Coulomb, SO, WSSurface, SolvingSE
+from boundState import BoundState
-
+    
-
+boundState = BoundState(16, 8, 1, 0, 0, 2, 2.5, -4.14)
-boundState = SolvingSE(16, 8, 1, 0, -4.14)
+boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65)
-boundState.SetRange(0, 0.1, 300)
+boundState.FindPotentialDepth(-70, -60, 0.5)
-boundState.PrintInput()
+# boundState.PrintWF()
-
+boundState.PlotBoundState()
-boundState.ClearPotential()
-boundState.AddPotential(WS(-40, 1.10, 0.65))
-
-boundState.SetLJ(0, 0.5)
-
-print(boundState.GetPotentialValue(0))

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})")
+    

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_A:3d}, {self.Z_A:3d})")
-    print(f"     a : ({self.a:3d}, {self.Za: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 = A
-    self.a = a
+    self.A_a = a
-    self.ZA = ZA
+    self.Z_A = ZA
-    self.Za = 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_A, self.Z_A)
-    self.mass_a = haha.GetMassFromAZ(self.a, self.Za)
+    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