snapshot. xsec scaling factor problem...

Raphael/IAEANuclearData.py

@@ -1 +0,0 @@
-../Cleopatra/IAEANuclearData.py

Raphael/dwba_zr.py

@@ -9,8 +9,8 @@ import time
 from sympy import S
 from sympy.physics.quantum.cg import wigner_9j
 
-# sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
+sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
-# from IAEANuclearData import IsotopeClass
+from IAEANuclearData import IsotopeClass
 
 from assLegendreP import associated_legendre_array
 from clebschGordan import clebsch_gordan, quantum_factorial, obeys_triangle_rule
@@ -21,9 +21,9 @@ import opticalPotentials as op
 from reactionData import approximate_to_half_integer, ReactionData
   
 class DWBA_ZR:
-  def __init__(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
+  def __init__(self, nu_A:str, nu_a:str, nu_b:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
 
-    self.reactDigest =  ReactionData(nu_A, nu_a, nu_b, nu_B, JB, orbital, ExB, ELabPerU)
+    self.reactDigest =  ReactionData(nu_A, nu_a, nu_b, JB, orbital, ExB, ELabPerU)
 
     if self.reactDigest.SpinBalanced == False :
       return
@@ -106,20 +106,37 @@ class DWBA_ZR:
     self.dwO.PrintPotentials()
 
     #---------------------------------------- other constants
+    print("========================================")
     D0 = 1.55e+4 # for (d,p)
 
     mass_I = self.dwI.mu
-    mass_O = self.dwO.mu
     k_I = self.dwI.k
+    mass_O = self.dwO.mu
     k_O = self.dwO.k # wave number of outgoing channel
 
+    # print(f" mu(I) : {mass_I}")
+    # print(f"  k(I) : {k_I}")
+    # print(f" mu(O) : {mass_O}")
+    # print(f"  k(O) : {k_O}")
+
     self.massBoverMassA = A_B/A_A
     self.ffactor = np.sqrt(4*np.pi)/k_I /k_O 
 
-    self.xsecScalingfactor = D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1)
+    # print(f"spin A : {self.spin_A}")
+    # print(f"spin a : {self.spin_a}")
+    # print(f"spin B : {self.spin_B}")
+
+    # self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.s +1)
+    self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1)
+
+    # print(f" spin factor : {self.spinFactor}")
+
+    self.xsecScalingfactor = D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * self.spinFactor
 
     self.radialInt = None
 
+    print(f"Xsec Scaling factor : {self.xsecScalingfactor:.6f}")
+
     self.PreCalNineJ()
     self.PreCalClebschGordan()
 
@@ -132,7 +149,7 @@ class DWBA_ZR:
       return f"{int(2*spin):+d}/2"
 
   def FindBoundState(self):
-    self.boundState.FindPotentialDepth(-70, -45, 0.5)
+    self.boundState.FindPotentialDepth(-80, -45, 0.5)
 
   def ConvertLJ2RadialIndex(self, L1:int, J1, L2:int, J2):
     index1 = int(J1 - L1 + self.spin_a)
@@ -200,6 +217,7 @@ class DWBA_ZR:
             pf2 = np.exp(1j*self.dwO.CoulombPhaseShift(L2))
             integral = simpson (bs*wf1*wf2, dx=self.boundState.dr)
             indexL2 = int(L2 - L1 + self.l)
+            # product = integral * pf1 * pf2 
             product = integral * pf1 * pf2 * self.massBoverMassA
             self.radialInt[L1][index1][indexL2][index2] = product
             # if J1 == L1 + self.spin_a and L2 == L1 + 1  and J2 == L2 - self.spin_b:
@@ -279,10 +297,10 @@ class DWBA_ZR:
     plt.show(block=False)
     input("Press Enter to continue...")
 
-  def PlotScatteringMatrix(self, isIncoming):
+  def PlotIncomingScatteringMatrix(self):
-    if isIncoming :
       self.dwI.PlotScatteringMatrix()
-    else:
+
+  def PlotOutgoingScatteringMatrix(self):
       self.dwO.PlotScatteringMatrix()
 
   def PlotIncomingDistortedWave(self,  L, J, maxR = None):
