diff --git a/docs/source/utilities.rst b/docs/source/utilities.rst
index 9c663cf..bc4e9d8 100644
--- a/docs/source/utilities.rst
+++ b/docs/source/utilities.rst
@@ -74,8 +74,8 @@ This function takes in quite a few parameters, due to its overlapping usage with
                                         chif_guess=0.5,
                                         t0=t0,
                                         tf=tf,
-                                        force_Omega0_optimization=False, #True if we want to perform a least squares fit for Omega0.
-                                        include_Omega0_as_parameter=True, #Use Omega0 as a third parameter in the optimization algorithm, along with the final mass and spin.
+                                        force_Omega0_optimization=False, #True if we want to perform a least squares fit for Omega_0.
+                                        include_Omega0_as_parameter=True, #Use Omega_0 as a third parameter in the optimization algorithm, along with the final mass and spin.
                                         start_with_wide_search=False, #If True we start with the scipy differential_evolution algorithm and finish with the scipy minimize algorithm
                                         t_shift_range = np.arange(-10,10,0.1),#Time shift array to be used in the mismatch calculation
                                         )
diff --git a/gwBOB/BOB_terms.py b/gwBOB/BOB_terms.py
index 9a31e0e..996b7fa 100644
--- a/gwBOB/BOB_terms.py
+++ b/gwBOB/BOB_terms.py
@@ -6,8 +6,8 @@
 
 def define_BOB_symbols():
     t, t0, tp, tau = sp.symbols('t t0 tp tau', real=True)
-    Omega0, Omega_QNM, Ap, Phi_0 = sp.symbols('Omega0 Omega_QNM Ap Phi_0', positive=True)
-    return t, t0, tp, tau, Omega0, Omega_QNM, Ap, Phi_0
+    Omega_0, Omega_QNM, Ap, Phi_0 = sp.symbols('Omega_0 Omega_QNM Ap Phi_0', positive=True)
+    return t, t0, tp, tau, Omega_0, Omega_QNM, Ap, Phi_0
 
 def BOB_amplitude_sym(t, tp, tau, Ap):
     '''
@@ -27,7 +27,7 @@ def BOB_amplitude_sym(t, tp, tau, Ap):
     x = (t - tp) / tau
     return Ap * sp.sech(x)
 
-def BOB_strain_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
+def BOB_strain_freq_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM):
     '''
     Frequency for the strain when assuming the BOB amplitude best models the strain (for finite t_0)
 
@@ -36,7 +36,7 @@ def BOB_strain_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
         t0 (sympy.Symbol): Initial Condition time
         tp (sympy.Symbol): Time of peak amplitude 
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
 
     Returns:
@@ -44,11 +44,11 @@ def BOB_strain_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
     '''
     x = (t - tp) / tau
     x0 = (t0 - tp) / tau
-    Omega_ratio = Omega0 / Omega_QNM
+    Omega_ratio = Omega_0 / Omega_QNM
     exponent = (sp.tanh(x) - 1) / (sp.tanh(x0) - 1)
     return Omega_QNM * Omega_ratio**exponent
 
-def BOB_strain_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
+def BOB_strain_phase_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM, Phi_0):
     '''
     Phase for the strain when the BOB amplitude models the strain (for finite t_0)
 
@@ -57,7 +57,7 @@ def BOB_strain_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
         t0 (sympy.Symbol): Initial Condition time
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
         Phi_0 (sympy.Symbol): Initial Condition Phase (phi)/(mode number)
 
@@ -66,10 +66,10 @@ def BOB_strain_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
 
         Omega - Strain frequency (waveform frequency)/(mode number)
     '''
-    Omega = BOB_strain_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM)
+    Omega = BOB_strain_freq_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM)
     x = (t - tp) / tau
     x0 = (t0 - tp) / tau
-    Omega_ratio = Omega0 / Omega_QNM
+    Omega_ratio = Omega_0 / Omega_QNM
     outer = Omega_QNM * tau / sp.Integer(2)
     tanh_tp_t0_tau_p1 = sp.tanh(-x0) + 1
     term1_exp = sp.Integer(2) / tanh_tp_t0_tau_p1
@@ -81,7 +81,7 @@ def BOB_strain_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
     Phi = outer * inner + Phi_0
     return Phi, Omega
 
-def BOB_news_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
+def BOB_news_freq_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM):
     '''
     Frequency for the news when the BOB amplitude models the news (for finite t_0)
 
@@ -90,7 +90,7 @@ def BOB_news_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
         t0 (sympy.Symbol): Initial Condition time
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
 
     Returns:
@@ -100,11 +100,11 @@ def BOB_news_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
     '''
     x = (t - tp) / tau
     x0 = (t0 - tp) / tau
-    F = (Omega_QNM**2 - Omega0**2) / (1 - sp.tanh(x0))
+    F = (Omega_QNM**2 - Omega_0**2) / (1 - sp.tanh(x0))
     Omega2 = Omega_QNM**2 + F * (sp.tanh(x) - 1)
     return sp.sqrt(Omega2)
 
-def BOB_news_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
+def BOB_news_phase_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM, Phi_0):
     '''
     Phase for the news when the BOB amplitude models the News (for finite t_0)
 
@@ -113,7 +113,7 @@ def BOB_news_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
         t0 (sympy.Symbol): Initial Condition time
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
         Phi_0 (sympy.Symbol): Initial Condition Phase (phi)/(mode number)
 
@@ -122,15 +122,15 @@ def BOB_news_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
 
         Omega - News frequency (waveform frequency)/(mode number)
     '''
-    Omega = BOB_news_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM)
+    Omega = BOB_news_freq_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM)
     x0 = (t0 - tp) / tau
-    F = (Omega_QNM**2 - Omega0**2) / (1 - sp.tanh(x0))
+    F = (Omega_QNM**2 - Omega_0**2) / (1 - sp.tanh(x0))
     delta = sp.Integer(2) * F - Omega_QNM**2
     term1 = Omega_QNM * tau / sp.Integer(2) * sp.log((Omega + Omega_QNM) / sp.Abs(Omega - Omega_QNM))
     term2 = -sp.sqrt(delta) * tau * sp.atan(Omega / sp.sqrt(delta))
     return term1 + term2 + Phi_0, Omega
 
-def BOB_psi4_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
+def BOB_psi4_freq_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM):
     '''
     Frequency for Psi4 when the BOB amplitude models Psi4 (for finite t_0)
 
@@ -139,7 +139,7 @@ def BOB_psi4_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
         t0 (sympy.Symbol): Initial Condition time
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
 
     Returns:
@@ -147,11 +147,11 @@ def BOB_psi4_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM):
     '''
     x = (t - tp) / tau
     x0 = (t0 - tp) / tau
-    k = (Omega_QNM**4 - Omega0**4) / (1 - sp.tanh(x0))
-    X = Omega0**4 + k * (sp.tanh(x) - sp.tanh(x0))
+    k = (Omega_QNM**4 - Omega_0**4) / (1 - sp.tanh(x0))
+    X = Omega_0**4 + k * (sp.tanh(x) - sp.tanh(x0))
     return X**sp.Rational(1, 4)
 
-def BOB_psi4_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
+def BOB_psi4_phase_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM, Phi_0):
     '''
     Phase for Psi4 when the BOB amplitude models Psi4 (for finite t_0)
 
