diff --git a/Cleopatra/IAEANuclearData.py b/Cleopatra/IAEANuclearData.py
index 874221c..59ae5e7 100644
--- a/Cleopatra/IAEANuclearData.py
+++ b/Cleopatra/IAEANuclearData.py
@@ -62,10 +62,26 @@ class IsotopeClass:
       return [1, 'H']
     if ASym == "d" :
       return [2, 'H']
+    if ASym == "t" :
+      return [3, 'H']
+    if ASym == "h" :
+      return [3, 'He']
+    if ASym == "a" :
+      return [4, 'He']
     match =  re.match(r'(\d+)(\D+)', ASym)
     return [int(match.group(1)), match.group(2) ]
 
   def GetAZ(self, ASym : str):
+    if ASym == "p" :
+      return [1, 1]
+    if ASym == "d" :
+      return [2, 1]
+    if ASym == "t" :
+      return [3, 1]
+    if ASym == "h" :
+      return [3, 2]
+    if ASym == "a" :
+      return [3, 2]
     [A, sym] = self.BreakDownName(ASym)
     try:
       dudu = self.data[(self.data['symbol']==sym) & (self.data['A']==A)]
@@ -119,6 +135,13 @@ class IsotopeClass:
       return dudu['jp'].iloc[0]
     except:
       return "unknown"
+  
+  def GetJpi(self, A : int, Z : int):
+    try:
+      dudu = self.data[(self.data['z']==Z) & (self.data['A']==A)]
+      return dudu['jp'].iloc[0]
+    except:
+      return "unknown"
 
   def GetHalfLife(self, ASym : str) -> float:
     [A, sym] = self.BreakDownName(ASym)
diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py
index edcfb1e..2d86212 100755
--- a/Raphael/DWBA_ZR.py
+++ b/Raphael/DWBA_ZR.py
@@ -1,9 +1,80 @@
 #!/usr/bin/env python3
 
 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
+
+from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, SolvingSE
+from mpmath import coulombf, coulombg
+
+import numpy as np
+from scipy.special import gamma
+
+# boundState = BoundState(16, 8, 1, 0, 1, 0, 0.5, -4.14)
+# boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.25)
+# boundState.FindPotentialDepth(-75, -60, 0.1)
+# # boundState.PrintWF()
+# boundState.PlotBoundState()
+
+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
+
+def FivePointsSlope(data, n):
+  return ( data[n + 2] - 8 * data[n + 1] + 8 * data[n - 1] -  data[n - 2] ) / 12
+
+
+dw = SolvingSE("60Ni", "p", 30)
+dw.SetRange(0, 0.1, 300)
+dw.SetLJ(0, 0.5)
+dw.dsolu0 = 1
+dw.CalCMConstants()
+dw.PrintInput()
+
+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)
+
+rpos = dw.rpos
+solU = dw.SolveByRK4()
+
+solU /= dw.maxSolU
+
+# for r, u in zip(rpos, solU):
+#     print(f"{r:.3f} {np.real(u):.6f} {np.imag(u):.6f}")
+
+def CoulombPhaseShift(L, eta):
+  return np.angle(gamma(L+1+1j*eta))
+
+sigma = CoulombPhaseShift(dw.L, dw.eta)
+# find pahse shift by using the asymptotic behavior of the wave function
+r1 = rpos[-2]
+f1 = float(coulombf(dw.L, dw.eta, dw.k*r1))
+g1 = float(coulombg(dw.L, dw.eta, dw.k*r1))
+u1 = solU[-2]
+
+r2 = rpos[-1]
+f2 = float(coulombf(dw.L, dw.eta, dw.k*r2))
+g2 = float(coulombg(dw.L, dw.eta, dw.k*r2))
+u2 = solU[-1]
+
+det = f2*g1 - f1*g2
+A = (f2*u1 - u2*f1) / det
+B = (u2*g1 - g2*u1) / det
+
+print(f"A = {np.real(A):.6f} + {np.imag(A):.6f} I")
+print(f"B = {np.real(B):.6f} + {np.imag(B):.6f} I")
+
+ScatMatrix = (B + A * 1j)/(B - A * 1j)
+print(f"Scat Matrix = {np.real(ScatMatrix):.6f} + {np.imag(ScatMatrix):.6f} I")
+
+solU = np.array(solU, dtype=np.complex128)
+solU *= np.exp(1j * sigma)/(B-A*1j)
+
+from matplotlib import pyplot as plt
+
+plt.plot(rpos, np.real(solU), label="Real")
+plt.plot(rpos, np.imag(solU), label="Imaginary")
+plt.legend()
+plt.show(block=False)
+
+input("Press Enter to continue...")
\ No newline at end of file
diff --git a/Raphael/boundState.py b/Raphael/boundState.py
index df20038..976d666 100755
--- a/Raphael/boundState.py
+++ b/Raphael/boundState.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python3
 