@@ -344,8 +362,8 @@ class DWBA_ZR:
 
     stop_time = time.time()
     print(f"Total time for pre-cal all CG {(stop_time - start_time) * 1000:.2f} msec")
-    print(f"self.maxL1, maxJ1, maxJ2, maxJ3")
+    print(f"max(L1 J1, L2, J2) = {self.maxL1}, {maxJ1}, {maxJ2}, {maxJ3}")
-    print(self.CG.shape)
+    print("CG shape : ",self.CG.shape)
 
   def GetPreCalCG(self, j1, m1, j2, m2, j3, m3):
     return self.CG[int(2*j1), int(2*m1 + 2*self.maxJ1+1), 
@@ -462,9 +480,14 @@ class DWBA_ZR:
     stop_time = time.time()
     print(f"\nTotal time {(stop_time - start_time) :.2f} sec")
 
-  def PrintAngDist(self):
+  def PrintAngDist(self, step:int = 1):
+    count = 0
     for th, xs in zip(self.angList, self.angDist):
-      print(f"{th:6.1f}, {xs:13.10f}")
+      if step > 1 and count % step != 0:
+        count += 1
+        continue
+      print(f"{{{th:6.1f}, {xs:13.10f}}},")
+      count += 1
 
   def PlotAngDist(self, angMin = None, angMax = None):
     plt.plot(self.angList, self.angDist)

Raphael/reactionData.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 import re
-# sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
+import sys, os
+sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
 from IAEANuclearData import IsotopeClass
 
 from clebschGordan import  obeys_triangle_rule
@@ -12,33 +13,39 @@ def approximate_to_half_integer(value):
   return round(value * 2) / 2
 
 class ReactionData:
-  def __init__(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
+  def __init__(self, nu_A:str, nu_a:str, nu_b:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
-    self.SpinBalanced = self.ReactionDigestion(nu_A, nu_a, nu_b, nu_B, JB, orbital, ExB, ELabPerU)
+    self.SpinBalanced = self.ReactionDigestion(nu_A, nu_a, nu_b, JB, orbital, ExB, ELabPerU)
     
-
+  def ReactionDigestion(self, nu_A:str, nu_a:str, nu_b:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
-  def ReactionDigestion(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
     iso = IsotopeClass()
       
     self.A_A, self.Z_A = iso.GetAZ(nu_A)
     self.A_a, self.Z_a = iso.GetAZ(nu_a)
     self.A_b, self.Z_b = iso.GetAZ(nu_b)
-    self.A_B, self.Z_B = iso.GetAZ(nu_B)
+
+    self.A_B = self.A_A + self.A_a - self.A_b
+    self.Z_B = self.Z_A + self.Z_a - self.Z_b
 
     self.ELab = self.A_a * ELabPerU
 
     mass_A = iso.GetMassFromSym(nu_A)
     mass_a = iso.GetMassFromSym(nu_a)
     mass_b = iso.GetMassFromSym(nu_b)
-    mass_B = iso.GetMassFromSym(nu_B)
+    mass_B = iso.GetMassFromAZ(self.A_B, self.Z_B)
     ExB = ExB
 
-    # sym_A = iso.GetSymbol(A_A, Z_A)
+    self.sym_A = iso.GetSymbol(self.A_A, self.Z_A)
-    # sym_B = iso.GetSymbol(A_B, Z_B)
+    self.sym_B = iso.GetSymbol(self.A_B, self.Z_B)
+
+    nu_B = f"{self.A_B}{self.sym_B}"
+    # print(nu_B)
 
     spin_A_str = iso.GetJpi(self.A_A, self.Z_A)
     self.spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str)))
     self.spin_B = float(eval(re.sub(r'[+-]', '', JB)))
 
+    print("-------- spin_B",self.spin_B)
+
     if self.A_a == 2 and self.Z_a == 1:
         self.spin_a = 1.0
         self.spin_b = 0.5
@@ -69,10 +76,10 @@ class ReactionData:
         index = match.start()  # Get position of the first letter
         
     self.node = int(orbital[:index])
-    l_sym = orbital[index:index+1]
+    self.l_sym = orbital[index:index+1]
     j_sym = orbital[index+1:]
     self.j = eval(j_sym)
-    self.l = op.ConvertLSym(l_sym)
+    self.l = op.ConvertLSym(self.l_sym)
 
     self.j = approximate_to_half_integer(self.j)
     self.s = approximate_to_half_integer(self.s)
@@ -109,11 +116,16 @@ class ReactionData:
     self.Q_value = mass_A + mass_a - mass_b - mass_B - ExB
     self.dwI = DistortedWave(nu_A, nu_a, self.ELab)
 
+    self.mass_I = self.dwI.mu
+    self.k_I = self.dwI.k
+
     Ecm_I = self.dwI.Ecm
     Ecm_O = Ecm_I + self.Q_value
     self.Eout = ((Ecm_O + mass_b + mass_B + ExB)**2 - (mass_b + mass_B + ExB)**2)/2/(mass_B + ExB)
     self.dwO = DistortedWave(nu_B, nu_b, self.Eout)
 
+    self.mass_O = self.dwO.mu
+
     Eout2 = self.ELab + self.Q_value #this is incorrec, but used in ptolmey infileCreator
 
     print("==================================================")

Raphael/solveSE.py

@@ -6,7 +6,7 @@ import re
 import sys, os
 import matplotlib.pyplot as plt
 
-# sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
+sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
 from IAEANuclearData import IsotopeClass
 
 class PotentialForm:

dwuck4/extractData.py

@@ -289,7 +289,8 @@ def plot_Xsec(data, isRuth = False):
     global plotID
 
     x_data, y_data = data
-    plt.figure(figsize=(8, 5))
+    # plt.figure(figsize=(8, 5))
+    plt.figure()
     plt.plot(x_data, y_data, linestyle="-", color="b", label="Extracted Data")
     plt.xlabel("Angle [deg]")
     if isRuth:
@@ -390,18 +391,19 @@ extract_LmaxSaSb() ## must be run first
 bs_data = extract_BoundState()
 # plot_BoundState(bs_data)
 
-# sAmpIn, sAmpOut = extract_ScatAmp()
+sAmpIn, sAmpOut = extract_ScatAmp()
-# plot_SMatrix(sAmpIn, sa)
+plot_SMatrix(sAmpIn, sa)
-# plot_SMatrix(sAmpOut, sb)
+plot_SMatrix(sAmpOut, sb)
 
 # elXsec_data = extract_ElasticXsec()
 # plot_Xsec(elXsec_data)
 
 xsec_data = extract_Xsec()
 plot_Xsec(xsec_data)
-
 x_data, y_data = xsec_data
 for i, r in enumerate(x_data):
+    if i % 5 != 0:
+        continue
     print(f"{{{r:7.3f}, {y_data[i]:10.7f}}},")
 
 def plot_RadialMatrix2(ma:float, mb:float, isPlot:bool=True):
@@ -425,8 +427,6 @@ def plot_RadialMatrix2(ma:float, mb:float, isPlot:bool=True):
 
     return radmat
 
-
-
 rList, dwIn, dwOut = extract_DistortedWave()
 def plot_DW(isIncoming:bool, L:int, m:float):
     if isIncoming :
@@ -444,7 +444,6 @@ def CoulombPS(L, eta):
 r_list, bsW = bs_data
 interp_radial = interp.interp1d(r_list, bsW, kind='cubic')
 
-
 def CalRadialIntgeral(L, ma, mb, isPlot:bool = True, verbose:int = 1):
 
     if isPlot :

dwuck4/inFileCreatorDW.py

@@ -2,6 +2,11 @@
 import sys
 import os
 
+if len(sys.argv) < 7:
+    print("Error: Not enough arguments provided.")
+    print("Usage: ./{sys.argv[0]} reaction target_gs-spin  orbital spin-pi  Ex  ELab[Mev/u]")
+    sys.exit(1)
+
 reaction = sys.argv[1]
 JA_pi = sys.argv[2]
 orbital = sys.argv[3]
@@ -9,176 +14,18 @@ JB_pi = sys.argv[4]
 Ex = float(sys.argv[5])
 ELab = float(sys.argv[6])
 
-if len(sys.argv) < 7:
-    print("Error: Not enough arguments provided.")
-    print("Usage: ./{sys.argv[0]} reaction target_gs-spin  orbital spin-pi  Ex  ELab[Mev/u]")
-    sys.exit(1)
-
 sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
+sys.path.append(os.path.join(os.path.dirname(__file__), '../Raphael'))
 from IAEANuclearData import IsotopeClass
+import opticalPotentials as op
+from reactionData import ReactionData
 
 #####################################################
 
 #  only for (d,p) or (p,d) using An & Cai, Kronning 
 
 #####################################################
-import numpy as np
 import re
-import matplotlib.pyplot as plt
-
-# Woods-Saxon
-v = 0
-r0 = 0
-a = 0
-vi = 0
-ri0 = 0
-ai = 0
-# Woods-Saxon Surface 
-vsi = 0
-rsi0 = 0
-asi = 0
-# Spin-orbit
-vso = 0
-rso0 = 0
-aso = 0
-vsoi = 0
-rsoi0 = 0
-asoi = 0
-# Coulomb 
-rc0 = 0
-
-def AnCai(A : int, Z : int, E : float):
-    global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0
-
-    A3 = A**(1./3.)
-    v  = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3
-    r0 = 1.152 - 0.00776 / A3
-    a  = 0.719 + 0.0126 * A3
-
-    vi  = 1.104 + 0.0622 * E
-    ri0 = 1.305 + 0.0997 / A3
-    ai  = 0.855 - 0.1 * A3
-
-    vsi  = 10.83 - 0.0306 * E
-    rsi0 = 1.334 + 0.152 / A3
-    asi  = 0.531 + 0.062 * A3
-
-    vso  = 3.557
-    rso0 = 0.972
-    aso  = 1.011
-
-    vsoi  = 0.0
-    rsoi0 = 0.0
-    asoi  = 0.0
-
-    rc0 = 1.303
-
-def Koning(A : int, Z : int, E : float, Zproj : float):
-    global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0
-
-    N   = A-Z
-    A3 = A**(1./3.)
-    
-    vp1 = 59.3 + 21.*(N-Z)/A - 0.024*A
-    vn1 = 59.3 - 21.*(N-Z)/A - 0.024*A
-    
-    vp2 = 0.007067 + 0.00000423*A
-    vn2 = 0.007228 - 0.00000148*A
-    
-    vp3 = 0.00001729 + 0.00000001136 * A
-    vn3 = 0.00001994 - 0.00000002 * A
-    
-    vp4 = 7e-9 # = vn4
-    vn4 = vp4
-    
-    wp1 = 14.667 + 0.009629*A
-    wn1 = 12.195 + 0.0167*A
-    
-    wp2 = 73.55 + 0.0795*A # = wn2
-    wn2 = wp2
-    
-    dp1 = 16 + 16.*(N-Z)/A
-    dn1 = 16 - 16.*(N-Z)/A
-    
-    dp2 = 0.018 + 0.003802/(1 + np.exp((A-156.)/8)) # = dn2
-    dn2 = dp2
-    
-    dp3 = 11.5  # = dn3
-    dn3 = dp3
-    
-    vso1 = 5.922 + 0.003 * A
-    vso2 = 0.004
-    
-    wso1 = -3.1
-    wso2 = 160
-    
-    epf = -8.4075 + 0.01378 *A
-    enf = -11.2814 + 0.02646 *A
-    
-    rc = 1.198 + 0.697/pow(A3,2) + 12.995/pow(A3,5)
-    vc = 1.73/rc * Z / A3
-
-    v  = vp1*(1 - vp2*(E-epf) + vp3*pow(E-epf,2) - vp4*pow(E-epf,3)) + vc * vp1 * (vp2 - 2*vp3*(E-epf) + 3*vp4*pow(E-epf,2))
-    #neutron
-    if  Zproj == 0 :
-        v  = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3))
-
-    r0 = 1.3039 - 0.4054 / A3
-    a  = 0.6778 - 0.000148 * A
-
-    vi  = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2))
-    if Zproj == 0 :
-        vi  = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2))
-    
-    ri0 = 1.3039 - 0.4054 / A3
-    ai  = 0.6778 - 0.000148 * A
-
-    vsi  = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf))
-    if Zproj == 0 :
-        vsi  = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf))
-
-    rsi0 = 1.3424 - 0.01585 * A3
-    asi  = 0.5187 + 0.0005205 * A
-    if Zproj == 0:
-        asi = 0.5446 - 0.0001656 * A
-
-    vso  = vso1 * np.exp(-vso2 * (E-epf))
-    if Zproj == 0:
-        vso = vso1 * np.exp(-vso2 * (E-enf))
-
-    rso0 = 1.1854 - 0.647/A3
-    aso  = 0.59
-
-    vsoi  = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2))
-    if Zproj == 0 :  
-        vsoi  = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2))
-
-    rsoi0 = 1.1854 - 0.647/A3
-    asoi  = 0.59
-
-    rc0 = rc
-
-def ConvertLSym(LSym :str) -> int:
-    if LSym == "s" :
-        return 0
-    elif LSym == "p" :
-        return 1
-    elif LSym == "d" :
-        return 2
-    elif LSym == "f" :
-        return 3
-    elif LSym == "g" :
-        return 4
-    elif LSym == "h" :
-        return 5
-    elif LSym == "i" :
-        return 6
-    elif LSym == "j" :
-        return 7
-    elif LSym == "k" :
-        return 8
-    else :
-        return -1
 
 #================== digest reaction
 
