Repository: Exo-TiC/ExoTiC-LD (language: python)

Issue:
Add in JWST mode sensitivity files
Add in ability to run limb-darkening for the following modes:

- [x] NIRCam F322W2
- [x] NIRCam F444
- [x] NIRISS SOSS order 1
- [x] NIRISS SOSS order 2
- [x] NIRSpec G140H/M f100lp
- [x] NIRSpec G140M/H f070lp
- [x] NIRSpec G235H/M f170lp
- [x] NIRSpec G395H/M f290lp
- [x] NIRSpec Prism
- [x] MIRI LRS
- [x] HST WFC3 UVIS G280 -1 spectral order
LD data files
The sensitivity files and kuruczlist files and template may need to be added to the zenodo repository. The way we call them in the script needs to be as universal as possible. 


@hrwakeford - Look into the setup of the repository as an installable package and if there is a best practice for used data files - for example adding their path to the .bashrc file


Discussion / hints:
At some point we will also want to convert the sensitivity files from IDL .sav files to something more appropriate for Python

- suggestions welcomed to the file format for these
Here is how you would do that with PandExo

```
import pandexo.engine.justdoit as jdi
import pandas as pd
#list of instruments
inst = jdi.print_instruments()
#omit jwst
jwst = list(inst)[1:]
thrus = {}
for i in jwst:
    if 'SOSS' in i: 
        name = i+'O1'
        thrus[name] = jdi.get_thruput(i,niriss=1)
        pd.DataFrame(thrus[name]).to_csv(name+'.csv')
        name = i+'O2'
        thrus[name] = jdi.get_thruput(i,niriss=2)
        pd.DataFrame(thrus[name]).to_csv(name+'.csv')
    elif 'NIRSpec G140' in i: 
        name = i+'f100lp'
        thrus[name] = jdi.get_thruput(i,nirspec='f100lp')
        pd.DataFrame(thrus[name]).to_csv(name+'.csv')
        name = i+'f070lp'
        thrus[name] = jdi.get_thruput(i,nirspec='f070lp')
        pd.DataFrame(thrus[name]).to_csv(name+'.csv')        
    else:
        thrus[i] = jdi.get_thruput(i)
        pd.DataFrame(thrus[i]).to_csv(i+'.csv')
    
        
    
```
Amazing thank you @natashabatalha I will get these incorporated ASAP
using [PandExo](https://github.com/natashabatalha/PandExo ) 
While working on this issue it has been found that the G750L and G750M sensitivity files did not read in correctly. 

In addition the following things will be done to improve the consistency of the code:


- [ ] Update all of the current sensitivity files to csv files 

- [ ] Make sure all files are in angstroms for wavelength and percentage for throughput/sensitivity

- [ ] Use consistent and clear file names

- [x] Make a figure showing all of the throughputs that are available in the code
- [ ] Add a gallery of each individual file


Figure showing all the throughputs covered in the code 
<img width="914" alt="Screenshot 2022-03-04 at 13 49 56" src="https://user-images.githubusercontent.com/7988823/156775028-9d1d56aa-c936-46cb-a6e9-bb77015bf526.png">



Fix patch:
diff --git a/README.md b/README.md
index 3ce68a0..e5541f5 100644
--- a/README.md
+++ b/README.md
@@ -12,18 +12,35 @@ This calculation is computed using 1D Kurucz stellar models or 3D stellar models
 This package was built from the original IDL code adapted by Hannah Wakeford and translated into python by Matthew Hill with improvements by Iva Laginja. The git history associated with these steps can be found in the ExoTiC-ISM package from which this is a spin-off repository.
 
 ## Data
-To run this package you will need to download the stellar model grids from [zenodo](https://zenodo.org/record/6052611)
+To run this package you will need to download the stellar model grids and supported instrument information from [ExoTiC-LD_data on zenodo](https://zenodo.org/record/6344946#.YistRy-l2ik)
 
 <!-- Once downloaded un-zip and alter the filepath indicated in limb_darkening.py -->
 
 ## Supported instruments and gratings
 Current supported instruments and gratings are:  
 
-**Hubble** *STIS*: G430L, G750L gratings
+### Spectroscopic:
+**Hubble** *STIS*: G750L, G750M, G430L gratings
 
-**Hubble** *WFC3*: IR/G102, IR/G141, UVIS/G280 grisms
+**Hubble** *WFC3*: UVIS/G280+1, UVIS/G280-1, IR/G102, IR/G141 grisms
+
+**Webb** *NIRSpec*: Prism, G395H, G395M, G235H, G235M, G140H-f100, G140M-f100, G140H-f070, G140M-f070
+
+**Webb** *NIRSpec*: SOSSo1, SOSSo2
+
+**Webb** *NIRSpec*: F322W2, F444
+
+**Webb** *NIRSpec*: LRS
+
+### Photometric:
+
+**TESS**
+
+**Spitzer** *IRAC*: Ch1 (3.6 microns), Ch2 (4.5 microns)
+
+<img src="Supported_spectroscopic_modes.png" width="80%" />  
+<img src="Supported_photometric_modes.png" width="80%" />  
 
-**Webb** *NIRSpec*: Prism, G395H
 
 ## About this repository
 
@@ -43,4 +60,4 @@ This project is licensed under the MIT License - see the [LICENSE.md](LICENSE.tx
 * Hannah Wakeford for IDL code and adaption to an independent Python package - [hrwakeford](https://github.com/hrwakeford)
 * Matthew Hill for a functional translation from IDL to Python - [@mattjhill](https://github.com/mattjhill)
 * Iva Laginja for implementing improvements to the script - [@ivalaginja](https://github.com/ivalaginja)
-* David Sing for providing up-to-date JWST sensitivity profiles
+* Natasha Batalha for providing the Webb throughput information from their PandExo package - [@natashabatalha](https://github.com/natashabatalha)
diff --git a/Sensitivity_files/G102.WFC3.sensitivity.sav b/Sensitivity_files/G102.WFC3.sensitivity.sav
deleted file mode 100755
index a8e859b..0000000
Binary files a/Sensitivity_files/G102.WFC3.sensitivity.sav and /dev/null differ
diff --git a/Sensitivity_files/G141.WFC3.sensitivity.sav b/Sensitivity_files/G141.WFC3.sensitivity.sav
deleted file mode 100755
index 1b5410e..0000000
Binary files a/Sensitivity_files/G141.WFC3.sensitivity.sav and /dev/null differ
diff --git a/Sensitivity_files/G280.WFC3.sensitivity.sav b/Sensitivity_files/G280.WFC3.sensitivity.sav
deleted file mode 100644
index 688e25f..0000000
Binary files a/Sensitivity_files/G280.WFC3.sensitivity.sav and /dev/null differ
diff --git a/Sensitivity_files/G430L.STIS.sensitivity.sav b/Sensitivity_files/G430L.STIS.sensitivity.sav
deleted file mode 100755
index 03d6942..0000000
Binary files a/Sensitivity_files/G430L.STIS.sensitivity.sav and /dev/null differ
diff --git a/Sensitivity_files/G750L.STIS.sensitivity.sav b/Sensitivity_files/G750L.STIS.sensitivity.sav
deleted file mode 100755
index 1baffcc..0000000
Binary files a/Sensitivity_files/G750L.STIS.sensitivity.sav and /dev/null differ
diff --git a/Sensitivity_files/G750M.STIS.sensitivity.sav b/Sensitivity_files/G750M.STIS.sensitivity.sav
deleted file mode 100755
index 06ec8e6..0000000
Binary files a/Sensitivity_files/G750M.STIS.sensitivity.sav and /dev/null differ
diff --git a/Sensitivity_files/NIRSpec.G395H.sensitivity.pandeia.sav b/Sensitivity_files/NIRSpec.G395H.sensitivity.pandeia.sav
deleted file mode 100644
index 23d285d..0000000
Binary files a/Sensitivity_files/NIRSpec.G395H.sensitivity.pandeia.sav and /dev/null differ
diff --git a/Sensitivity_files/NIRSpec.prism.sensitivity.pandeia.sav b/Sensitivity_files/NIRSpec.prism.sensitivity.pandeia.sav
deleted file mode 100644
index ef887ba..0000000
Binary files a/Sensitivity_files/NIRSpec.prism.sensitivity.pandeia.sav and /dev/null differ
diff --git a/Sensitivity_files/kuruczlist.sav b/Sensitivity_files/kuruczlist.sav
deleted file mode 100755
index f2521b2..0000000
Binary files a/Sensitivity_files/kuruczlist.sav and /dev/null differ
diff --git a/Sensitivity_files/templates.sav b/Sensitivity_files/templates.sav
deleted file mode 100755
index 650821c..0000000
Binary files a/Sensitivity_files/templates.sav and /dev/null differ
diff --git a/Supported_photometric_modes.png b/Supported_photometric_modes.png
new file mode 100644
index 0000000..76642ff
Binary files /dev/null and b/Supported_photometric_modes.png differ
diff --git a/Supported_spectroscopic_modes.png b/Supported_spectroscopic_modes.png
new file mode 100644
index 0000000..a1a5392
Binary files /dev/null and b/Supported_spectroscopic_modes.png differ
diff --git a/limb_darkening.py b/limb_darkening.py
index aa80483..f371ace 100755
--- a/limb_darkening.py
+++ b/limb_darkening.py
@@ -17,19 +17,26 @@ def limb_dark_fit(mode, wsdata, M_H, Teff, logg, dirsen, ld_model='1D'):
     Calculates stellar limb-darkening coefficients for a given wavelength bin.
 
     Modes Currently Supported:
+    Spectroscopic:
     HST STIS G750L, G750M, G430L gratings
-    HST WFC3 UVIS/G280, IR/G102, IR/G141 grisms
-    JWST NIRSpec Prism, G395H
+    HST WFC3 UVIS/G280+1, UVIS/G280-1, IR/G102, IR/G141 grisms
+    JWST NIRSpec Prism, G395H, G395M, G235H, G235M, G140H-f100, G140M-f100, G140H-f070, G140M-f070
+    JWST NIRISS SOSSo1, SOSSo2
+    JWST NIRCam F322W2, F444
+    JWST MIRI LRS
+    Photometric:
+    TESS
+    Spitzer IRAC Ch1 (3.6 microns), Ch2 (4.5 microns)
 
     Procedure from Sing et al. (2010, A&A, 510, A21).
     Uses 3D limb darkening from Magic et al. (2015, A&A, 573, 90).
     Uses photon FLUX Sum over (lambda*dlamba).
-    :param mode: string; mode to use ('G430L','G750L','G750M', 'G280', 'G102', 'G141', 'NIRSpec_Prism', 'NIRSpec_G395H')
-    :param wsdata: array; data wavelength solution
+    :param mode: string; mode to use Spectroscopic: ('STIS_G430L','STIS_G750L', 'WFC3_G280p1', 'WFC3_G280n1', 'WFC3_G102', 'WFC3_G141', 'NIRSpec_Prism', 'NIRSpec_G395H', 'NIRSpec_G395M', 'NIRSpec_G235H', 'NIRSpec_G235M', 'NIRSpec_G140Hf100', 'NIRSpec_G140Mf100', 'NIRSpec_G140Hf070', 'NIRSpec_G140Mf070', 'NIRISS_SOSSo1', 'NIRISS_SOSSo2', 'NIRCam_F322W2', 'NIRCam_F444', 'MIRI_LRS'), Photometric: ('IRAC_Ch1', 'IRAC_Ch2', 'TESS')
+    :param wsdata: array; data wavelength solution for range required
     :param M_H: float; stellar metallicity
     :param Teff: float; stellar effective temperature (K)
     :param logg: float; stellar gravity
-    :param dirsen: string; path to main limb darkening directory
+    :param dirsen: string; path to main limb darkening directory downloaded from Zenodo V2.1
     :param ld_model: string; '1D' or '3D', makes choice between limb darkening models; default is 1D
     :return: uLD: float; linear limb darkening coefficient
     aLD, bLD: float; quadratic limb darkening coefficients
@@ -54,9 +61,8 @@ def limb_dark_fit(mode, wsdata, M_H, Teff, logg, dirsen, ld_model='1D'):
         MH_ind = M_H_Grid_load[optM]
 
         # Determine which model is to be used, by using the input metallicity M_H to figure out the file name we need
-        direc = 'Kurucz'
         file_list = 'kuruczlist.sav'
-        sav1 = readsav(os.path.join(dirsen, file_list))
+        sav1 = readsav(os.path.join(direc, file_list))
         model = bytes.decode(sav1['li'][MH_ind])  # Convert object of type "byte" to "string"
 
         # Select Teff and subsequently logg
@@ -230,57 +236,154 @@ def limb_dark_fit(mode, wsdata, M_H, Teff, logg, dirsen, ld_model='1D'):
 
     ### Load response function and interpolate onto kurucz model grid
 
+    ## Spectroscopic Modes
     # FOR Hubble STIS
-    if mode == 'G430L':
-        sav = readsav(os.path.join(dirsen, 'G430L.STIS.sensitivity.sav'))  # wssens,sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
-        wdel = 3
-
-    if mode == 'G750M':
-        sav = readsav(os.path.join(dirsen, 'G750M.STIS.sensitivity.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
-        wdel = 0.554
-
-    if mode == 'G750L':
-        sav = readsav(os.path.join(dirsen, 'G750L.STIS.sensitivity.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
-        wdel = 4.882
+    if mode == 'STIS_G430L':
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/HST_STIS_G430L_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
 
-    # FOR Hubble WFC3
-    if mode == 'G141':  # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
-        sav = readsav(os.path.join(dirsen, 'G141.WFC3.sensitivity.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
+    if mode == 'STIS_G750L':
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/HST_STIS_G750L_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
         wdel = 1
 
-    if mode == 'G102':  # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
-        sav = readsav(os.path.join(dirsen, 'G102.WFC3.sensitivity.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
+    # FOR Hubble WFC3
+    if mode == 'WFC3_G141': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/HST_WFC3_G141_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+        
+    if mode == 'WFC3_G102': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/HST_WFC3_G102_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
         wdel = 1
 
-    if mode == 'G280':  # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
-        sav = readsav(os.path.join(dirsen, 'G280.WFC3.sensitivity.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
+    if mode == 'WFC3_G280p1': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/HST_WFC3_G280p1_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
         wdel = 1
 
+    if mode == 'WFC3_G280n1':  
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/HST_WFC3_G280n1_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1        
+
     # FOR Webb NIRSpec
-    if mode == 'NIRSpec_Prism':  # Provided by D.K. Sing (JHU)
-        sav = readsav(os.path.join(dirsen, 'NIRSpec.prism.sensitivity.pandeia.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
+    if mode == 'NIRSpec_Prism':  
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_prism_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+
+    if mode == 'NIRSpec_G395H': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G395H_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
         wdel = 1
 
-    if mode == 'NIRSpec_G395H':  # Provided by D.K. Sing (JHU)
-        sav = readsav(os.path.join(dirsen, 'NIRSpec.G395H.sensitivity.pandeia.sav'))  # wssens, sensitivity
-        wssens = sav['wssens']
-        sensitivity = sav['sensitivity']
+    if mode == 'NIRSpec_G395M': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G395M_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
         wdel = 1        
 
+    if mode == 'NIRSpec_G235H': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G235H_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+        
+    if mode == 'NIRSpec_G235M': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G235M_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1        
+        
+    if mode == 'NIRSpec_G140Hf100': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G140H-f100_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+        
+    if mode == 'NIRSpec_G140Mf100': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G140M-f100_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+        
+    if mode == 'NIRSpec_G140Hf070': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G140H-f070_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+        
+    if mode == 'NIRSpec_G140Mf070': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRSpec_G140M-f070_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1
+        
+    # For Webb NIRISS SOSS
+    if mode == 'NIRISS_SOSSo1': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRISS_SOSSo1_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+        
+    if mode == 'NIRISS_SOSSo2': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRISS_SOSSo2_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+
+    # For Webb NIRCam
+    if mode == 'NIRCam_F322W2': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRCam_F322W2_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+        
+    if mode == 'NIRCam_F444': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_NIRCam_F444_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+        
+    # For Webb MIRI
+    if mode == 'MIRI_LRS': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/JWST_MIRI_LRS_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+
+    ## Photometric Modes    
+    # Spitzer IRAC
+    if mode == 'IRAC_Ch1': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/Spitzer_IRAC_Ch1_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+
+    if mode == 'IRAC_Ch2': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/Spitzer_IRAC_Ch2_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1            
+
+    # TESS    
+    if mode == 'TESS': 
+        sav = pd.read_csv(os.path.join(dirsen,'sensitivity_files/TESS_throughput.csv'))
+        wssens = sav['wave'].values
+        sensitivity = sav['tp'].values
+        wdel = 1    
+        
     widek = np.arange(len(wsdata))
     wsmode = wssens
     wsmode = np.concatenate((np.array([wsmode[0] - wdel - wdel, wsmode[0] - wdel]),
@@ -355,8 +458,6 @@ def limb_dark_fit(mode, wsdata, M_H, Teff, logg, dirsen, ld_model='1D'):
 
     return uLD, c1, c2, c3, c4, cp1, cp2, cp3, cp4, aLD, bLD
 
-
-
 def int_tabulated(X, F, sort=False):
     Xsegments = len(X) - 1
 
@@ -402,15 +503,16 @@ def quadratic_limb_darkening(x, aLD=0.0, bLD=0.0):
     model = 1. - aLD * (1. - x) - bLD * (1. - x) ** (4. / 2.)
     return model
 
+
 if __name__ == '__main__': 
     M_H = 0.1
     Teff = 6545
     logg = 4.2
 
-    grating = 'G141'
-    wsdata = 
+    mode = 'WFC3_G141'
+    wsdata = np.arange(11100,11200,0.5)
 
-    dirsen = ''
+    dirsen = '/Users/iz19726/Downloads/LD_data/' #Where the Zenodo Download has been placed
     ld_model = '1D'
 
-    limb_dark_fit(grating, wsdata, M_H, Teff, logg, dirsen, ld_model='1D')
\ No newline at end of file
+    limb_dark_fit(mode, wsdata, M_H, Teff, logg, dirsen, ld_model='1D')
\ No newline at end of file


Test patch:
diff --git a/Notebooks/ld_tests.ipynb b/Notebooks/ld_tests.ipynb
deleted file mode 100644
index 250e002..0000000
--- a/Notebooks/ld_tests.ipynb
+++ /dev/null
@@ -1,529 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Limb-darkening functions"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import os\n",
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "from scipy.io import readsav\n",
-    "from scipy.interpolate import interp1d, splev, splrep\n",
-    "from astropy.modeling.models import custom_model\n",
-    "from astropy.modeling.fitting import LevMarLSQFitter"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### What information we need to get the coefficients\n",
-    "- Stellar parameters\n",
-    "    - Teff\n",
-    "    - \\[Fe/H\\]\n",
-    "    - log(g)\n",
-    "- mode used in the observation\n",
-    "- wavelength array of the wanted coefficients\n",
-    "- Which model grid\n",
-    "    - 1D\n",
-    "    - 3D"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def limb_dark_fit(mode, wsdata, M_H, Teff, logg, dirsen, ld_model='1D'):\n",
-    "    \"\"\"\n",
-    "    Calculates stellar limb-darkening coefficients for a given wavelength bin.\n",
-    "\n",
-    "    Modes Currently Supported:\n",
-    "    HST STIS G750L, G750M, G430L gratings\n",
-    "    HST WFC3 UVIS/G280, IR/G102, IR/G141 grisms\n",
-    "    JWST NIRSpec Prism, G395H\n",
-    "\n",
-    "    Procedure from Sing et al. (2010, A&A, 510, A21).\n",
-    "    Uses 3D limb darkening from Magic et al. (2015, A&A, 573, 90).\n",
-    "    Uses photon FLUX Sum over (lambda*dlamba).\n",
-    "    :param mode: string; mode to use ('G430L','G750L','G750M', 'G280', 'G102', 'G141', 'NIRSpec_Prism', 'NIRSpec_G395H')\n",
-    "    :param wsdata: array; data wavelength solution\n",
-    "    :param M_H: float; stellar metallicity\n",
-    "    :param Teff: float; stellar effective temperature (K)\n",
-    "    :param logg: float; stellar gravity\n",
-    "    :param dirsen: string; path to main limb darkening directory\n",
-    "    :param ld_model: string; '1D' or '3D', makes choice between limb darkening models; default is 1D\n",
-    "    :return: uLD: float; linear limb darkening coefficient\n",
-    "    aLD, bLD: float; quadratic limb darkening coefficients\n",
-    "    cp1, cp2, cp3, cp4: float; three-parameter limb darkening coefficients\n",
-    "    c1, c2, c3, c4: float; non-linear limb-darkening coefficients\n",
-    "    \"\"\"\n",
-    "\n",
-    "    print('You are using the', str(ld_model), 'limb darkening models.')\n",
-    "\n",
-    "    if ld_model == '1D':\n",
-    "\n",
-    "        direc = os.path.join(dirsen, 'Kurucz')\n",
-    "\n",
-    "        print('Current Directories Entered:')\n",
-    "        print('  ' + dirsen)\n",
-    "        print('  ' + direc)\n",
-    "\n",
-    "        # Select metallicity\n",
-    "        M_H_Grid = np.array([-0.1, -0.2, -0.3, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0, -3.5, -4.0, -4.5, -5.0, 0.0, 0.1, 0.2, 0.3, 0.5, 1.0])\n",
-    "        M_H_Grid_load = np.array([0, 1, 2, 3, 5, 7, 8, 9, 10, 11, 12, 13, 14, 17, 20, 21, 22, 23, 24])\n",
-    "        optM = (abs(M_H - M_H_Grid)).argmin()\n",
-    "        MH_ind = M_H_Grid_load[optM]\n",
-    "\n",
-    "        # Determine which model is to be used, by using the input metallicity M_H to figure out the file name we need\n",
-    "        direc = 'Kurucz'\n",
-    "        file_list = 'kuruczlist.sav'\n",
-    "        sav1 = readsav(os.path.join(dirsen, file_list))\n",
-    "        model = bytes.decode(sav1['li'][MH_ind])  # Convert object of type \"byte\" to \"string\"\n",
-    "\n",
-    "        # Select Teff and subsequently logg\n",
-    "        Teff_Grid = np.array([3500, 3750, 4000, 4250, 4500, 4750, 5000, 5250, 5500, 5750, 6000, 6250, 6500])\n",
-    "        optT = (abs(Teff - Teff_Grid)).argmin()\n",
-    "\n",
-    "        logg_Grid = np.array([4.0, 4.5, 5.0])\n",
-    "        optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        if logg_Grid[optG] == 4.0:\n",
-    "            Teff_Grid_load = np.array([8, 19, 30, 41, 52, 63, 74, 85, 96, 107, 118, 129, 138])\n",
-    "\n",
-    "        elif logg_Grid[optG] == 4.5:\n",
-    "            Teff_Grid_load = np.array([9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 119, 129, 139])\n",
-    "\n",
-    "        elif logg_Grid[optG] == 5.0:\n",
-    "            Teff_Grid_load = np.array([10, 21, 32, 43, 54, 65, 76, 87, 98, 109, 120, 130, 140])\n",
-    "\n",
-    "        # Where in the model file is the section for the Teff we want? Index T_ind tells us that.\n",
-    "        T_ind = Teff_Grid_load[optT]\n",
-    "        header_rows = 3    #  How many rows in each section we ignore for the data reading\n",
-    "        data_rows = 1221   # How  many rows of data we read\n",
-    "        line_skip_data = (T_ind + 1) * header_rows + T_ind * data_rows   # Calculate how many lines in the model file we need to skip in order to get to the part we need (for the Teff we want).\n",
-    "        line_skip_header = T_ind * (data_rows + header_rows)\n",
-    "\n",
-    "        # Read the header, in case we want to have the actual Teff, logg and M_H info.\n",
-    "        # headerinfo is a pandas object.\n",
-    "        headerinfo = pd.read_csv(os.path.join(dirsen, direc, model), delim_whitespace=True, header=None,\n",
-    "                                 skiprows=line_skip_header, nrows=1)\n",
-    "\n",
-    "        Teff_model = headerinfo[1].values[0]\n",
-    "        logg_model = headerinfo[3].values[0]\n",
-    "        MH_model = headerinfo[6].values[0]\n",
-    "        MH_model = float(MH_model[1:-1])\n",
-    "\n",
-    "        print('\\nClosest values to your inputs:')\n",
-    "        print('Teff: ', Teff_model)\n",
-    "        print('M_H: ', MH_model)\n",
-    "        print('log(g): ', logg_model)\n",
-    "\n",
-    "        # Read the data; data is a pandas object.\n",
-    "        data = pd.read_csv(os.path.join(dirsen, direc, model), delim_whitespace=True, header=None,\n",
-    "                              skiprows=line_skip_data, nrows=data_rows)\n",
-    "\n",
-    "        # Unpack the data\n",
-    "        ws = data[0].values * 10   # Import wavelength data\n",
-    "        f0 = data[1].values / (ws * ws)\n",
-    "        f1 = data[2].values * f0 / 100000.\n",
-    "        f2 = data[3].values * f0 / 100000.\n",
-    "        f3 = data[4].values * f0 / 100000.\n",
-    "        f4 = data[5].values * f0 / 100000.\n",
-    "        f5 = data[6].values * f0 / 100000.\n",
-    "        f6 = data[7].values * f0 / 100000.\n",
-    "        f7 = data[8].values * f0 / 100000.\n",
-    "        f8 = data[9].values * f0 / 100000.\n",
-    "        f9 = data[10].values * f0 / 100000.\n",
-    "        f10 = data[11].values * f0 / 100000.\n",
-    "        f11 = data[12].values * f0 / 100000.\n",
-    "        f12 = data[13].values * f0 / 100000.\n",
-    "        f13 = data[14].values * f0 / 100000.\n",
-    "        f14 = data[15].values * f0 / 100000.\n",
-    "        f15 = data[16].values * f0 / 100000.\n",
-    "        f16 = data[17].values * f0 / 100000.\n",
-    "\n",
-    "        # Make single big array of them\n",
-    "        fcalc = np.array([f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16])\n",
-    "        phot1 = np.zeros(fcalc.shape[0])\n",
-    "\n",
-    "        # Define mu\n",
-    "        mu = np.array([1.000, .900, .800, .700, .600, .500, .400, .300, .250, .200, .150, .125, .100, .075, .050, .025, .010])\n",
-    "\n",
-    "        # Passed on to main body of function are: ws, fcalc, phot1, mu\n",
-    "\n",
-    "    elif ld_model == '3D':\n",
-    "\n",
-    "        direc = os.path.join(dirsen, '3DGrid')\n",
-    "\n",
-    "        print('Current Directories Entered:')\n",
-    "        print('  ' + dirsen)\n",
-    "        print('  ' + direc)\n",
-    "\n",
-    "        # Select metallicity\n",
-    "        M_H_Grid = np.array([-3.0, -2.0, -1.0, 0.0])  # Available metallicity values in 3D models\n",
-    "        M_H_Grid_load = ['30', '20', '10', '00']  # The according identifiers to individual available M_H values\n",
-    "        optM = (abs(M_H - M_H_Grid)).argmin()  # Find index at which the closes M_H values from available values is to the input M_H.\n",
-    "\n",
-    "        # Select Teff\n",
-    "        Teff_Grid = np.array([4000, 4500, 5000, 5500, 5777, 6000, 6500, 7000])  # Available Teff values in 3D models\n",
-    "        optT = (abs(Teff - Teff_Grid)).argmin()  # Find index at which the Teff values is, that is closest to input Teff.\n",
-    "\n",
-    "        # Select logg, depending on Teff. If several logg possibilities are given for one Teff, pick the one that is\n",
-    "        # closest to user input (logg).\n",
-    "\n",
-    "        if Teff_Grid[optT] == 4000:\n",
-    "            logg_Grid = np.array([1.5, 2.0, 2.5])\n",
-    "            optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 4500:\n",
-    "            logg_Grid = np.array([2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])\n",
-    "            optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 5000:\n",
-    "            logg_Grid = np.array([2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])\n",
-    "            optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 5500:\n",
-    "            logg_Grid = np.array([3.0, 3.5, 4.0, 4.5, 5.0])\n",
-    "            optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 5777:\n",
-    "            logg_Grid = np.array([4.4])\n",
-    "            optG = 0\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 6000:\n",
-    "            logg_Grid = np.array([3.5, 4.0, 4.5])\n",
-    "            optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 6500:\n",
-    "            logg_Grid = np.array([4.0, 4.5])\n",
-    "            optG = (abs(logg - logg_Grid)).argmin()\n",
-    "\n",
-    "        elif Teff_Grid[optT] == 7000:\n",
-    "            logg_Grid = np.array([4.5])\n",
-    "            optG = 0\n",
-    "\n",
-    "        # Select Teff and Log g. Mtxt, Ttxt and Gtxt are then put together as string to load correct files.\n",
-    "        Mtxt = M_H_Grid_load[optM]\n",
-    "        Ttxt = \"{:2.0f}\".format(Teff_Grid[optT] / 100)\n",
-    "        if Teff_Grid[optT] == 5777:\n",
-    "            Ttxt = \"{:4.0f}\".format(Teff_Grid[optT])\n",
-    "        Gtxt = \"{:2.0f}\".format(logg_Grid[optG] * 10)\n",
-    "\n",
-    "        #\n",
-    "        file = 'mmu_t' + Ttxt + 'g' + Gtxt + 'm' + Mtxt + 'v05.flx'\n",
-    "        print('Filename:', file)\n",
-    "\n",
-    "        # Read data from IDL .sav file\n",
-    "        sav = readsav(os.path.join(direc, file))  # readsav reads an IDL .sav file\n",
-    "        ws = sav['mmd'].lam[0]  # read in wavelength\n",
-    "        flux = sav['mmd'].flx  # read in flux\n",
-    "        Teff_model = Teff_Grid[optT]\n",
-    "        logg_model = logg_Grid[optG]\n",
-    "        MH_model = str(M_H_Grid[optM])\n",
-    "\n",
-    "        print('\\nClosest values to your inputs:')\n",
-    "        print('Teff: ', Teff_model)\n",
-    "        print('M_H: ', MH_model)\n",
-    "        print('log(g): ', logg_model)\n",
-    "\n",
-    "        f0 = flux[0]\n",
-    "        f1 = flux[1]\n",
-    "        f2 = flux[2]\n",
-    "        f3 = flux[3]\n",
-    "        f4 = flux[4]\n",
-    "        f5 = flux[5]\n",
-    "        f6 = flux[6]\n",
-    "        f7 = flux[7]\n",
-    "        f8 = flux[8]\n",
-    "        f9 = flux[9]\n",
-    "        f10 = flux[10]\n",
-    "\n",
-    "        # Make single big array of them\n",
-    "        fcalc = np.array([f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10])\n",
-    "        phot1 = np.zeros(fcalc.shape[0])\n",
-    "\n",
-    "        # Mu from grid\n",
-    "        # 0.00000    0.0100000    0.0500000     0.100000     0.200000     0.300000   0.500000     0.700000     0.800000     0.900000      1.00000\n",
-    "        mu = sav['mmd'].mu\n",
-    "\n",
-    "        # Passed on to main body of function are: ws, fcalc, phot1, mu\n",
-    "\n",
-    "    ### Load response function and interpolate onto kurucz model grid\n",
-    "\n",
-    "    # FOR Hubble STIS\n",
-    "    if mode == 'G430L':\n",
-    "        sav = readsav(os.path.join(dirsen, 'G430L.STIS.sensitivity.sav'))  # wssens,sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 3\n",
-    "\n",
-    "    if mode == 'G750M':\n",
-    "        sav = readsav(os.path.join(dirsen, 'G750M.STIS.sensitivity.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 0.554\n",
-    "\n",
-    "    if mode == 'G750L':\n",
-    "        sav = readsav(os.path.join(dirsen, 'G750L.STIS.sensitivity.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 4.882\n",
-    "\n",
-    "    # FOR Hubble WFC3\n",
-    "    if mode == 'G141':  # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html\n",
-    "        sav = readsav(os.path.join(dirsen, 'G141.WFC3.sensitivity.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 1\n",
-    "\n",
-    "    if mode == 'G102':  # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html\n",
-    "        sav = readsav(os.path.join(dirsen, 'G102.WFC3.sensitivity.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 1\n",
-    "\n",
-    "    if mode == 'G280':  # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html\n",
-    "        sav = readsav(os.path.join(dirsen, 'G280.WFC3.sensitivity.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 1\n",
-    "\n",
-    "    # FOR Webb NIRSpec\n",
-    "    if mode == 'NIRSpec_Prism':  # Provided by D.K. Sing (JHU)\n",
-    "        sav = readsav(os.path.join(dirsen, 'NIRSpec.prism.sensitivity.pandeia.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 1\n",
-    "\n",
-    "    if mode == 'NIRSpec_G395H':  # Provided by D.K. Sing (JHU)\n",
-    "        sav = readsav(os.path.join(dirsen, 'NIRSpec.G395H.sensitivity.pandeia.sav'))  # wssens, sensitivity\n",
-    "        wssens = sav['wssens']\n",
-    "        sensitivity = sav['sensitivity']\n",
-    "        wdel = 1        \n",
-    "\n",
-    "    widek = np.arange(len(wsdata))\n",
-    "    wsmode = wssens\n",
-    "    wsmode = np.concatenate((np.array([wsmode[0] - wdel - wdel, wsmode[0] - wdel]),\n",
-    "                            wsmode,\n",
-    "                            np.array([wsmode[len(wsmode) - 1] + wdel,\n",
-    "                                      wsmode[len(wsmode) - 1] + wdel + wdel])))\n",
-    "\n",
-    "    respoutmode = sensitivity / np.max(sensitivity)\n",
-    "    respoutmode = np.concatenate((np.zeros(2), respoutmode, np.zeros(2)))\n",
-    "\n",
-    "    inter_resp = interp1d(wsmode, respoutmode, bounds_error=False, fill_value=0)\n",
-    "    respout = inter_resp(ws)  # interpolate sensitivity curve onto model wavelength grid\n",
-    "\n",
-    "    wsdata = np.concatenate((np.array([wsdata[0] - wdel - wdel, wsdata[0] - wdel]), wsdata,\n",
-    "                             np.array([wsdata[len(wsdata) - 1] + wdel, wsdata[len(wsdata) - 1] + wdel + wdel])))\n",
-    "    respwavebin = wsdata / wsdata * 0.0\n",
-    "    widek = widek + 2  # need to add two indicies to compensate for padding with 2 zeros\n",
-    "    respwavebin[widek] = 1.0\n",
-    "    data_resp = interp1d(wsdata, respwavebin, bounds_error=False, fill_value=0)\n",
-    "    reswavebinout = data_resp(ws)  # interpolate data onto model wavelength grid\n",
-    "\n",
-    "    # Integrate over the spectra to make synthetic photometric points.\n",
-    "    for i in range(fcalc.shape[0]):  # Loop over spectra at diff angles\n",
-    "        fcal = fcalc[i, :]\n",
-    "        Tot = int_tabulated(ws, ws * respout * reswavebinout)\n",
-    "        phot1[i] = (int_tabulated(ws, ws * respout * reswavebinout * fcal, sort=True)) / Tot\n",
-    "\n",
-    "    if ld_model == '1D':\n",
-    "        yall = phot1 / phot1[0]\n",
-    "    elif ld_model == '3D':\n",
-    "        yall = phot1 / phot1[10]\n",
-    "\n",
-    "    Co = np.zeros((6, 4))   # NOT-REUSED\n",
-    "\n",
-    "    A = [0.0, 0.0, 0.0, 0.0]  # c1, c2, c3, c4      # NOT-REUSED\n",
-    "    x = mu[1:]     # wavelength\n",
-    "    y = yall[1:]   # flux\n",
-    "    weights = x / x   # NOT-REUSED\n",
-    "\n",
-    "    # Start fitting the different models\n",
-    "    fitter = LevMarLSQFitter()\n",
-    "\n",
-    "    # Fit a four parameter non-linear limb darkening model and get fitted variables, c1, c2, c3, c4.\n",
-    "    corot_4_param = nonlinear_limb_darkening()\n",
-    "    corot_4_param = fitter(corot_4_param, x, y)\n",
-    "    c1, c2, c3, c4 = corot_4_param.parameters\n",
-    "\n",
-    "    # Fit a three parameter non-linear limb darkening model and get fitted variables, cp2, cp3, cp4 (cp1 = 0).\n",
-    "    corot_3_param = nonlinear_limb_darkening()\n",
-    "    corot_3_param.c0.fixed = True  # 3 param is just 4 param with c0 = 0.0\n",
-    "    corot_3_param = fitter(corot_3_param, x, y)\n",
-    "    cp1, cp2, cp3, cp4 = corot_3_param.parameters\n",
-    "\n",
-    "    # Fit a quadratic limb darkening model and get fitted parameters aLD and bLD.\n",
-    "    quadratic = quadratic_limb_darkening()\n",
-    "    quadratic = fitter(quadratic, x, y)\n",
-    "    aLD, bLD = quadratic.parameters\n",
-    "\n",
-    "    # Fit a linear limb darkening model and get fitted variable uLD.\n",
-    "    linear = nonlinear_limb_darkening()\n",
-    "    linear.c0.fixed = True\n",
-    "    linear.c2.fixed = True\n",
-    "    linear.c3.fixed = True\n",
-    "    linear = fitter(linear, x, y)\n",
-    "    uLD = linear.c1.value\n",
-    "\n",
-    "    print('\\nLimb darkening parameters:')\n",
-    "    print(\"4param \\t{:0.8f}\\t{:0.8f}\\t{:0.8f}\\t{:0.8f}\".format(c1, c2, c3, c4))\n",
-    "    print(\"3param \\t{:0.8f}\\t{:0.8f}\\t{:0.8f}\".format(cp2, cp3, cp4))\n",
-    "    print(\"Quad \\t{:0.8f}\\t{:0.8f}\".format(aLD, bLD))\n",
-    "    print(\"Linear \\t{:0.8f}\".format(uLD))\n",
-    "\n",
-    "    return uLD, c1, c2, c3, c4, cp1, cp2, cp3, cp4, aLD, bLD\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def int_tabulated(X, F, sort=False):\n",
-    "    Xsegments = len(X) - 1\n",
-    "\n",
-    "    # Sort vectors into ascending order.\n",
-    "    if not sort:\n",
-    "        ii = np.argsort(X)\n",
-    "        X = X[ii]\n",
-    "        F = F[ii]\n",
-    "\n",
-    "    while (Xsegments % 4) != 0:\n",
-    "        Xsegments = Xsegments + 1\n",
-    "\n",
-    "    Xmin = np.min(X)\n",
-    "    Xmax = np.max(X)\n",
-    "\n",
-    "    # Uniform step size.\n",
-    "    h = (Xmax + 0.0 - Xmin) / Xsegments\n",
-    "    # Compute the interpolates at Xgrid.\n",
-    "    # x values of interpolates >> Xgrid = h * FINDGEN(Xsegments + 1L) + Xmin\n",
-    "    z = splev(h * np.arange(Xsegments + 1) + Xmin, splrep(X, F))\n",
-    "\n",
-    "    # Compute the integral using the 5-point Newton-Cotes formula.\n",
-    "    ii = (np.arange((len(z) - 1) / 4, dtype=int) + 1) * 4\n",
-    "\n",
-    "    return np.sum(2.0 * h * (7.0 * (z[ii - 4] + z[ii]) + 32.0 * (z[ii - 3] + z[ii - 1]) + 12.0 * z[ii - 2]) / 45.0)\n",
-    "\n",
-    "\n",
-    "@custom_model\n",
-    "def nonlinear_limb_darkening(x, c0=0.0, c1=0.0, c2=0.0, c3=0.0):\n",
-    "    \"\"\"\n",
-    "    Define non-linear limb darkening model with four parameters c0, c1, c2, c3.\n",
-    "    \"\"\"\n",
-    "    model = (1. - (c0 * (1. - x ** (1. / 2)) + c1 * (1. - x ** (2. / 2)) + c2 * (1. - x ** (3. / 2)) + c3 *\n",
-    "                   (1. - x ** (4. / 2))))\n",
-    "    return model\n",
-    "\n",
-    "\n",
-    "@custom_model\n",
-    "def quadratic_limb_darkening(x, aLD=0.0, bLD=0.0):\n",
-    "    \"\"\"\n",
-    "    Define linear limb darkening model with parameters aLD and bLD.\n",
-    "    \"\"\"\n",
-    "    model = 1. - aLD * (1. - x) - bLD * (1. - x) ** (4. / 2.)\n",
-    "    return model"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "_uLD, c1, c2, c3, c4, _cp1, _cp2, _cp3, _cp4, _aLD, _bLD = limb_dark_fit(grat, wavelength, M_H, Teff,\n",
-    "                                                                         logg, limbDir, ld_model)\n",
-    "\n",
-    "print('\\nc1 = {}'.format(c1))\n",
-    "print('c2 = {}'.format(c2))\n",
-    "print('c3 = {}'.format(c3))\n",
-    "print('c4 = {}'.format(c4))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.7.3"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}