# This fits indivually the cases where we don't get a good results
# `9 March 2024 JEC edited to fit hydrocarbon bands in 3 mu region
#names_miri = ['2MASSJ085747','2MASSJ130152','2MASSJ150958','2MASSJ170756','2MASSJ173628','2MASSJ181129',\
# '2MASSJ182302','2MASSJ203110','2MASSJ203234','2MASSJ203311','2MASSJ203326','2MASSJ204521']
#names_nircam = ['GSC-08152-02121','TYC-8989-436-1','CPD-59-5831','CD-40-11169','TYC-7380-1046-1','TYC-6272-339-1',\
# 'LS-4992','VI-Cyg-1','CPR2002-A38','ALS-15181','GSC-03157-00327','2MASS-J20452110+42235132']
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from matplotlib.ticker import FormatStrFormatter
import astropy.units as u
import astropy.constants as const
from astropy.modeling import models, fitting
from numpy import linspace,exp
from models_mcmc_extension import EmceeFitter
from spectres import spectres
from astropy.modeling import models, fitting
from astropy import units as u
from astropy.convolution import convolve_models
from astropy.modeling.models import Polynomial1D
from astropy.modeling.physical_models import Drude1D
from scipy.interpolate import UnivariateSpline
from astropy.table import Table
from astropy.modeling.models import custom_model
from scipy import special
from scipy.integrate import quad
import math
# Define model
@custom_model
def skewed_gaussian(x, amplitude=0.2, mean=9.8, stddev=2., gamma=0.01):
"""
One dimensional Gaussian model.
Parameters
----------
amplitude : float or `~astropy.units.Quantity`.
Amplitude (peak value) of the Gaussian - for a normalized profile
(integrating to 1), set amplitude = 1 / (stddev * np.sqrt(2 * np.pi))
mean : float or `~astropy.units.Quantity`.
Mean of the Gaussian.
stddev : float or `~astropy.units.Quantity`.
Standard deviation of the Gaussian with FWHM = 2 * stddev * np.sqrt(2 * np.log(2)).
"""
#return amplitude * np.exp(-0.5 * (x - mean) ** 2 / stddev**2)*(1-special.erf((gamma*(x-mean))/(stddev*math.sqrt(2))))
return amplitude * 2/stddev/np.sqrt(2 * np.pi) * np.exp(-0.5 * ((x-mean)/stddev)**2)*0.5 * (1+special.erf(gamma * (x-mean)/stddev/np.sqrt(2)))
def selecting_tyc7380():
direc = "/home/zeegers/wisci_first_shot/data/jwst/nircam/TYC7380_new/"
dataset = "TYC-7380-1046-1_F322W2_fullSED.dat"
wave_nircam_F322W2, flux_nircam_F322W2, uncs_nircam_F322W2, aaa, bbb = np.loadtxt(direc + dataset, unpack=True, skiprows=5)
data_sel=((wave_nircam_F322W2 > 2.43) & (wave_nircam_F322W2 <= 3.955))
wave_nircam_F322W2_sel = wave_nircam_F322W2[data_sel]
flux_nircam_F322W2_sel = flux_nircam_F322W2[data_sel]
uncs_nircam_F322W2_sel = uncs_nircam_F322W2[data_sel]
return(wave_nircam_F322W2_sel, flux_nircam_F322W2_sel, uncs_nircam_F322W2_sel)
def optical_depthfit(continuum_wav, continuum_flux, sel):
linfitter = fitting.LinearLSQFitter()
fit = fitting.SimplexLSQFitter()
poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel])
model_power = models.PowerLaw1D(amplitude = 0.1, x_0=3., alpha=1.6)
gl_init_gaus1=models.Gaussian1D(amplitude=0.1, stddev=0.01, mean=2.7586,bounds={'mean':(2.7586-0.001,2.7586+0.001)})
gl_init_lor1=models.Lorentz1D(amplitude = -0.14, x_0 = 2.7586, fwhm = 0.01,bounds={'x_0':(2.7586-0.001,2.7586+0.001)})
gl_init_lor2=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6257, fwhm = 0.01,bounds={'x_0':(2.6257-0.001,2.6257+0.001)})
gl_init_lor3=models.Lorentz1D(amplitude = -0.14, x_0 = 2.67250, fwhm = 0.01,bounds={'x_0':(2.67250-0.001,2.67250+0.001)})
gl_init_lor4=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6126, fwhm = 0.01,bounds={'x_0':(2.6126-0.001,2.6126+0.001)})
gl_init_lor5=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7400, fwhm = 0.01,bounds={'x_0':(3.7400-0.001,3.7400+0.001)})
#gl_init_lor6=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7035, fwhm = 0.01,bounds={'x_0':(3.7035-0.001,3.7035+0.001)})
power_cont = fit(model_power+gl_init_lor1+gl_init_lor2+gl_init_lor3+gl_init_lor4+gl_init_lor5, continuum_wav[sel], continuum_flux[sel], maxiter=10000)
print(power_cont)
tau_data = np.log(power_cont(continuum_wav)/continuum_flux)
#tau_data = np.log(poly_cont(continuum_wav)/continuum_flux)
plt.plot(continuum_wav, continuum_flux*continuum_wav**2.)
plt.plot(continuum_wav[sel], continuum_flux[sel]*continuum_wav[sel]**2.)
plt.plot(continuum_wav,power_cont(continuum_wav)*continuum_wav**2.)
plt.show()
plt.plot(continuum_wav, continuum_flux)
plt.plot(continuum_wav,power_cont(continuum_wav))
plt.show()
tuple_return2 = [tau_data, power_cont(continuum_wav)]
return(tuple_return2)
def fit_carbonyl_gauss(settings_array1, settings_array2, settings_array3, settings_array4, settings_array5, continuum_wavs, continuum_flux, sel, sel_lines):
fit2 = fitting.SimplexLSQFitter()
fit = fitting.LevMarLSQFitter()
linfitter = fitting.LinearLSQFitter()
poly_cont = linfitter(models.Polynomial1D(3), continuum_wavs[sel], continuum_flux[sel])
model_power = models.PowerLaw1D(amplitude = 0.1, x_0=3., alpha=1.6)
gl_init_gaus1=models.Gaussian1D(amplitude=0.1, stddev=0.01, mean=2.7586,bounds={'mean':(2.7586-0.001,2.7586+0.001)})
gl_init_lor1=models.Lorentz1D(amplitude = -0.14, x_0 = 2.7586, fwhm = 0.01,bounds={'x_0':(2.7586-0.001,2.7586+0.001)})
gl_init_lor2=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6257, fwhm = 0.01,bounds={'x_0':(2.6257-0.001,2.6257+0.001)})
gl_init_lor3=models.Lorentz1D(amplitude = -0.14, x_0 = 2.67250, fwhm = 0.01,bounds={'x_0':(2.67250-0.001,2.67250+0.001)})
gl_init_lor4=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6126, fwhm = 0.01,bounds={'x_0':(2.6126-0.001,2.6126+0.001)})
gl_init_lor5=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7400, fwhm = 0.01,bounds={'x_0':(3.7400-0.001,3.7400+0.001)})
#gl_init_lor6=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7035, fwhm = 0.01,bounds={'x_0':(3.7035-0.001,3.7035+0.001)})
power_cont = fit2(model_power+gl_init_lor1+gl_init_lor2+gl_init_lor3+gl_init_lor4+gl_init_lor5, continuum_wav[sel], continuum_flux[sel], maxiter=10000)
tau_data = np.log(power_cont(continuum_wavs)/continuum_flux)
stddevs = np.std(tau_data[sel])
weights = 1.0 / (stddevs)
# load the trapped water ice 200 K model
a,b = np.loadtxt('/home/zeegers/wisci_first_shot/models/ice_models/MgSiO3+H2O_7.5_200K.txt', unpack=True)
wav_model_ice = 10000./a
tau_model_ice = -1.*np.log(b)
gl_init_skewgaus=skewed_gaussian(amplitude=0.02, stddev=0.5, mean=2.9,gamma=0.2)
line = models.Linear1D(slope=0.01, intercept=0.0001)
model_sel=((wav_model_ice >= 2.5) & (wav_model_ice <= 3.9))
ice_fit2 = fit(line + gl_init_skewgaus, wav_model_ice[model_sel], tau_model_ice[model_sel], maxiter=10000)
icemodel_line = models.Linear1D(slope=ice_fit2.slope_0.value, intercept=ice_fit2.intercept_0.value, fixed={'slope': True, 'intercept': True})
icemodel_skewedgaus = skewed_gaussian(amplitude=ice_fit2.amplitude_1.value, stddev=ice_fit2.stddev_1.value, mean=ice_fit2.mean_1.value,gamma=ice_fit2.gamma_1.value, fixed={ 'stddev': True, 'mean': True, 'gamma': True})
#gl_init_gaus=models.Gaussian1D(amplitude=0.1, stddev=0.05, mean=5.81,bounds={'mean':(5.81-0.02,5.81+0.02), 'stddev':(1e-4,0.1)})
#gl_init_gaus1=models.Gaussian1D(amplitude=settings_array1[0], stddev=settings_array1[1], mean=settings_array1[2],bounds={'mean':(settings_array1[2]-settings_array1[3],settings_array1[2]+settings_array1[4]), 'stddev':(settings_array1[5],settings_array1[6])})
#gl_init_gaus2=models.Gaussian1D(amplitude=settings_array2[0], stddev=settings_array2[1], mean=settings_array2[2],bounds={'mean':(settings_array2[2]-settings_array2[3],settings_array2[2]+settings_array2[4]), 'stddev':(settings_array2[5],settings_array2[6])})
#gl_init_gaus3=models.Gaussian1D(amplitude=settings_array3[0], stddev=settings_array3[1], mean=settings_array3[2],bounds={'mean':(settings_array3[2]-settings_array3[3],settings_array3[2]+settings_array3[4]), 'stddev':(settings_array3[5],settings_array3[6])})
#gl_init_gaus4=models.Gaussian1D(amplitude=settings_array4[0], stddev=settings_array4[1], mean=settings_array4[2],bounds={'mean':(settings_array4[2]-settings_array4[3],settings_array4[2]+settings_array4[4]), 'stddev':(settings_array4[5],settings_array4[6])})
#gl_init_gaus5=models.Gaussian1D(amplitude=settings_array5[0], stddev=settings_array5[1], mean=settings_array5[2],bounds={'mean':(settings_array5[2]-settings_array5[3],settings_array5[2]+settings_array5[4]), 'stddev':(settings_array5[5],settings_array5[6])})
gl_init_gaus1=models.Gaussian1D(amplitude=settings_array1[0], stddev=settings_array1[1], mean=settings_array1[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus2=models.Gaussian1D(amplitude=settings_array2[0], stddev=settings_array2[1], mean=settings_array2[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus3=models.Gaussian1D(amplitude=settings_array3[0], stddev=settings_array3[1], mean=settings_array3[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus4=models.Gaussian1D(amplitude=settings_array4[0], stddev=settings_array4[1], mean=settings_array4[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus5=models.Gaussian1D(amplitude=settings_array5[0], stddev=settings_array5[1], mean=settings_array5[2], fixed={'stddev': True, 'mean': True})
#stellar_line_init = models.Lorentz1D(amplitude = 0.10, x_0 = 3.2966, fwhm = 0.005,bounds={'x_0':(3.2966-0.001,3.2966+0.002)}, fixed={'fwhm': True})
# 2.8728 Pf 11 3.03919 Pf eps
stellar_line_init1 = models.Gaussian1D(amplitude=0.15, stddev=0.003, mean=2.8728, bounds={'mean':(2.8728-0.001,2.8728+0.002),'stddev':(0.001,0.01)})
stellar_line_init2 = models.Gaussian1D(amplitude=0.15, stddev=0.003, mean=3.03919, bounds={'mean':(3.03919-0.001,3.03919+0.002),'stddev':(0.001,0.01)})
stellar_line_init3 = models.Gaussian1D(amplitude=0.15, stddev=0.003, mean=3.2966, bounds={'mean':(3.2966-0.001,3.2966+0.002),'stddev':(0.001,0.01)})
#, fixed={'stddev': True, 'mean': True}
print(gl_init_gaus1,gl_init_gaus2,gl_init_gaus3, gl_init_gaus4, gl_init_gaus5)
print(continuum_wavs, tau_data)
gl_fit = fit(gl_init_gaus1+gl_init_gaus2+gl_init_gaus3+gl_init_gaus4+gl_init_gaus5+stellar_line_init1+stellar_line_init2+stellar_line_init3+icemodel_line + icemodel_skewedgaus, continuum_wavs[sel_lines], tau_data[sel_lines], weights=weights, maxiter=10000) # this is where something goes wrong
return gl_fit
def fit_carbonyl_drude(settings_array, continuum_wav, continuum_flux, sel, sel_lines):
fit2 = fitting.SimplexLSQFitter()
fit = fitting.LevMarLSQFitter()
linfitter = fitting.LinearLSQFitter()
gl_init_gaus1=models.Gaussian1D(amplitude=0.1, stddev=0.01, mean=2.7586,bounds={'mean':(2.7586-0.001,2.7586+0.001)})
gl_init_lor1=models.Lorentz1D(amplitude = -0.14, x_0 = 2.7586, fwhm = 0.01,bounds={'x_0':(2.7586-0.001,2.7586+0.001)})
gl_init_lor2=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6257, fwhm = 0.01,bounds={'x_0':(2.6257-0.001,2.6257+0.001)})
gl_init_lor3=models.Lorentz1D(amplitude = -0.14, x_0 = 2.67250, fwhm = 0.01,bounds={'x_0':(2.67250-0.001,2.67250+0.001)})
gl_init_lor4=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6126, fwhm = 0.01,bounds={'x_0':(2.6126-0.001,2.6126+0.001)})
gl_init_lor5=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7400, fwhm = 0.01,bounds={'x_0':(3.7400-0.001,3.7400+0.001)})
#gl_init_lor6=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7035, fwhm = 0.01,bounds={'x_0':(3.7035-0.001,3.7035+0.001)})
poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel])
model_power = models.PowerLaw1D(amplitude = 0.1, x_0=3., alpha=1.6)
power_cont = fit2(model_power+gl_init_lor1+gl_init_lor2+gl_init_lor3+gl_init_lor4+gl_init_lor5, continuum_wav[sel], continuum_flux[sel], maxiter=10000)
tau_data = np.log(power_cont(continuum_wav)/continuum_flux)
stddev = np.std(tau_data[sel])
weights = 1.0 / (stddev)
fwhm = 2 * settings_array[1] * np.sqrt(2 * np.log(2))
fwhm_lower = 2 * settings_array[5] * np.sqrt(2 * np.log(2))
print('fwhm_lower',fwhm_lower)
fwhm_upper = 2 * settings_array[6] * np.sqrt(2 * np.log(2))
gl_init_gaus=models.Drude1D(amplitude=settings_array[0], fwhm=fwhm, x_0=settings_array[2],bounds={'x_0':(settings_array[2]-settings_array[3],settings_array[2]-settings_array[4]), 'fwhm':(fwhm_lower,fwhm_upper)})
gl_fit = fit(gl_init_gaus, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights, maxiter=10000)
return gl_fit
def fit_carbonyl_errors_gauss(name, settings_array1, settings_array2, settings_array3, settings_array4, settings_array5, continuum_wav, continuum_flux, sel, sel_lines):
#define water ice model
a,b = np.loadtxt('/home/zeegers/wisci_first_shot/models/ice_models/MgSiO3+H2O_7.5_200K.txt', unpack=True)
wav_model_ice = 10000./a
tau_model_ice = -1.*np.log(b)
gl_init_skewgaus=skewed_gaussian(amplitude=0.02, stddev=0.5, mean=2.9,gamma=0.2)
model_sel=((wav_model_ice >= 2.5) & (wav_model_ice <= 3.9))
gl_fit2 = fit(line + gl_init_skewgaus, wav_model_ice[model_sel], tau_model_ice[model_sel], maxiter=10000)
# calculates errors
# and fitting with Emcee
print("settings_array", settings_array1[0],settings_array1[1],settings_array1[2])
fit2 = fitting.SimplexLSQFitter()
fit = fitting.LevMarLSQFitter()
linfitter = fitting.LinearLSQFitter()
gl_init_gaus1=models.Gaussian1D(amplitude=0.1, stddev=0.01, mean=2.7586,bounds={'mean':(2.7586-0.001,2.7586+0.001)})
gl_init_lor1=models.Lorentz1D(amplitude = -0.14, x_0 = 2.7586, fwhm = 0.01,bounds={'x_0':(2.7586-0.001,2.7586+0.001)})
gl_init_lor2=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6257, fwhm = 0.01,bounds={'x_0':(2.6257-0.001,2.6257+0.001)})
gl_init_lor3=models.Lorentz1D(amplitude = -0.14, x_0 = 2.67250, fwhm = 0.01,bounds={'x_0':(2.67250-0.001,2.67250+0.001)})
gl_init_lor4=models.Lorentz1D(amplitude = -0.14, x_0 = 2.6126, fwhm = 0.01,bounds={'x_0':(2.6126-0.001,2.6126+0.001)})
gl_init_lor5=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7400, fwhm = 0.01,bounds={'x_0':(3.7400-0.001,3.7400+0.001)})
#gl_init_lor6=models.Lorentz1D(amplitude = -0.01, x_0 = 3.7035, fwhm = 0.01,bounds={'x_0':(3.7035-0.001,3.7035+0.001)})
poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel])
model_power = models.PowerLaw1D(amplitude = 0.1, x_0=3., alpha=1.6)
power_cont = fit2(model_power+gl_init_lor1+gl_init_lor2+gl_init_lor3+gl_init_lor4+gl_init_lor5, continuum_wav[sel], continuum_flux[sel], maxiter=10000)
tau_data = np.log(power_cont(continuum_wav)/continuum_flux)
stddevs = np.std(tau_data[sel])
weights = 1.0 / (stddevs)
#gl_init_gaus=models.Gaussian1D(amplitude=0.1, stddev=0.05, mean=5.81,bounds={'mean':(5.81-0.03,5.81+0.03), 'stddev':(1e-4,0.15)})
#gl_init_gaus1=models.Gaussian1D(amplitude=settings_array1[0], stddev=settings_array1[1], mean=settings_array1[2],bounds={'mean':(settings_array1[2]-settings_array1[3],settings_array1[2]+settings_array1[4]), 'stddev':(settings_array1[5],settings_array1[6])})
#gl_init_gaus2=models.Gaussian1D(amplitude=settings_array2[0], stddev=settings_array2[1], mean=settings_array2[2],bounds={'mean':(settings_array2[2]-settings_array2[3],settings_array2[2]+settings_array2[4]), 'stddev':(settings_array2[5],settings_array2[6])})
#gl_init_gaus3=models.Gaussian1D(amplitude=settings_array3[0], stddev=settings_array3[1], mean=settings_array3[2],bounds={'mean':(settings_array3[2]-settings_array3[3],settings_array3[2]+settings_array3[4]), 'stddev':(settings_array3[5],settings_array3[6])})
#gl_init_gaus4=models.Gaussian1D(amplitude=settings_array4[0], stddev=settings_array4[1], mean=settings_array4[2],bounds={'mean':(settings_array4[2]-settings_array4[3],settings_array4[2]+settings_array4[4]), 'stddev':(settings_array4[5],settings_array4[6])})
#gl_init_gaus5=models.Gaussian1D(amplitude=settings_array5[0], stddev=settings_array5[1], mean=settings_array5[2],bounds={'mean':(settings_array5[2]-settings_array5[3],settings_array5[2]+settings_array5[4]), 'stddev':(settings_array5[5],settings_array5[6])})
gl_init_gaus1=models.Gaussian1D(amplitude=settings_array1[0], stddev=settings_array1[1], mean=settings_array1[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus2=models.Gaussian1D(amplitude=settings_array2[0], stddev=settings_array2[1], mean=settings_array2[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus3=models.Gaussian1D(amplitude=settings_array3[0], stddev=settings_array3[1], mean=settings_array3[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus4=models.Gaussian1D(amplitude=settings_array4[0], stddev=settings_array4[1], mean=settings_array4[2], fixed={'stddev': True, 'mean': True})
gl_init_gaus5=models.Gaussian1D(amplitude=settings_array5[0], stddev=settings_array5[1], mean=settings_array5[2], fixed={'stddev': True, 'mean': True})
#stellar_line_init = models.Lorentz1D(amplitude = 0.08, x_0 = 3.2966, fwhm = 0.01,bounds={'x_0':(3.2966-0.001,3.2966+0.001)}, fixed={'fwhm': True})
stellar_line_init = models.Gaussian1D(amplitude=0.15, stddev=0.003, mean=3.2966, bounds={'mean':(3.2966-0.001,3.2966+0.002),'stddev':(0.001,0.01)})
gl_fit = fit(gl_init_gaus1+gl_init_gaus2+gl_init_gaus3+gl_init_gaus4+gl_init_gaus5+stellar_line_init, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights, maxiter=100000)
# plotting the initial models
init_model = gl_init_gaus1+gl_init_gaus2+gl_init_gaus3+gl_init_gaus4+gl_init_gaus5
plt.plot(continuum_wav, tau_data)
plt.errorbar(continuum_wav,tau_data, yerr=continuum_error, linestyle="-",marker='')
plt.plot(continuum_wav, init_model(continuum_wav))
plt.ylabel('optical depth mod', fontsize=14)
plt.xlabel(r' Wavelength [$\mu$m]',fontsize=18)
plt.title('initial model')
plt.show()
plt.plot(continuum_wav, tau_data)
plt.errorbar(continuum_wav,tau_data, yerr=continuum_error, linestyle="-",marker='')
plt.plot(continuum_wav, gl_fit(continuum_wav))
plt.ylabel('optical depth mod', fontsize=14)
plt.xlabel(r' Wavelength [$\mu$m]',fontsize=18)
plt.title('first fit model')
plt.show()
fit2 = EmceeFitter(nsteps=1000, burnfrac=0.1) #, save_samples=emcee_samples_file
fit_mcmc_result = fit2(gl_fit, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights)
# fit2.plot_emcee_results(fit_mcmc_result, filebase="/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_rem_github/carbonyl/carbonyl_indv_fits/plots/emcee_res")
fit2.plot_emcee_results(fit_mcmc_result, filebase="/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/water_ice/plots/emcee_res")
plt.show()
print(fit_mcmc_result.parameters)
print(fit_mcmc_result.uncs)
chains = fit2.fit_info['sampler'].get_chain(flat=True,discard=np.int32(0.1*1000))
tuple_return = [fit_mcmc_result, chains,stddevs]
print('tuple',tuple_return)
return tuple_return
#return fitparams, self.fit_info
def fit_carbonyl_errors_drude(name, settings_array, continuum_wav, continuum_flux, sel, sel_lines):
# calculates errors
# and fitting with Emcee for the Drude profile
fit2 = fitting.SimplexLSQFitter()
fit = fitting.LevMarLSQFitter()
linfitter = fitting.LinearLSQFitter()
model_power = models.PowerLaw1D(amplitude = 0.1, x_0=3., alpha=1.6)
power_cont = fit2(model_power, continuum_wav[sel], continuum_flux[sel], maxiter=10000)
poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel])
tau_data = np.log(power_cont(continuum_wav)/continuum_flux)
stddevs = np.std(tau_data[sel])
weights = 1.0 / (stddevs)
gl_init_gaus=models.Drude1D(amplitude=settings_array[0], fwhm=fwhm, x_0=settings_array[2],bounds={'x_0':(settings_array[2]-settings_array[3],settings_array[2]-settings_array[4]), 'fwhm':(fwhm_lower,fwhm_upper)})
gl_fit = fit(gl_init_gaus, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights, maxiter=100000)
fit2 = EmceeFitter(nsteps=100, burnfrac=0.1) #, save_samples=emcee_samples_file
fit_mcmc_result = fit2(gl_fit, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights)
# Get the chains:
print(fit2.fit_info['sampler'].get_chain(flat=True,discard=np.int32(0.1*100)))
fit2.plot_emcee_results(fit_mcmc_result, filebase="/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/water_ice/emcee_res")
plt.show()
print(fit_mcmc_result.parameters)
print(fit_mcmc_result.uncs)
chains = fit2.fit_info['sampler'].get_chain(flat=True,discard=np.int32(0.1*100))
tuple_return = tuple(fit_mcmc_result, chains)
print('tuple',tuple_return)
return tuple_return
#return fitparams, self.fit_inf
def plot_fits(filename,continuum_wav,tau_data,continuum_error):
# plots the fits and writes it to a file
fig,axs=plt.subplots(2,1)
axs[0].plot(continuum_wav, tau_data)
axs[0].errorbar(continuum_wav,tau_data, yerr=continuum_error, linestyle="-",marker='')
axs[0].plot(continuum_wav, tau_model)
axs[0].set_ylabel('optical depth', fontsize=14)
axs[0].set_xlabel(r' Wavelength [$\mu$m]',fontsize=18)
axs[1].plot(continuum_wav, tau_data-tau_model)
axs[1].set_ylabel('residuals', fontsize=14)
axs[1].set_xlabel(r'Wavelength ($\mu$m)', fontsize=14)
plt.subplots_adjust(hspace=0,wspace=0.17)
plt.show()
def print_results(filename, source_name, Av, wav_mean, fwhm, fwhm_error, surface_area, surface_area_error):
# in the table we need: name of the source, mean wavelength, FWHM + error, surface area in per cm^-1
# create both a table with a text file and with a tex format for the paper
# create empty tables
# table txt
table_path = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/water_ice/water_ice_res/'
table_txt = Table(
names=(
"Source",
"AV",
"wavelength[micron]",
"FWHM[micron]",
"FWHM error",
"Integrated Area",
"Integrated Area error"
),
dtype=(
"str",
"float64",
"float64",
"float64",
"float64",
"float64",
"float64"
),
)
# table latex
table_latex = Table(
names=(
"Source",
r"$A_{V}$"
r"\lambda (\micron)",
r"$\Delta \lambda$",
"Amplitude",
"Integrated Area"
),
dtype=("str", "str", "str", "str", "str"),
)
# add the fitting results to the tables
for i in range(0,len(source_name)):
table_txt.add_row(
(
source_name[i],
Av[i],
wav_mean[i],
fwhm[i],
fwhm_error[i],
surface_area[i],
surface_area_error[i]
)
)
table_latex.add_row(
(
source_name[i],
f"${Av[i]:.3f}$",
f"${wav_mean[i]:.3f}$",
f"${fwhm[i]:.3f}\pm{fwhm_error[i]:.3f}$",
f"${surface_area[i]:.3f}\pm {surface_area_error[i]:.3f}$"
)
)
tabname = filename
# write the tables to files
table_txt.write(
table_path + f"{tabname}.txt",
format="ascii.commented_header",
overwrite=True,
)
table_latex.write(
table_path + f"{tabname}.tex",
format="aastex",
col_align="lccc",
latexdict={
"caption": r"Fitting results. \label{tab:fit_results}",
},
overwrite=True,
)
# to do list:
# Get 10 Lac on the list of sources and fit it
# put names of sources on the plots
#
# 10lac_nircam_mrs_merged.fits
# names sources
names_miri = ['2MASSJ085747','2MASSJ150958']
names_miri_2mass = ['2MASSJ08574757-4609145','2MASSJ15095841-5958463']
# main program
direc = '/home/zeegers/wisci_first_shot/'
datafile = '_nircam_mrs_merged.fits'
# creating a table for the results
output_directory = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/water_ice/water_ice_res/'
filename_results = 'output_res_hydrocarbon.txt'
# output directory for plots
output_directory_plots = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/water_ice/plots/'
output_directory_plots_emcee = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/water_ice/plots/emcee_res'
Av_array = np.array([ 5.1, 4.7]) # old estimates
Av_array_new = np.array([4.98, 4.35])
Rv_array_new = ([3.32, 3.13])
fwhm1 = 0.05
fwhm2 = 0.05
fwhm3 = 0.05
fwhm4 = 0.05
fwhm5 = 0.09
stddev1 = fwhm1/(2.0*math.sqrt(2 * math.log(2)))
stddev2 = fwhm2/(2.0*math.sqrt(2 * math.log(2)))
stddev3 = fwhm3/(2.0*math.sqrt(2 * math.log(2)))
stddev4 = fwhm4/(2.0*math.sqrt(2 * math.log(2)))
stddev5 = fwhm5/(2.0*math.sqrt(2 * math.log(2)))
# fix standard deviation and the mean
# settings array gives the amplitude, stddev, mean, mean_bounds and stddev bounds
#gl_init_gaus1=models.Gaussian1D(amplitude=-0.0005, stddev=stddev1, mean=3.376, fixed={'stddev': True, 'mean': True})
#gl_init_gaus2=models.Gaussian1D(amplitude=-0.001, stddev=stddev2, mean=3.420, fixed={'stddev': True, 'mean': True})
#gl_init_gaus3=models.Gaussian1D(amplitude=-0.001, stddev=stddev3, mean=3.474, fixed={'stddev': True, 'mean': True})
#gl_init_gaus4=models.Gaussian1D(amplitude=-0.0005, stddev=stddev4, mean=3.520, fixed={'stddev': True, 'mean': True})
#gl_init_gaus5=models.Gaussian1D(amplitude=-0.001, stddev=stddev5, mean=3.289, fixed={'stddev': True, 'mean': True})
settings_array_gauss1 = np.array([[0.01, stddev1, 3.376,0.03,0.03,1e-4,0.2],
[0.01,stddev1,3.376,0.03,0.03,1e-4,0.2]
])
settings_array_gauss2 = np.array([[0.01,stddev2,3.420,0.03,0.03,1e-4,0.2],
[0.01,stddev2,3.420,0.03,0.03,1e-4,0.2]
])
settings_array_gauss3 = np.array([[0.01, stddev3, 3.474,0.03,0.03,1e-4,0.2],
[0.01,stddev3,3.474,0.03,0.03,1e-4,0.2]
])
settings_array_gauss4 = np.array([[0.01, stddev4, 3.520,0.03,0.03,1e-4,0.2],
[0.01,stddev4,3.520,0.03,0.03,1e-4,0.2]
])
settings_array_gauss5 = np.array([[0.01, stddev5, 3.289,0.03,0.03,1e-4,0.2],
[0.01,stddev5,3.289,0.03,0.03,1e-4,0.2]
])
# output arrays
amplitude_array = np.zeros((2,6))
amplitude_array_error = np.zeros((2,6))
amplitude_array_unc_plus = np.zeros((2,6))
amplitude_array_unc_minus = np.zeros((2,6))
mean_array = np.zeros((2,6))
mean_array_error = np.zeros((2,6))
mean_array_unc_plus = np.zeros((2,6))
mean_array_unc_minus = np.zeros((2,6))
stddev_array = np.zeros((2,6))
stddev_array_error = np.zeros((2,6))
stddev_array_unc_plus = np.zeros((2,6))
stddev_array_unc_minus = np.zeros((2,6))
fwhm_array = np.zeros((2,6))
fwhm_array_error = np.zeros((2,6))
fwhm_array_unc_plus = np.zeros((2,6))
fwhm_array_unc_minus = np.zeros((2,6))
surface_area_cm = np.zeros((2,6))
surface_area_error = np.zeros((2,6))
surface_area_unc_plus = np.zeros((2,6))
surface_area_unc_minus = np.zeros((2,6))
surface_area_cm_total = np.zeros((2))
surface_area_error_total = np.zeros((2))
surface_area_unc_plus_total = np.zeros((2))
surface_area_unc_minus_total = np.zeros((2))
surface_wavenm_total = np.empty(shape=[2,576000])
max_tau_values = np.zeros((2))
settings_array_drude = np.array([[0.2, 0.1,5.81,0.001,0.001,0.01,0.3],
[0.2, 0.1,5.81,0.001,0.001,0.01,0.3]
])
for i in range(0,len(names_miri)):
data_miri_merged = fits.getdata(direc + 'data/jwst/delivery_v6/' + names_miri[i]+datafile)
# changing this to an if else statement to select the dataset of TYC 7380 from KM new data reduction
# i = 4 or 2MASSJ173628
if i == 4:
wavelength_merged_new, flux_merged_new, uncs_merged_new = selecting_tyc7380()
else:
wavelength_merged = data_miri_merged['WAVELENGTH']
flux_merged = data_miri_merged['FLUX']
uncs_merged = data_miri_merged['UNC']
#Get rid of NaN values from spectrum
good = np.where(np.isfinite(flux_merged)&(flux_merged > 0.))
flux_merged_new = flux_merged[good]
wavelength_merged_new = wavelength_merged[good]
uncs_merged_new = uncs_merged[good]
continuum_sel=((wavelength_merged_new > 2.6) & (wavelength_merged_new <= 3.9)) # aliphatic hydrocarbon NIRCAM
continuum_wav = wavelength_merged_new[continuum_sel]
continuum_flux = flux_merged_new[continuum_sel]
# We need to get the signal to noise from the ETC calculation, probably the best we've got so far?
continuum_error = uncs_merged_new[continuum_sel]
continuum_error_large = np.argwhere(continuum_error>0.5)
continuum_error[continuum_error_large] = 0.
# Check if in some cases you need to improve the line selection
# We get strong lines for ...
# Creating a line list from the stellar models:
# weak line at 2.6126 (Pf_14)
# line at 2.6257 (Br_beta)
# line at 2.67250
# line at 2.7586
# lines at 2.8731 3.0391 3.2696 3.2966 3.7035 3.7400
# If we want to fit the waterice feature we need to include the lines inside the waterice feature in the fit
# The model will be Alexey's ice model + gaussians + carbon
# Alexey's model will need a scaling parameter A to fit it to the feature, beause otherwise there's not much to fit
# Think about how to rewrite
feature_carbonyl = ((continuum_wav > 2.8) & (continuum_wav <= 3.65))
#stellar_line1 = ((continuum_wav > 2.6163) & (continuum_wav <= 2.6323))
#stellar_line2 = ((continuum_wav > 2.6687) & (continuum_wav <= 2.6836))
#stellar_line3 = ((continuum_wav > 2.7494) & (continuum_wav <= 2.7696))
#stellar_line4 = ((continuum_wav > 3.0271) & (continuum_wav <= 3.0494))
#stellar_line5 = ((continuum_wav > 3.2897) & (continuum_wav <= 3.3049))
#stellar_line6 = ((continuum_wav > 3.2715) & (continuum_wav <= 3.3071))
stellar_line7 = ((continuum_wav > 3.7001) & (continuum_wav <= 3.7107))
stellar_line8 = ((continuum_wav > 3.72) & (continuum_wav <= 3.76))
stellar_line1_extreme = ((continuum_wav > 3.058) & (continuum_wav <= 3.11)) # 3.058 – 3.11
sel =~ (feature_carbonyl|stellar_line7) # can be expanded if more features will be added
sel_lines =~ (stellar_line7) # can be expanded if more features will be added
optical_depth_return = optical_depthfit(continuum_wav, continuum_flux, sel)
tau_data = optical_depth_return[0]
power_cont_array = optical_depth_return[1]
gl_fit_gauss = fit_carbonyl_gauss(settings_array_gauss1[i], settings_array_gauss2[i], settings_array_gauss3[i],settings_array_gauss4[i],settings_array_gauss5[i], continuum_wav, continuum_flux, sel, sel_lines) # something wrong here?
#gl_fit_drude = fit_carbonyl_drude(settings_array_drude[i], continuum_wav, continuum_flux, sel, sel_lines)
#tau_model = gl_fit_drude(continuum_wav)
tau_model = gl_fit_gauss(continuum_wav)
#print(gl_fit_drude)
print(gl_fit_gauss)
fig,axs=plt.subplots(2,1)
axs[0].plot(continuum_wav, tau_data)
axs[0].errorbar(continuum_wav,tau_data, yerr=continuum_error, linestyle="-",marker='')
axs[0].plot(continuum_wav, tau_model, zorder =20)
axs[0].set_ylabel('optical depth', fontsize=14)
axs[0].set_xlabel(r' Wavelength [$\mu$m]',fontsize=18)
axs[1].plot(continuum_wav, tau_data-tau_model)
axs[1].set_ylabel('residuals', fontsize=14)
axs[1].set_xlabel(r'Wavelength ($\mu$m)', fontsize=14)
plt.subplots_adjust(hspace=0,wspace=0.17)
plt.show()
fit_result = fit_carbonyl_errors_gauss(names_miri[i], settings_array_gauss1[i], settings_array_gauss2[i], settings_array_gauss3[i],settings_array_gauss4[i],settings_array_gauss5[i], continuum_wav, continuum_flux, sel, sel_lines)
print(fit_result)
fit_mcmc_result = fit_result[0]
tau_model2 = fit_mcmc_result(continuum_wav)
stddev = fit_result[2]
fig,axs=plt.subplots(3,1, figsize=(6, 10))
axs[0].plot(continuum_wav, continuum_flux*continuum_wav**2., color = "black")
axs[0].plot(continuum_wav[sel], continuum_flux[sel]*continuum_wav[sel]**2., marker='o',linestyle=" ", markersize = 1)
axs[0].plot(continuum_wav,power_cont_array*continuum_wav**2., color="green")
axs[0].set_ylabel(r'$\lambda^2$*F($\lambda$) [$\mu\mathrm{m}^2$ Jy]', fontsize=10)
axs[1].errorbar(continuum_wav, tau_data, yerr=stddev, marker='',linestyle="-", color="grey")#, errorevery=60)
axs[1].plot(continuum_wav, tau_data, marker='',linestyle="-", color = "black", zorder=10)
axs[1].plot(continuum_wav, tau_model2, color="orange", lw=2, zorder=100)
axs[1].axhline(y=0, color='k', ls=":", c="k")
axs[1].set_ylabel('optical depth', fontsize=11)
axs[1].set_xlabel(r' Wavelength [$\mu$m]',fontsize=18)
axs[2].plot(continuum_wav, (tau_data-tau_model2)/stddev, color="grey")
axs[2].set_ylabel('residuals', fontsize=14)
axs[2].set_xlabel(r'Wavelength ($\mu$m)', fontsize=11)
axs[2].axhline(y=0, color='k', ls=":", c="k")
plt.subplots_adjust(hspace=0,wspace=0.17)
plt.show()
fig.savefig(output_directory_plots+names_miri[i]+"plot_res_carbonyl.pdf",dpi=200)
plt.close()
fig2,ax=plt.subplots()
ax.plot(continuum_wav, tau_data)
ax.errorbar(continuum_wav, tau_data, yerr=continuum_error, marker='',linestyle="-", color="black")
ax.plot(continuum_wav, tau_model2, color="crimson", lw=2, zorder=10)
ax.axhline(y=0, color='k', ls=":", c="k")
ax.set_ylabel('optical depth', fontsize=14)
ax.set_xlabel(r'Wavelength ($\mu$m)', fontsize=14)
plt.show()
fig2.savefig(output_directory_plots+names_miri[i]+"plot_noresiduals_carbon34.pdf",dpi=200, format="pdf")
plt.close()
aaa=(np.where((continuum_wav>3.3) & (continuum_wav<3.5)))
near_34 = tau_model2[aaa]
max_tau_values[i] = np.max(near_34)
######
# Calculating all values and uncertainties for the five lines
######
amp_unc_ave = ([fit_mcmc_result.amplitude_0.unc,fit_mcmc_result.amplitude_1.unc,fit_mcmc_result.amplitude_2.unc,fit_mcmc_result.amplitude_3.unc,fit_mcmc_result.amplitude_4.unc])
stddev_unc_ave =([fit_mcmc_result.stddev_0.unc, fit_mcmc_result.stddev_1.unc,fit_mcmc_result.stddev_2.unc,fit_mcmc_result.stddev_3.unc,fit_mcmc_result.stddev_4.unc])
mean_unc_ave = ([fit_mcmc_result.mean_0.unc, fit_mcmc_result.mean_1.unc,fit_mcmc_result.mean_2.unc,fit_mcmc_result.mean_3.unc,fit_mcmc_result.mean_4.unc])
amplitudes = ([fit_mcmc_result.amplitude_0.value, fit_mcmc_result.amplitude_1.value, fit_mcmc_result.amplitude_2.value,fit_mcmc_result.amplitude_3.value,fit_mcmc_result.amplitude_4.value])
means = ([fit_mcmc_result.mean_0.value,fit_mcmc_result.mean_1.value,fit_mcmc_result.mean_2.value,fit_mcmc_result.mean_3.value,fit_mcmc_result.mean_4.value])
stddevs = ([fit_mcmc_result.stddev_0.value,fit_mcmc_result.stddev_1.value,fit_mcmc_result.stddev_2.value,fit_mcmc_result.stddev_3.value,fit_mcmc_result.stddev_4.value])
# we need to save the results for each model and for the total
# below we store the results for each of the 5 models
# needs to be rewritten for the case where we don't evaluate the mean and the standard deviation, but only the amplitude
for j in range(0,5):
amplitude_array[i,j]=amplitudes[j]
amplitude_array_error[i,j]=amp_unc_ave[j]
mean_array[i,j]=means[j]
mean_array_error[i,j]=mean_unc_ave[j]
stddev_array[i,j] = stddevs[j]
stddev_array_error[i,j]=stddev_unc_ave[j]
# calculate uncertainties using chains, this still needs to be adapted for the for loop
chains = fit_result[1]
# FWHM and uncertainties
fwhm = 2. * stddevs[j] * np.sqrt(2. * np.log(2.))
fwhm_unc_array = np.percentile(fwhm, [16, 50, 84])
fwhm_array[i,j] = fwhm_unc_array[1]
fwhm_array_error[i,j] = 0.5 * (fwhm_unc_array[2] - fwhm_unc_array[0])
fwhm_array_unc_plus[i,j] = fwhm_unc_array[2] - fwhm_unc_array[1]
fwhm_array_unc_minus[i,j] = fwhm_unc_array[1] - fwhm_unc_array[0]
stddev_wavenm = 1./(means[j]**2)*stddevs[j]*10000.
surface_wavenm = chains[:,0+j]*stddev_wavenm * np.sqrt(2. * np.pi)
surface_wavenm_unc_array = np.percentile(surface_wavenm, [16, 50, 84])
surface_area_cm[i,j] = surface_wavenm_unc_array[1]
surface_area_error[i,j] = 0.5 * (surface_wavenm_unc_array[2] - surface_wavenm_unc_array[0])
surface_area_unc_plus[i,j] = surface_wavenm_unc_array[2] - surface_wavenm_unc_array[1]
surface_area_unc_minus[i,j] = surface_wavenm_unc_array[1] - surface_wavenm_unc_array[0]
# adding chains for total surface area
surface_wavenm_total += surface_wavenm
fwhm_unc_array = 0
stddev_wavenm = 0
surface_wavenm = 0
surface_wavenm_unc_array = 0
chains = 0
fwhm = 0
# below the forloop which is used when all the parameters are free or fitted with bounds, but not fixed
#for j in range(0,5):
#amplitude_array[i,j]=amplitudes[j]
#amplitude_array_error[i,j]=amp_unc_ave[j]
#mean_array[i,j]=means[j]
#mean_array_error[i,j]=mean_unc_ave[j]
#stddev_array[i,j] = stddevs[j]
#stddev_array_error[i,j]=stddev_unc_ave[j]
## calculate uncertainties using chains, this still needs to be adapted for the for loop
#chains = fit_result[1]
## FWHM and uncertainties
#fwhm = 2. * chains[:,2+j*3] * np.sqrt(2. * np.log(2.))
#fwhm_unc_array = np.percentile(fwhm, [16, 50, 84])
#fwhm_array[i,j] = fwhm_unc_array[1]
#fwhm_array_error[i,j] = 0.5 * (fwhm_unc_array[2] - fwhm_unc_array[0])
#fwhm_array_unc_plus[i,j] = fwhm_unc_array[2] - fwhm_unc_array[1]
#fwhm_array_unc_minus[i,j] = fwhm_unc_array[1] - fwhm_unc_array[0]
#stddeviation = chains[:,2+j*3]
#stddev_wavenm = 1./(chains[:,1+j*3]**2)*stddeviation*10000.
#surface_wavenm = chains[:,0+j*3]*stddev_wavenm * np.sqrt(2. * np.pi)
#surface_wavenm_unc_array = np.percentile(surface_wavenm, [16, 50, 84])
#surface_area_cm[i,j] = surface_wavenm_unc_array[1]
#surface_area_error[i,j] = 0.5 * (surface_wavenm_unc_array[2] - surface_wavenm_unc_array[0])
#surface_area_unc_plus[i,j] = surface_wavenm_unc_array[2] - surface_wavenm_unc_array[1]
#surface_area_unc_minus[i,j] = surface_wavenm_unc_array[1] - surface_wavenm_unc_array[0]
## adding chains for total surface area
#surface_wavenm_total += surface_wavenm[i]
#fwhm_unc_array = 0
#stddev_wavenm = 0
#surface_wavenm = 0
#surface_wavenm_unc_array = 0
#chains = 0
#fwhm = 0
######
# Calculating surface area of all the features added
######
surface_area_cm_total[i] = np.sum(surface_area_cm[i,:]) # collapse and add in one dimension
surface_wavenm_unc_array_total = np.percentile(surface_wavenm_total, [16, 50, 84])
surface_area_error_total[i] = 0.5 * (surface_wavenm_unc_array_total[2] - surface_wavenm_unc_array_total[0])
surface_area_unc_plus_total[i] = surface_wavenm_unc_array_total[2] - surface_wavenm_unc_array_total[1]
surface_area_unc_minus_total[i] = surface_wavenm_unc_array_total[1] - surface_wavenm_unc_array_total[0]
surface_wavenm_unc_array_total = 0
surface_wavenm_total = 0
#break
# plotting for each of the 5 individual lines and the total
names_features = (["feature_1", "feature_2", "feature_3", "feature_4", "feature_5"])
for k in range(0,4):
plt.plot(Av_array,amplitude_array[:,k],linestyle="",marker="o")
plt.errorbar(Av_array,amplitude_array[:,k], yerr=amplitude_array_error[:,k], marker='o',linestyle="")
plt.ylabel(r'Amplitude', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"amplAv"+names_features[k]+".pdf",dpi=200)
plt.show()
plt.plot(Av_array,mean_array[:,k],linestyle="",marker="o")
plt.errorbar(Av_array,mean_array[:,k], yerr=mean_array_error[:,k], marker='o',linestyle="")
plt.ylabel(r'Mean (micron)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"meanAv"+names_features[k]+".pdf",dpi=200)
plt.show()
plt.plot(Av_array,stddev_array[:,k],linestyle="",marker="o")
plt.errorbar(Av_array,stddev_array[:,k], yerr=stddev_array_error[:,k], marker='o',linestyle="")
plt.ylabel(r'stddev (micron)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"stddevAv"+names_features[k]+".pdf",dpi=200)
plt.show()
plt.plot(Av_array,fwhm_array[:,k],linestyle="",marker="o")
plt.errorbar(Av_array,fwhm_array[:,k], yerr=fwhm_array_error[:,k], marker='o',linestyle="")
plt.ylabel(r'FWHM (micron)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"fwhmAv"+names_features[k]+".pdf",dpi=200)
plt.show()
plt.plot(Av_array,surface_area_cm[:,k],linestyle="",marker="o")
plt.errorbar(Av_array,surface_area_cm[:,k], yerr=surface_area_error[:,k], marker='o',linestyle="")
plt.ylabel(r'surface area (cm^-1)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaAv"+names_features[k]+".pdf",dpi=200)
plt.show()
surface_asymmetric_error = [surface_area_unc_minus[:,k], surface_area_unc_plus[:,k]]
plt.plot(Av_array_new,surface_area_cm[:,k],linestyle="",marker="o")
plt.errorbar(Av_array_new,surface_area_cm[:,k], yerr=surface_asymmetric_error, marker='o',linestyle="")
plt.ylabel(r'surface area', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaAvnew"+names_features[k]+".pdf",dpi=200)
plt.show()
plt.plot(Rv_array_new,surface_area_cm[:,k], linestyle="",marker="o")
plt.errorbar(Rv_array_new,surface_area_cm[:,k], yerr=surface_asymmetric_error, marker='o',linestyle="")
plt.ylabel(r'surface area', fontsize=14)
plt.xlabel(r'Rv', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaRv"+names_features[k]+".pdf",dpi=200)
plt.show()
print_results(names_features[k], names_miri_2mass, Av_array, mean_array[:,k], fwhm_array[:,k], fwhm_array_error[:,k], surface_area_cm[:,k], surface_area_error[:,k])
surface_asymmetric_error_total = [surface_area_unc_minus_total, surface_area_unc_plus_total]
#print_results('', names_miri_2mass, Av_array, mean_array, fwhm_array, fwhm_array_error, surface_area_cm, surface_area_error)
plt.plot(Av_array_new,surface_area_cm_total,linestyle="",marker="o")
plt.errorbar(Av_array_new,surface_area_cm_total, yerr=surface_asymmetric_error_total, marker='o',linestyle="")
plt.ylabel(r'surface area', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaAvnew_total.pdf",dpi=200)
plt.show()
plt.plot(Rv_array_new,surface_area_cm_total,linestyle="",marker="o")
plt.errorbar(Rv_array_new,surface_area_cm_total, yerr=surface_asymmetric_error_total, marker='o',linestyle="")
plt.ylabel(r'surface area', fontsize=14)
plt.xlabel(r'Rv', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaRv_total.pdf",dpi=200)
plt.show()
plt.plot(Av_array_new, max_tau_values,linestyle="",marker="o")
plt.ylabel(r'feature_depth', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.xlim(xmin=0)
plt.savefig(output_directory_plots+"featuredepthAvnew_total.pdf",dpi=200)
plt.show()