@@ -189,82 +36,46 @@ nu_a = nuclei[1]
 nu_b = nuclei[2]
 nu_B = nuclei[3]
 
-iso = IsotopeClass()
+reactionData = ReactionData(nu_A, nu_a, nu_b, JB_pi, orbital, Ex, ELab)
 
-A_A, Z_A = iso.GetAZ(nu_A)
+sym_A = reactionData.sym_A
-A_a, Z_a = iso.GetAZ(nu_a)
+A_A = reactionData.A_A
-A_b, Z_b = iso.GetAZ(nu_b)
+Z_A = reactionData.Z_A
-A_B, Z_B = iso.GetAZ(nu_B)
+A_a = reactionData.A_a
+Z_a = reactionData.Z_a
+A_B = reactionData.A_B
+Z_B = reactionData.Z_B
+A_b = reactionData.A_b
+Z_b = reactionData.Z_b
+A_x = reactionData.A_x
+Z_x = reactionData.Z_x
+A_c = reactionData.A_c
+Z_c = reactionData.Z_c
+node = reactionData.node
+l_sym = reactionData.l_sym
 
-A_x = abs(A_a - A_b)
+spin_a = reactionData.spin_a
-Z_x = abs(Z_a - Z_b)
+spin_b = reactionData.spin_b
 
