improve the g(x) in RK4 method

Raphael/DWBA.py

@@ -30,38 +30,64 @@ from dwba_zr import DWBA_ZR
 
 
 ####################################### Simple distorted wave calculation
-kaka = DistortedWave("11C", "p", 60)
+'''
-kaka.SetRange(0, 0.1, 1000)
+kaka = DistortedWave("148Sm", "a", 50)
-kaka.maxL = 14
 kaka.PrintInput()
 kaka.ClearPotential()
-kaka.AddPotential(WoodsSaxonPot(-34.714-6.749j, 1.122, 0.676), False) # False = only use 11C for radius calculation
+kaka.AddPotential(WoodsSaxonPot(-65.500 -29.800j, 1.427, 0.671), False) # False = only use 11C for radius calculation
-kaka.AddPotential(WS_SurfacePot(       -3.194j, 1.307, 0.524), False)
+kaka.AddPotential(CoulombPotential(1.4), False)
-kaka.AddPotential(SpinOrbit_Pot( -4.532+0.477j, 0.894, 0.590), False)
-kaka.AddPotential(CoulombPotential(1.578), False)
-kaka.PrintPotentials()
-
-kaka.CalScatteringMatrix(True)
-kaka.PrintScatteringMatrix()
-kaka.PlotScatteringMatrix()
-
-# kaka.PlotDistortedWave(1, 1.5)
-
-# kaka.PlotDCSUnpolarized(180, 1, None, True)
-
-exit()
-
-kaka = DistortedWave("60Ni", "p", 30)
-kaka.PrintInput()
-kaka.ClearPotential()
-kaka.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.120, 0.669), False) # False = only use 11C for radius calculation
-kaka.AddPotential(WS_SurfacePot(       -6.878j, 1.280, 0.550), False)
-kaka.AddPotential(SpinOrbit_Pot( -5.250+0.162j, 1.020, 0.590), False)
-kaka.AddPotential(CoulombPotential(1.258), False)
 kaka.PrintPotentials()
 
 kaka.CalScatteringMatrix()
 kaka.PrintScatteringMatrix()
 # kaka.PlotScatteringMatrix()
 
-kaka.PlotDCSUnpolarized(180, 1, None, True)
+kaka.CalAngDistribution(180, 0.5, None, False)
+kaka.PlotDCSUnpolarized()
+'''
+
+'''
+kaka = DistortedWave("60Ni", "p", 30)
+# kaka.SetRange(0, 0.02, 1000)
+kaka.PrintInput()
+kaka.ClearPotential()
+kaka.AddPotential(WoodsSaxonPot(-47.937 -2.853j, 1.200, 0.669), False) # False = only use 11C for radius calculation
+kaka.AddPotential(WS_SurfacePot(        -6.878j, 1.280, 0.550), False)
+kaka.AddPotential(SpinOrbit_Pot( -5.250 +0.162j, 1.020, 0.590), False)
+kaka.AddPotential(CoulombPotential(1.258), False)
+kaka.PrintPotentials()
+
+# kaka.PlotPotential(0, 0.5, 10)
+
+kaka.CalScatteringMatrix(True)
+kaka.PrintScatteringMatrix()
+kaka.PlotScatteringMatrix()
+# kaka.PlotDistortedWave(1, 1.5)
+
+kaka.CalAngDistribution(180, 0.5, None, False)
+kaka.PlotDCSUnpolarized()
+'''
+
+kaka = DistortedWave("11C", "p", 60)
+# kaka.SetRange(0, 0.1, 1000)
+kaka.maxL = 14
+kaka.PrintInput()
+kaka.ClearPotential()
+kaka.AddPotential(WoodsSaxonPot(-34.714 -6.749j, 1.122, 0.676), False) # False = only use 11C for radius calculation
+kaka.AddPotential(WS_SurfacePot(        -3.194j, 1.307, 0.524), False)
+kaka.AddPotential(SpinOrbit_Pot( -4.532 +0.477j, 0.894, 0.590), False)
+kaka.AddPotential(CoulombPotential(1.578), False)
+kaka.PrintPotentials()
+
+# kaka.PlotPotential(10)
+
+kaka.CalScatteringMatrix(True)
+kaka.PrintScatteringMatrix()
+# kaka.PlotScatteringMatrix()
+
+# kaka.PlotDistortedWave(1, 1.5)
+kaka.CalAngDistribution(180, 0.5, None, False)
+kaka.PlotDCSUnpolarized()
+
+
+

