diff --git a/Raphael/DWBA.py b/Raphael/DWBA.py
index 348675f..b3b9dec 100755
--- a/Raphael/DWBA.py
+++ b/Raphael/DWBA.py
@@ -2,7 +2,8 @@
 
 import time
 import matplotlib.pyplot as plt
-from distortedWave import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, DistortedWave
+from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot
+from distortedWave import  DistortedWave
 from dwba_zr import DWBA_ZR
 
 # haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10)
diff --git a/Raphael/dwba_zr.py b/Raphael/dwba_zr.py
index b9961ad..9232a01 100755
--- a/Raphael/dwba_zr.py
+++ b/Raphael/dwba_zr.py
@@ -66,29 +66,29 @@ class DWBA_ZR:
     
     print("====================== Incoming wave function ")
     if A_a == 2 and Z_a == 1:
-      op.AnCai(A_A, Z_A, self.ELab)
+      pot = op.AnCai(A_A, Z_A, self.ELab)
     if A_a == 1 and Z_a == 1:
-      op.Koning(A_A, Z_A, self.ELab , Z_a)
+      pot = op.Koning(A_A, Z_A, self.ELab , Z_a)
 
     self.dwI = self.reactDigest.dwI
     self.dwI.spin_A = self.spin_A
     self.dwI.PrintInput()
     self.dwI.ClearPotential()
-    self.dwI.AddPotential(WoodsSaxonPot(   -op.v,    op.r0,    op.a), False)
-    self.dwI.AddPotential(WoodsSaxonPot(-1j*op.vi,   op.ri0,   op.ai), False)
-    self.dwI.AddPotential(WS_SurfacePot(-1j*op.vsi,  op.rsi0,  op.asi), False)
-    self.dwI.AddPotential(SpinOrbit_Pot(   -op.vso,  op.rso0,  op.aso), False)
-    self.dwI.AddPotential(SpinOrbit_Pot(-1j*op.vsoi, op.rsoi0, op.asoi), False)
-    self.dwI.AddPotential(CoulombPotential( op.rc0), False)
+    self.dwI.AddPotential(WoodsSaxonPot(   -pot.v,    pot.r0,    pot.a), False)
+    self.dwI.AddPotential(WoodsSaxonPot(-1j*pot.vi,   pot.ri0,   pot.ai), False)
+    self.dwI.AddPotential(WS_SurfacePot(-1j*pot.vsi,  pot.rsi0,  pot.asi), False)
+    self.dwI.AddPotential(SpinOrbit_Pot(   -pot.vso,  pot.rso0,  pot.aso), False)
+    self.dwI.AddPotential(SpinOrbit_Pot(-1j*pot.vsoi, pot.rsoi0, pot.asoi), False)
+    self.dwI.AddPotential(CoulombPotential( pot.rc0), False)
     self.dwI.PrintPotentials()
 
     self.maxL1 = self.dwI.maxL
 
     print("====================== Outgoing wave function ")
     if A_b == 1 and Z_b == 1:
-      op.Koning(A_B, Z_B, self.Eout, Z_b)
+      pot = op.Koning(A_B, Z_B, self.Eout, Z_b)
     if A_b == 2 and Z_b == 1:
-      op.AnCai(A_B, Z_B, self.Eout)
+      pot = op.AnCai(A_B, Z_B, self.Eout)
 
     self.maxL2 = self.maxL1 + self.l
 
@@ -97,12 +97,12 @@ class DWBA_ZR:
     self.dwO.maxL = self.maxL2
     self.dwO.PrintInput()
     self.dwO.ClearPotential()
