diff --git a/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb b/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb
index cff661f3a..2682b28b1 100644
--- a/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb
+++ b/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb
@@ -5,8 +5,8 @@
    "metadata": {},
    "source": [
     "# BOTS Time Series Observations\n",
-    "\n ",
-    "**Use case:** Bright Object Time Series; extracting exoplanet spectra.<br>\n",
+    "\n",
+    " **Use case:** Bright Object Time Series; extracting exoplanet spectra.<br>\n",
     "**Data:** JWST simulated NIRSpec data from ground-based campaign; GJ436b spectra from the Goyal et al. (2018).<br>\n",
     "**Tools:**  scikit, lmfit, scipy, matplotlip, astropy, pandas.<br>\n",
     "**Cross-intrument:** . <br>\n",
@@ -61,22 +61,20 @@
     "import matplotlib.pyplot as plt\n",
     "import matplotlib.ticker as ticker\n",
     "from matplotlib.backends.backend_pdf import PdfPages\n",
-    "from astropy.utils.data import get_pkg_data_filename, download_file\n",
-    "from astropy.table import Table, Column, MaskedColumn\n",
+    "from astropy.utils.data import download_file\n",
+    "from astropy.table import Table\n",
     "from astropy.io import fits, ascii\n",
     "from astropy.modeling.models import custom_model\n",
     "from astropy.modeling.fitting import LevMarLSQFitter\n",
     "import astropy.units as u\n",
     "from scipy.interpolate import interp1d, splev, splrep\n",
-    "import scipy.optimize as opt\n",
     "from scipy.io import readsav\n",
     "from scipy import stats\n",
     "import glob\n",
     "import lmfit\n",
     "import pickle\n",
-    "from os import path,mkdir\n",
+    "from os import path, mkdir\n",
     "from sklearn.linear_model import LinearRegression\n",
-    "import warnings\n",
     "import pandas as pd\n",
     "import os\n",
     "import shutil"
@@ -97,39 +95,37 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#---------------------------------------------------------\n",
     "# SETUP  ----------------------------------------------\n",
     "# Setup directories\n",
-    "save_directory      ='./notebookrun2/'          #local directory to save files to\n",
-    "data_directory      ='./'                       #local data to work with fits files if desired\n",
+    "save_directory = './notebookrun2/'          # Local directory to save files to\n",
+    "data_directory = './'                       # Local data to work with fits files if desired\n",
     "\n",
     "# Setup Detector Properties & Rednoise measurement timescale \n",
-    "gain               = 1.0  # 2D spectra has already converted to counts, gain of detector is 1.0\n",
-    "binmeasure         = 256  # Binning technique to measure rednoise, choose bin size to evaluate sigma_r\n",
-    "number_of_images   = 8192 # Number of images in the dataset\n",
-    "\n",
-    "# Setup Planet Properies\n",
-    "grating    = 'NIRSpecPrism'\n",
-    "ld_model   = '3D'         # 3D/1D stellar model choise (transit was injected with the 3D model)\n",
-    "\n",
-    "# Setup Stellar Properies for Limb-Darkening Calculation\n",
-    "Teff     = 4500           # Effective Temperature (K)\n",
-    "logg     = 4.5            # Surface Gravity\n",
-    "M_H      = 0.0            # Stellar Metallicity log_10[M/H]\n",
-    "Rstar    = 0.455          # planet radius (in units of solar radii Run)\n",
-    "\n",
-    "#Setup Transit parameters (can get from Nasa exoplanet archive)\n",
-    "t0       = 2454865.084034              # bjd time of inferior conjunction \n",
-    "per      = 2.64389803                  # orbital period (days) BJD_TDB\n",
-    "rp       = 0.0804                      # planet radius (in units of stellar radii)\n",
-    "a_Rs     = 14.54                       # semi-major axis (input a/Rstar so units of stellar radii)\n",
-    "inc      = 86.858  *(2*np.pi/360)      # orbital inclination (in degrees->radians)\n",
-    "ecc      = 0.0                         # eccentricity\n",
-    "omega    = 0.0  *(2*np.pi/360)         # longitude of periastron (in degrees->radians)\n",
+    "gain = 1.0   # 2D spectra has already converted to counts, gain of detector is 1.0\n",
+    "binmeasure = 256   # Binning technique to measure rednoise, choose bin size to evaluate sigma_r\n",
+    "number_of_images = 8192  # Number of images in the dataset\n",
+    "\n",
+    "# Setup Planet Properties\n",
+    "grating = 'NIRSpecPrism'\n",
+    "ld_model = '3D'      # 3D/1D stellar model choice (transit was injected with the 3D model)\n",
+    "\n",
+    "# Setup Stellar Properties for Limb-Darkening Calculation\n",
+    "Teff = 4500           # Effective Temperature (K)\n",
+    "logg = 4.5            # Surface Gravity\n",
+    "M_H = 0.0            # Stellar Metallicity log_10[M/H]\n",
+    "Rstar = 0.455          # Planet radius (in units of solar radii Run)\n",
+    "\n",
+    "# Setup Transit parameters (can get from NASA exoplanet archive)\n",
+    "t0 = 2454865.084034              # bjd time of inferior conjunction \n",
+    "per = 2.64389803                  # orbital period (days) BJD_TDB\n",
+    "rp = 0.0804                      # Planet radius (in units of stellar radii)\n",
+    "a_Rs = 14.54                       # Semi-major axis (input a/Rstar so units of stellar radii)\n",
+    "inc = 86.858 * (2*np.pi/360)      # Orbital inclination (in degrees->radians)\n",
+    "ecc = 0.0                         # Eccentricity\n",
+    "omega = 0.0 * (2*np.pi/360)         # Longitude of periastron (in degrees->radians)\n",
     "\n",
     "rho_star = (3*np.pi)/(6.67259E-8*(per*86400)**2)*(a_Rs)**3     # Stellar Density (g/cm^3) from a/Rs\n",
-    "# a_Rs=(rho_star*6.67259E-8*per_sec*per_sec/(3*np.pi))**(1/3)  # a/Rs from Stellar Density (g/cm^3)\n",
-    "#---------------------------------------------------------\n"
+    "# a_Rs=(rho_star*6.67259E-8*per_sec*per_sec/(3*np.pi))**(1/3)  # a/Rs from Stellar Density (g/cm^3)"
    ]
   },
   {
@@ -139,9 +135,11 @@
    "outputs": [],
    "source": [
     "# Create local directories\n",
-    "if path.exists(save_directory) == False: mkdir(save_directory)      #Create new directory to save outputs to if needed\n",
-    "if path.exists(save_directory+'3DGrid') == False: mkdir(save_directory+'3DGrid')      #Create new directory to save\n",
-    "limb_dark_directory = save_directory               # point to limb darkeing directory contaning 3DGrid/ directory"
+    "if not path.exists(save_directory):\n",
+    "    mkdir(save_directory)      # Create a new directory to save outputs to if needed\n",
+    "if not path.exists(save_directory+'3DGrid'): \n",
+    "    mkdir(save_directory+'3DGrid')      # Create new directory to save\n",
+    "limb_dark_directory = save_directory    # Point to limb darkeing directory contaning 3DGrid/ directory"
    ]
   },
   {
@@ -168,19 +166,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Download 1GB NIRSpec Data\n",
-    "fn_jw=save_directory+'jwst_data.pickle'\n",
-    "if path.exists(fn_jw) == False:\n",
-    "    fn   = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/jwst_data.pickle')\n",
+    "# Download 1GB NIRSpec Data\n",
+    "fn_jw = save_directory+'jwst_data.pickle'\n",
+    "if not path.exists(fn_jw):\n",
+    "    fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/jwst_data.pickle')\n",
     "    dest = shutil.move(fn,save_directory+'jwst_data.pickle')  \n",
     "    print('JWST Data Download Complete')\n",
     "\n",
-    "#Download further files needed, move to local directory for easier repeated access\n",
-    "fn_sens= download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/NIRSpec.prism.sensitivity.sav')\n",
-    "dest   = shutil.move(fn_sens,save_directory+'NIRSpec.prism.sensitivity.sav')  #Move files to save_directory\n",
+    "# Download further files needed, move to local directory for easier repeated access\n",
+    "fn_sens = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/NIRSpec.prism.sensitivity.sav')\n",
+    "dest = shutil.move(fn_sens, save_directory+'NIRSpec.prism.sensitivity.sav')  # Move files to save_directory\n",
     "\n",
-    "fn_ld  = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/3DGrid/mmu_t45g45m00v05.flx')\n",
-    "destld = shutil.move(fn_ld,save_directory+'3DGrid/mmu_t45g45m00v05.flx')        "
+    "fn_ld = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/3DGrid/mmu_t45g45m00v05.flx')\n",
+    "destld = shutil.move(fn_ld, save_directory+'3DGrid/mmu_t45g45m00v05.flx')        "
    ]
   },
   {
@@ -196,79 +194,74 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "if path.exists(fn_jw) == True:\n",
+    "if path.exists(fn_jw):\n",
     "    dbfile = open(fn_jw, 'rb') # for reading also binary mode is important\n",
     "    jwst_data = pickle.load(dbfile)\n",
     "    print('Loading JWST data from Pickle File')\n",
-    "    bjd            =jwst_data['bjd']\n",
-    "    wsdata_all     =jwst_data['wsdata_all']\n",
-    "    shx            =jwst_data['shx']\n",
-    "    shy            =jwst_data['shy']\n",
-    "    common_mode    =jwst_data['common_mode']\n",
-    "    all_spec       =jwst_data['all_spec']\n",
-    "    exposure_length=jwst_data['exposure_length']\n",
+    "    bjd = jwst_data['bjd']\n",
+    "    wsdata_all =jwst_data['wsdata_all']\n",
+    "    shx = jwst_data['shx']\n",
+    "    shy = jwst_data['shy']\n",
+    "    common_mode = jwst_data['common_mode']\n",
+    "    all_spec = jwst_data['all_spec']\n",
+    "    exposure_length = jwst_data['exposure_length']\n",
     "    dbfile.close()    \n",
     "    print('Done')\n",
-    "elif path.exists(fn_jw) == False:\n",
-    "    #---------------------------------------------------------\n",
-    "    # load all fits images\n",
+    "elif not path.exists(fn_jw):\n",
+    "    # Load all fits images\n",
     "    # Arrays created for BJD time, and the white light curve total_counts\n",
-    "    list=glob.glob(data_directory+\"*.fits\")\n",
-    "    index_of_images=np.arange(number_of_images) #\n",
-    "    bjd=np.zeros((number_of_images))\n",
-    "    exposure_length=np.zeros((number_of_images))\n",
-    "    all_spec=np.zeros((32,512,number_of_images))\n",
+    "    list = glob.glob(data_directory+\"*.fits\")\n",
+    "    index_of_images = np.arange(number_of_images) \n",
+    "    bjd = np.zeros((number_of_images))\n",
+    "    exposure_length = np.zeros((number_of_images))\n",
+    "    all_spec = np.zeros((32, 512, number_of_images))\n",
     "    for i in index_of_images:\n",
-    "        img=list[i]\n",
+    "        img = list[i]\n",
     "        print(img)\n",
-    "        hdul=fits.open(img)\n",
-    "        #hdul.info()\n",
-    "        bjd_image=hdul[0].header['BJD_TDB']\n",
-    "        BZERO=hdul[0].header['BZERO']        #flux value offset\n",
-    "        bjd[i]=bjd_image\n",
-    "        expleng=hdul[0].header['INTTIME']    #Total integration time for one MULTIACCUM (seconds)\n",
-    "        exposure_length[i]=expleng/86400.    #Total integration time for one MULTIACCUM (days)\n",
+    "        hdul = fits.open(img)\n",
+    "        # hdul.info()\n",
+    "        bjd_image = hdul[0].header['BJD_TDB']\n",
+    "        BZERO = hdul[0].header['BZERO']        # flux value offset\n",
+    "        bjd[i] = bjd_image\n",
+    "        expleng = hdul[0].header['INTTIME']    # Total integration time for one MULTIACCUM (seconds)\n",
+    "        exposure_length[i] = expleng/86400.    # Total integration time for one MULTIACCUM (days)\n",
     "        print(bjd[i])\n",
     "        data = hdul[0].data\n",
-    "        #total counts in image\n",
-    "        #total_counts[i]=gain*np.sum(data[11:18,170:200]-BZERO) #total counts in 12 pix wide aperature around pixel 60 image\n",
-    "        all_spec[:,:,i]=gain*(data-BZERO)     #load all spectra into an array subtract flux value offset\n",
+    "        # total counts in image\n",
+    "        # total_counts[i]=gain*np.sum(data[11:18,170:200]-BZERO) # total counts in 12 pix wide aperture around pixel 60 image\n",
+    "        all_spec[:, :, i] = gain * (data-BZERO)     # Load all spectra into an array subtract flux value offset\n",
     "        hdul.close()\n",
-    "    #Sort data\n",
-    "    srt=np.argsort(bjd) #index to sort\n",
-    "    bjd=bjd[srt]\n",
-    "    #total_counts=total_counts[srt]\n",
-    "    exposure_length=exposure_length[srt]\n",
-    "    all_spec[:,:,:]=all_spec[:,:,srt]\n",
-    "\n",
-    "    # Get Wavelegnth of Data\n",
+    "    # Sort data\n",
+    "    srt = np.argsort(bjd) # index to sort\n",
+    "    bjd = bjd[srt]\n",
+    "    # total_counts=total_counts[srt]\n",
+    "    exposure_length = exposure_length[srt]\n",
+    "    all_spec[:, :, :] = all_spec[:, :, srt]\n",
+    "\n",
+    "    # Get Wavelength of Data\n",
     "    file_wave = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/JWST_NIRSpec_wavelength_microns.txt')\n",
     "    f = open(file_wave, 'r')\n",
     "    wsdata_all = np.genfromtxt(f)\n",
     "    \n",
-    "    print('wsdata size :',wsdata_all.shape)\n",
-    "    print('Data wavelength Loaded :',wsdata_all)\n",
-    "    print('wsdata new size :',wsdata_all.shape)\n",
+    "    print('wsdata size :', wsdata_all.shape)\n",
+    "    print('Data wavelength Loaded :', wsdata_all)\n",
+    "    print('wsdata new size :', wsdata_all.shape)\n",
     "    \n",
-    "    #---------------------------------------------------------\n",
     "    # Read in Detrending parameters\n",
     "    # Mean of parameter must be 0.0 to be properly normalized\n",
-    "    # Idealy standard deviation of paramter = 1.0\n",
-    "    file_xy=download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/JWST_NIRSpec_Xposs_Yposs_CM_detrending.txt')\n",
+    "    # Idealy standard deviation of parameter = 1.0\n",
+    "    file_xy = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/JWST_NIRSpec_Xposs_Yposs_CM_detrending.txt')\n",
     "    f = open(file_xy, 'r')\n",
     "    data = np.genfromtxt(f, delimiter=',')\n",
-    "    shx        = data[:,0]\n",
-    "    shy        = data[:,1]\n",
-    "    common_mode= data[:,2]\n",
+    "    shx = data[:, 0]\n",
+    "    shy = data[:, 1]\n",
+    "    common_mode = data[:, 2]\n",
     "    \n",
-    "    #Store Data into a pickle file\n",
-    "    jwst_data={'bjd':bjd, 'wsdata_all':wsdata_all, 'shx':shx , 'shy':shy , 'common_mode':common_mode, 'all_spec':all_spec, 'exposure_length':exposure_length}\n",
+    "    # Store Data in a pickle file\n",
+    "    jwst_data = {'bjd':bjd, 'wsdata_all':wsdata_all, 'shx':shx , 'shy':shy , 'common_mode':common_mode, 'all_spec':all_spec, 'exposure_length':exposure_length}\n",
     "    dbfile = open('jwst_data.pickle', 'ab') # Its important to use binary mode\n",
-    "    pickle.dump(jwst_data,dbfile)\n",
-    "    dbfile.close()\n",
-    "    #---------------------------------------------------------\n",
-    "\n",
-    "    "
+    "    pickle.dump(jwst_data, dbfile)\n",
+    "    dbfile.close()"
    ]
   },
   {
@@ -284,34 +277,34 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "expnum=2                                           #Choose Exposure number to view\n",
+    "expnum = 2                                           # Choose Exposure number to view\n",
     "\n",
     "plt.rcParams['figure.figsize'] = [10.0, 3.0]           # Figure dimensions\n",
-    "plt.rcParams['figure.dpi']     = 200                   # Resolution\n",
-    "plt.rcParams['savefig.dpi']    = 200\n",
-    "plt.rcParams['image.aspect']   = 5                     # Aspect ratio (the CCD is quite long!!!)\n",
+    "plt.rcParams['figure.dpi'] = 200                   # Resolution\n",
+    "plt.rcParams['savefig.dpi'] = 200\n",
+    "plt.rcParams['image.aspect'] = 5                     # Aspect ratio (the CCD is quite long!!!)\n",
     "plt.cmap = plt.cm.magma\n",
     "plt.cmap.set_bad('k',1.)\n",
-    "plt.rcParams['image.cmap'] = 'magma'                   # Colormap.\n",
+    "plt.rcParams['image.cmap'] = 'magma'                   # Colormap\n",
     "plt.rcParams['image.interpolation'] = 'none'\n",
     "plt.rcParams['image.origin'] = 'lower'\n",
     "plt.rcParams['font.family'] = \"monospace\"\n",
     "plt.rcParams['font.monospace'] = 'DejaVu Sans Mono'\n",
     "\n",
-    "img=all_spec[:,:,expnum]\n",
-    "zeros=np.where(img <= 0)     #Plot on a log scale, so set zero or negative values to a small number \n",
-    "img[zeros]=1E-10\n",
-    "fig,axs = plt.subplots()\n",
-    "f=axs.imshow(np.log10(img),vmin=0) #Plot image\n",
+    "img = all_spec[:, :, expnum]\n",
+    "zeros = np.where(img <= 0)     # Plot on a log scale, so set zero or negative values to a small number \n",
+    "img[zeros] = 1E-10\n",
+    "fig, axs = plt.subplots()\n",
+    "f = axs.imshow(np.log10(img), vmin=0) # Plot image\n",
     "plt.xlabel('x-pixel')\n",
     "plt.ylabel('y-pixel')\n",
     "axs.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
     "axs.yaxis.set_minor_locator(ticker.MultipleLocator(1))\n",
     "axs.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
     "axs.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
-    "plt.title('2D NIRSpec Image of Exposure '+str(expnum))\n",
-    "fig.colorbar(f,label='Log$_{10}$ Electron counts',ax=axs)\n",
-    "plt.show()\n"
+    "plt.title('2D NIRSpec Image of Exposure ' + str(expnum))\n",
+    "fig.colorbar(f, label='Log$_{10}$ Electron counts', ax=axs)\n",
+    "plt.show()"
    ]
   },
   {
@@ -344,41 +337,39 @@
    "outputs": [],
    "source": [
     "all_spec.shape\n",
-    "y_lower = 11                                             # Lower extraction aperature\n",
-    "y_upper = 18                                             # Upper extraction aperature\n",
+    "y_lower = 11                                             # Lower extraction aperture\n",
+    "y_upper = 18                                             # Upper extraction aperture\n",
     "all_spec_1D=np.sum(all_spec[y_lower:y_upper,:,:],axis=0) # Sum along Y-axis from pixels 11 to 18\n",
-    "\n",
-    "#Plot \n",
+    "# Plot \n",
     "\n",
     "plt.rcParams['figure.figsize'] = [10.0, 3.0]           # Figure dimensions\n",
-    "plt.rcParams['figure.dpi']     = 200                   # Resolution\n",
-    "plt.rcParams['savefig.dpi']    = 200\n",
-    "plt.rcParams['image.aspect']   = 5                     # Aspect ratio (the CCD is quite long!!!)\n",
+    "plt.rcParams['figure.dpi'] = 200                   # Resolution\n",
+    "plt.rcParams['savefig.dpi'] = 200\n",
+    "plt.rcParams['image.aspect'] = 5                     # Aspect ratio (the CCD is quite long!!!)\n",
     "plt.cmap = plt.cm.magma\n",
-    "plt.cmap.set_bad('k',1.)\n",
-    "plt.rcParams['image.cmap'] = 'magma'                   # Colormap.\n",
+    "plt.cmap.set_bad('k', 1.)\n",
+    "plt.rcParams['image.cmap'] = 'magma'                   # Colormap\n",
     "plt.rcParams['image.interpolation'] = 'none'\n",
     "plt.rcParams['image.origin'] = 'lower'\n",
     "plt.rcParams['font.family'] = \"monospace\"\n",
     "plt.rcParams['font.monospace'] = 'DejaVu Sans Mono'\n",
     "\n",
-    "img=all_spec[:,:,expnum]\n",
-    "zeros=np.where(img <= 0)     #Plot on a log scale, so set zero or negitive values to a small number \n",
-    "img[zeros]=1E-10\n",
-    "fig,axs = plt.subplots()\n",
-    "f=axs.imshow(np.log10(img),vmin=0) #Plot image\n",
+    "img = all_spec[:, :, expnum]\n",
+    "zeros = np.where(img <= 0)     # Plot on a log scale, so set zero or negitive values to a small number \n",
+    "img[zeros] = 1E-10\n",
+    "fig, axs = plt.subplots()\n",
+    "f = axs.imshow(np.log10(img), vmin=0) #Plot image\n",
     "plt.xlabel('x-pixel')\n",
     "plt.ylabel('y-pixel')\n",
     "axs.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
     "axs.yaxis.set_minor_locator(ticker.MultipleLocator(1))\n",
     "axs.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
     "axs.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
-    "plt.axhline(y_lower, color = 'w', ls = 'dashed')\n",
-    "plt.axhline(y_upper, color = 'w', ls = 'dashed')\n",
+    "plt.axhline(y_lower, color='w', ls='dashed')\n",
+    "plt.axhline(y_upper, color='w', ls='dashed')\n",
     "plt.title('2D NIRSpec Image of Exposure '+str(expnum))\n",
-    "fig.colorbar(f,label='Log$_{10}$ Electron counts',ax=axs)\n",
-    "plt.show()\n",
-    "\n"
+    "fig.colorbar(f,label='Log$_{10}$ Electron counts', ax=axs)\n",
+    "plt.show()"
    ]
   },
   {
@@ -395,7 +386,7 @@
    "outputs": [],
    "source": [
     "fig,axs = plt.subplots()\n",
-    "f=plt.plot(wsdata_all,all_spec_1D[:,0], linewidth=2,zorder=0)  #overplot Transit model at data\n",
+    "f = plt.plot(wsdata_all,all_spec_1D[:, 0], linewidth=2, zorder=0)  # overplot Transit model at data\n",
     "plt.xlabel('Wavelength ($\\mu$m)')\n",
     "plt.ylabel('Flux (e-)')\n",
     "axs.xaxis.set_major_locator(ticker.MultipleLocator(0.5))\n",
@@ -429,18 +420,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#---------------------------------------------------------\n",
     "#Calculate Orbital Phase\n",
-    "phase=(bjd-t0)/(per)  #phase in days relative to T0 ephemeris\n",
-    "phase=phase-np.fix(phase[number_of_images-1]) # Have current phase occur at value 0.0\n",
+    "phase = (bjd-t0) / (per)  # phase in days relative to T0 ephemeris\n",
+    "phase = phase - np.fix(phase[number_of_images-1]) # Have current phase occur at value 0.0\n",
     "\n",
-    "t_fine = np.linspace(np.min(bjd), np.max(bjd), 1000) #times at which to calculate light curve\n",
-    "phase_fine=(t_fine-t0)/(per)  #phase in days relative to T0 ephemeris\n",
-    "phase_fine=phase_fine-np.fix(phase[number_of_images-1]) # Have current phase occur at value 0.0\n",
+    "t_fine = np.linspace(np.min(bjd), np.max(bjd), 1000) # times at which to calculate light curve\n",
+    "phase_fine = (t_fine-t0)/(per)  # phase in days relative to T0 ephemeris\n",
+    "phase_fine = phase_fine-np.fix(phase[number_of_images-1]) # Have current phase occur at value 0.0\n",
     "\n",
-    "b0=a_Rs * np.sqrt((np.sin(phase * 2* np.pi)) ** 2 + (np.cos(inc) * np.cos(phase * 2 * np.pi)) ** 2)\n",
-    "intransit=(b0-rp < 1.0E0).nonzero()  #Select indicies between first and fourth contact\n",
-    "outtransit=(b0-rp > 1.0E0).nonzero() #Select indicies out of transit\n"
+    "b0 = a_Rs * np.sqrt((np.sin(phase * 2* np.pi)) ** 2 + (np.cos(inc) * np.cos(phase * 2 * np.pi)) ** 2)\n",
+    "intransit = (b0-rp < 1.0E0).nonzero()  # Select indicies between first and fourth contact\n",
+    "outtransit = (b0-rp > 1.0E0).nonzero() # Select indicies out of transit\n"
    ]
   },
   {
@@ -464,17 +454,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "shx_tmp=shx/np.mean(shx)-1.0E0       #Set Mean around 0.0\n",
-    "shx_detrend=shx_tmp/np.std(shx_tmp)  #Set standard deviation to 1.0\n",
-    "shy_tmp=shy/np.mean(shy)-1.0E0       #Set Mean around 0.0\n",
-    "shy_detrend=shy_tmp/np.std(shy_tmp)  #Set standard deviation to 1.0\n",
+    "shx_tmp = shx / np.mean(shx) - 1.0E0       # Set Mean around 0.0\n",
+    "shx_detrend = shx_tmp/np.std(shx_tmp)  # Set standard deviation to 1.0\n",
+    "shy_tmp = shy / np.mean(shy) - 1.0E0       # Set Mean around 0.0\n",
+    "shy_detrend = shy_tmp/np.std(shy_tmp)  # Set standard deviation to 1.0\n",
     "\n",
-    "cm=common_mode/np.mean(common_mode)-1.0E0\n",
-    "cm_detrend=cm/np.std(cm)\n",
+    "cm = common_mode / np.mean(common_mode) - 1.0E0\n",
+    "cm_detrend = cm/np.std(cm)\n",
     "\n",
     "fig,axs = plt.subplots()\n",
-    "plt.plot(shx_detrend,label='X-possition')\n",
-    "plt.plot(shy_detrend,label='Y-possition')\n",
+    "plt.plot(shx_detrend, label='X-possition')\n",
+    "plt.plot(shy_detrend, label='Y-possition')\n",
     "plt.xlabel('Image Sequence Number')\n",
     "plt.ylabel('Relative Detector Possition')\n",
     "plt.title('Time-series Detrending Vectors')\n",
@@ -511,15 +501,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "shx=shx_detrend\n",
-    "shy=shy_detrend\n",
+    "shx = shx_detrend\n",
+    "shy = shy_detrend\n",
     "common_mode=cm_detrend\n",
     "\n",
-    "XX=np.array([shx,shy,shx**2,shy**2,shx*shy,common_mode,np.ones(number_of_images)])  #Detrending array without linear time trend\n",
+    "XX=np.array([shx, shy, shx**2, shy**2, shx*shy, common_mode, np.ones(number_of_images)])  # Detrending array without linear time trend\n",
     "XX=np.transpose(XX)\n",
-    "XXX=np.array([shx,shy,shx**2,shy**2,shx*shy,common_mode,phase,np.ones(number_of_images)])  #Detrending array with with linear time trend\n",
-    "XXX=np.transpose(XXX)\n",
-    "\n"
+    "XXX=np.array([shx,shy, shx**2, shy**2, shx*shy, common_mode, phase, np.ones(number_of_images)])  # Detrending array with with linear time trend\n",
+    "XXX=np.transpose(XXX)"
    ]
   },
   {
@@ -537,27 +526,27 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pix1=170       # wavelength bin lower range\n",
-    "pix2=200       # wavelength bin upper range\n",
-    "y=np.sum(all_spec_1D[pix1:pix2,:],axis=0)    # flux over a selected wavelength bin\n",
+    "pix1 = 170       # wavelength bin lower range\n",
+    "pix2 = 200       # wavelength bin upper range\n",
+    "y = np.sum(all_spec_1D[pix1:pix2, :], axis=0)    # flux over a selected wavelength bin\n",
     "\n",
-    "msize=plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
+    "msize = plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
     "plt.rcParams['figure.figsize'] = [10.0, 3.0]           # Figure dimensions\n",
     "\n",
     "fig,axs = plt.subplots()\n",
-    "f=plt.plot(wsdata_all,all_spec_1D[:,0], linewidth=2,zorder=0)  #Plot Region of wavelength bin\n",
-    "plt.fill_between(wsdata_all[pix1:pix2],0,all_spec_1D[pix1:pix2,0],alpha=0.5)\n",
+    "f = plt.plot(wsdata_all,all_spec_1D[:, 0], linewidth=2, zorder=0)  #P lot Region of wavelength bin\n",
+    "plt.fill_between(wsdata_all[pix1:pix2], 0, all_spec_1D[pix1:pix2,0], alpha=0.5)\n",
     "plt.xlabel('Wavelength ($\\mu$m)')\n",
     "plt.ylabel('Flux (e-)')\n",
     "plt.title('1D Extracted Spectrum')\n",
     "axs.xaxis.set_major_locator(ticker.MultipleLocator(0.5))\n",
     "axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.1))\n",
-    "plt.annotate('H$_2$O',xy=(3.0,42000))\n",
-    "plt.annotate('CO$_2$',xy=(4.2,42000))\n",
+    "plt.annotate('H$_2$O', xy=(3.0, 42000))\n",
+    "plt.annotate('CO$_2$', xy=(4.2, 42000))\n",
     "plt.show()\n",
     "\n",
-    "fig,axs = plt.subplots()\n",
-    "plt.scatter(bjd,y/np.mean(y[outtransit]),label='$f(t)$ Data',zorder=1,s=msize*0.75,linewidth=1 ,alpha=0.4, marker='+',edgecolors='blue')\n",
+    "fig, axs = plt.subplots()\n",
+    "plt.scatter(bjd, y/np.mean(y[outtransit]), label='$f(t)$ Data', zorder=1, s=msize*0.75, linewidth=1 , alpha=0.4, marker='+',edgecolors='blue')\n",
     "plt.xlabel('Barycentric Julian Date (days)')\n",
     "plt.ylabel('Relative Flux')\n",
     "plt.title('Time-series Transit Light Curve  $\\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2])+'] $\\mu$m')\n",
@@ -585,12 +574,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "yfit=regressor.predict(XX)                         # Project the fit over the whole time series\n",
+    "yfit = regressor.predict(XX)                         # Project the fit over the whole time series\n",
     "\n",
     "plt.rcParams['figure.figsize'] = [10.0, 7.0]           # Figure dimensions\n",
-    "msize=plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
-    "plt.scatter(bjd,y/np.mean(y[outtransit]),label='$f(t)$ Data',zorder=1,s=msize*0.75,linewidth=1 ,alpha=0.5, marker='+',edgecolors='blue')\n",
-    "f=plt.plot(bjd,yfit,label='$S(x)$ Regression fit ', linewidth=2,color='orange',zorder=2,alpha=0.85)\n",
+    "msize = plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
+    "plt.scatter(bjd, y/np.mean(y[outtransit]),label='$f(t)$ Data',zorder=1, s=msize*0.75, linewidth=1 , alpha=0.5, marker='+', edgecolors='blue')\n",
+    "f = plt.plot(bjd, yfit, label='$S(x)$ Regression fit ', linewidth=2, color='orange', zorder=2, alpha=0.85)\n",
     "plt.xlabel('Barycentric Julian Date (days)')\n",
     "plt.ylabel('Relative Flux')\n",
     "plt.title('Time-series Transit Light Curve  $\\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2])+'] $\\mu$m')\n",
@@ -598,8 +587,8 @@
     "axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.005))\n",
     "axs.yaxis.set_major_locator(ticker.MultipleLocator(0.002))\n",
     "axs.yaxis.set_minor_locator(ticker.MultipleLocator(0.001))\n",
-    "yplot=y/np.mean(y[outtransit])\n",
-    "plt.ylim(yplot.min()*0.999, yplot.max()*1.001)\n",
+    "yplot = y / np.mean(y[outtransit])\n",
+    "plt.ylim(yplot.min() * 0.999, yplot.max()*1.001)\n",
     "plt.xlim(bjd.min()-0.001, bjd.max()+0.001)\n",
     "plt.legend(loc='lower left')\n",
     "plt.show()"
@@ -1167,22 +1156,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Functions to call and calculate models\n",
+    "# Functions to call and calculate models\n",
     "def residual(p,phase,x,y,err,c1, c2, c3, c4):\n",
-    "    #calculate new orbit\n",
-    "    b0=p['a_Rs'].value * np.sqrt((np.sin(phase * 2* np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase * 2 * np.pi)) ** 2)\n",
-    "    #Select indicies between first and fourth contact\n",
-    "    intransit=(b0-p['rprs'].value < 1.0E0).nonzero()\n",
-    "    #Make light curve model, set all values initially to 1.0\n",
-    "    light_curve=b0/b0\n",
-    "    mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0[intransit])  #Madel and Agol\n",
+    "    # calculate new orbit\n",
+    "    b0 = p['a_Rs'].value * np.sqrt((np.sin(phase * 2* np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase * 2 * np.pi)) ** 2)\n",
+    "    # Select indicies between first and fourth contact\n",
+    "    intransit = (b0-p['rprs'].value < 1.0E0).nonzero()\n",
+    "    # Make light curve model, set all values initially to 1.0\n",
+    "    light_curve = b0/b0\n",
+    "    mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0[intransit])  # Madel and Agol\n",
     "    light_curve[intransit]=mulimb0\n",
     "    model=(light_curve)*p['f0'].value * (p['Fslope'].value*phase + p['xsh'].value*shx + p['x2sh'].value*shx**2. + p['ysh'].value*shy + p['y2sh'].value*shy**2. + p['xysh'].value*shy*shx +  p['comm'].value*common_mode + 1.0) # transit model is baseline flux X transit model X systematics model\n",
     "    chi2now=np.sum((y-model)**2/err**2)\n",
     "    res=np.std((y-model)/p['f0'].value)\n",
-    "    print(\"rprs: \",p['rprs'].value,\"current chi^2=\",chi2now,' scatter ',res,end=\"\\r\")\n",
+    "    print(\"rprs: \", p['rprs'].value, \"current chi^2=\", chi2now,' scatter ', res, end=\"\\r\")\n",
     "    return (y-model)/err\n",
-    "    #return np.sum((y-model)**2/err**2)"
+    "    # return np.sum((y-model)**2/err**2)"
    ]
   },
   {
@@ -1198,10 +1187,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def model_fine(p):  #Make Transit model with a fine grid for plotting purposes\n",
-    "    b0=p['a_Rs'].value * np.sqrt((np.sin(phase_fine * 2* np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase_fine * 2 * np.pi)) ** 2)\n",
-    "    mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0)  #Madel and Agol\n",
-    "    model_fine=mulimb0\n",
+    "def model_fine(p):  # Make Transit model with a fine grid for plotting purposes\n",
+    "    b0 = p['a_Rs'].value * np.sqrt((np.sin(phase_fine * 2* np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase_fine * 2 * np.pi)) ** 2)\n",
+    "    mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0)  # Madel and Agol\n",
+    "    model_fine = mulimb0\n",
     "    return model_fine"
    ]
   },
