diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py
index 3584c16..5dc5f49 100755
--- a/Raphael/DWBA_ZR.py
+++ b/Raphael/DWBA_ZR.py
@@ -2,6 +2,7 @@
 
 from boundState import BoundState
 from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot
+import matplotlib.pyplot as plt
 
 # 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)
@@ -11,18 +12,37 @@ from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePo
 
 from distortedWave import DistortedWave
 
-dw = DistortedWave("60Ni", "p", 30)
+# dw = DistortedWave("60Ni", "p", 30)
+# 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)
 
+dw = DistortedWave("60Ni", "d", 60)
+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)
+dw.AddPotential(WoodsSaxonPot(-81.919, 1.15, 0.768), False)
+dw.AddPotential(WoodsSaxonPot(-4.836j, 1.33, 0.464), False)
+dw.AddPotential(WS_SurfacePot(-8.994j, 1.373, 0.774), False)
+dw.AddPotential(SpinOrbit_Pot(-3.557, 0.972, 1.011), False)
+dw.AddPotential(CoulombPotential(1.303), False)
+
+
+
 
 dw.CalScatteringMatrix()
-
 dw.PrintScatteringMatrix()
 
-# dw.PlotScatteringMatrix()
+dw.PlotDCSUnpolarized(180, 1)
 
-dw.PlotDCSUnpolarized()
\ No newline at end of file
+# for i in range(1, 19):
+#   theta = 10*i
+#   # ruth = dw.RutherFord(theta)
+#   # coulAmp = dw.CoulombScatterintAmp(theta)
+#   dw.CalLegendre(theta)
+#   nuAmp1 = dw.NuclearScatteringAmp(-0.5, 0.5, 14)
+#   nuAmp2 = dw.NuclearScatteringAmp(0.5, -0.5, 14)
+#   # dsc = dw.DCSUnpolarized(theta, 14)
+#   # print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}, {dsc:10.6f}, {ruth:10.6f}")
+#   print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}")
diff --git a/Raphael/assLegendreP.py b/Raphael/assLegendreP.py
index 6da6a7e..46ef562 100755
--- a/Raphael/assLegendreP.py
+++ b/Raphael/assLegendreP.py
@@ -2,36 +2,36 @@
 
 import numpy as np
 
-def associated_legendre_polynomials_single_angle(L, M, theta_deg):
+def associated_legendre_array(maxL, maxM, theta_deg):
     # Convert theta from degrees to radians
     theta = np.radians(theta_deg)
     x = np.cos(theta)
     
-    P = np.zeros((L + 1, M + 1))
+    P = np.zeros((maxL + 1, maxM + 1))
     
     # P^m_l for m = 0, l = 0
     P[0, 0] = 1.0
     
     # P^m_l for m = 0, l > 0
-    for l in range(1, L + 1):
+    for l in range(1, maxL + 1):
         P[l, 0] = ((2*l - 1) * x * P[l-1, 0] - (l - 1) * P[l-2, 0]) / l
     
     # P^m_l for m > 0 (using recursion)
-    for m in range(1, M + 1):
+    for m in range(1, maxM + 1):
         # P^m_m
         P[m, m] = (1 - 2*m) * P[m-1, m-1]
         P[m, m] *= np.sqrt(1 - x**2)
         
-        for l in range(m + 1, L + 1):
+        for l in range(m + 1, maxL + 1):
             # P^m_l for m < l
             P[l, m] = ((2*l - 1) * x * P[l-1, m] - (l + m - 1) * P[l-2, m]) / (l - m)
     
     return P
 
 # Example usage
-#L = 15  # Maximum l degree
-#M = 5  # Maximum m order
+# L = 15  # Maximum l degree
+# M = 3  # Maximum m order
 
-#legendre_polynomials = associated_legendre_polynomials_single_angle(L, M, 45)
+# legendre_polynomials = associated_legendre_array(L, M, 10)
 
-#print(legendre_polynomials)
\ No newline at end of file
+# print(legendre_polynomials)
\ No newline at end of file
diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py
index c90bd4b..979ecac 100755
--- a/Raphael/distortedWave.py
+++ b/Raphael/distortedWave.py
@@ -9,6 +9,7 @@ import numpy as np
 from scipy.special import gamma, sph_harm, factorial
 import matplotlib.pyplot as plt
 
+from assLegendreP import associated_legendre_array
 
 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
@@ -29,6 +30,7 @@ def KroneckerDelta(i, j):
   else:
     return 0
 
+############################################################
 class DistortedWave(SolvingSE):
   def __init__(self, target, projectile, ELab):
     super().__init__(target, projectile, ELab)
@@ -38,6 +40,8 @@ class DistortedWave(SolvingSE):
     self.ScatMatrix = []
     self.distortedWaveU = []
 
+    self.legendreArray = []
+
   def SetLJ(self, L, J):
     self.L = L
     self.J = J
@@ -105,9 +109,15 @@ class DistortedWave(SolvingSE):
 
   def PrintScatteringMatrix(self):
     for L in range(0, len(self.ScatMatrix)):
+      print("{", end="")
       for i in range(0, len(self.ScatMatrix[L])):