-    self.dwO.AddPotential(WoodsSaxonPot(   -op.v,    op.r0,    op.a), False)
-    self.dwO.AddPotential(WoodsSaxonPot(-1j*op.vi,   op.ri0,   op.ai), False)
-    self.dwO.AddPotential(WS_SurfacePot(-1j*op.vsi,  op.rsi0,  op.asi), False)
-    self.dwO.AddPotential(SpinOrbit_Pot(   -op.vso,  op.rso0,  op.aso), False)
-    self.dwO.AddPotential(SpinOrbit_Pot(-1j*op.vsoi, op.rsoi0, op.asoi), False)
-    self.dwO.AddPotential(CoulombPotential( op.rc0), False)
+    self.dwO.AddPotential(WoodsSaxonPot(   -pot.v,    pot.r0,    pot.a), False)
+    self.dwO.AddPotential(WoodsSaxonPot(-1j*pot.vi,   pot.ri0,   pot.ai), False)
+    self.dwO.AddPotential(WS_SurfacePot(-1j*pot.vsi,  pot.rsi0,  pot.asi), False)
+    self.dwO.AddPotential(SpinOrbit_Pot(   -pot.vso,  pot.rso0,  pot.aso), False)
+    self.dwO.AddPotential(SpinOrbit_Pot(-1j*pot.vsoi, pot.rsoi0, pot.asoi), False)
+    self.dwO.AddPotential(CoulombPotential( pot.rc0), False)
 
     self.dwO.PrintPotentials()
 
diff --git a/Raphael/opticalPotentials.py b/Raphael/opticalPotentials.py
index d159c9c..5dec4bd 100755
--- a/Raphael/opticalPotentials.py
+++ b/Raphael/opticalPotentials.py
@@ -4,137 +4,180 @@
 
 import numpy as np
 
-# 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
+class OpticalPotential:
+    # 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
+    def Print(self):
+        print(f"Wood-Saxon (Re) : {self.v:7.3f}, {self.r0:5.3f}, {self.a:5.3f}")
+        print(f"           (Im) : {self.vi:7.3f}, {self.ri0:5.3f}, {self.ai:5.3f}")
+        print(f"WS Surface (Im) : {self.vsi:7.3f}, {self.rsi0:5.3f}, {self.asi:5.3f}")
+        print(f"Spin-Orbit (Re) : {self.vso:7.3f}, {self.rso0:5.3f}, {self.aso:5.3f}")
+        print(f"           (Im) : {self.vsoi:7.3f}, {self.rsoi0:5.3f}, {self.asoi:5.3f}")
+        print(f"Coulomb    (Re) : {'':7s}, {self.rc0:5.3f}")
 
-    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
+class SuAndHan(OpticalPotential):
+    def __init__(self,  A : int, Z :int , E : float):
+        N = A - Z
+        A3 = A**(1./3.)
 
-    vsi  = 10.83 - 0.0306 * E
-    rsi0 = 1.334 + 0.152 / A3
-    asi  = 0.531 + 0.062 * A3
+        vsiCOND = 27.5816 - 0.0797 * E + 48.0*(N-Z)/A
+        viCOND = -4.0174 + 0.1409 * E 
+        
+        self.v  = 175.0881 - 0.6236 * E + 0.0006*E*E + 30.*(N-Z)/A - 0.236 * Z/A3
+        self.r0 = 1.3421
+        self.a  = 0.6578
 
-    vso  = 3.557
-    rso0 = 0.972
-    aso  = 1.011
+        self.vi  = viCOND
+        if  viCOND < 0 :
+            self.vi = 0.0
+        self.ri0 = 1.4259
+        self.ai  = 0.6578
 
-    vsoi  = 0.0
-    rsoi0 = 0.0
-    asoi  = 0.0
+        self.vsi  = vsiCOND
+        if vsiCOND < 0 :
+            self.vsi = 0.0
+        self.rsi0 = 1.2928
+        self.asi  = 0.6359
 
-    rc0 = 1.303
+        self.vso  = 0.0
+        self.rso0 = 1.2686
+        self.aso  = 0.85
 
-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
+        self.vsoi  = 0.0
+        self.rsoi0 = 0.0
+        self.asoi  = 0.0
 
-    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
+        self.rc0 = 1.350
 