-#---- check mass number and charge number is balnaced
+l = reactionData.l
-if A_A + A_a - A_b - A_B != 0 or Z_A + Z_a - Z_b - Z_B != 0 :
+j = reactionData.j
-    print("reaction is incorrect, mass or charge not balanced.")
-    exit()
 
-#---- check is (d,p) or (p, d)
+Q_value = reactionData.Q_value
-if (Z_a !=1 or Z_b != 1) or (A_a + A_b != 3) :
+BindingEnergy = reactionData.BindingEnergy
-    print("not (d,p) or (p,d) reaction. stop.")
-    exit()
-
-mass_A = iso.GetMassFromSym(nu_A)
-mass_a = iso.GetMassFromSym(nu_a)
-mass_b = iso.GetMassFromSym(nu_b)
-mass_B = iso.GetMassFromSym(nu_B)
-
-mass_x = iso.GetMassFromAZ( A_x, Z_x)
-
-#.... core 
-if A_A < A_B : # (d,p)
-    A_c = A_A
-    Z_c = Z_A
-    BindingEnergy = mass_B - mass_A - mass_x + Ex
-else:  #(p,d)
-    A_c = A_B
-    Z_c = Z_B
-    BindingEnergy = mass_A - mass_B - mass_x
-
-sym_A = iso.GetSymbol(A_A, Z_A)
-sym_B = iso.GetSymbol(A_B, Z_B)
-
-if A_a == 2 and Z_a == 1:
-    spin_a = 1.0
-    spin_b = 0.5   
-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")
-
-#=================== digest orbital
-match = re.search(r'[a-zA-Z]', orbital)  # Find first letter
-if match:
-    index = match.start()  # Get position of the first letter
-    
-node = int(orbital[:index])
-l_sym = orbital[index:index+1]
-j_sym = orbital[index+1:]
-j = eval(j_sym)
-l = ConvertLSym(l_sym)
 
 #=================== outfile name
 fileOutName = str(sym_A) + str(A_A) + "_" + str(nu_a) + str(nu_b) + "_" \
            + str(node) + str(l_sym) + str(int(2*j)) + "_" + str(Ex) + "_" + str(ELab) + ".in"
 
