diff --git a/cubespecfit/__init__.py b/cubespecfit/__init__.py
new file mode 100644
index 0000000..d6fc94d
--- /dev/null
+++ b/cubespecfit/__init__.py
@@ -0,0 +1,2 @@
+from .fitting import *
+from .models import *
\ No newline at end of file
diff --git a/cubespecfit/fitting.py b/cubespecfit/fitting.py
index 0aa3fb4..91da496 100644
--- a/cubespecfit/fitting.py
+++ b/cubespecfit/fitting.py
@@ -1,3 +1,5 @@
+__all__ = ['fitspec_lmfit','bic','fitcube','fitcube_two_model']
+
 import numpy as np
 import matplotlib
 from matplotlib import pyplot as plt
@@ -7,7 +9,7 @@
 import logging
 import warnings
 
-from cubespecfit.wcs_helpers import _save_results
+from .wcs_helpers import _save_results
 
 # Suppress specific numpy warnings globally
 warnings.filterwarnings('ignore', message='Mean of empty slice')
@@ -119,7 +121,6 @@ def constrained_residual(params, x, data, weights):
         # Reconstruct model result
         param_vals = [result.params[p].value for p in params_dict.keys()]
         bestfit = model(wave, *param_vals)
-        chi2 = np.sum((flux - bestfit) ** 2 / err ** 2)
 
         # Extract results - only varying parameters
         varying_params = [p for p in params_dict.keys() if result.params[p].vary]
@@ -127,7 +128,6 @@ def constrained_residual(params, x, data, weights):
         param_errors = np.array([result.params[p].stderr if result.params[p].stderr is not None
                                  else np.nan for p in varying_params])
 
-        result_obj = result
 
     else:
         # Standard unconstrained fit
@@ -140,8 +140,11 @@ def constrained_residual(params, x, data, weights):
                                  else np.nan for p in varying_params])
 
         bestfit = result.best_fit
-        chi2 = result.chisqr
-        result_obj = result
+
+    # NOTE lmfit doesn't appear to return what we think of as the "standard" chisquared
+    #chi2 = result.chisqr
+    chi2 = np.sum((flux - bestfit) ** 2 / err ** 2)
+    result_obj = result
 
     return param_values, param_errors, bestfit, chi2, result_obj
 
diff --git a/cubespecfit/models.py b/cubespecfit/models.py
index 3772356..46308d7 100644
--- a/cubespecfit/models.py
+++ b/cubespecfit/models.py
@@ -1,3 +1,6 @@
+__all__ = ['ha_nii_model','ha_nii_broad_model','ha_nii_sii_model','sii_model','oiii_model',
+    'oiii_broad_model','oiii_hbeta_model','oiii_hbeta_broad_model','hbeta_model']
+
 import numpy as np
 
 # ------------------------------------------------------------------------------------------------------------
diff --git a/cubespecfit/wcs_helpers.py b/cubespecfit/wcs_helpers.py
index ab95d09..a43729b 100644
--- a/cubespecfit/wcs_helpers.py
+++ b/cubespecfit/wcs_helpers.py
@@ -1,3 +1,5 @@
+__all__ = ['_propagate_celestial_wcs','_save_results']
+
 import os
 import numpy as np
 from astropy.io import fits