@@ -1218,13 +1207,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "wave1=wsdata_all[pix1]\n",
-    "wave2=wsdata_all[pix2]\n",
+    "wave1 = wsdata_all[pix1]\n",
+    "wave2 = wsdata_all[pix2]\n",
     "bin_wave_index = ((wsdata_all > wave1) & (wsdata_all <= wave2)).nonzero()\n",
-    "wsdata=wsdata_all[bin_wave_index]*1E4 # Select wavelength bin values (um=> angstroms)\n",
+    "wsdata = wsdata_all[bin_wave_index]*1E4 # Select wavelength bin values (um=> angstroms)\n",
     "\n",
-    "_uLD, c1, c2, c3, c4, _cp1, _cp2, _cp3, _cp4, aLD, bLD = limb_dark_fit(grating,wsdata, M_H, Teff,logg, limb_dark_directory, ld_model)\n",
-    "\n"
+    "_uLD, c1, c2, c3, c4, _cp1, _cp2, _cp3, _cp4, aLD, bLD = limb_dark_fit(grating,wsdata, M_H, Teff,logg, limb_dark_directory, ld_model)"
    ]
   },
   {
@@ -1242,23 +1230,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Run the Transit Model\n",
+    "# Run the Transit Model\n",
     "rl = 0.0825     # Planet-to-star Radius Ratio\n",
     "\n",
-    "b0=a_Rs * np.sqrt((np.sin(phase * 2* np.pi)) ** 2 + (np.cos(inc) * np.cos(phase * 2 * np.pi)) ** 2)\n",
-    "intransit=(b0-rl < 1.0E0).nonzero()  #Select indicies between first and fourth contact\n",
+    "b0 = a_Rs * np.sqrt((np.sin(phase * 2* np.pi)) ** 2 + (np.cos(inc) * np.cos(phase * 2 * np.pi)) ** 2)\n",
+    "intransit = (b0-rl < 1.0E0).nonzero()  #Select indicies between first and fourth contact\n",
     "\n",
     "mulimb0, mulimbf = occultnl(rl, c1, c2, c3, c4, b0)  #Mandel & Agol non-linear limb darkened transit model\n",
-    "model=mulimb0*yfit \n",
+    "model = mulimb0*yfit \n",
     "\n",
-    "#plot\n",
+    "# plot\n",
     "plt.rcParams['figure.figsize'] = [10.0, 7.0]           # Figure dimensions\n",
-    "msize=plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
-    "fig=plt.figure(constrained_layout=True)\n",
-    "gs = fig.add_gridspec(3, 1,hspace=0.00, wspace=0.00)\n",
-    "ax1=fig.add_subplot(gs[0:2,:])\n",
-    "ax1.scatter(bjd,y/np.mean(y[outtransit]),label='$f(t)$ Data',zorder=1,s=msize*0.75,linewidth=1 ,alpha=0.5, marker='+',edgecolors='blue')\n",
-    "ax1.plot(bjd,model,label='$S(x)$ Regression fit ', linewidth=2,color='orange',zorder=2,alpha=0.85)\n",
+    "msize = plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
+    "fig = plt.figure(constrained_layout=True)\n",
+    "gs = fig.add_gridspec(3, 1, hspace=0.00, wspace=0.00)\n",
+    "ax1 = fig.add_subplot(gs[0:2, :])\n",
+    "ax1.scatter(bjd, y/np.mean(y[outtransit]), label='$f(t)$ Data', zorder=1, s=msize*0.75, linewidth=1 , alpha=0.5, marker='+', edgecolors='blue')\n",
+    "ax1.plot(bjd,model, label='$S(x)$ Regression fit ', linewidth=2, color='orange', zorder=2, alpha=0.85)\n",
     "ax1.xaxis.set_ticklabels([])\n",
     "plt.ylabel('Relative Flux')\n",
     "plt.title('Time-series Transit Light Curve  $\\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2])+'] $\\mu$m')\n",
@@ -1266,31 +1254,31 @@
     "ax1.xaxis.set_minor_locator(ticker.MultipleLocator(0.005))\n",
     "ax1.yaxis.set_major_locator(ticker.MultipleLocator(0.002))\n",
     "ax1.yaxis.set_minor_locator(ticker.MultipleLocator(0.001))\n",
-    "yplot=y/np.mean(y[outtransit])\n",
+    "yplot = y/np.mean(y[outtransit])\n",
     "plt.ylim(yplot.min()*0.999, yplot.max()*1.001)\n",
     "plt.xlim(bjd.min()-0.001, bjd.max()+0.001)\n",
     "plt.legend()\n",
     "fig.add_subplot(ax1)\n",
-    "#Residual\n",
-    "ax2=fig.add_subplot(gs[2,:])\n",
-    "ax2.scatter(bjd,1E6*(y/np.mean(y[outtransit])-model),label='$f(t)$ Data',zorder=1,s=msize*0.75,linewidth=1 ,alpha=0.5, marker='+',edgecolors='blue')\n",
+    "# Residual\n",
+    "ax2=fig.add_subplot(gs[2, :])\n",
+    "ax2.scatter(bjd, 1E6*(y/np.mean(y[outtransit])-model), label='$f(t)$ Data', zorder=1, s=msize*0.75, linewidth=1, alpha=0.5, marker='+', edgecolors='blue')\n",
     "wsb, wsb_bin_edges,binnumber = stats.binned_statistic(bjd,1E6*(y/np.mean(y[outtransit])-model), bins=256)\n",
     "plt.scatter(wsb_bin_edges[1:],wsb, linewidth=2,alpha=0.75,facecolors='orange',edgecolors='none', marker='o',zorder=25)\n",
     "plt.xlabel('Barycentric Julian Date (days)')\n",
     "plt.ylabel('Residual (ppm)')\n",
     "ax2.xaxis.set_major_locator(ticker.MultipleLocator(0.01))\n",
     "ax2.xaxis.set_minor_locator(ticker.MultipleLocator(0.005))\n",
-    "yplot=y/np.mean(y[outtransit])\n",
+    "yplot = y / np.mean(y[outtransit])\n",
     "plt.xlim(bjd.min()-0.001, bjd.max()+0.001)\n",
     "fig.add_subplot(ax2)\n",
     "plt.show()\n",
     "\n",
     "#print chi^2 value\n",
-    "err=np.sqrt(y)/np.mean(y[outtransit])\n",
+    "err = np.sqrt(y) / np.mean(y[outtransit])\n",
     "print('Chi^2 = '+str(np.sum((y/np.mean(y[outtransit])-model)**2/err**2)))\n",
     "\n",
     "print('Residual Standard Deviation : '+str(1E6*np.std((y/np.mean(y[outtransit])-model)))+' ppm')\n",
-    "print('256 Bin Standard Deviation  :'+str(np.std(wsb))+' ppm')\n"
+    "print('256 Bin Standard Deviation  :'+str(np.std(wsb))+' ppm')"
    ]
   },
   {
@@ -1348,17 +1336,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "k0   = 113 #98   #100\n",
+    "k0 = 113 #98   #100\n",
     "kend = 392 #422\n",
-    "wk   = 15\n",
+    "wk = 15\n",
     "number_of_bins = int((kend-k0)/wk)\n",
-    "wsd  = np.zeros((number_of_bins))\n",
+    "wsd = np.zeros((number_of_bins))\n",
     "werr = np.zeros((number_of_bins))\n",
     "rprs = np.zeros((number_of_bins))\n",
     "rerr = np.zeros((number_of_bins))\n",
     "sig_r = np.zeros((number_of_bins))\n",
     "sig_w = np.zeros((number_of_bins))\n",
-    "beta  = np.zeros((number_of_bins))\n",
+    "beta = np.zeros((number_of_bins))\n",
     "depth = np.zeros((number_of_bins))\n",
     "depth_err = np.zeros((number_of_bins))"
    ]