-print(fileOutName)
 
 #=================== find the maximum L for partial wave
-mass_I = mass_A * mass_a / (mass_A + mass_a) # reduced mass of incoming channel
+mass_I = reactionData.mass_I # reduced mass of incoming channel
-hbarc = 197.3269788 # MeV.fm
+k_I = reactionData.k_I # wave number of incoming channel
-k_I = np.sqrt(2*mass_I * A_a * ELab)/hbarc # wave number of incoming channel
 touching_Radius = 1.25*(A_A**(1./3) + A_a**(1./3)) + 10 # add 10 fm 
 maxL = int(touching_Radius * k_I) # maximum partial wave
-print(f"max L : {maxL}")
+
+print(f"file out : {fileOutName}")
+print(f"   max L : {maxL}")
 
 #=================== create outfile
 with open(fileOutName, "w") as file:
@@ -274,87 +85,87 @@ with open(fileOutName, "w") as file:
     file.write(f"{0.1:+08.4f}{15:+08.4f}\n")
 #===== Block 5
     if A_a == 2 :
-        AnCai(A_A, Z_A, A_a*ELab)
+        op.AnCai(A_A, Z_A, A_a*ELab)
     else:
-        Koning(A_A, Z_A, A_a*ELab, Z_a)
+        op.Koning(A_A, Z_A, A_a*ELab, Z_a)
 
     file.write(f"{A_a*ELab:+08.4f}")
     file.write(f"{A_a:+08.4f}")
     file.write(f"{Z_a:+08.4f}")
     file.write(f"{A_A:+08.4f}")
     file.write(f"{Z_A:+08.4f}")
