diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py
index 992d6c8..2e11ec3 100755
--- a/Raphael/DWBA_ZR.py
+++ b/Raphael/DWBA_ZR.py
@@ -5,22 +5,22 @@ 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, -3.273)
-boundState.SetPotential(1.25, 0.65, -6, 1.10, 0.65, 1.30)
-boundState.FindPotentialDepth(-75, -40, 0.1)
-# boundState.PrintWF()
-boundState.PlotBoundState()
+# boundState = BoundState(16, 8, 1, 0, 1, 0, 0.5, -3.273)
+# boundState.SetPotential(1.25, 0.65, -6, 1.10, 0.65, 1.30)
+# boundState.FindPotentialDepth(-75, -40, 0.1)
+# # boundState.PrintWF()
+# boundState.PlotBoundState()
 
-exit()
+# exit()
 
 from distortedWave import DistortedWave
 
-# 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", "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()
@@ -32,10 +32,12 @@ from distortedWave import DistortedWave
 # dw.AddPotential(CoulombPotential(1.303), False)
 
 
-# dw.CalScatteringMatrix()
+dw.CalScatteringMatrix()
 # dw.PrintScatteringMatrix()
 
-# dw.PlotDCSUnpolarized(180, 1)
+dw.PlotDCSUnpolarized(180, 1)
+
+exit()
 
 # for i in range(1, 19):
 #   theta = 10*i
diff --git a/Raphael/clebschGordan.py b/Raphael/clebschGordan.py
new file mode 100755
index 0000000..8b000e2
--- /dev/null
+++ b/Raphael/clebschGordan.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from scipy.special import gamma
+
+# from sympy.physics.quantum.cg import CG
+# from sympy import S
+# def clebsch_gordan(j1, m1, j2, m2, j, m):
+#     cg = CG(S(j1), S(m1), S(j2), S(m2), S(j), S(m))
+#     result = cg.doit()
+#     return np.complex128(result)
+
+import numpy as np
+from math import sqrt
+
+def quantum_factorial(n):
+    """
+    Calculate factorial for integer or half-integer numbers using gamma function.
+    For integer n: n! = n * (n-1) * ... * 1
+    For half-integer n: n! = Γ(n + 1)
+    """
+    if n < 0:
+        return 0.0
+    return gamma(n + 1)
+
+def clebsch_gordan(j1, m1,j2, m2, j, m):
+    """
+    Calculate Clebsch-Gordan coefficient <j1 m1 j2 m2 | j m>
+    
+    Parameters:
+    j1, j2: angular momentum quantum numbers
+    m1, m2: magnetic quantum numbers
+    j: total angular momentum quantum number
+    m: total magnetic quantum number
+    
+    Returns:
+    float: Clebsch-Gordan coefficient value
+    """
+    # Check validity of inputs using triangular inequalities and conservation
+    if not np.isclose(m, m1 + m2, atol=1e-10):
+        return 0.0
+    if abs(m1) > j1 or abs(m2) > j2 or abs(m) > j:
+        return 0.0
+    if not (abs(j1 - j2) <= j <= j1 + j2):
+        return 0.0
+    if j1 < 0 or j2 < 0 or j < 0:
+        return 0.0
+    
+    # Ensure all quantum numbers are either integer or half-integer
+    if not (np.mod(2*j1, 1) < 1e-10 or np.isclose(np.mod(2*j1, 1), 1, atol=1e-10)):
+        return 0.0
+    if not (np.mod(2*j2, 1) < 1e-10 or np.isclose(np.mod(2*j2, 1), 1, atol=1e-10)):
+        return 0.0
+    if not (np.mod(2*j, 1) < 1e-10 or np.isclose(np.mod(2*j, 1), 1, atol=1e-10)):
+        return 0.0
+
+    # Calculate the coefficient
+    prefactor = sqrt((2*j + 1) * quantum_factorial(j1 + j2 - j) * 
+                    quantum_factorial(j1 - j2 + j) * 
+                    quantum_factorial(-j1 + j2 + j) / 
+                    quantum_factorial(j1 + j2 + j + 1))
+    
+    prefactor *= sqrt(quantum_factorial(j + m) * quantum_factorial(j - m) * 
+                     quantum_factorial(j1 - m1) * quantum_factorial(j1 + m1) * 
+                     quantum_factorial(j2 - m2) * quantum_factorial(j2 + m2))
+    
+    # Sum over k
+    sum_result = 0.0
+    k_min = max(0, max(j2 - j - m1, j1 + m2 - j))
+    k_max = min(j1 + j2 - j, min(j1 - m1, j2 + m2))
+    
+    # Ensure k_min and k_max are integers for the range
+    k_min = int(np.ceil(k_min))
+    k_max = int(np.floor(k_max))
+    
+    for k in range(k_min, k_max + 1):
+        denominator = (quantum_factorial(k) * 
+                      quantum_factorial(j1 + j2 - j - k) * 
+                      quantum_factorial(j1 - m1 - k) * 
+                      quantum_factorial(j2 + m2 - k) * 
+                      quantum_factorial(j - j2 + m1 + k) * 
+                      quantum_factorial(j - j1 - m2 + k))
+        if np.isclose(denominator, 0, atol=1e-10):
+            continue
+        term = (-1)**k / denominator
+        sum_result += term
+    
+    return prefactor * sum_result
+
diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py
index 979ecac..44a74ed 100755
--- a/Raphael/distortedWave.py
+++ b/Raphael/distortedWave.py
@@ -17,18 +17,20 @@ def SevenPointsSlope(data, n):
 def FivePointsSlope(data, n):
   return ( data[n + 2] - 8 * data[n + 1] + 8 * data[n - 1] -  data[n - 2] ) / 12
 