@@ -160,7 +160,7 @@ def BOB_psi4_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
         t0 (sympy.Symbol): Initial Condition time
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
         Phi_0 (sympy.Symbol): Initial Condition Phase (phi)/(mode number)
 
@@ -169,19 +169,19 @@ def BOB_psi4_phase_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM, Phi_0):
 
         Omega - Psi4 frequency (waveform frequency)/(mode number)
     '''
-    Omega = BOB_psi4_freq_finite_t0_sym(t, t0, tp, tau, Omega0, Omega_QNM)
+    Omega = BOB_psi4_freq_finite_t0_sym(t, t0, tp, tau, Omega_0, Omega_QNM)
     x0 = (t0 - tp) / tau
-    k = (Omega_QNM**4 - Omega0**4) / (1 - sp.tanh(x0))
-    KappaP = (Omega0**4 + k * (1 - sp.tanh(x0)))**sp.Rational(1, 4)
-    KappaM = (Omega0**4 - k * (1 + sp.tanh(x0)))**sp.Rational(1, 4)
-    term_P_tanh = KappaP * tau * (sp.atanh(Omega / KappaP) - sp.atanh(Omega0 / KappaP))
-    term_M_tanh = KappaM * tau * (sp.atanh(Omega / KappaM) - sp.atanh(Omega0 / KappaM))
-    term_P_tan  = KappaP * tau * (sp.atan(Omega / KappaP) - sp.atan(Omega0 / KappaP))
-    term_M_tan  = KappaM * tau * (sp.atan(Omega / KappaM) - sp.atan(Omega0 / KappaM))
+    k = (Omega_QNM**4 - Omega_0**4) / (1 - sp.tanh(x0))
+    KappaP = (Omega_0**4 + k * (1 - sp.tanh(x0)))**sp.Rational(1, 4)
+    KappaM = (Omega_0**4 - k * (1 + sp.tanh(x0)))**sp.Rational(1, 4)
+    term_P_tanh = KappaP * tau * (sp.atanh(Omega / KappaP) - sp.atanh(Omega_0 / KappaP))
+    term_M_tanh = KappaM * tau * (sp.atanh(Omega / KappaM) - sp.atanh(Omega_0 / KappaM))
+    term_P_tan  = KappaP * tau * (sp.atan(Omega / KappaP) - sp.atan(Omega_0 / KappaP))
+    term_M_tan  = KappaM * tau * (sp.atan(Omega / KappaM) - sp.atan(Omega_0 / KappaM))
     return term_P_tanh + term_P_tan - term_M_tanh - term_M_tan + Phi_0, Omega
 
 # --- Asymptotic (t0 -> -inf) ---
-def BOB_strain_freq_sym(t, tp, Omega0, Omega_QNM, tau):
+def BOB_strain_freq_sym(t, tp, Omega_0, Omega_QNM, tau):
     '''
     Frequency for the strain when the BOB amplitude models the strain (for t_0 = -inf)
 
@@ -189,16 +189,16 @@ def BOB_strain_freq_sym(t, tp, Omega0, Omega_QNM, tau):
         t (sympy.Symbol): Time 
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
 
     Returns:
         Omega - Strain frequency (waveform frequency)/(mode number)
     '''
     x = (t - tp)/tau
-    return Omega_QNM * (Omega0 / Omega_QNM)**((sp.tanh(x) - 1) / sp.Integer(-2))
+    return Omega_QNM * (Omega_0 / Omega_QNM)**((sp.tanh(x) - 1) / sp.Integer(-2))
 
-def BOB_strain_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
+def BOB_strain_phase_sym(t, tp, Omega_0, Omega_QNM, tau, Phi_0):
     '''
 
     Phase for the strain when the BOB amplitude models the strain (for t_0 = -inf)
@@ -207,7 +207,7 @@ def BOB_strain_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
         t (sympy.Symbol): Time 
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
         Phi_0 (sympy.Symbol): Initial Condition Phase (phi)/(mode number)
 
@@ -216,15 +216,15 @@ def BOB_strain_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
         
         Omega - Strain frequency (waveform frequency)/(mode number)
     '''
-    Omega = BOB_strain_freq_sym(t, tp, Omega0, Omega_QNM, tau)
+    Omega = BOB_strain_freq_sym(t, tp, Omega_0, Omega_QNM, tau)
     x = (t - tp)/tau
     outer = tau / sp.Integer(2)
-    log_sqrt_ratio = sp.Rational(1,2) * sp.log(Omega_QNM/Omega0)
+    log_sqrt_ratio = sp.Rational(1,2) * sp.log(Omega_QNM/Omega_0)
     term1_arg = log_sqrt_ratio * (sp.tanh(x) + 1)
     term2_arg = log_sqrt_ratio * (sp.tanh(x) - 1)
-    return outer * (Omega0 * Ei_sp(term1_arg) - Omega_QNM * Ei_sp(term2_arg)) + Phi_0, Omega
+    return outer * (Omega_0 * Ei_sp(term1_arg) - Omega_QNM * Ei_sp(term2_arg)) + Phi_0, Omega
 
-def BOB_news_freq_sym(t, tp, Omega0, Omega_QNM, tau):
+def BOB_news_freq_sym(t, tp, Omega_0, Omega_QNM, tau):
     '''
     
     Frequency for the news when the BOB amplitude models the news (for t_0 = -inf)
@@ -233,7 +233,7 @@ def BOB_news_freq_sym(t, tp, Omega0, Omega_QNM, tau):
         t (sympy.Symbol): Time 
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
 
     Returns:
@@ -241,11 +241,11 @@ def BOB_news_freq_sym(t, tp, Omega0, Omega_QNM, tau):
     '''
 
     x = (t - tp)/tau
-    Omega_minus = Omega_QNM**2 - Omega0**2
-    Omega_plus = Omega_QNM**2 + Omega0**2
+    Omega_minus = Omega_QNM**2 - Omega_0**2
+    Omega_plus = Omega_QNM**2 + Omega_0**2
     return sp.sqrt(Omega_minus * sp.tanh(x) / sp.Integer(2) + Omega_plus / sp.Integer(2))
 
-def BOB_news_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
+def BOB_news_phase_sym(t, tp, Omega_0, Omega_QNM, tau, Phi_0):
     '''
     Phase for the strain when the BOB amplitude models the strain (for t_0 = -inf)
 
@@ -253,7 +253,7 @@ def BOB_news_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
         t (sympy.Symbol): Time 
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
         Phi_0 (sympy.Symbol): Initial Condition Phase (phi)/(mode number)
 
@@ -262,13 +262,13 @@ def BOB_news_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
 
         Omega - Strain frequency (waveform frequency)/(mode number)
         '''
-    Omega = BOB_news_freq_sym(t, tp, Omega0, Omega_QNM, tau)
+    Omega = BOB_news_freq_sym(t, tp, Omega_0, Omega_QNM, tau)
     outer = tau / sp.Integer(2)
     inner1 = sp.log((Omega + Omega_QNM) / sp.Abs(Omega - Omega_QNM))
-    inner2 = sp.log((Omega + Omega0) / sp.Abs(Omega - Omega0))
-    return outer * (Omega_QNM * inner1 - Omega0 * inner2) + Phi_0, Omega
+    inner2 = sp.log((Omega + Omega_0) / sp.Abs(Omega - Omega_0))
+    return outer * (Omega_QNM * inner1 - Omega_0 * inner2) + Phi_0, Omega
 
-def BOB_psi4_freq_sym(t, tp, Omega0, Omega_QNM, tau):
+def BOB_psi4_freq_sym(t, tp, Omega_0, Omega_QNM, tau):
     '''
     Frequency for Psi4 when the BOB amplitude models Psi4 (for t_0 = -inf)
 