-    file.write(f"{rc0:+08.4f}")
+    file.write(f"{op.rc0:+08.4f}")
     file.write(f"{"":8s}")
     file.write(f"{"":8s}")
     file.write(f"{2*spin_a:+08.4f}\n")
     # Woods-Saxon
     file.write(f"{1:+08.4f}") 
-    file.write(f"{-v:+08.4f}") # real
+    file.write(f"{-op.v:+08.4f}") # real
-    file.write(f"{r0:+08.4f}") # 
+    file.write(f"{op.r0:+08.4f}") # 
-    file.write(f"{a:+08.4f}") # 
+    file.write(f"{op.a:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-vi:+08.4f}") # imag
+    file.write(f"{-op.vi:+08.4f}") # imag
-    file.write(f"{ri0:+08.4f}") # 
+    file.write(f"{op.ri0:+08.4f}") # 
-    file.write(f"{ai:+08.4f}\n") # 
+    file.write(f"{op.ai:+08.4f}\n") # 
     # Woods-Saxon surface
     file.write(f"{2:+08.4f}") 
     file.write(f"{"":8s}") # real
     file.write(f"{"":8s}") # 
     file.write(f"{"":8s}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{4*vsi:+08.4f}") # imag
+    file.write(f"{4*op.vsi:+08.4f}") # imag
-    file.write(f"{rsi0:+08.4f}") # 
+    file.write(f"{op.rsi0:+08.4f}") # 
-    file.write(f"{asi:+08.4f}\n") # 
+    file.write(f"{op.asi:+08.4f}\n") # 
     # Spin-Orbit 
     file.write(f"{-4:+08.4f}") 
