diff --git a/calculateDipoleTrapPotential.py b/calculateDipoleTrapPotential.py
index b829611..9b3d51e 100644
--- a/calculateDipoleTrapPotential.py
+++ b/calculateDipoleTrapPotential.py
@@ -54,7 +54,6 @@ def particleDensity(w_x, w_z, Power, Polarizability, N, T, m = 164*u.u): # For a
     v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x')
     v_y = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'y')
     v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z')
-
     return N * (2 * np.pi)**3 * (v_x * v_y * v_z) *  (m / (2 * np.pi * ac.k_B * T))**(3/2)
 
 def thermaldeBroglieWavelength(T, m = 164*u.u):
@@ -69,7 +68,7 @@ def scatteringLength(B):
     return a_bkg
 
 def dipolarLength(mu = 9.93 * ac.muB, m = 164*u.u):
-    return (m * ac.mu0 * mu**2) / (8 * np.pi * ac.hbar**2)
+    return (m * ac.mu0 * mu**2) / (12 * np.pi * ac.hbar**2)
 
 def scatteringCrossSection(B):
     return 8 * np.pi * scatteringLength(B)**2 + ((32*np.pi)/45) * dipolarLength()**2
@@ -230,9 +229,6 @@ def computeTrapPotential(w_x, w_z, Power, Polarizability, options):
     TrapFrequency = [v, dv]
     ExtractedTrapFrequencies = [IdealTrapFrequency, TrapFrequency]
 
-    # Gamma_elastic = calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N = 1.13 * 1e7, T = 22 * u.uK, B = 0 * u.G)
-    # PSD = calculatePSD(w_x, w_z, Power, Polarizability, N = 1.13 * 1e7, T = 22 * u.uK).decompose()
-
     return Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies
 
 #####################################################################
@@ -324,6 +320,9 @@ if __name__ == '__main__':
     Polarizability = 184.4 # in a.u, most precise measured value of Dy polarizability
     w_x, w_z = 27.5*u.um, 33.8*u.um # Beam Waists in the x and y directions
     
+    AspectRatio = 4.6
+    w_x = AspectRatio * w_x
+
     options = {
         'axis': 1, # axis referenced to the beam along which you want the dipole trap potential
         'extent': 1e4, # range of spatial coordinates in one direction to calculate trap potential over
@@ -338,13 +337,33 @@ if __name__ == '__main__':
     ComputedPotentials = [] 
     Params = []  
 
-    Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies = computeTrapPotential(w_x, w_z, Power, Polarizability, options)
-    ComputedPotentials.append(IdealTrappingPotential)
-    ComputedPotentials.append(TrappingPotential)
-    Params.append([TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies])
+    # Positions, IdealTrappingPotential, TrappingPotential, TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies = computeTrapPotential(w_x, w_z, Power, Polarizability, options)
+    # ComputedPotentials.append(IdealTrappingPotential)
+    # ComputedPotentials.append(TrappingPotential)
+    # Params.append([TrapDepthsInKelvin, CalculatedTrapFrequencies, ExtractedTrapFrequencies])
 
-    ComputedPotentials = np.asarray(ComputedPotentials)
-    plotPotential(Positions, ComputedPotentials, options['axis'], Params)
+    # ComputedPotentials = np.asarray(ComputedPotentials)
+    # plotPotential(Positions, ComputedPotentials, options['axis'], Params)
+
+    AtomNumber = 1.13 * 1e7
+    Temperature = 30 * u.uK
+    BField = 0 * u.G
+
+    n = particleDensity(w_x, w_z, Power, Polarizability, N = AtomNumber, T = Temperature, m = 164*u.u).decompose().to(u.cm**(-3))
+    Gamma_elastic = calculateElasticCollisionRate(w_x, w_z, Power, Polarizability, N = AtomNumber, T = Temperature, B = BField)
+    PSD = calculatePSD(w_x, w_z, Power, Polarizability, N = AtomNumber, T = Temperature).decompose()
+    
+    print('%.2E' % (n.value)* (n.unit))
+    print('%.2f' % (Gamma_elastic.value) * (Gamma_elastic.unit))
+    print('%.2E' % (PSD.value))
+
+    v_x = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'x')
+    v_y = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'y')
+    v_z = calculateTrapFrequency(w_x, w_z, Power, Polarizability, dir = 'z')
+
+    print(v_x)
+    print(v_y)
+    print(v_z)
 
     # v, dv, popt, pcov = extractTrapFrequency(Positions, TrappingPotential, options['axis'])
     # plotHarmonicFit(Positions, TrappingPotential, TrapDepthsInKelvin, options['axis'], popt, pcov)