@@ -276,17 +276,17 @@ def BOB_psi4_freq_sym(t, tp, Omega0, Omega_QNM, tau):
         t (sympy.Symbol): Time 
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
 
     Returns:
         Omega - Psi4 frequency (waveform frequency) / (mode number)
     '''
     x = (t - tp)/tau
-    k = (Omega_QNM**4 - Omega0**4) / sp.Integer(2)
-    return (Omega0**4 + k * (sp.tanh(x) + 1))**sp.Rational(1, 4)
+    k = (Omega_QNM**4 - Omega_0**4) / sp.Integer(2)
+    return (Omega_0**4 + k * (sp.tanh(x) + 1))**sp.Rational(1, 4)
 
-def BOB_psi4_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
+def BOB_psi4_phase_sym(t, tp, Omega_0, Omega_QNM, tau, Phi_0):
     '''
     Phase for Psi4 when the BOB amplitude models Psi4 (for t_0 = -inf)
 
@@ -294,7 +294,7 @@ def BOB_psi4_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
         t (sympy.Symbol): Time 
         tp (sympy.Symbol): Time of peak amplitude
         tau (sympy.Symbol): Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
-        Omega0 (sympy.Symbol): Initial Condition Frequency
+        Omega_0 (sympy.Symbol): Initial Condition Frequency
         Omega_QNM (sympy.Symbol): Real part of Quasinormal mode (QNM) frequency (little omega)/(mode number)
         Phi_0 (sympy.Symbol): Initial Condition Phase (phi)/(mode number)
 
@@ -303,16 +303,16 @@ def BOB_psi4_phase_sym(t, tp, Omega0, Omega_QNM, tau, Phi_0):
 
         Omega - Psi4 frequency (waveform frequency) / (mode number)
     '''
-    Omega = BOB_psi4_freq_sym(t, tp, Omega0, Omega_QNM, tau)
-    Omega_minus_q0 = Omega_QNM - Omega0
-    Omega_plus_q0 = Omega_QNM + Omega0
+    Omega = BOB_psi4_freq_sym(t, tp, Omega_0, Omega_QNM, tau)
+    Omega_minus_q0 = Omega_QNM - Omega_0
+    Omega_plus_q0 = Omega_QNM + Omega_0
     outer_num = sp.sqrt(Omega_minus_q0 * Omega_plus_q0) * tau
     outer_den = 2 * sp.sqrt(sp.Abs(Omega_minus_q0)) * sp.sqrt(sp.Abs(Omega_plus_q0))
     outer = outer_num / outer_den
     inner1 = Omega_QNM * (sp.log(sp.Abs(Omega + Omega_QNM)) - sp.log(sp.Abs(Omega - Omega_QNM)))
-    inner2 = -Omega0 * (sp.log(sp.Abs(Omega + Omega0)) - sp.log(sp.Abs(Omega - Omega0)))
+    inner2 = -Omega_0 * (sp.log(sp.Abs(Omega + Omega_0)) - sp.log(sp.Abs(Omega - Omega_0)))
     inner3 = 2 * Omega_QNM * sp.atan(Omega / Omega_QNM)
-    inner4 = -2 * Omega0 * sp.atan(Omega / Omega0)
+    inner4 = -2 * Omega_0 * sp.atan(Omega / Omega_0)
     return outer * (inner1 + inner2 + inner3 + inner4) + Phi_0, Omega
 
 
diff --git a/gwBOB/BOB_terms_jax.py b/gwBOB/BOB_terms_jax.py
index 9a5220e..bf84c74 100644
--- a/gwBOB/BOB_terms_jax.py
+++ b/gwBOB/BOB_terms_jax.py
@@ -33,7 +33,7 @@ def convert_BOB_to_JAXBOB(BOB):
     temp =  JAXBOB(BOB.t, BOB.Omega_0, BOB.Omega_QNM, BOB.tau, BOB.Ap,BOB.tp,BOB.m)
     return temp
 
-def BOB_amplitude_jax(t, tau, Ap, t_p):
+def BOB_amplitude_jax(t, tau, Ap, tp):
     '''
     BOB amplitude evolution
 
@@ -48,16 +48,16 @@ def BOB_amplitude_jax(t, tau, Ap, t_p):
     Returns:
         A(t) : Waveform amplitude at time t
     '''
-    tt = (t - t_p) / tau
+    tt = (t - tp) / tau
     return Ap / jnp.cosh(tt)
 
-def BOB_news_freq_jax(t, Omega_0, Omega_QNM, tau, t_p, m):
+def BOB_news_freq_jax(t, Omega_0, Omega_QNM, tau, tp, m):
     '''
-    Waveform frequency for the news when the BOB amplitude models the news (taking t_0 = -inf)
+    Waveform frequency for the news when the BOB amplitude models the news (taking t0 = -inf)
 
     Args:
         t : Time 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -66,19 +66,19 @@ def BOB_news_freq_jax(t, Omega_0, Omega_QNM, tau, t_p, m):
     Returns:
         omega - News waveform frequency
     '''
-    tt = (t - t_p) / tau
+    tt = (t - tp) / tau
     Omega_minus = Omega_QNM**2 - Omega_0**2
     Omega_plus  = Omega_QNM**2 + Omega_0**2
     Omega2 = Omega_minus * jnp.tanh(tt) / 2. + Omega_plus / 2.
     return m*jnp.sqrt(jnp.maximum(Omega2, 1e-12)) 
 
-def BOB_news_phase_jax(t, Omega_0, Omega_QNM, tau, t_p, Phi_0, m=2):
+def BOB_news_phase_jax(t, Omega_0, Omega_QNM, tau, tp, Phi_0, m=2):
     '''
-    Waveform phase for the news when the BOB amplitude models the news (taking t_0 = -inf)
+    Waveform phase for the news when the BOB amplitude models the news (taking t0 = -inf)
 
     Args:
         t : Time 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -90,7 +90,7 @@ def BOB_news_phase_jax(t, Omega_0, Omega_QNM, tau, t_p, Phi_0, m=2):
 
         omega - News waveform frequency
     '''
-    omega = BOB_news_freq_jax(t, Omega_0, Omega_QNM, tau, t_p, m) #news_freq_jax returns little omega
+    omega = BOB_news_freq_jax(t, Omega_0, Omega_QNM, tau, tp, m) #news_freq_jax returns little omega
     Omega = omega/m
 
     Omega_minus_Q = jnp.abs(Omega - Omega_QNM) 
@@ -109,13 +109,13 @@ def BOB_news_phase_jax(t, Omega_0, Omega_QNM, tau, t_p, Phi_0, m=2):
     phase = (outer * (Omega_QNM * inner1 - Omega_0 * inner2) + Phi_0)*m
     
     return phase,omega
-def BOB_psi4_freq_jax(t, Omega_0, Omega_QNM, tau, t_p,m):
+def BOB_psi4_freq_jax(t, Omega_0, Omega_QNM, tau, tp,m):
     '''