-    file.write(f"{-4*vso:+08.4f}") # real
+    file.write(f"{-4*op.vso:+08.4f}") # real
-    file.write(f"{rso0:+08.4f}") # 
+    file.write(f"{op.rso0:+08.4f}") # 
-    file.write(f"{aso:+08.4f}") # 
+    file.write(f"{op.aso:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-4*vsoi:+08.4f}") # imag
+    file.write(f"{-4*op.vsoi:+08.4f}") # imag
-    file.write(f"{rsoi0:+08.4f}") # 
+    file.write(f"{op.rsoi0:+08.4f}") # 
-    file.write(f"{asoi:+08.4f}\n") # 
+    file.write(f"{op.asoi:+08.4f}\n") # 
 #===== Block 6
     if A_a == 2 :
-        Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b)
+        op.Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b)
     else:
-        AnCai(A_B, Z_B, A_a*ELab + Q_value - Ex)
+        op.AnCai(A_B, Z_B, A_a*ELab + Q_value - Ex)
     file.write(f"{Q_value:+08.4f}")
     file.write(f"{A_b:+08.4f}")
     file.write(f"{Z_b:+08.4f}")
     file.write(f"{A_B:+08.4f}")
     file.write(f"{Z_B:+08.4f}")
-    file.write(f"{rc0:+08.4f}")
+    file.write(f"{op.rc0:+08.4f}")
     file.write(f"{"":8s}")
     file.write(f"{"":8s}")
     file.write(f"{2*spin_b:+08.4f}\n")
     # Woods-Saxon
     file.write(f"{1:+08.4f}") 
-    file.write(f"{-v:+08.4f}") # real
+    file.write(f"{-op.v:+08.4f}") # real
-    file.write(f"{r0:+08.4f}") # 
+    file.write(f"{op.r0:+08.4f}") # 
-    file.write(f"{a:+08.4f}") # 
+    file.write(f"{op.a:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-vi:+08.4f}") # imag
+    file.write(f"{-op.vi:+08.4f}") # imag
-    file.write(f"{ri0:+08.4f}") # 
+    file.write(f"{op.ri0:+08.4f}") # 
-    file.write(f"{ai:+08.4f}\n") # 
+    file.write(f"{op.ai:+08.4f}\n") # 
     # Woods-Saxon surface
     file.write(f"{2:+08.4f}") 
     file.write(f"{"":8s}") # real
     file.write(f"{"":8s}") # 
     file.write(f"{"":8s}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{4*vsi:+08.4f}") # imag
+    file.write(f"{4*op.vsi:+08.4f}") # imag
-    file.write(f"{rsi0:+08.4f}") # 
+    file.write(f"{op.rsi0:+08.4f}") # 
-    file.write(f"{asi:+08.4f}\n") # 
+    file.write(f"{op.asi:+08.4f}\n") # 
     # Spin-Orbit 
     file.write(f"{-4:+08.4f}") 
-    file.write(f"{-4*vso:+08.4f}") # real
+    file.write(f"{-4*op.vso:+08.4f}") # real
-    file.write(f"{rso0:+08.4f}") # 
+    file.write(f"{op.rso0:+08.4f}") # 
-    file.write(f"{aso:+08.4f}") # 
+    file.write(f"{op.aso:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-4*vsoi:+08.4f}") # imag
+    file.write(f"{-4*op.vsoi:+08.4f}") # imag
-    file.write(f"{rsoi0:+08.4f}") # 
+    file.write(f"{op.rsoi0:+08.4f}") # 
-    file.write(f"{asoi:+08.4f}\n") # 
+    file.write(f"{op.asoi:+08.4f}\n") # 
 #====== bound state
     file.write(f"{BindingEnergy:+08.4f}")
     file.write(f"{A_x:+08.4f}")