@@ -1380,27 +1368,24 @@
    "outputs": [],
    "source": [
     "k=k0   #wavelength to start\n",
-    "#--------------------------------------------------------------------------\n",
-    "#Loop over wavelength bins and fit for each one\n",
+    "# --------------------------------------------------------------------------\n",
+    "# Loop over wavelength bins and fit for each one\n",
     "for bin in range(0,number_of_bins):\n",
     "    \n",
-    "    #---------------------------------------------------------\n",
     "    # Select wavelength bin\n",
-    "    wave1=wsdata_all[k]\n",
-    "    wave2=wsdata_all[k+wk]\n",
+    "    wave1 = wsdata_all[k]\n",
+    "    wave2 = wsdata_all[k+wk]\n",
     "\n",
-    "    #Indicies to select for wavelgth bin\n",
+    "    # Indicies to select for wavelgth bin\n",
     "    bin_wave_index = ((wsdata_all > wave1) & (wsdata_all <= wave2)).nonzero()\n",
     "\n",
-    "    #make light curve bin\n",
-    "    wave_bin_counts=np.sum(all_spec_1D[k+1:k+wk,:],axis=0) #Sum Wavelength pixels\n",
-    "    wave_bin_counts_err=np.sqrt(wave_bin_counts)           #adopt photon noise for errors\n",
-    "    #---------------------------------------------------------\n",
+    "    # make light curve bin\n",
+    "    wave_bin_counts = np.sum(all_spec_1D[k+1:k+wk, :], axis=0)  # Sum Wavelength pixels\n",
+    "    wave_bin_counts_err = np.sqrt(wave_bin_counts)            # adopt photon noise for errors\n",
     "\n",
-    "    #---------------------------------------------------------\n",
     "    # Calculate Limb Darkening\n",
     "\n",
-    "    wsdata=wsdata_all[bin_wave_index]*1E4 # Select wavelength bin values (um=> angstroms)\n",
+    "    wsdata = wsdata_all[bin_wave_index]*1E4 # Select wavelength bin values (um=> angstroms)\n",
     "    _uLD, c1, c2, c3, c4, _cp1, _cp2, _cp3, _cp4, aLD, bLD = limb_dark_fit(grating,wsdata, M_H, Teff,logg, limb_dark_directory, ld_model)\n",
     "\n",
     "    print('\\nc1 = {}'.format(c1))\n",
@@ -1408,58 +1393,55 @@
     "    print('c3 = {}'.format(c3))\n",
     "    print('c4 = {}'.format(c4))\n",
     "    print('')\n",
-    "    #u   = [c1,c2,c3,c4]      #limb darkening coefficients\n",
+    "    # u   = [c1,c2,c3,c4]      # limb darkening coefficients\n",
     "    u   = [aLD,bLD]\n",
-    "    #---------------------------------------------------------\n",
     "\n",
-    "    #---------------------------------------------------------\n",
     "    # Make initial model\n",
     "    \n",
     "    #Setup LMFIT\n",
-    "    x=bjd                    # X data\n",
-    "    y=wave_bin_counts        # Y data\n",
-    "    err=wave_bin_counts_err  # Y Error\n",
+    "    x = bjd                    # X data\n",
+    "    y = wave_bin_counts        # Y data\n",
+    "    err = wave_bin_counts_err  # Y Error\n",
     "\n",
-    "    #Perform Quick Linear regression on out-of-transit data to obtain accurate starting Detector fit values\n",
+    "    # Perform Quick Linear regression on out-of-transit data to obtain accurate starting Detector fit values\n",
     "    if wave1 > 2.7 and wave1 < 3.45:\n",
     "        regressor.fit(XXX[outtransit], y[outtransit]/np.mean(y[outtransit]))\n",
     "    else:\n",
     "        regressor.fit(XX[outtransit], y[outtransit]/np.mean(y[outtransit]))\n",
     "\n",
     "    # create a set of Parameters for LMFIT https://lmfit.github.io/lmfit-py/parameters.html\n",
-    "    #class Parameter(name, value=None, vary=True, min=- inf, max=inf, expr=None, brute_step=None, user_data=None)¶\n",
-    "    #Set vary=0 to fix\n",
-    "    #Set vary=1 to fit\n",
+    "    # class Parameter(name, value=None, vary=True, min=- inf, max=inf, expr=None, brute_step=None, user_data=None)¶\n",
+    "    # Set vary=0 to fix\n",
+    "    # Set vary=1 to fit\n",
     "    p = lmfit.Parameters()  #object to store L-M fit Parameters           # Parameter Name\n",
-    "    p.add('cosinc'  , value=np.cos(inc)           ,vary=0)                # inclination, vary cos(inclin)\n",
-    "    p.add('rho_star', value=rho_star              ,vary=0)                # stellar density\n",
-    "    p.add('a_Rs'    , value=a_Rs                  ,vary=0)                # a/Rstar\n",
-    "    p.add('rprs'    , value=rp                     ,vary=1, min=0, max=1)  # planet-to-star radius ratio\n",
-    "    p.add('t0'      , value=t0                    ,vary=0)                # Transit T0\n",
-    "    p.add('f0'      , value=np.mean(y[outtransit]),vary=1, min=0)         # Baseline Flux\n",
-    "    p.add('ecc'     , value=ecc                   ,vary=0, min=0 , max=1) # eccentricity\n",
-    "    p.add('omega'   , value=omega                 ,vary=0)                # arguments of periatron\n",
-    "    #Turn on a linear slope in water feature to account for presumably changing H2O ice builtup on widow during cryogenic test\n",
+    "    p.add('cosinc', value=np.cos(inc), vary=0)                # inclination, vary cos(inclin)\n",
+    "    p.add('rho_star', value=rho_star, vary=0)                # stellar density\n",
+    "    p.add('a_Rs', value=a_Rs, vary=0)                # a/Rstar\n",
+    "    p.add('rprs', value=rp, vary=1, min=0, max=1)  # planet-to-star radius ratio\n",
+    "    p.add('t0', value=t0, vary=0)                # Transit T0\n",
+    "    p.add('f0', value=np.mean(y[outtransit]), vary=1, min=0)         # Baseline Flux\n",
+    "    p.add('ecc', value=ecc, vary=0, min=0, max=1) # eccentricity\n",
+    "    p.add('omega', value=omega, vary=0)                # arguments of periatron\n",
+    "    # Turn on a linear slope in water feature to account for presumably changing H2O ice builtup on widow during cryogenic test\n",
     "    if wave1 > 2.7 and wave1 < 3.45:\n",
-    "        p.add('Fslope', value=regressor.coef_[6]  ,vary=1)                # Orbital phase\n",
+    "        p.add('Fslope', value=regressor.coef_[6], vary=1)                # Orbital phase\n",
     "    else:\n",
-    "        p.add('Fslope', value=0                   ,vary=0)                # Orbital phase\n",
-    "    p.add('xsh'     , value=regressor.coef_[0]    ,vary=1)                # Detector X-shift detrending\n",
-    "    p.add('ysh'     , value=regressor.coef_[1]    ,vary=1)                # Detector X-shift detrending\n",
-    "    p.add('x2sh'    , value=regressor.coef_[2]    ,vary=1)                # Detector X^2-shift detrending\n",
-    "    p.add('y2sh'    , value=regressor.coef_[3]    ,vary=1)                # Detector Y^2-shift detrending\n",
-    "    p.add('xysh'    , value=regressor.coef_[4]    ,vary=1)                # Detector X*Y detrending\n",
-    "    p.add('comm'    , value=regressor.coef_[5]    ,vary=1)                # Common-Mode detrending\n",
+    "        p.add('Fslope', value=0, vary=0)                          # Orbital phase\n",
+    "    p.add('xsh', value=regressor.coef_[0], vary=1)                 # Detector X-shift detrending\n",
+    "    p.add('ysh', value=regressor.coef_[1], vary=1)                 # Detector X-shift detrending\n",
+    "    p.add('x2sh', value=regressor.coef_[2], vary=1)                # Detector X^2-shift detrending\n",
+    "    p.add('y2sh', value=regressor.coef_[3], vary=1)                # Detector Y^2-shift detrending\n",
+    "    p.add('xysh', value=regressor.coef_[4], vary=1)                # Detector X*Y detrending\n",
+    "    p.add('comm', value=regressor.coef_[5], vary=1)                # Common-Mode detrending\n",
     "    \n",
-    "    #--------------------------------------------------------------------------\n",
     "    # Perform Minimization https://lmfit.github.io/lmfit-py/fitting.html\n",
     "    # create Minimizer\n",
     "    # mini = lmfit.Minimizer(residual, p, nan_policy='omit',fcn_args=(phase,x,y,err)\n",
     "    print('')\n",
-    "    print('Fitting Bin',bin,' Wavelength =',np.mean(wsdata)/1E4, '  Range= [',wave1,':',wave2,']')\n",
+    "    print('Fitting Bin', bin,' Wavelength =', np.mean(wsdata)/1E4, '  Range= [',wave1, ':', wave2,']')\n",
     "\n",
     "    #  solve with Levenberg-Marquardt using the\n",
-    "    result = lmfit.minimize(residual,params=p,args=(phase,x,y,err,c1, c2, c3, c4))\n",
+    "    result = lmfit.minimize(residual, params=p, args=(phase, x, y, err, c1, c2, c3, c4))\n",
     "    #result = mini.minimize(method='emcee')\n",
     "\n",
     "    print('')\n",
@@ -1484,54 +1466,56 @@
     "\n",
     "    #Update with best-fit parameters\n",
     "    p['rho_star'].value = result.params['rho_star'].value\n",
-    "    p['cosinc'].value   = result.params['cosinc'].value\n",
-    "    p['rprs'].value     = result.params['rprs'].value\n",
-    "    p['t0'].value       = result.params['t0'].value\n",
-    "    p['f0'].value       = result.params['f0'].value\n",
-    "    p['Fslope'].value   = result.params['Fslope'].value\n",
-    "    p['xsh'].value      = result.params['xsh'].value\n",
-    "    p['ysh'].value      = result.params['ysh'].value\n",
-    "    p['x2sh'].value     = result.params['x2sh'].value\n",
-    "    p['y2sh'].value     = result.params['y2sh'].value\n",
-    "    p['xysh'].value     = result.params['xysh'].value\n",
-    "    p['comm'].value     = result.params['comm'].value\n",
+    "    p['cosinc'].value = result.params['cosinc'].value\n",
+    "    p['rprs'].value = result.params['rprs'].value\n",
+    "    p['t0'].value = result.params['t0'].value\n",
+    "    p['f0'].value = result.params['f0'].value\n",
+    "    p['Fslope'].value = result.params['Fslope'].value\n",
+    "    p['xsh'].value = result.params['xsh'].value\n",
+    "    p['ysh'].value = result.params['ysh'].value\n",
+    "    p['x2sh'].value = result.params['x2sh'].value\n",
+    "    p['y2sh'].value = result.params['y2sh'].value\n",
+    "    p['xysh'].value = result.params['xysh'].value\n",
+    "    p['comm'].value = result.params['comm'].value\n",
     "    # Update Fit Spectra arrays\n",
-    "    wsd[bin]=np.mean(wsdata)/1E4\n",
-    "    werr[bin]=(wsdata.max()-wsdata.min())/2E4\n",
-    "    rprs[bin]=result.params['rprs'].value\n",
-    "    rerr[bin]=result.params['rprs'].stderr\n",
+    "    wsd[bin] = np.mean(wsdata)/1E4\n",
+    "    werr[bin] = (wsdata.max()-wsdata.min())/2E4\n",
+    "    rprs[bin] = result.params['rprs'].value\n",
+    "    rerr[bin] = result.params['rprs'].stderr\n",
     "\n",
     "    # Calculate Bestfit Model\n",
-    "    final_model=y-result.residual*err\n",
-    "    final_model_fine=model_fine(p)\n",
-    "\n",
-    "    #More Stats\n",
-    "    resid=(y-final_model)/p['f0'].value\n",
-    "    residppm=1E6*(y-final_model)/p['f0'].value\n",
-    "    residerr=err/p['f0'].value\n",
-    "    sigma=np.std((y-final_model)/p['f0'].value)*1E6\n",
-    "    print(\"Residual standard deviation  (ppm) : \",1E6*np.std((y-final_model)/p['f0'].value))\n",
+    "    final_model = y-result.residual*err\n",
+    "    final_model_fine = model_fine(p)\n",
+    "\n",
+    "    # More Stats\n",
+    "    resid = (y-final_model)/p['f0'].value\n",
+    "    residppm = 1E6*(y-final_model)/p['f0'].value\n",
+    "    residerr = err/p['f0'].value\n",
+    "    sigma = np.std((y-final_model)/p['f0'].value)*1E6\n",
+    "    print(\"Residual standard deviation  (ppm) : \", 1E6*np.std((y-final_model)/p['f0'].value))\n",
     "    print(\"Photon noise                 (ppm) : \", (1/np.sqrt(p['f0'].value))*1E6     )\n",
     "    print(\"Photon noise performance       (%) : \", (1/np.sqrt(p['f0'].value))*1E6 / (sigma) *100 )\n",
     "    n = text_file.write(\"\\nResidual standard deviation  (ppm) : \"+str(1E6*np.std((y-final_model)/p['f0'].value)))\n",
     "    n = text_file.write(\"\\nPhoton noise                 (ppm) : \"+str((1/np.sqrt(p['f0'].value))*1E6))\n",
     "    n = text_file.write(\"\\nPhoton noise performance       (%) : \"+str((1/np.sqrt(p['f0'].value))*1E6 / (sigma) *100 ))\n",
     " \n",
-    "    #Measure Rednoise with Binning Technique\n",
-    "    sig0=np.std(resid)\n",
-    "    bins=number_of_images/binmeasure\n",
+    "    # Measure Rednoise with Binning Technique\n",
+    "    sig0 = np.std(resid)\n",
+    "    bins = number_of_images / binmeasure\n",
     "    wsb, wsb_bin_edges,binnumber = stats.binned_statistic(bjd,resid, bins=bins)\n",
-    "    sig_binned=np.std(wsb)\n",
-    "    sigrednoise=np.sqrt(sig_binned**2-sig0**2/binmeasure)\n",
-    "    if np.isnan(sigrednoise) == True : sigrednoise=0   #if no rednoise detected, set to zero\n",
-    "    sigwhite   =np.sqrt(sig0**2-sigrednoise**2)\n",
-    "    sigrednoise=np.sqrt(sig_binned**2-sigwhite**2/binmeasure)\n",
-    "    if np.isnan(sigrednoise) == True : sigrednoise=0   #if no rednoise detected, set to zero\n",
-    "    beta[bin]=np.sqrt(sig0**2+binmeasure*sigrednoise**2)/sig0\n",
+    "    sig_binned = np.std(wsb)\n",
+    "    sigrednoise = np.sqrt(sig_binned**2-sig0**2/binmeasure)\n",
+    "    if np.isnan(sigrednoise):\n",
+    "        sigrednoise=0   # if no rednoise detected, set to zero\n",
+    "    sigwhite =np.sqrt(sig0**2-sigrednoise**2)\n",
+    "    sigrednoise = np.sqrt(sig_binned**2-sigwhite**2/binmeasure)\n",
+    "    if np.isnan(sigrednoise): \n",
+    "        sigrednoise=0   #if no rednoise detected, set to zero\n",
+    "    beta[bin] = np.sqrt(sig0**2+binmeasure*sigrednoise**2)/sig0\n",
     "    \n",
-    "    print(\"White noise                  (ppm) : \",1E6*sigwhite)\n",
-    "    print(\"Red noise                    (ppm) : \",1E6*sigrednoise)\n",
-    "    print(\"Transit depth measured error (ppm) : \",2E6*result.params['rprs'].value*result.params['rprs'].stderr)\n",
+    "    print(\"White noise                  (ppm) : \", 1E6*sigwhite)\n",
+    "    print(\"Red noise                    (ppm) : \", 1E6*sigrednoise)\n",
+    "    print(\"Transit depth measured error (ppm) : \", 2E6*result.params['rprs'].value*result.params['rprs'].stderr)\n",
     "    \n",
     "    n = text_file.write(\"\\nWhite noise                  (ppm) : \"+str(1E6*sigwhite))\n",
     "    n = text_file.write(\"\\nRed noise                    (ppm) : \"+str(1E6*sigrednoise))\n",
@@ -1549,9 +1533,9 @@
     "    ascii.write(ascii_data, save_directory+'JWST_NIRSpec_Prism_fit_transmission_spectra.csv', format='csv',overwrite=True)\n",
     "    #---------------------------------------------------------\n",
     "    msize=plt.rcParams['lines.markersize'] ** 2. #default marker size\n",
-    "    #Plot data models\n",
+    "    # Plot data models\n",
     "     \n",
-    "    #plot\n",
+    "    # plot\n",
     "    plt.rcParams['figure.figsize'] = [10.0, 7.0]           # Figure dimensions\n",
     "    msize=plt.rcParams['lines.markersize'] ** 2.           # default marker size\n",
     "    fig=plt.figure(constrained_layout=True)\n",
@@ -1696,7 +1680,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.11.6"
   }
  },
  "nbformat": 4,