-from sympy.physics.quantum.cg import CG
-from sympy import S
-def clebsch_gordan(j1, m1, j2, m2, j, m):
-    cg = CG(S(j1), S(m1), S(j2), S(m2), S(j), S(m))
-    result = cg.doit()
-    return np.complex128(result)
+# from sympy.physics.quantum.cg import CG
+# from sympy import S
+# def clebsch_gordan(j1, m1, j2, m2, j, m):
+#     cg = CG(S(j1), S(m1), S(j2), S(m2), S(j), S(m))
+#     result = cg.doit()
+#     return np.complex128(result)
 
 def KroneckerDelta(i, j):
   if i == j:
     return 1
   else:
     return 0
+  
+from clebschGordan import clebsch_gordan
 
 ############################################################
 class DistortedWave(SolvingSE):
@@ -181,8 +183,9 @@ class DistortedWave(SolvingSE):
       value = 0
       for J in  np.arange(Jmin, Jmax + 1, 1):
         index = int(J - Jmin)
-        value += clebsch_gordan(l, 0, self.S, v0, J, v0) * clebsch_gordan(l, v0-v, self.S, v, J, v0) * self.ScatMatrix[l][index]
-
+        cg1 = clebsch_gordan(l, 0, self.S, v0, J, v0)
+        cg2 = clebsch_gordan(l, v0 - v, self.S, v, J, v0)
+        value += cg1 * cg2 * self.ScatMatrix[l][index]
       return value - KroneckerDelta(v, v0)
 
   def CalLegendre(self, theta_deg, maxL = None, maxM = None):
@@ -224,7 +227,7 @@ class DistortedWave(SolvingSE):
     value = value / (2 * self.S + 1) 
     return value
 
-  def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None):
+  def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None, verbose = False):
     theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1)
 
     thetaTick = 30
@@ -237,7 +240,8 @@ class DistortedWave(SolvingSE):
         y_values.append(1)
       else:
         y_values.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta))
-        print(f"{theta:6.2f}, {y_values[-1]:10.6f}")
+        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')