-    Waveform frequency for psi4 when assuming the BOB amplitude best models psi4 (taking t_0 = -inf)
+    Waveform frequency for psi4 when assuming the BOB amplitude best models psi4 (taking t0 = -inf)
 
     Args:
         t : Time 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -124,18 +124,18 @@ def BOB_psi4_freq_jax(t, Omega_0, Omega_QNM, tau, t_p,m):
     Returns:
         omega - Psi4 waveform frequency
     '''
-    tt = (t - t_p) / tau
+    tt = (t - tp) / tau
     k = (Omega_QNM**4 - Omega_0**4) / 2.0
     X = Omega_0**4 + k * (jnp.tanh(tt) + 1.0)
     return m*jnp.sqrt(jnp.sqrt(jnp.maximum(X, 1e-12)))
 
-def BOB_strain_freq(t, Omega_0, Omega_QNM, tau, t_p,m):
+def BOB_strain_freq(t, Omega_0, Omega_QNM, tau, tp,m):
     '''
-    Waveform frequency for strain when assuming the BOB amplitude best models the strain (taking t_0 = -inf)
+    Waveform frequency for strain when assuming the BOB amplitude best models the strain (taking t0 = -inf)
 
     Args:
         t : Time 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -144,19 +144,19 @@ def BOB_strain_freq(t, Omega_0, Omega_QNM, tau, t_p,m):
     Returns:
         omega - Strain waveform frequency
     '''
-    tt = (t - t_p) / tau
+    tt = (t - tp) / tau
     Omega_ratio = Omega_0/Omega_QNM
     tanh_tt_m1 = jnp.tanh(tt)-1
     return m*Omega_QNM*(Omega_ratio**(tanh_tt_m1/(-2.)))
 
-def BOB_psi4_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t_0, t_p,m):
+def BOB_psi4_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t0, tp,m):
     '''
-    Waveform frequency for psi4 when assuming the BOB amplitude best models psi4 (for finite t_0)
+    Waveform frequency for psi4 when assuming the BOB amplitude best models psi4 (for finite t0)
 
     Args:
         t : Time 
-        t_0 : Initial Condition time
-        t_p : Time of peak amplitude
+        t0 : Initial Condition time
+        tp : Time of peak amplitude
         tau : Damping term; can also be described as 1/gamma (gamma is imaginry QNM fre)
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -165,21 +165,21 @@ def BOB_psi4_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t_0, t_p,m):
     Returns:
         omega - Psi4 waveform frequency
     '''
-    tt = (t - t_p) / tau
-    t0p = (t_0-t_p) / tau
+    tt = (t - tp) / tau
+    t0p = (t0-tp) / tau
     k_denom = 1 - jnp.tanh(t0p)
     k = (Omega_QNM**4 - Omega_0**4) / k_denom
     X = Omega_0**4 + k * (jnp.tanh(tt) - jnp.tanh(t0p))
     return m*(jnp.sqrt(jnp.sqrt(jnp.maximum(X, 1e-12))))
 
-def BOB_news_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t_0, t_p,m):
+def BOB_news_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t0, tp,m):
     '''
-    Waveform frequency for the news when assuming the BOB amplitude best models the news (for finite t_0)
+    Waveform frequency for the news when assuming the BOB amplitude best models the news (for finite t0)
 
     Args:
         t : Time 
-        t_0 : Initial Condition time
-        t_p : Time of peak amplitude
+        t0 : Initial Condition time
+        tp : Time of peak amplitude
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -188,21 +188,21 @@ def BOB_news_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t_0, t_p,m):
     Returns:
         omega - News waveform frequency
     '''
-    tt = (t - t_p) / tau
-    t0p = (t_0-t_p) / tau
+    tt = (t - tp) / tau
+    t0p = (t0-tp) / tau
     F_denom = 1 - jnp.tanh(t0p)
     F = (Omega_QNM**2 - Omega_0**2) / F_denom
     Omega2 = Omega_QNM**2 + F * (jnp.tanh(tt) - 1)
     return m*jnp.sqrt(jnp.maximum(Omega2, 1e-12))
 
-def BOB_strain_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t_0, t_p,m):
+def BOB_strain_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t0, tp,m):
     '''
-    Waveform frequency for the strain when assuming the BOB amplitude best models the strain (for finite t_0)
+    Waveform frequency for the strain when assuming the BOB amplitude best models the strain (for finite t0)
     
     Args:
         t : Time 
-        t_0 : Initial Condition time
-        t_p : Time of peak amplitude
+        t0 : Initial Condition time
+        tp : Time of peak amplitude
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Omega_0 : Initial Condition Frequency
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
@@ -211,8 +211,8 @@ def BOB_strain_freq_finite_t0(t, Omega_0, Omega_QNM, tau, t_0, t_p,m):
     Returns:
         omega - Strain waveform frequency  
     '''
-    tt = (t - t_p) / tau
-    t0p = (t_0-t_p) / tau
+    tt = (t - tp) / tau
+    t0p = (t0-tp) / tau
     Omega_ratio = Omega_0/Omega_QNM
     tanh_tt_m1 = jnp.tanh(tt)-1
     tanh_t0p_m1 = jnp.tanh(t0p)-1
@@ -228,10 +228,10 @@ def deriv_g(t):
         return g_prime
     return deriv_g
 @partial(jit, static_argnames=('omega_func', 'A_func','N'))
-def get_series_terms_ad(t, Omega_0, Omega_QNM, tau, Ap, t_p, omega_func, A_func, m, N):
+def get_series_terms_ad(t, Omega_0, Omega_QNM, tau, Ap, tp, omega_func, A_func, m, N):
     """
     Generates the raw, unsigned series terms [f₀, Df₀, D²f₀, ..., Dⁿf₀]
-    using JAX's automatic differentiation for t_0 = -inf scenarios
+    using JAX's automatic differentiation for t0 = -inf scenarios
 
     Args:
         t : Time 
@@ -239,7 +239,7 @@ def get_series_terms_ad(t, Omega_0, Omega_QNM, tau, Ap, t_p, omega_func, A_func,
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Ap : Peak waveform amplitude 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         omega_func: frequency function
         A_func: amplitude function  
         m : Mode number
@@ -250,13 +250,13 @@ def get_series_terms_ad(t, Omega_0, Omega_QNM, tau, Ap, t_p, omega_func, A_func,
     """
     # Define the base function f₀(t) = A(t) / (i * ω(t))
     def f0_func(time):
-        A = A_func(time, tau, Ap, t_p)
-        omega = omega_func(time, Omega_0, Omega_QNM, tau, t_p, m)
+        A = A_func(time, tau, Ap, tp)
+        omega = omega_func(time, Omega_0, Omega_QNM, tau, tp, m)
         return A / (1j * omega)
 
     # Define the operator D's pre-factor g(t) = 1 / (i * ω(t))
     def g_func(time):
-        omega = omega_func(time, Omega_0, Omega_QNM, tau, t_p, m)
+        omega = omega_func(time, Omega_0, Omega_QNM, tau, tp, m)
         return 1.0 / (1j * omega)
 
     # List to hold the functions that compute [f₀, Df₀, D²f₀, ...]
@@ -280,7 +280,7 @@ def g_func(time):
     
     return all_terms
 @partial(jit, static_argnames=('omega_func', 'A_func','N'))
-def get_series_terms_ad_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_0, omega_func, A_func, m, N):
+def get_series_terms_ad_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, tp, t0, omega_func, A_func, m, N):
     """
     Generates the raw, unsigned series terms [f₀, Df₀, D²f₀, ..., Dⁿf₀]
     using JAX's automatic differentiation for finite t0 scenarios.
@@ -291,7 +291,7 @@ def get_series_terms_ad_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_0, omeg
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Ap : Peak waveform amplitude 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         omega_func: frequency function
         A_func: amplitude function  
         m : Mode number
@@ -304,13 +304,13 @@ def get_series_terms_ad_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_0, omeg
     """
     # Define the base function f₀(t) = A(t) / (i * ω(t))
     def f0_func(time):
-        A = A_func(time, tau, Ap, t_p)
-        omega = omega_func(time, Omega_0, Omega_QNM, tau, t_0, t_p, m)
+        A = A_func(time, tau, Ap, tp)
+        omega = omega_func(time, Omega_0, Omega_QNM, tau, t0, tp, m)
         return A / (1j * omega)
 
     # Define the operator D's pre-factor g(t) = 1 / (i * ω(t))
     def g_func(time):
-        omega = omega_func(time, Omega_0, Omega_QNM, tau, t_0, t_p, m)
+        omega = omega_func(time, Omega_0, Omega_QNM, tau, t0, tp, m)
         return 1.0 / (1j * omega)
 
     # List to hold the functions that compute [f₀, Df₀, D²f₀, ...]
@@ -355,12 +355,12 @@ def fast_truncated_sum(all_raw_terms):
     
     return series_sum
 @partial(jit, static_argnames=('omega_func', 'A_func','N'))
-def calculate_strain_from_news(t, Omega_0, Omega_QNM, tau, Ap, t_p, 
+def calculate_strain_from_news(t, Omega_0, Omega_QNM, tau, Ap, tp, 
                                  omega_func, A_func, m, N):
     '''
     Calculate the strain from the news using the series aapproximation
     Generates the raw, unsigned series terms [f₀, Df₀, D²f₀, ..., Dⁿf₀]