-from solveSE import WSPotential, CoulombPotential, SpinOrbitPotential,  SolvingSE
+from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot,  SolvingSE
 
 from scipy.optimize import curve_fit
 from scipy.interpolate import interp1d
@@ -12,22 +12,14 @@ 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.CalCMConstants( True)
     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.PrintInput()
+    self.node = node # number of nodes of the wave function r > 0
     self.FoundBounfState = False
 
   def SetPotential(self, r0, a0, Vso, rso, aso, rc = 0.0):
@@ -38,8 +30,8 @@ class BoundState(SolvingSE):
     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
+    self.AddPotential(WoodsSaxonPot(-60, r0, a0), False) # not use mass number of a
+    self.AddPotential(SpinOrbit_Pot(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
 
@@ -50,8 +42,8 @@ class BoundState(SolvingSE):
     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
+      self.AddPotential(WoodsSaxonPot(V0, self.r0, self.a0), False)
+      self.AddPotential(SpinOrbit_Pot(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()
@@ -62,16 +54,16 @@ class BoundState(SolvingSE):
         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.")
+      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
+      return
 
     f = interp1d(lastSolU, V0List, kind='cubic')
     self.V_BS = f(0)
@@ -79,8 +71,8 @@ class BoundState(SolvingSE):
 
     # 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
+    self.AddPotential(WoodsSaxonPot(self.V_BS, self.r0, self.a0), False)
+    self.AddPotential(SpinOrbit_Pot(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)
 
@@ -109,7 +101,7 @@ class BoundState(SolvingSE):
     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
+    #extrapolate wf with quadrotic 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])
@@ -117,18 +109,15 @@ class BoundState(SolvingSE):
     
     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]
+        W_values = [np.real(whitw(-1j*self.eta, self.L + 0.5,  2*self.k * r))/(self.k * r) for r in self.rpos[1:]]
+        W_values.insert(0, 0)
 
-        self.ANC = self.wf[rIndex] / W_values[rIndex]
+        self.ANC = float(self.wf[rIndex] / W_values[rIndex])
         print(f"ANC : {self.ANC:10.6e}")
         self.wf[rIndex:] = self.ANC * np.array(W_values[rIndex:])
 
diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py
index 1b8aca0..cd943ab 100755
--- a/Raphael/solveSE.py
+++ b/Raphael/solveSE.py
@@ -2,7 +2,7 @@
 
 import math
 import numpy as np
-
+import re
 import sys, os
 
 sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
@@ -29,7 +29,7 @@ class CoulombPotential:
     else:
         return (self.Charge * self.ee) / (2 * self.Rc) * (3 - (x / self.Rc)**2)
 
-class WSPotential:
+class WoodsSaxonPot:
   def __init__(self, V0, r0, a0) :
     self.V0 = V0
     self.r0 = r0
@@ -45,7 +45,7 @@ class WSPotential:
   def output(self, x):
     return self.V0/(1 + math.exp((x-self.R0)/self.a0))
 
-class SpinOrbitPotential:
+class SpinOrbit_Pot:
   def __init__(self, VSO, rSO, aSO) :
     # the LS factor is put in the SolvingSE Class
     self.VSO = VSO
@@ -65,7 +65,7 @@ class SpinOrbitPotential:
     else :
       return 4*1e+19
 
-class WSSurface:
+class WS_SurfacePot:
   def __init__(self, V0, r0, a0):
     self.V0 = V0
     self.r0 = r0
@@ -80,7 +80,7 @@ class WSSurface:
 
   def output(self, x):
     exponent = (x - self.R0) / self.a0
-    return self.V0 * math.exp(exponent) / (1 + math.exp(exponent))**2
+    return 4* self.V0 * math.exp(exponent) / (1 + math.exp(exponent))**2
 
 #========================================
 class SolvingSE:
@@ -109,52 +109,78 @@ class SolvingSE:
   potential_List = []
 
   def PrintInput(self):
-    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")
-#    print(f"     k : {self.k: 10.3f} MeV/c")
-#    print(f"   eta : {self.eta: 10.3f}")
-#    print(f"     L : {self.L},  maxL : {self.maxL}")
+    print(f"     A : ({self.A_A:3d}, {self.Z_A:3d}), spin : {self.spin_A},")
+    print(f"     a : ({self.A_a:3d}, {self.Z_a:3d}), spin : {self.spin_a},")
+    print(f"  Elab : {self.Energy : 10.5f} MeV")
+    print(f"    mu : {self.mu: 10.5f} MeV/c2")
+    print(f"   Ecm : {self.Ecm: 10.5f} MeV")
+    print(f"     k : {self.k: 10.5f} MeV/c")
+    print(f"   eta : {self.eta: 10.5f}")
+    print(f"     L : {self.L},  maxL : {self.maxL}")
     print(f"    dr : {self.dr} fm, nStep : {self.nStep}")
     print(f"rStart : {self.rStart} fm,  rMax : {self.nStep * self.dr} fm")
-    print(f"spin-A : {self.sA}, spin-a : {self.sa} ")
 
+  def __init__(self, A_or_SymA = None, ZA_or_Syma = None \
+               , a_or_ELabPerA = None, Za_or_none = None, ELabPerA_or_none = None):
+    if Za_or_none is None :
+      self.ConstructUsingSymbol(A_or_SymA, ZA_or_Syma, a_or_ELabPerA)
+    else:
+      self.ConstructUsingAZ(A_or_SymA, ZA_or_Syma, a_or_ELabPerA, Za_or_none, ELabPerA_or_none)
 
-  def __init__(self, A, ZA, a, Za, Energy):
+    haha = IsotopeClass()
+    self.mass_A = haha.GetMassFromAZ(self.A_A, self.Z_A)
+    self.mass_a = haha.GetMassFromAZ(self.A_a, self.Z_a)
+    self.spin_A = float(eval(re.sub(r'[+-]', '', haha.GetJpi(self.A_A, self.Z_A))))
+    self.spin_a = float(eval(re.sub(r'[+-]', '', haha.GetJpi(self.A_a, self.Z_a))))
+    self.S = self.spin_a
+
+    self.mu = (self.mass_A * self.mass_a)/(self.mass_A + self.mass_a)
+    self.Ecm = 0.0
+
+  def ConstructUsingAZ(self, A, ZA, a, Za, ELabPerA):
+    print(f"ConstructUsingAZ : {A}, {ZA}, {a}, {Za}, {ELabPerA}")
     self.A_A = A
     self.A_a = a
     self.Z_A = ZA
     self.Z_a = Za
     self.Z = ZA * Za
 
-    self.sA = 0
-    self.sa = 0
+    self.L = 0
+    self.S = 0
+    self.J = 0
+    self.Energy = ELabPerA
 
+  def ConstructUsingSymbol(self, Sym_A, Sym_a, ELabPerA):
+    print(f"ConstructUsingSymbol : {Sym_A}, {Sym_a}, {ELabPerA}")
     self.L = 0
     self.S = 0
     self.J = 0
 
     haha = IsotopeClass()
-    self.mass_A = haha.GetMassFromAZ(self.A_A, self.Z_A)
-    self.mass_a = haha.GetMassFromAZ(self.A_a, self.Z_a)
+    self.A_A, self.Z_A = haha.GetAZ(Sym_A)
+    self.A_a, self.Z_a = haha.GetAZ(Sym_a)
+    self.Z = self.Z_A * self.Z_a
+    self.Energy = ELabPerA
 
-    #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.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
-#    self.eta = self.Z * self.ee * math.sqrt( self.mu/2/self.Ecm ) / self.hbarc
+  def CalCMConstants(self, useELabAsEcm = False):
+    if useELabAsEcm:
+      self.E_tot = self.Energy # total energy in CM
+      self.Ecm = self.Energy # KE in cm
+    else:
+      self.E_tot = math.sqrt(math.pow((self.mass_a+self.mass_A),2) + 2 * self.mass_A * self.Energy)
+      self.Ecm = self.E_tot - (self.mass_a + self.mass_A) 
+    
+    self.k = math.sqrt(self.mu * 2 * abs(self.Ecm)) / self.hbarc
+    if self.Z == 0 :
+      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**(1/3) + self.a**(1/3)) + 3))
-
-  def SetSpin(self, sA, sa):
-    self.sA = sA
-    self.sa = sa
-    self.S = self.sa
+  # def SetA_ExSpin(self, ExA, sA ):
+  #   self.ExA = ExA
+  #   self.spin_A = sA
 
   def SetLJ(self, L, J):
     self.L = L
@@ -234,7 +260,7 @@ class SolvingSE:
       self.SolU.append(y + dy)
       dSolU.append(z + dz)
 
-      if abs(self.SolU[-1]) > self.maxSolU:
+      if np.real(self.SolU[-1]) > self.maxSolU:
         self.maxSolU = abs(self.SolU[-1])
 
     return self.SolU