diff --git a/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb b/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb
index 29ae67c92..8862c721b 100644
--- a/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb
+++ b/notebooks/transit_spectroscopy_notebook/Exoplanet_Transmission_Spectra_JWST.ipynb
@@ -6,7 +6,7 @@
"source": [
"# BOTS Time Series Observations\n",
"\n",
- "**Use case:** Bright Object Time Series; extracting exoplanet spectra.
\n",
+ " **Use case:** Bright Object Time Series; extracting exoplanet spectra.
\n",
"**Data:** JWST simulated NIRSpec data from ground-based campaign; GJ436b spectra from the Goyal et al. (2018).
\n",
"**Tools:** scikit, lmfit, scipy, matplotlip, astropy, pandas.
\n",
"**Cross-intrument:** .
\n",
@@ -61,22 +61,19 @@
"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 +94,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 +134,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 +165,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",
- " dest = shutil.move(fn,save_directory+'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 +193,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 +276,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.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 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 +336,39 @@
"outputs": [],
"source": [
"all_spec.shape\n",
- "y_lower = 11 # Lower extraction aperature\n",
- "y_upper = 18 # Upper extraction aperature\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",
+ "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",
+ "# 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 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.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()"
]
},
{
@@ -394,14 +384,14 @@
"metadata": {},
"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",
- "plt.xlabel('Wavelength ($\\mu$m)')\n",
+ "fig, axs = plt.subplots()\n",
+ "f = plt.plot(wsdata_all, all_spec_1D[:, 0], linewidth=2, zorder=0) # overplot Transit model at data\n",
+ "plt.xlabel(r'Wavelength ($\\mu$m)')\n",
"plt.ylabel('Flux (e-)')\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()"
]
},
@@ -429,18 +419,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",
- "\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"
+ "# 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",
+ "\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"
]
},
{
@@ -464,17 +453,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",
+ "fig, axs = plt.subplots()\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 +500,14 @@
"metadata": {},
"outputs": [],
"source": [
- "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.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"
+ "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.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)"
]
},
{
@@ -537,30 +525,30 @@
"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",
- "plt.xlabel('Wavelength ($\\mu$m)')\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",
+ "plt.xlabel(r'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",
+ "plt.title(r'Time-series Transit Light Curve $\\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2]) + r'] $\\mu$m')\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
@@ -585,21 +573,21 @@
"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",
+ "plt.title(r'Time-series Transit Light Curve $\\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2]) + r'] $\\mu$m')\n",
"axs.xaxis.set_major_locator(ticker.MultipleLocator(0.01))\n",
"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()"
@@ -722,7 +710,7 @@
"\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",
+ " 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",
@@ -744,7 +732,7 @@
"\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",
+ " skiprows=line_skip_data, nrows=data_rows)\n",
"\n",
" # Unpack the data\n",
" ws = data[0].values * 10 # Import wavelength data\n",
@@ -873,7 +861,7 @@
"\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",
+ " # Load response function and interpolate onto kurucz model grid\n",
"\n",
" # FOR STIS\n",
" if grating == 'G430L':\n",
@@ -920,7 +908,6 @@
" sensitivity = sav['sensitivity']\n",
" wdel = 12\n",
"\n",
- "\n",
" widek = np.arange(len(wsdata))\n",
" wsHST = wssens\n",
" wsHST = np.concatenate((np.array([wsHST[0] - wdel - wdel, wsHST[0] - wdel]),\n",
@@ -952,12 +939,12 @@
" elif ld_model == '3D':\n",
" yall = phot1 / phot1[10]\n",
"\n",
- " Co = np.zeros((6, 4)) # NOT-REUSED\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",
+ " # 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",
+ " # weights = x / x # NOT-REUSED\n",
"\n",
" # Start fitting the different models\n",
" fitter = LevMarLSQFitter()\n",
@@ -1019,8 +1006,7 @@
" # 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"
+ " 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)"
]
},
{
@@ -1069,7 +1055,7 @@
" dmumax = 1.0\n",
"\n",
" while (dmumax > fac * 1.e-3) and (nr <= 131072):\n",
- " #print(nr)\n",
+ " # print(nr)\n",
" mulimbp = mulimb\n",
" nr = nr * 2\n",
" dt = 0.5 * np.pi / nr\n",
@@ -1096,7 +1082,7 @@
" if len(ix1) == 0:\n",
" ix1 = -1\n",
"\n",
- " #print(ix1)\n",
+ " # print(ix1)\n",
" # python cannot index on single values so you need to use atlest_1d for the below to work when mulimb is a single value\n",
" dmumax = np.max(np.abs(np.atleast_1d(mulimb)[ix1] - np.atleast_1d(mulimbp)[ix1]) / (\n",
" np.atleast_1d(mulimb)[ix1] + np.atleast_1d(mulimbp)[ix1]))\n",
@@ -1127,11 +1113,10 @@
" nb = len(np.atleast_1d(b0))\n",
" muo1 = np.zeros(nb)\n",
"\n",
- "\n",
" for i in range(nb):\n",
" # substitute z=b0(i) to shorten expressions\n",
" z = np.atleast_1d(b0)[i]\n",
- " #z = z.value # stripping it of astropy units\n",
+ " # z = z.value # stripping it of astropy units\n",
" if z >= 1+w:\n",
" muo1[i] = 1.0\n",
" continue\n",
@@ -1151,7 +1136,7 @@
" muo1[i] = 1 - w ** 2\n",
" continue\n",
"\n",
- " return muo1\n"
+ " return muo1"
]
},
{
@@ -1167,22 +1152,22 @@
"metadata": {},
"outputs": [],
"source": [
- "#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",
- " 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",
+ "# 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",
+ " 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",
" return (y-model)/err\n",
- " #return np.sum((y-model)**2/err**2)"
+ " # return np.sum((y-model)**2/err**2)"
]
},
{
@@ -1198,10 +1183,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 +1203,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,55 +1226,55 @@
"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",
+ "mulimb0, mulimbf = occultnl(rl, c1, c2, c3, c4, b0) # Mandel & Agol non-linear limb darkened transit model\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",
+ "plt.title(r'Time-series Transit Light Curve $\\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2]) + r'] $\\mu$m')\n",
"ax1.xaxis.set_major_locator(ticker.MultipleLocator(0.01))\n",
"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",
- "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",
+ "# 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",
+ "# print chi^2 value\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 +1332,17 @@
"metadata": {},
"outputs": [],
"source": [
- "k0 = 113 #98 #100\n",
- "kend = 392 #422\n",
- "wk = 15\n",
+ "k0 = 113 # 98 #100\n",
+ "kend = 392 # 422\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))"
]
@@ -1379,98 +1363,92 @@
"metadata": {},
"outputs": [],
"source": [
- "k=k0 #wavelength to start\n",
- "#--------------------------------------------------------------------------\n",
- "#Loop over wavelength bins and fit for each one\n",
- "for bin in range(0,number_of_bins):\n",
+ "k = k0 # wavelength to start\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",
- " _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",
+ " 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",
" print('c2 = {}'.format(c2))\n",
" print('c3 = {}'.format(c3))\n",
" print('c4 = {}'.format(c4))\n",
" print('')\n",
- " #u = [c1,c2,c3,c4] #limb darkening coefficients\n",
- " u = [aLD,bLD]\n",
- " #---------------------------------------------------------\n",
+ " # u = [c1,c2,c3,c4] # limb darkening coefficients\n",
+ " # u = [aLD, bLD]\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",
+ " # Setup LMFIT\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",
- " 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",
+ " # 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",
" 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 = mini.minimize(method='emcee')\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",
- " print('Re-Fitting Bin',bin,' Wavelength =',np.mean(wsdata)/1E4, ' Range= [',wave1,':',wave2,']')\n",
- " #--------------------------------------------------------------------------\n",
+ " print('Re-Fitting Bin', bin, ' Wavelength =', np.mean(wsdata)/1E4, ' Range= [', wave1, ':', wave2, ']')\n",
+ " # --------------------------------------------------------------------------\n",
" print(\"\")\n",
- " print(\"redchi\",result.redchi)\n",
- " print(\"chi2\",result.chisqr)\n",
- " print(\"nfree\",result.nfree)\n",
- " print(\"bic\",result.bic)\n",
- " print(\"aic\",result.aic)\n",
+ " print(\"redchi\", result.redchi)\n",
+ " print(\"chi2\", result.chisqr)\n",
+ " print(\"nfree\", result.nfree)\n",
+ " print(\"bic\", result.bic)\n",
+ " print(\"aic\", result.aic)\n",
" print(\"L-M FIT Variable\")\n",
" print(lmfit.fit_report(result.params))\n",
" text_file = open(save_directory+'JWST_NIRSpec_Prism_fit_light_curve_bin'+str(bin)+'_statistics.txt', \"w\")\n",
@@ -1482,135 +1460,137 @@
" n = text_file.write(lmfit.fit_report(result.params))\n",
" # file-output.py\n",
"\n",
- " #Update with best-fit parameters\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",
- " 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",
+ " 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 = 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",
- " 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",
+ " # 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):\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",
" n = text_file.write(\"\\nTransit depth measured error (ppm) : \"+str(2E6*result.params['rprs'].value*result.params['rprs'].stderr))\n",
" text_file.close()\n",
- " depth[bin]=1E6*result.params['rprs'].value**2\n",
- " depth_err[bin]=2E6*result.params['rprs'].value*result.params['rprs'].stderr\n",
- "\n",
- " sig_r[bin]=sigrednoise*1E6\n",
- " sig_w[bin]=sigwhite*1E6\n",
- " #--------------------------------------------------------------------------\n",
- " #---------------------------------------------------------\n",
- " #Write Fit Spectra to ascii file\n",
- " ascii_data = Table([wsd, werr, rprs, rerr,depth,depth_err,sig_w,sig_r,beta], names=['Wavelength Center (um)', 'Wavelength half-width (um)','Rp/Rs','Rp/Rs 1-sigma error','Transit Depth (ppm)','Transit Depth error','Sigma_white (ppm)','Sigma_red (ppm)','Beta Rednoise Inflation factor'])\n",
- " 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",
+ " depth[bin] = 1E6*result.params['rprs'].value**2\n",
+ " depth_err[bin] = 2E6*result.params['rprs'].value*result.params['rprs'].stderr\n",
+ "\n",
+ " sig_r[bin] = sigrednoise*1E6\n",
+ " sig_w[bin] = sigwhite*1E6\n",
+ " # --------------------------------------------------------------------------\n",
+ " # ---------------------------------------------------------\n",
+ " # Write Fit Spectra to ascii file\n",
+ " ascii_data = Table([wsd, werr, rprs, rerr, depth, depth_err, sig_w, sig_r, beta], names=['Wavelength Center (um)', 'Wavelength half-width (um)', 'Rp/Rs', 'Rp/Rs 1-sigma error', 'Transit Depth (ppm)', 'Transit Depth error', 'Sigma_white (ppm)', 'Sigma_red (ppm)', 'Beta Rednoise Inflation factor'])\n",
+ " 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",
" \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(x,y/p['f0'].value,s=msize*0.75, linewidth=1,zorder=0,alpha=0.5, marker='+',edgecolors='blue')\n",
- " ax1.plot(x,final_model/p['f0'].value, linewidth=1,color='orange',alpha=0.8,zorder=15) #overplot Transit model at data\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(x, y/p['f0'].value, s=msize*0.75, linewidth=1, zorder=0, alpha=0.5, marker='+', edgecolors='blue')\n",
+ " ax1.plot(x, final_model/p['f0'].value, linewidth=1, color='orange', alpha=0.8, zorder=15) # overplot Transit model at data\n",
" ax1.xaxis.set_ticklabels([])\n",
" plt.ylabel('Relative Flux')\n",
- " plt.title('Time-series Transit Light Curve $\\lambda=$['+str(wave1)+':'+str(wave2)+'] $\\mu$m')\n",
+ " plt.title(r'Time-series Transit Light Curve $\\lambda=$['+str(wave1)+':'+str(wave2) + r'] $\\mu$m')\n",
" ax1.xaxis.set_major_locator(ticker.MultipleLocator(0.01))\n",
" 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",
" fig.add_subplot(ax1)\n",
- " #Residual\n",
- " ax2=fig.add_subplot(gs[2,:])\n",
- " ax2.scatter(x,residppm, s=msize*0.75,linewidth=1 ,alpha=0.5, marker='+',edgecolors='blue',zorder=0) #overplot Transit model at data\n",
- " wsb, wsb_bin_edges,binnumber = stats.binned_statistic(bjd,residppm, bins=256)\n",
- " plt.scatter(wsb_bin_edges[1:],wsb, linewidth=2,alpha=0.75,facecolors='orange',edgecolors='none', marker='o',zorder=25)\n",
+ " # Residual\n",
+ " ax2 = fig.add_subplot(gs[2, :])\n",
+ " ax2.scatter(x, residppm, s=msize*0.75, linewidth=1, alpha=0.5, marker='+', edgecolors='blue', zorder=0) # overplot Transit model at data\n",
+ " wsb, wsb_bin_edges, binnumber = stats.binned_statistic(bjd, residppm, 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",
- " plt.plot([bjd.min(),bjd.max()],[0,0],color='black',zorder=10)\n",
- " plt.plot([bjd.min(),bjd.max()],[sigma,sigma],linestyle='--',color='black',zorder=15)\n",
- " plt.plot([bjd.min(),bjd.max()],[-sigma,-sigma],linestyle='--',color='black',zorder=20)\n",
+ " plt.plot([bjd.min(), bjd.max()], [0, 0], color='black', zorder=10)\n",
+ " plt.plot([bjd.min(), bjd.max()], [sigma, sigma], linestyle='--', color='black', zorder=15)\n",
+ " plt.plot([bjd.min(), bjd.max()], [-sigma, -sigma], linestyle='--', color='black', zorder=20)\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",
- " #save\n",
+ " # save\n",
" pp = PdfPages(save_directory+'JWST_NIRSpec_Prism_fit_light_curve_bin'+str(bin)+'_lightcurve.pdf')\n",
- " plt.savefig(pp,format='pdf')\n",
+ " plt.savefig(pp, format='pdf')\n",
" pp.close()\n",
" plt.clf() \n",
- " #--------------------------------------------------------------------------\n",
- " #plot systematic corrected light curve\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",
- " intransit=(b0-p['rprs'].value < 1.0E0).nonzero()\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",
- " fig,axs = plt.subplots()\n",
- " plt.scatter(x,light_curve+resid,s=msize*0.75, linewidth=1,zorder=0,alpha=0.5, marker='+',edgecolors='blue')\n",
+ " \n",
+ " # plot systematic corrected light curve\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",
+ " intransit = (b0-p['rprs'].value < 1.0E0).nonzero()\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",
+ " fig, axs = plt.subplots()\n",
+ " plt.scatter(x, light_curve+resid, s=msize*0.75, linewidth=1, zorder=0, alpha=0.5, marker='+', edgecolors='blue')\n",
" plt.xlabel('BJD')\n",
" plt.ylabel('Relative Flux')\n",
- " plt.plot(x,light_curve, linewidth=2,color='orange',alpha=0.8,zorder=15) #overplot Transit model at data\n",
+ " plt.plot(x, light_curve, linewidth=2, color='orange', alpha=0.8, zorder=15) # overplot Transit model at data\n",
" pp = PdfPages(save_directory+'JWST_NIRSpec_Prism_fit_light_curve_bin'+str(bin)+'_corrected.pdf')\n",
- " plt.savefig(pp,format='pdf')\n",
+ " plt.savefig(pp, format='pdf')\n",
" pp.close()\n",
" plt.clf()\n",
- " plt.close('all') #close all figures\n",
- " #--------------------------------------------------------------------------\n",
- " k=k+wk #step wavelength index to next bin\n",
+ " plt.close('all') # close all figures\n",
+ " # --------------------------------------------------------------------------\n",
+ " k = k + wk # step wavelength index to next bin\n",
" \n",
- " print('** Can Now View Output PDFs in ',save_directory)\n"
+ " print('** Can Now View Output PDFs in ', save_directory)"
]
},
{
@@ -1626,46 +1606,46 @@
"metadata": {},
"outputs": [],
"source": [
- "#--------------------------------------------------------------------------\n",
- "#Load Injected Transmission spectra to compare with recovered value\n",
+ "# --------------------------------------------------------------------------\n",
+ "# Load Injected Transmission spectra to compare with recovered value\n",
"\n",
- "#Download Injected Spectra\n",
- "fn_tm = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt')\n",
- "destld = shutil.move(fn_tm,save_directory+'trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt') \n",
+ "# Download Injected Spectra\n",
+ "fn_tm = download_file(r'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt')\n",
+ "destld = shutil.move(fn_tm, save_directory+'trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt') \n",
"\n",
"f = open(save_directory+'trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt', 'r')\n",
"data = np.genfromtxt(f, delimiter=' ')\n",
- "model_ws = data[:,0]\n",
- "model_spec = data[:,1]\n",
+ "model_ws = data[:, 0]\n",
+ "model_spec = data[:, 1]\n",
"\n",
- "#Read fit transit depths\n",
+ "# Read fit transit depths\n",
"data = ascii.read(save_directory+'JWST_NIRSpec_Prism_fit_transmission_spectra.csv', format='csv')\n",
- "wsd = data['Wavelength Center (um)']\n",
+ "wsd = data['Wavelength Center (um)']\n",
"werr = data['Wavelength half-width (um)']\n",
"rprs = data['Rp/Rs']\n",
"rerr = data['Rp/Rs 1-sigma error']\n",
"beta = data['Beta Rednoise Inflation factor']\n",
"\n",
- "#plot\n",
- "fig,axs = plt.subplots()\n",
- "plt.plot(model_ws,model_spec**2*1E6, linewidth=2,zorder=0,color='blue',label='Injected Spectra') #overplot Transit model at data\n",
- "plt.errorbar(wsd,rprs**2*1E6,xerr=werr,yerr=2*rerr*rprs*1E6*beta, fmt='o',zorder=5,alpha=0.4,color='orange',label='Recovered Spectra with $\\sigma_r$')\n",
- "plt.errorbar(wsd,rprs**2*1E6,xerr=werr,yerr=2*rerr*rprs*1E6, fmt='o',zorder=10,color='orange',label='Recovered Spectra')\n",
- "plt.xlabel('Wavelength ($\\mu$m)')\n",
+ "# plot\n",
+ "fig, axs = plt.subplots()\n",
+ "plt.plot(model_ws, model_spec**2*1E6, linewidth=2, zorder=0, color='blue', label='Injected Spectra') # overplot Transit model at data\n",
+ "plt.errorbar(wsd, rprs**2*1E6, xerr=werr, yerr=2*rerr*rprs*1E6*beta, fmt='o', zorder=5, alpha=0.4, color='orange', label=r'Recovered Spectra with $\\sigma_r$')\n",
+ "plt.errorbar(wsd, rprs**2*1E6, xerr=werr, yerr=2*rerr*rprs*1E6, fmt='o', zorder=10, color='orange', label='Recovered Spectra')\n",
+ "plt.xlabel(r'Wavelength ($\\mu$m)')\n",
"plt.ylabel('Transit Depth ($R_p/R_s$)$^2$ (ppm)')\n",
"axs.yaxis.set_major_locator(ticker.MultipleLocator(200))\n",
"axs.yaxis.set_minor_locator(ticker.MultipleLocator(100))\n",
"axs.xaxis.set_major_locator(ticker.MultipleLocator(0.5))\n",
"axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.1))\n",
- "axs.text(3.3,6850,'CH$_4$')\n",
- "axs.text(4.25,6700,'CO')\n",
- "axs.text(2.3,6750,'CH$_4$')\n",
- "axs.text(2.75,6550,'H$_2$O')\n",
- "plt.ylim(5700,7000)\n",
- "plt.xlim(0.9,5.25)\n",
+ "axs.text(3.3, 6850, 'CH$_4$')\n",
+ "axs.text(4.25, 6700, 'CO')\n",
+ "axs.text(2.3, 6750, 'CH$_4$')\n",
+ "axs.text(2.75, 6550, 'H$_2$O')\n",
+ "plt.ylim(5700, 7000)\n",
+ "plt.xlim(0.9, 5.25)\n",
"plt.legend(loc='lower right') \n",
"plt.show()\n",
- "plt.clf()\n"
+ "plt.clf()"
]
},
{
@@ -1696,7 +1676,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.11.5"
+ "version": "3.11.6"
}
},
"nbformat": 4,
diff --git a/notebooks/transit_spectroscopy_notebook/requirements.txt b/notebooks/transit_spectroscopy_notebook/requirements.txt
index 243d0e235..f44d2de9e 100644
--- a/notebooks/transit_spectroscopy_notebook/requirements.txt
+++ b/notebooks/transit_spectroscopy_notebook/requirements.txt
@@ -1,8 +1,8 @@
-lmfit==1.2.2
-scikit-learn==1.3.0
-numpy==1.25.2
-scipy==1.11.2
-joblib==1.3.2
-matplotlib==3.7.3
-astropy==5.3.3
-pandas==2.1.0
+lmfit>=1.2.2
+scikit-learn>=1.3.0
+numpy>=1.25.2
+scipy>=1.11.2
+joblib>=1.3.2
+matplotlib>=3.7.3
+astropy>=5.3.3
+pandas>=2.1.0