-    using JAX's automatic differentiation for scenarios with t_0 = -inf
+    using JAX's automatic differentiation for scenarios with t0 = -inf
 
     Args:
         t : Time 
@@ -368,7 +368,7 @@ def calculate_strain_from_news(t, Omega_0, Omega_QNM, tau, Ap, t_p,
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Ap : Peak waveform amplitude 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         omega_func: frequency function
         A_func: amplitude function  
         m : Mode number
@@ -380,18 +380,18 @@ def calculate_strain_from_news(t, Omega_0, Omega_QNM, tau, Ap, t_p,
         series sum
     '''
     # 1. Generate the raw, unsigned derivative terms
-    all_raw_terms = get_series_terms_ad(t, Omega_0, Omega_QNM, tau, Ap, t_p,
+    all_raw_terms = get_series_terms_ad(t, Omega_0, Omega_QNM, tau, Ap, tp,
                                         omega_func, A_func, m, N)
     sum = fast_truncated_sum(all_raw_terms)
     return all_raw_terms,sum
 
 @partial(jit, static_argnames=('omega_func', 'A_func','N'))
-def calculate_strain_from_news_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_0,
+def calculate_strain_from_news_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, tp, t0,
                                  omega_func, A_func, m, N):
     '''
     Calculate the strain from the news using the series aapproximation
     Generates the raw, unsigned series terms [f₀, Df₀, D²f₀, ..., Dⁿf₀]
-    using JAX's automatic differentiation for scenarios invoving finite t_0 values.
+    using JAX's automatic differentiation for scenarios invoving finite t0 values.
 
     Args:
         t : Time 
@@ -399,7 +399,7 @@ def calculate_strain_from_news_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Ap : Peak waveform amplitude 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         omega_func: frequency function
         A_func: amplitude function  
         m : Mode number
@@ -411,7 +411,7 @@ def calculate_strain_from_news_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_
         series sum
     '''
     # 1. Generate the raw, unsigned derivative terms
-    all_raw_terms = get_series_terms_ad_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p, t_0,
+    all_raw_terms = get_series_terms_ad_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, tp, t0,
                                         omega_func, A_func, m, N)
     sum = fast_truncated_sum(all_raw_terms)
     return all_raw_terms,sum
@@ -438,12 +438,12 @@ def _build_symbolic_series(base_func, g_func, N_order):
     return term_funcs
 
 @partial(jit, static_argnames=('omega_func', 'A_func', 'N'))
-def calculate_strain_from_psi4(t, Omega_0, Omega_QNM, tau, Ap, t_p,
+def calculate_strain_from_psi4(t, Omega_0, Omega_QNM, tau, Ap, tp,
                                             omega_func, A_func, m, N):
     '''
     Calculate the strain from psi4 using the series aapproximation
     Generates the raw, unsigned series terms [f₀, Df₀, D²f₀, ..., Dⁿf₀]
-    using JAX's automatic differentiation for scenarios with t_0 = -inf.
+    using JAX's automatic differentiation for scenarios with t0 = -inf.
 
     Args:
         t : Time 
@@ -451,7 +451,7 @@ def calculate_strain_from_psi4(t, Omega_0, Omega_QNM, tau, Ap, t_p,
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Ap : Peak waveform amplitude 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         omega_func: frequency function
         A_func: amplitude function  
         m : Mode number
@@ -465,13 +465,13 @@ def calculate_strain_from_psi4(t, Omega_0, Omega_QNM, tau, Ap, t_p,
     M = N
     # --- Define the single D operator pre-factor ONCE ---
     def g_func(time):
-        omega = omega_func(time, Omega_0, Omega_QNM, tau, t_p, m)
+        omega = omega_func(time, Omega_0, Omega_QNM, tau, tp, m)
         return 1.0 / (1j * omega)
 
     # --- Stage 1: Build the symbolic function for the News sum ---
     def f0_psi4_func(time):
-        A = A_func(time, tau, Ap, t_p)
-        return A / (1j * omega_func(time, Omega_0, Omega_QNM, tau, t_p, m))
+        A = A_func(time, tau, Ap, tp)
+        return A / (1j * omega_func(time, Omega_0, Omega_QNM, tau, tp, m))
 
     news_series_term_funcs = _build_symbolic_series(f0_psi4_func, g_func, N)
     
@@ -484,7 +484,7 @@ def news_sum_func(time):
     def f0_strain_func(time):
         # The "amplitude" is the full sum from the previous stage
         A_news = news_sum_func(time)
-        return A_news / (1j * omega_func(time, Omega_0, Omega_QNM, tau, t_p, m))
+        return A_news / (1j * omega_func(time, Omega_0, Omega_QNM, tau, tp, m))
 
     strain_series_term_funcs = _build_symbolic_series(f0_strain_func, g_func, M)
 
@@ -500,12 +500,12 @@ def f0_strain_func(time):
     return all_raw_terms_for_strain,strain_sum
 
 @partial(jit, static_argnames=('omega_func', 'A_func', 'N'))
-def calculate_strain_from_psi4_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p,
-                                            t_0,omega_func, A_func, m, N):
+def calculate_strain_from_psi4_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, tp,
+                                            t0,omega_func, A_func, m, N):
     '''
     Calculate the strain from psi4 using the series aapproximation
     Generates the raw, unsigned series terms [f₀, Df₀, D²f₀, ..., Dⁿf₀]
-    using JAX's automatic differentiation for scenarios with finite t_0.
+    using JAX's automatic differentiation for scenarios with finite t0.
 
     Args:
         t : Time 
@@ -513,7 +513,7 @@ def calculate_strain_from_psi4_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p,
         Omega_QNM : Real part of Quasinormal mode (QNM) frequency/(mode number)
         tau : Damping time; inverse of the imaginary part of the QNM frequency
         Ap : Peak waveform amplitude 
-        t_p : Time of peak amplitude
+        tp : Time of peak amplitude
         omega_func: frequency function
         A_func: amplitude function  
         m : Mode number
@@ -527,13 +527,13 @@ def calculate_strain_from_psi4_finite_t0(t, Omega_0, Omega_QNM, tau, Ap, t_p,
     M = N
     # --- Define the single D operator pre-factor ONCE ---
     def g_func(time):
-        omega = omega_func(time, Omega_0, Omega_QNM, tau, t_0, t_p, m)
+        omega = omega_func(time, Omega_0, Omega_QNM, tau, t0, tp, m)
         return 1.0 / (1j * omega)
 
     # --- Stage 1: Build the symbolic function for the News sum ---
     def f0_psi4_func(time):
-        A = A_func(time, tau, Ap, t_p)
-        return A / (1j * omega_func(time, Omega_0, Omega_QNM, tau, t_0, t_p, m))
+        A = A_func(time, tau, Ap, tp)
+        return A / (1j * omega_func(time, Omega_0, Omega_QNM, tau, t0, tp, m))
 
     news_series_term_funcs = _build_symbolic_series(f0_psi4_func, g_func, N)
     
@@ -546,7 +546,7 @@ def news_sum_func(time):
     def f0_strain_func(time):
         # The "amplitude" is the full sum from the previous stage
         A_news = news_sum_func(time)
-        return A_news / (1j * omega_func(time, Omega_0, Omega_QNM, tau, t_0, t_p, m))
+        return A_news / (1j * omega_func(time, Omega_0, Omega_QNM, tau, t0, tp, m))
 
     strain_series_term_funcs = _build_symbolic_series(f0_strain_func, g_func, M)
 
diff --git a/gwBOB/BOB_utils.py b/gwBOB/BOB_utils.py
index f0fdfbe..54dc3ca 100644
--- a/gwBOB/BOB_utils.py
+++ b/gwBOB/BOB_utils.py
@@ -627,9 +627,9 @@ def get_correct_Phi_and_Omega(self):
             ValueError: If ``what_should_BOB_create`` is not one of the
                 supported flavors.
         '''
-        #Even in the cases of strain_using_news, we still want to use the news frequency in all of the Omega0 optimizations because the analytical news frequency term
+        #Even in the cases of strain_using_news, we still want to use the news frequency in all of the Omega_0 optimizations because the analytical news frequency term
         #is built assuming the BOB amplitude best describes the news. While in principle, the accuracy could be improved for strain_using_news (and all X_using_Y cases)
-        #by optimizing Omega0 against the NR strain frequency, this would be unphysical.
+        #by optimizing Omega_0 against the NR strain frequency, this would be unphysical.
         if('psi4' in self._wf_config.what_to_create):
             if(self.minf_t0 is True):
                 Phi,Omega = BOB_terms.BOB_psi4_phase(self)
@@ -829,7 +829,7 @@ def fit_Omega0(self):
             ValueError: If ``minf_t0`` is False (use ``fit_t0`` for finite t0).
         '''
         if(self.minf_t0 is False):
-            raise ValueError("You are setup for a finite t0 right now. Omega0 fitting is only defined for t0 = infinity.")
+            raise ValueError("You are setup for a finite t0 right now. Omega_0 fitting is only defined for t0 = infinity.")
         if(self._wf_config.end_after_tpeak<self.end_fit_after_tpeak):
             logger.warning("end_after_tpeak is less than end_fit_after_tpeak. Setting end_fit_after_tpeak to end_after_tpeak")
             self.end_fit_after_tpeak = self._wf_config.end_after_tpeak
@@ -921,7 +921,7 @@ def construct_BOB_finite_t0(self,N):
         '''
         #Perform parameter sanity checks
         if(self.optimize_Omega0):
-            raise ValueError("Cannot optimize Omega0 for finite t0 values.")
+            raise ValueError("Cannot optimize Omega_0 for finite t0 values.")
 
         if(self._fit_config.optimize_t0_and_Omega0):
             self.fit_t0_and_Omega0()
@@ -1364,11 +1364,11 @@ def initialize_with_sxs_data(self,sxs_id,l=2,m=2,download=True,resample_dt = 0.1
 
         self.chif = np.linalg.norm(self.chif)
         self.chif_with_sign = self.chif*sign
-        self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.chif_with_sign,self.mf))
+        self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.mf, self.chif_with_sign))
         self.Omega_0 = self.Omega_ISCO
         self.l = l
         self.m = m