Raphael/distortedWave.py

@@ -44,6 +44,13 @@ class DistortedWave(SolvingSE):
       eta = self.eta
     return np.angle(gamma(L+1+1j*eta))
   
+  def CoulombHankel(self, sign, rho, L = None, eta = None) -> complex:
+    if L is None:
+      L = self.L
+    if eta is None:
+      eta = self.eta
+    return float(coulombg(L, eta, rho)) + sign * 1j * float(coulombf(L, eta, rho)) 
+
   def CalScatteringMatrix(self, normTo1 = False, maxL = None, verbose = False):
     start_time = time.time()  # Start the timer
 
@@ -56,7 +63,6 @@ class DistortedWave(SolvingSE):
     tempZeroList = np.zeros(len(self.rpos), dtype=np.complex128)
 
     for L in range(0, maxL+1):
-      # sigma = self.CoulombPhaseShift()
 
       temp_ScatMatrix = []
       temp_distortedWaveU = []
@@ -70,22 +76,67 @@ class DistortedWave(SolvingSE):
         self.SetLJ(L, J)
         self.SolveByRK4()
 
+        #============ 2-point using Hankel
+        # r1 = self.rpos[-2]
+        # u1 = self.solU[-2]
+        # Hp1 = self.CoulombHankel( 1, self.k * r1)
+        # Hm1 = self.CoulombHankel(-1, self.k * r1)
+
+        # r2 = self.rpos[-1]
+        # u2 = self.solU[-1]
+        # Hp2 = self.CoulombHankel( 1, self.k * r2)
+        # Hm2 = self.CoulombHankel(-1, self.k * r2)
+
+        # ScatMatrix = (u1 * Hm2 - Hm1 * u2) / (u1 * Hp2 - Hp1 * u2)
+        # norm = 2j * (u1 * Hp2 - Hp1 * u2 ) / (Hp1 * Hm2 - Hm1 * Hp2)
+
+        #=========== 2 point G, F, this is fastest
         r1 = self.rpos[-2]
+        u1 = self.solU[-2]
         f1 = float(coulombf(self.L, self.eta, self.k*r1))
         g1 = float(coulombg(self.L, self.eta, self.k*r1))
-        u1 = self.solU[-2]
 
         r2 = self.rpos[-1]
+        u2 = self.solU[-1]
         f2 = float(coulombf(self.L, self.eta, self.k*r2))
         g2 = float(coulombg(self.L, self.eta, self.k*r2))
-        u2 = self.solU[-1]
+
+        print(f"{L:2d}, {J:4.1f} | {r1:.3f}, {f1:10.6f}, {g1:10.6f} | {r2:.3f}, {f2:10.6f}, {g2:10.6f}")
         
         det = f2*g1 - f1*g2
-
         A = (f2*u1 - u2*f1) / det
         B = (u2*g1 - g2*u1) / det
-
         ScatMatrix = (B + A * 1j)/(B - A * 1j)