-    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
+class AnCai(OpticalPotential):
+    def __init__(self,  A : int, Z :int , E : float):
 
-    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
+        A3 = A**(1./3.)
+        self.v  = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3
+        self.r0 = 1.152 - 0.00776 / A3
+        self.a  = 0.719 + 0.0126 * A3
 
-    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))
+        self.vi  = 1.104 + 0.0622 * E
+        self.ri0 = 1.305 + 0.0997 / A3
+        self.ai  = 0.855 - 0.1 * A3
 
-    rsi0 = 1.3424 - 0.01585 * A3
-    asi  = 0.5187 + 0.0005205 * A
-    if Zproj == 0:
-        asi = 0.5446 - 0.0001656 * A
+        self.vsi  = 10.83 - 0.0306 * E
+        self.rsi0 = 1.334 + 0.152 / A3
+        self.asi  = 0.531 + 0.062 * A3
 
-    vso  = vso1 * np.exp(-vso2 * (E-epf))
-    if Zproj == 0:
-        vso = vso1 * np.exp(-vso2 * (E-enf))
+        self.vso  = 3.557
+        self.rso0 = 0.972
+        self.aso  = 1.011
 
-    rso0 = 1.1854 - 0.647/A3
-    aso  = 0.59
+        self.vsoi  = 0.0
+        self.rsoi0 = 0.0
+        self.asoi  = 0.0
 
-    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))
+        self.rc0 = 1.303
 
-    rsoi0 = 1.1854 - 0.647/A3
-    asoi  = 0.59
+class Koning(OpticalPotential):
+    def __init__(self,  A : int, Z :int , E : float, Zproj : float):
 
-    rc0 = rc
+        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
+
+        self.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 :
+            self.v  = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3))
+
+        self.r0 = 1.3039 - 0.4054 / A3
+        self.a  = 0.6778 - 0.000148 * A
+
+        self.vi  = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2))
+        if Zproj == 0 :
+            self.vi  = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2))
+        
+        self.ri0 = 1.3039 - 0.4054 / A3
+        self.ai  = 0.6778 - 0.000148 * A
+
+        self.vsi  = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf))
+        if Zproj == 0 :
+            self.vsi  = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf))
+
+        self.rsi0 = 1.3424 - 0.01585 * A3
+        self.asi  = 0.5187 + 0.0005205 * A
+        if Zproj == 0:
+            self.asi = 0.5446 - 0.0001656 * A
+
+        self.vso  = vso1 * np.exp(-vso2 * (E-epf))
+        if Zproj == 0:
+            self.vso = vso1 * np.exp(-vso2 * (E-enf))
+
+        self.rso0 = 1.1854 - 0.647/A3
+        self.aso  = 0.59
+
+        self.vsoi  = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2))
+        if Zproj == 0 :  
+            self.vsoi  = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2))
+
+        self.rsoi0 = 1.1854 - 0.647/A3
+        self.asoi  = 0.59
 
 def ConvertLSym(LSym :str) -> int:
     if LSym == "s" :
diff --git a/dwuck4/inFileCreatorDW.py b/dwuck4/inFileCreatorDW.py
index 69635e6..f0af1f3 100755
--- a/dwuck4/inFileCreatorDW.py
+++ b/dwuck4/inFileCreatorDW.py
@@ -85,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 :
-        op.AnCai(A_A, Z_A, A_a*ELab)
+        pot = op.AnCai(A_A, Z_A, A_a*ELab)
     else:
-        op.Koning(A_A, Z_A, A_a*ELab, Z_a)
+        pot = 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"{op.rc0:+08.4f}")
+    file.write(f"{pot.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"{-op.v:+08.4f}") # real
-    file.write(f"{op.r0:+08.4f}") # 
-    file.write(f"{op.a:+08.4f}") # 
+    file.write(f"{-pot.v:+08.4f}") # real
+    file.write(f"{pot.r0:+08.4f}") # 
+    file.write(f"{pot.a:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-op.vi:+08.4f}") # imag
-    file.write(f"{op.ri0:+08.4f}") # 
-    file.write(f"{op.ai:+08.4f}\n") # 
+    file.write(f"{-pot.vi:+08.4f}") # imag
+    file.write(f"{pot.ri0:+08.4f}") # 
+    file.write(f"{pot.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*op.vsi:+08.4f}") # imag
-    file.write(f"{op.rsi0:+08.4f}") # 
-    file.write(f"{op.asi:+08.4f}\n") # 
+    file.write(f"{4*pot.vsi:+08.4f}") # imag
+    file.write(f"{pot.rsi0:+08.4f}") # 
+    file.write(f"{pot.asi:+08.4f}\n") # 
     # Spin-Orbit 
     file.write(f"{-4:+08.4f}") 
-    file.write(f"{-4*op.vso:+08.4f}") # real
-    file.write(f"{op.rso0:+08.4f}") # 
-    file.write(f"{op.aso:+08.4f}") # 
+    file.write(f"{-4*pot.vso:+08.4f}") # real
+    file.write(f"{pot.rso0:+08.4f}") # 
+    file.write(f"{pot.aso:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-4*op.vsoi:+08.4f}") # imag
-    file.write(f"{op.rsoi0:+08.4f}") # 
-    file.write(f"{op.asoi:+08.4f}\n") # 
+    file.write(f"{-4*pot.vsoi:+08.4f}") # imag
+    file.write(f"{pot.rsoi0:+08.4f}") # 
+    file.write(f"{pot.asoi:+08.4f}\n") # 
 #===== Block 6
     if A_a == 2 :
-        op.Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b)
+        pot = op.Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b)
     else:
-        op.AnCai(A_B, Z_B, A_a*ELab + Q_value - Ex)
+        pot = 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"{op.rc0:+08.4f}")
+    file.write(f"{pot.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"{-op.v:+08.4f}") # real
-    file.write(f"{op.r0:+08.4f}") # 
-    file.write(f"{op.a:+08.4f}") # 
+    file.write(f"{-pot.v:+08.4f}") # real
+    file.write(f"{pot.r0:+08.4f}") # 
+    file.write(f"{pot.a:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-op.vi:+08.4f}") # imag
-    file.write(f"{op.ri0:+08.4f}") # 
-    file.write(f"{op.ai:+08.4f}\n") # 
+    file.write(f"{-pot.vi:+08.4f}") # imag
+    file.write(f"{pot.ri0:+08.4f}") # 
+    file.write(f"{pot.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*op.vsi:+08.4f}") # imag
-    file.write(f"{op.rsi0:+08.4f}") # 
-    file.write(f"{op.asi:+08.4f}\n") # 
+    file.write(f"{4*pot.vsi:+08.4f}") # imag
+    file.write(f"{pot.rsi0:+08.4f}") # 
+    file.write(f"{pot.asi:+08.4f}\n") # 
     # Spin-Orbit 
     file.write(f"{-4:+08.4f}") 
-    file.write(f"{-4*op.vso:+08.4f}") # real
-    file.write(f"{op.rso0:+08.4f}") # 
-    file.write(f"{op.aso:+08.4f}") # 
+    file.write(f"{-4*pot.vso:+08.4f}") # real
+    file.write(f"{pot.rso0:+08.4f}") # 
+    file.write(f"{pot.aso:+08.4f}") # 
     file.write(f"{"":8s}") # spin-orbit skipped
-    file.write(f"{-4*op.vsoi:+08.4f}") # imag
-    file.write(f"{op.rsoi0:+08.4f}") # 
-    file.write(f"{op.asoi:+08.4f}\n") # 
+    file.write(f"{-4*pot.vsoi:+08.4f}") # imag
+    file.write(f"{pot.rsoi0:+08.4f}") # 
+    file.write(f"{pot.asoi:+08.4f}\n") # 
 #====== bound state
     file.write(f"{BindingEnergy:+08.4f}")
     file.write(f"{A_x:+08.4f}")