-        w_r,tau = gen_utils.get_qnm(self.chif,self.mf,self.l,np.abs(self.m),n=0,sign=sign)
+        w_r,tau = gen_utils.get_qnm(self.mf, self.chif, self.l, np.abs(self.m), n=0, sign=sign)
         self.w_r = np.abs(w_r)
         self.tau = np.abs(tau)
         self.Omega_QNM = self.w_r/np.abs(self.m)
@@ -1524,12 +1524,12 @@ def initialize_with_cce_data(self,cce_id,l=2,m=2,perform_superrest_transformatio
         sign = np.sign(self.chif[2])
         self.chif = np.linalg.norm(self.chif)
         self.chif_with_sign = self.chif*sign
-        self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.chif_with_sign,self.mf))
+        self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.mf, self.chif_with_sign))
         self.Omega_0 = self.Omega_ISCO
         self.l = l
         self.m = m
         
-        w_r,tau = gen_utils.get_qnm(self.chif,self.mf,self.l,np.abs(self.m),n=0,sign=sign)
+        w_r,tau = gen_utils.get_qnm(self.mf, self.chif, self.l, np.abs(self.m), n=0, sign=sign)
         self.w_r = np.abs(w_r)
         self.tau = np.abs(tau)
         self.Omega_QNM = self.w_r/np.abs(self.m)
@@ -1665,12 +1665,12 @@ def initialize_with_NR_mode(self,t_NR,y_psi4,y_strain,mf,chif,l=2,m=2,w_r = -1,t
         self.l = l
         self.m = m
         
-        self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.chif_with_sign,self.mf))
+        self.Omega_ISCO = np.abs(gen_utils.get_Omega_isco(self.mf, self.chif_with_sign))
         self.Omega_0 = self.Omega_ISCO
         
         if(w_r<0 or tau<0):
             logger.info("Calculating Kerr QNM parameters from provided Mf and chif")
-            w_r,tau = gen_utils.get_qnm(self.chif,self.mf,self.l,np.abs(self.m),n=0,sign=sign)
+            w_r,tau = gen_utils.get_qnm(self.mf, self.chif, self.l, np.abs(self.m), n=0, sign=sign)
             self.w_r = np.abs(w_r)
             self.tau = np.abs(tau)
         else:
diff --git a/gwBOB/gen_utils.py b/gwBOB/gen_utils.py
index 39f248e..19e332d 100644
--- a/gwBOB/gen_utils.py
+++ b/gwBOB/gen_utils.py
@@ -149,16 +149,16 @@ def get_frequency(ts):
     if(freq.y[find_nearest_index(freq.t,tp)]<0):
         freq.y = -freq.y
     return kuibit_ts(ts.t,freq.y)