+        norm = B - A * 1j
+
+        #============ 5 points slopes
+        # r_list = self.rpos[-5:]
+        # u_list = self.solU[-5:]
+        # f_list = [float(coulombf(self.L, self.eta, self.k * r)) for r in r_list]
+        # g_list = [float(coulombg(self.L, self.eta, self.k * r)) for r in r_list]
+        # u0 = u_list[2]
+        # f0 = f_list[2]
+        # g0 = g_list[2]
+        # du = FivePointsSlope(u_list, 2)
+        # df = FivePointsSlope(f_list, 2)
+        # dg = FivePointsSlope(g_list, 2)
+
+        # det = df*g0 - dg*f0
+        # A = (df*u0 - du*f0) /det
+        # B = (du*g0 - dg*u0) /det
+        # ScatMatrix = (B + A * 1j)/(B - A * 1j)
+        # norm = B - A * 1j
+
+        # #============ 100 points fitting
+        # r_list = self.rpos[-5:]
+        # u_list = self.solU[-5:]
+        # f_list = [float(coulombf(self.L, self.eta, self.k * r)) for r in r_list]
+        # g_list = [float(coulombg(self.L, self.eta, self.k * r)) for r in r_list]
+
+        # # Fit u_list to A * g_list + B * f_list
+        # A, B = np.linalg.lstsq(np.vstack([g_list, f_list]).T, u_list, rcond=None)[0]
+        # ScatMatrix = (B + A * 1j) / (B - A * 1j)
+        # norm = B - A * 1j
 
         if verbose:
           print(f"{{{L},{J}, {np.real(ScatMatrix):10.6f} +  {np.imag(ScatMatrix):10.6f}I}}")
@@ -96,8 +147,7 @@ class DistortedWave(SolvingSE):
         if normTo1 :
           dwU /= self.maxSolU
         else:
-          #dwU *= np.exp(-1j*sigma)/(B-A*1j)
+          dwU *= 1./norm
-          dwU *= 1./(B-A*1j)
         temp_distortedWaveU.append(dwU)
       
       self.ScatMatrix.append(temp_ScatMatrix)
@@ -239,25 +289,31 @@ class DistortedWave(SolvingSE):
     value = value / (2 * self.S + 1) 
     return value
 
-  def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None, verbose = False):
+  def CalAngDistribution(self,  thetaRange = 180, thetaStepDeg = 0.2, maxL = None, verbose = False):
-    theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1)
+    self.theta_list = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1)
+    self.angDist = []
+    for theta in self.theta_list:
+      if theta == 0:
+        self.angDist.append(1)
+      else:
+        self.angDist.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta))
+        if verbose :
+          print(f"{theta:6.2f}, {self.angDist[-1]:10.6f}")
 
+  def PrintAngDistribution(self):
+    print("======================= Angular Distribution")
+    for theta, dist in zip(self.theta_list, self.angDist):
+      print(f"{theta:5.1f}, {dist:10.3f}")
+
+  def PlotDCSUnpolarized(self):
     thetaTick = 30
+    thetaRange = max(self.theta_list)
     if thetaRange < 180:
       thetaTick = 10
 
-    y_values = []
-    for theta in theta_values:
-      if theta == 0:
-        y_values.append(1)
-      else:
-        y_values.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta))
-        if verbose :
-          print(f"{theta:6.2f}, {y_values[-1]:10.6f}")
-
     plt.figure(figsize=(8, 6))
     # plt.plot(theta_values, y_values, marker='o', linestyle='-', color='blue')    
-    plt.plot(theta_values, y_values, linestyle='-', color='blue')
+    plt.plot(self.theta_list, self.angDist, linestyle='-', color='blue')
     plt.title("Differential Cross Section (Unpolarized)")
     plt.xlabel("Angle [deg]")
     plt.ylabel("D.C.S / Ruth")

Raphael/solveSE.py

@@ -4,6 +4,7 @@ import math
 import numpy as np
 import re
 import sys, os
+import matplotlib.pyplot as plt
 
 sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
 from IAEANuclearData import IsotopeClass
@@ -22,7 +23,7 @@ class PotentialForm:
   #   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))
+    self.R0= self.r0 * (A**(1./3.) +a**(1./3.))
 
   def output(self, x):
     return 0
@@ -58,7 +59,7 @@ class WoodsSaxonPot(PotentialForm):
     self.id = 1
 
   def output(self, x):
-    if self.V0 == 0.0:
+    if abs(self.V0) == 0.0:
       return 0
     else:
       return self.V0/(1 + math.exp((x-self.R0)/self.a0))