-        print(f"{{{L:2d},{i-self.S:4.1f}, {np.real(self.ScatMatrix[L][i]):10.6f} +  {np.imag(self.ScatMatrix[L][i]):10.6f}I}}", end=" ")
-      print()
+        print("{", end="")
+        print(f"{L:2d},{L+i-self.S:4.1f}, {np.real(self.ScatMatrix[L][i]):10.6f} +  {np.imag(self.ScatMatrix[L][i]):10.6f}I", end="")
+        if i < len(self.ScatMatrix[L])-1 :
+          print("}, ", end="")
+        else:
+          print("}", end="")
+      print("},")
 
   def GetScatteringMatrix(self, L, J):
     return self.ScatMatrix[L][J-L+self.S]
@@ -152,15 +162,15 @@ class DistortedWave(SolvingSE):
     plt.show(block=False)
     input("Press Enter to continue...")
 
-  def RutherFord(self, theta):
-    sin_half_theta = np.sin(theta / 2)
+  def RutherFord(self, theta_deg):
+    sin_half_theta = np.sin(np.radians(theta_deg + 1e-20) / 2)
     result = self.eta**2 / (4 * (self.k**2) * (sin_half_theta**4))
     return result
 
   def CoulombScatterintAmp(self, theta_deg):
-    sin_sq = pow(np.sin(np.radians(theta_deg)/2), 2)
+    sin_sq = pow(np.sin(np.radians(theta_deg + 1e-20)/2), 2)
     coulPS = self.CoulombPhaseShift(0)
-    return - self.eta/ (2*self.k*sin_sq) * np.exp(2j*(coulPS - self.eta*np.log(sin_sq)))
+    return - self.eta / (2 * self.k * sin_sq) * np.exp(1j * (2*coulPS - self.eta * np.log(sin_sq)))
   
   def GMatrix1Spin(self, v, v0, l ) -> complex:
     if self.S == 0 :
@@ -175,7 +185,21 @@ class DistortedWave(SolvingSE):
 
       return value - KroneckerDelta(v, v0)
 
-  def NuclearScatteringAmp(self, v, v0, theta, phi, maxL = None ) -> complex:
+  def CalLegendre(self, theta_deg, maxL = None, maxM = None):
+    if maxL is None:
+      maxL = self.maxL
+    if maxM is None:
+      maxM = int(2*self.S)
+    self.legendreArray = associated_legendre_array(maxL, maxM, theta_deg)
+    return self.legendreArray
+
+  def GetPreCalLegendre(self, L, M):
+    if abs (M) <= int(2*self.S):
+      return self.legendreArray[L][int(abs(M))]
+    else :
+      return 0
+
+  def NuclearScatteringAmp(self, v, v0, maxL = None ) -> complex:
     value = 0
     if maxL is None:
       maxL = self.maxL
@@ -184,21 +208,24 @@ class DistortedWave(SolvingSE):
         value += 0
       else:
         coulPS = self.CoulombPhaseShift(l)
-        value += np.sqrt(2*l+1) * sph_harm(v0 - v, l, phi, theta) * np.exp(2j * coulPS)* self.GMatrix1Spin(v, v0, l)
+        fact = pow(-1, v0-v) * np.sqrt(factorial(l - abs(v0-v))/factorial(l + abs(v0-v)))
+        value += (2*l+1) * fact * self.GetPreCalLegendre(l, v0-v) * np.exp(2j * coulPS)* self.GMatrix1Spin(v, v0, l)
 
-    return value * np.sqrt(4*np.pi)/ 2 / 1j / self.k
-
-  def DCSUnpolarized(self, theta, phi, maxL = None):
+    return value / 2j / self.k
+  
+  def DCSUnpolarized(self, theta_deg, maxL = None):
     value = 0
+    self.CalLegendre(theta_deg)
+    jaja = self.CoulombScatterintAmp(theta_deg)
     for v in np.arange(-self.S, self.S + 1, 1):
       for v0 in np.arange(-self.S, self.S + 1, 1):
-        value += abs(self.CoulombScatterintAmp(theta) * KroneckerDelta(v, v0) +  self.NuclearScatteringAmp(v, v0, theta, phi, maxL))**2
+        value += abs( jaja * KroneckerDelta(v, v0) +  self.NuclearScatteringAmp(v, v0, maxL))**2
 
-    value = value / (2 * self.S + 1)
+    value = value / (2 * self.S + 1) 
     return value
 
   def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None):
-    theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg))
+    theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1)
 
     thetaTick = 30
     if thetaRange < 180:
@@ -206,22 +233,21 @@ class DistortedWave(SolvingSE):
 
     y_values = []
     for theta in theta_values:
-      theta = np.deg2rad(theta)
       if theta == 0:
         y_values.append(1)
       else:
-        y_values.append(self.DCSUnpolarized(theta , 0, maxL)/ self.RutherFord(theta))
-      print(f"{np.rad2deg(theta):6.2f}, {y_values[-1]:10.6f}")
+        y_values.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta))
+        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', label='Real Part')    
+    # plt.plot(theta_values, y_values, marker='o', linestyle='-', color='blue')    
+    plt.plot(theta_values, y_values, linestyle='-', color='blue')
     plt.title("Differential Cross Section (Unpolarized)")
-    plt.xlabel("Theta (degrees)")
-    plt.ylabel("Value")
+    plt.xlabel("Angle [deg]")
+    plt.ylabel("D.C.S / Ruth")
     plt.yscale("log")
     plt.xticks(np.arange(0, thetaRange + 1, thetaTick))
     plt.xlim(0, thetaRange)
-    plt.legend()
     plt.grid()
     plt.show(block=False)
     input("Press Enter to continue...")
\ No newline at end of file