-def get_r_isco(chi_with_sign,M):
+def get_r_isco(mf, chif_with_sign):
     '''
-    Compute the prograde ISCO radius for a Kerr black hole. For negative chi_with_sign values, we return the retrograde value
-    
+    Compute the prograde ISCO radius for a Kerr black hole. For negative chif_with_sign values, we return the retrograde value
+
     Parameters
     ----------
-    chi_with_sign : float
-        Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.
-    M : float
+    mf : float
         Mass of the remnant (in geometric units).
+    chif_with_sign : float
+        Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.
 
     Returns
     -------
@@ -167,29 +167,29 @@ def get_r_isco(chi_with_sign,M):
     '''
     #Bardeen Press Teukolskly eq 2.21
     #defined for prograde orbits
-    if(chi_with_sign>=0):
+    if(chif_with_sign>=0):
         sign = 1.0
     else:
         sign = -1.0
-    chi = np.abs(chi_with_sign)
-    a = chi*M
-    a_M = a/M
+    chi = np.abs(chif_with_sign)
+    a = chi*mf
+    a_M = a/mf
 
     z1 = 1 + (((1-a_M**2)**(1./3.)) * ((1+a_M)**(1./3.) + (1-a_M)**(1./3.))) #good
     z2 = (3*(a_M**2) + z1**2)**0.5 #good
-    r_isco = M * (3 + z2 - (sign)*((3-z1)*(3+z1+2*z2))**0.5) #good
+    r_isco = mf * (3 + z2 - (sign)*((3-z1)*(3+z1+2*z2))**0.5) #good
     return r_isco
-def get_Omega_isco(chi_with_sign,M):
+def get_Omega_isco(mf, chif_with_sign):
     '''
     Compute the orbital angular velocity at the ISCO for a Kerr black hole.
-    For negative chi_with_sign values, we return the retrograde value
+    For negative chif_with_sign values, we return the retrograde value
 
     Parameters
     ----------
-    chi_with_sign : float
-        Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.
-    M : float
+    mf : float
         Mass of the remnant (in geometric units).
+    chif_with_sign : float
+        Dimensionless spin of the remnant. Negative value indicates a spin pointed opposite the direction of the initial anular momentum of the binary.
 
     Returns
     -------
@@ -198,26 +198,26 @@ def get_Omega_isco(chi_with_sign,M):
     '''
     #Bardeen Press Teukolskly eq 2.16
     #defined for prograde orbits
-    if(chi_with_sign>=0):
+    if(chif_with_sign>=0):
         sign = 1.0
     else:
         sign = -1.0
-    chi = np.abs(chi_with_sign)
-    r_isco = get_r_isco(chi_with_sign,M)
-    a = chi*M
-    Omega = sign*np.sqrt(M)/(r_isco**1.5 + sign*a*np.sqrt(M)) # = dphi/dt
+    chi = np.abs(chif_with_sign)
+    r_isco = get_r_isco(mf, chif_with_sign)
+    a = chi*mf
+    Omega = sign*np.sqrt(mf)/(r_isco**1.5 + sign*a*np.sqrt(mf)) # = dphi/dt
     return Omega
 
-def get_qnm(chif,Mf,l,m,n=0,sign=1):
+def get_qnm(mf, chif, l, m, n=0, sign=1):
     '''
     Get the fundamental quasinormal mode frequency components for a Kerr black hole.
 
     Parameters
     ----------
+    mf : float
+        Final black-hole mass (geometric units).
     chif : float
         Dimensionless final spin magnitude (chi).
-    Mf : float
-        Final black-hole mass (geometric units).
     l : int
         Spherical-harmonic index ``l``.
     m : int
@@ -232,7 +232,6 @@ def get_qnm(chif,Mf,l,m,n=0,sign=1):
     tuple[float, float]
         ``(w_r, tau)`` where ``w_r`` is the real angular frequency and ``tau`` is the damping time.
     '''
-    #omega_qnm, all_C, ells = qnmfits.read_qnms.qnm_from_tuple((l,m,n,1),chif,M=M)
     if(sign==-1):
         grav_lmn = qnm.modes_cache(s=-2,l=l,m=-m,n=n)
         omega_qnm, A, C = grav_lmn(a=chif)#qnm package uses M = 1 so a = chi here
@@ -240,7 +239,7 @@ def get_qnm(chif,Mf,l,m,n=0,sign=1):
     else:
         grav_lmn = qnm.modes_cache(s=-2,l=l,m=m,n=n)
         omega_qnm, A, C = grav_lmn(a=chif)#qnm package uses M = 1 so a = chi here
-    omega_qnm /= Mf #rescale to remnant black hole mass
+    omega_qnm /= mf #rescale to remnant black hole mass
     w_r = np.abs(omega_qnm.real)
     imag_qnm = np.abs(omega_qnm.imag)
     tau = 1./imag_qnm