@@ -75,13 +76,15 @@ class SpinOrbit_Pot(PotentialForm):
     self.id = 2
 
   def output(self, x):
-    if self.V0 == 0.0 :
+    if abs(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
+        exponent = (x - self.R0) / self.a0
+        return 4*(self.V0 * math.exp(exponent))/(self.a0*math.pow(1+math.exp(exponent),2))/x
       else :
-        return 4*1e+19
+        # return 4*1e+19
+        return 0
 
   def printPot(self):
     return super().printPot("Spin-Orbit")
@@ -94,7 +97,7 @@ class WS_SurfacePot(PotentialForm):
     self.id = 3
 
   def output(self, x):
-    if self.V0 == 0 :
+    if abs(self.V0) == 0 :
       return 0
     else:
       exponent = (x - self.R0) / self.a0
@@ -158,6 +161,8 @@ class SolvingSE:
     self.mu = (self.mass_A * self.mass_a)/(self.mass_A + self.mass_a)
     self.Ecm = 0.0
 
+    self.Factor = 2*self.mu/math.pow(self.hbarc,2)
+
   def ConstructUsingAZ(self, A, ZA, a, Za, ELabPerA):
     # print(f"ConstructUsingAZ : {A}, {ZA}, {a}, {Za}, {ELabPerA:.3f}")
     self.A_A = A
@@ -183,7 +188,6 @@ class SolvingSE:
     self.Z = self.Z_A * self.Z_a
     self.Energy = ELabPerA
 
-
   def CalCMConstants(self, useELabAsEcm = False):
     if useELabAsEcm:
       self.E_tot = self.Energy # total energy in CM
@@ -240,10 +244,17 @@ class SolvingSE:
     value = 0
     for pot in self.potential_List:
       if isinstance(pot, PotentialForm):
-        if pot.id == 2 and self.L > 0:
+        if pot.id == 2 :
-          value = value + self.LS() * pot.output(x)
+          if self.L > 0:
+            value = value + self.LS() * pot.output(x)
         else:
           value = value + pot.output(x)
+
+    if x > 1e-6:
+      value = self.Factor*(value) + self.L*(1+self.L)/x/x
+    else:
+      value = 0
+
     return value
   
   def PrintPotentials(self):
@@ -251,14 +262,38 @@ class SolvingSE:
       if isinstance(pot, PotentialForm):
         pot.printPot()
 
+  def PlotPotential(self, L, J, maxR = None):
+    self.SetLJ(L, J)
+    pot = []
+    e_list = []
+    for r in self.rpos:
+      pot.append(self.__PotentialValue(r))
+
+    pot_real = [np.real(self.__PotentialValue(r)) for r in self.rpos]
+    pot_imag = [np.imag(self.__PotentialValue(r)) for r in self.rpos]
+
+    plt.figure()
+    plt.plot(self.rpos, pot_real, label='Real Part')
+    plt.plot(self.rpos, pot_imag, label='Imaginary Part')
+    plt.xlabel('r (fm)')
+    plt.ylabel('Potential (MeV)')
+    plt.title('Potential vs. r')
+    if maxR != None:
+      plt.xlim(-1, maxR)
+    plt.ylim([min(pot_real)*1.1, abs(min(pot_real)*1.1)])
+    plt.legend()
+    plt.grid(True)
+    plt.show(block=False)
+    input("press anykey to continous....")
+
   def GetPotentialValue(self, x):
     return self.__PotentialValue(x)
 
   # The G-function, u''[r] = G[r, u[r], u'[r]]
   def __G(self, x, y, dy):
     #return  -2*x*dy -2*y  # solution of gaussian
-    if x > 0 :
+    if x > 1e-6 :
-      return 2*self.mu/math.pow(self.hbarc,2)*(self.__PotentialValue(x) - self.Ecm)*y + self.L*(1+self.L)/x/x*y
+      return self.__PotentialValue(x)*y - self.Factor*self.Ecm*y
     else:
       return 0