@@ -454,7 +453,7 @@ def estimate_parameters(BOB,
     tf : float, optional
         End time after peak for mismatch window, by default 75.
     force_Omega0_optimization : bool, optional
-        If True, enforce Omega0 optimization during construction, by default False.
+        If True, enforce Omega_0 optimization during construction, by default False.
     NR_data : kuibit_ts, optional
         Reference data to compare against. If None, inferred from ``BOB`` selection.
     make_current_naturally : bool, optional
@@ -462,9 +461,9 @@ def estimate_parameters(BOB,
     make_mass_naturally : bool, optional
         If True, construct mass quadrupole via natural method, by default False.
     include_Omega0_as_parameter : bool, optional
-        Include Omega0 as an optimization parameter, by default False.
+        Include Omega_0 as an optimization parameter, by default False.
     include_2Omega0_as_parameters : bool, optional
-        Include both lm and lmm Omega0 as parameters (quadrupole builds), by default False.
+        Include both lm and lmm Omega_0 as parameters (quadrupole builds), by default False.
     perform_phase_alignment_first : bool, optional
         If True, perform phase alignment before quadrupole combination, by default False.
     start_with_wide_search : bool, optional
@@ -545,7 +544,7 @@ def create_guess(x):
         if(include_Omega0_as_parameter):
             #keep this for ordinary (l,m) &(l,-m) modes
             BOB.Omega_0 = lm_Omega0_guess
-        w_r,tau = get_qnm(BOB.chif,BOB.mf,BOB.l,np.abs(BOB.m),sign=np.sign(BOB.chif_with_sign))
+        w_r,tau = get_qnm(BOB.mf, BOB.chif, BOB.l, np.abs(BOB.m), sign=np.sign(BOB.chif_with_sign))
         BOB.Omega_QNM = w_r/np.abs(BOB.m)
         BOB.Phi_0 = 0
         BOB.tau = tau
@@ -922,13 +921,13 @@ def load_lower_lev_SXS(sim):
     else:
         raise ValueError("only one Level found")
     return sim_lower
-def Omega_0_fit_psi4(Mf,chif_with_sign):
+def Omega_0_fit_psi4(mf, chif_with_sign):
     '''
-    Omega0 for psi4 using the fit from Kankani and McWilliams (2025)
-    
+    Omega_0 for psi4 using the fit from Kankani and McWilliams (2025)
+
     Parameters
     ----------
-    Mf : float
+    mf : float
         Remnant mass.
     chif_with_sign : float
         Remnant spin. (negative values indicate a final spin pointing opposite to the initial orbital angular momentum)
@@ -937,19 +936,18 @@ def Omega_0_fit_psi4(Mf,chif_with_sign):
     -------
     float
         Omega_0 fit.
-    
     '''
     A = 1.42968337
     B = 0.08424419
     C = -1.22848524
-    return A*Mf + B*chif_with_sign + C
-def Omega_0_fit_news(Mf,chif_with_sign):
+    return A*mf + B*chif_with_sign + C
+def Omega_0_fit_news(mf, chif_with_sign):
     '''
-    Omega0 for news using the fit from Kankani and McWilliams (2025)
+    Omega_0 for news using the fit from Kankani and McWilliams (2025)
 
     Parameters
     ----------
-    Mf : float
+    mf : float
         Remnant mass.
     chif_with_sign : float
         Remnant spin. (negative values indicate a final spin pointing opposite to the initial orbital angular momentum)
@@ -958,19 +956,18 @@ def Omega_0_fit_news(Mf,chif_with_sign):
     -------
     float
         Omega_0 fit.
-    
     '''
     A = 0.33568227
     B = 0.03450997
     C = -0.18763176
-    return A*Mf + B*chif_with_sign + C
-def Omega_0_fit_strain(Mf,chif_with_sign):
+    return A*mf + B*chif_with_sign + C
+def Omega_0_fit_strain(mf, chif_with_sign):
     '''
-    Omega0 for strain using the fit from Kankani and McWilliams (2025)
+    Omega_0 for strain using the fit from Kankani and McWilliams (2025)
 
     Parameters
     ----------
-    Mf : float
+    mf : float
         Remnant mass.
     chif_with_sign : float
         Remnant spin. (negative values indicate a final spin pointing opposite to the initial orbital angular momentum)
@@ -979,12 +976,11 @@ def Omega_0_fit_strain(Mf,chif_with_sign):
     -------
     float
         Omega_0 fit.
-    
     '''
     A = 0.01663248
     B = 0.01798275
     C = 0.07882578
-    return A*Mf + B*chif_with_sign + C
+    return A*mf + B*chif_with_sign + C
     
 
 
diff --git a/tests/unit/test_BOB_orchestrator.py b/tests/unit/test_BOB_orchestrator.py
index 140975f..71b447b 100644
--- a/tests/unit/test_BOB_orchestrator.py
+++ b/tests/unit/test_BOB_orchestrator.py
@@ -187,7 +187,7 @@ def test_construct_BOB_finite_t0_rejects_optimize_Omega0(self):
         bob = BOB_utils.BOB()
         bob.optimize_Omega0 = True
         # Force the finite-t0 branch by direct call (skipping the setter)
-        with pytest.raises(ValueError, match="Cannot optimize Omega0 for finite t0"):
+        with pytest.raises(ValueError, match="Cannot optimize Omega_0 for finite t0"):
             bob.construct_BOB_finite_t0(N=1)
 
     def test_setting_set_initial_time_before_what_to_create_raises(self):
diff --git a/tests/unit/test_gen_utils_math.py b/tests/unit/test_gen_utils_math.py
index be17484..d290162 100644
--- a/tests/unit/test_gen_utils_math.py
+++ b/tests/unit/test_gen_utils_math.py
@@ -16,34 +16,34 @@
 
 
 def test_get_r_isco_values():
-    """ISCO radius for representative (chi, M) pairs.
+    """ISCO radius for representative (mf, chif) pairs.
 
-    Reference values: chi=0 is Schwarzschild (r_isco = 6M); other values
+    Reference values: chif=0 is Schwarzschild (r_isco = 6M); other values
     were computed once and treated as the trusted reference.
     """
-    chi_vals = np.array([0.0, 0.5, 0.9])
-    M_vals = np.array([1.0, 2.0, 5.0])
+    mf_vals   = np.array([1.0, 2.0, 5.0])
+    chif_vals = np.array([0.0, 0.5, 0.9])
     expected = [
-        6.0,                  # (chi=0,   M=1) -> Schwarzschild ISCO = 6M
-        8.466005059061652,    # (chi=0.5, M=2)
-        11.604415208809435,   # (chi=0.9, M=5)
+        6.0,                  # (mf=1, chif=0)   -> Schwarzschild ISCO = 6M
+        8.466005059061652,    # (mf=2, chif=0.5)
+        11.604415208809435,   # (mf=5, chif=0.9)
     ]
-    for chi, M, exp in zip(chi_vals, M_vals, expected):
-        result = gen_utils.get_r_isco(chi, M)
+    for mf, chif, exp in zip(mf_vals, chif_vals, expected):
+        result = gen_utils.get_r_isco(mf, chif)
         assert np.isclose(result, exp, rtol=1e-8)
 
 
 def test_get_Omega_isco_values():
-    """ISCO orbital frequency for representative (chi, M) pairs."""
-    chi_vals = np.array([0.0, 0.5, 0.9])
-    M_vals = np.array([1.0, 2.0, 5.0])
+    """ISCO orbital frequency for representative (mf, chif) pairs."""
+    mf_vals   = np.array([1.0, 2.0, 5.0])
+    chif_vals = np.array([0.0, 0.5, 0.9])
     expected = [
         0.06804138174397717,
         0.05429417949013838,
         0.0450883417670616,
     ]
-    for chi, M, exp in zip(chi_vals, M_vals, expected):
-        result = gen_utils.get_Omega_isco(chi, M)
+    for mf, chif, exp in zip(mf_vals, chif_vals, expected):
+        result = gen_utils.get_Omega_isco(mf, chif)
         assert np.isclose(result, exp, rtol=1e-8)
 
 
@@ -52,8 +52,8 @@ def test_get_qnm():
 
     Spans ell=2,3, n=0,1, retrograde and prograde branches.
     """
-    chi_vals  = np.array([0.0, 0.0, 0.0, 0.5, 0.5, 0.5])
-    M_vals    = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 2.0])
+    mf_vals   = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 2.0])
+    chif_vals = np.array([0.0, 0.0, 0.0, 0.5, 0.5, 0.5])
     l_vals    = np.array([2,   3,   2,   2,   2,   2])
     m_vals    = np.array([2,   2,   2,   2,   2,   2])
     n_vals    = np.array([0,   0,   1,   0,   0,   0])
@@ -68,10 +68,10 @@ def test_get_qnm():
         11.231973996651769, 11.676945396785948, 23.353890793571896,
     ])
 
-    for chi, M, l, m, n, sgn, exp_w, exp_tau in zip(
-        chi_vals, M_vals, l_vals, m_vals, n_vals, sign_vals,
+    for mf, chif, l, m, n, sgn, exp_w, exp_tau in zip(
+        mf_vals, chif_vals, l_vals, m_vals, n_vals, sign_vals,
         expected_w_r_vals, expected_tau_vals,
     ):
-        result_w, result_tau = gen_utils.get_qnm(chi, M, l, m, n=n, sign=sgn)
+        result_w, result_tau = gen_utils.get_qnm(mf, chif, l, m, n=n, sign=sgn)
         assert np.isclose(result_w,   exp_w,   rtol=1e-8)
         assert np.isclose(result_tau, exp_tau, rtol=1e-8)