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 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 optical_depthfit(continuum_wav, continuum_flux, sel):
#linfitter = fitting.LinearLSQFitter()
#fit = fitting.LevMarLSQFitter()
fit = fitting.SimplexLSQFitter()
#poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel]*continuum_wav[sel]**2.)
model_power = models.PowerLaw1D(amplitude = 0.5, x_0=1., alpha=2.)
power_cont = fit(model_power, continuum_wav[sel], continuum_flux[sel], maxiter=10000)
print(power_cont)
tau_data = np.log(power_cont(continuum_wav)/continuum_flux)
plt.plot(continuum_wav, continuum_flux*continuum_wav**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()
return(tau_data)
def fit_carbonyl_gauss_skewed(settings_array, continuum_wav, continuum_flux, sel, sel_lines):
fit_cont = fitting.SimplexLSQFitter()
fit = fitting.LevMarLSQFitter()
#linfitter = fitting.LinearLSQFitter()
#poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel])
power_cont = fit_cont(models.PowerLaw1D(amplitude = 0.5, x_0=1., alpha=2.), continuum_wav[sel], continuum_flux[sel], maxiter=10000)
#tau_data = np.log(poly_cont(continuum_wav)/continuum_flux)
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_skewgaus=skewed_gaussian(amplitude=settings_array[0], stddev=settings_array[1], mean=settings_array[2],gamma=settings_array[3])
print()
#gl_init_gaus=models.Gaussian1D(amplitude=settings_array[0], stddev=settings_array[3], mean=settings_array[6],bounds={'amplitude':(settings_array[1],settings_array[2]), 'stddev':(settings_array[4],settings_array[5]), 'mean':(settings_array[6]-settings_array[7],settings_array[6]+settings_array[8])})
# plotting the initial models
init_model = gl_init_skewgaus
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')
gl_fit = fit(gl_init_skewgaus, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights, maxiter=10000)
return gl_fit
def fit_carbonyl_errors_gauss_skewed(name, settings_array, continuum_wav, continuum_flux, sel, sel_lines):
# calculates errors
# and fitting with Emcee
fit_cont = fitting.SimplexLSQFitter()
fit = fitting.LevMarLSQFitter()
#linfitter = fitting.LinearLSQFitter()
#poly_cont = linfitter(models.Polynomial1D(3), continuum_wav[sel], continuum_flux[sel])
power_cont = fit_cont(models.PowerLaw1D(amplitude = 0.5, x_0=1., alpha=2.), continuum_wav[sel], continuum_flux[sel], maxiter=10000)
#tau_data = np.log(poly_cont(continuum_wav)/continuum_flux)
tau_data = np.log(power_cont(continuum_wav)/continuum_flux)
stddevs = np.std(tau_data[sel])
weights = 1.0 / (stddevs)
gl_init_skewgaus=skewed_gaussian(amplitude=settings_array[0], stddev=settings_array[1], mean=settings_array[2],gamma=settings_array[3])
gl_fit = fit(gl_init_skewgaus, continuum_wav[sel_lines], tau_data[sel_lines], weights=weights, maxiter=100000)
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_24_3/carbonyl/silicon_feature/plots_fshot/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))
#log_probs = fit2.fit_info['sampler'].get_log_prob(flat=True,discard=np.int32(0.1*1000)) # Also extract the log_probs
tuple_return = [fit_mcmc_result, chains]
print('tuple',tuple_return)
return tuple_return
#return fitparams, self.fit_info
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_mode, wav_mode_error, optical_depths, optical_depths_error, fwhm, fwhm_error, surface_area, surface_area_error):
# in the table we need: name of the source, mean wavelength, FWHM + error, integrated 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/silicon_feature/sil_res_fshot/'
table_txt = Table(
names=(
"Source",
"AV",
"Mode [micron]",
"Mode error [micron]",
"FWHM[micron]",
"FWHM error [micron]",
"Optical depth [micron]",
"Optical depth error [micron]",
"Integrated Area [micron]",
"Integrated Area error [micron]"
),
dtype=(
"str",
"float64",
"float64",
"float64",
"float64",
"float64",
"float64",
"float64",
"float64",
"float64"
),
)
# table latex
table_latex = Table(
names=(
"Source",
r"$A_{V}$",
r"Mode (\micron)",
r"Optical depth (\micron)",
"FWHM",
"Integrated Area"
),
dtype=("str", "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_mode[i],
wav_mode_error[i],
optical_depths[i],
optical_depths_error[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_mode[i]:.3f}\pm{wav_mode_error[i]:.3f}$',
f'${optical_depths[i]:.3f}\pm{optical_depths_error[i]:.4f}$',
f'${fwhm[i]:.3f}\pm{fwhm_error[i]:.3f}$',
f'${surface_area[i]:.3f}\pm {surface_area_error[i]:.3f}$'
)
)
tabname = "fit_fshots_micron"
# 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 = ['10lac', '2MASSJ085747','2MASSJ150958']
names_miri_2mass = ['10lac','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/silicon_feature'
filename_results = 'output_res_carbonyl.txt'
# output directory for plots
output_directory_plots = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/silicon_feature/plots_fshot/'
output_directory_plots_emcee = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/silicon_feature/plots_fshot/emcee_res'
Av_array = np.array([0.21, 5.1, 4.7]) # old estimates
Av_array_new = np.array([0.21, 4.98, 4.35])
Rv_array_new = ([3.1, 3.32, 3.13])
logg = ([36000, 24380, 15330])
temperature =([4.03, 3.25, 3.4])
# trying to fit with skewed gaussian:
# amplitude, ampl range, stdev, stdev range, wavelength, wav range + and -
settings_array_gauss = np.array([[0.2, 0.1, 1.0,0.8, 0.3, 3.0, 9.8, 0.2, 0.2],
[0.15, 0.05, 1.0, 0.8, 0.3, 3.0, 9.8, 0.05, 0.05],
[0.3, 0.05, 1.0, 0.8, 0.3, 3.0, 9.8, 0.05, 0.05]
])
# output arrays
amplitude_array = np.zeros((3))
amplitude_array_error = np.zeros((3))
amplitude_array_unc_plus = np.zeros((3))
amplitude_array_unc_minus = np.zeros((3))
mode_array = np.zeros((3))
mode_array_error = np.zeros((3))
mode_array_unc_plus = np.zeros((3))
mode_array_unc_minus = np.zeros((3))
mean_array = np.zeros((3))
mean_array_error = np.zeros((3))
mean_array_unc_plus = np.zeros((3))
mean_array_unc_minus = np.zeros((3))
stddev_array = np.zeros((3))
stddev_array_error = np.zeros((3))
stddev_array_unc_plus = np.zeros((3))
stddev_array_unc_minus = np.zeros((3))
fwhm_array = np.zeros((3))
fwhm_array_error = np.zeros((3))
fwhm_array_unc_plus = np.zeros((3))
fwhm_array_unc_minus = np.zeros((3))
surface_area_cm = np.zeros((3))
surface_area_error = np.zeros((3))
surface_area_unc_plus = np.zeros((3))
surface_area_unc_minus = np.zeros((3))
optical_depths = np.zeros((3))
optical_depths_error = np.zeros((3))
optical_depths_unc_plus = np.zeros((3))
optical_depths_unc_minus = np.zeros((3))
# try to fit with Drude profile
settings_array_drude = np.array([[0.2, 0.1,9.8,0.001,0.001,0.01,0.3],
[0.2, 2.0,9.8,0.001,0.001,0.01,0.3],
[0.2, 2.0,9.8,0.001,0.001,0.01,0.3]
])
# scale (central amplitude), x_o (central wavelength), gamma_o(full-width-half-maximum of profile), asym(asymmetry where a value of 0 results in a standard Drude profile)
settings_array_drude_modified = np.array([[0.2, 0.1,9.8,0.001],
[0.2, 9.8,2.0,0.01],
[0.2, 9.8,2.0,0.01]
])
#amplitude, stdev, wavelength, gamma
#settings_array_skewed_gauss = np.array([[0.02, 0.1,9.8,0.0001],
#[0.2, 0.1,9.8,0.001],
#[0.2, 0.1,9.8,0.001]
#])
settings_array_skewed_gauss = np.array([[0.001, 1.0,9.8,0.01],
[0.1, 0.1,9.8,0.01],
[0.1, 0.1,9.8,0.01]
])
#filename_plots = names_miri+''
for i in range(0,len(names_miri)):
data_miri_merged = fits.getdata(direc + 'data/jwst/delivery_v6/' + names_miri[i]+datafile)
wavelength_merged = data_miri_merged['WAVELENGTH']
flux_merged = data_miri_merged['FLUX']
uncs_merged = data_miri_merged['UNC']
plt.plot(wavelength_merged,flux_merged)
plt.show()
#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 > 6.6) & (wavelength_merged_new <= 13.0))
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]
feature_sil = ((continuum_wav > 8.0) & (continuum_wav <= 12.3))
stellar_line1 = ((continuum_wav > 7.45) & (continuum_wav <= 7.466))
stellar_line2 = ((continuum_wav > 7.49) & (continuum_wav <=7.51))
stellar_line3 = ((continuum_wav > 11.25) & (continuum_wav <= 11.35))
stellar_line4 = ((continuum_wav > 12.30) & (continuum_wav <= 12.40))
stellar_line5 = ((continuum_wav > 8.7) & (continuum_wav <= 8.8))
stellar_line6 = ((continuum_wav > 9.6) & (continuum_wav <= 9.8))
stellar_line7 = ((continuum_wav > 13.0) & (continuum_wav <= 13.2))
# i = 0 lines at 7 micron?
# i = 1
# i = 2 7.46 7.503 11.30 12.37
# i = 3 has above lines and more less strong
# i = 4 not many lines doesn't have to ignored
# i = 5 7.46 7.503 11.30 12.37
# i = 6 no lines that are strong
# i = 7 many lines ignore, the strongest are: two lines 7.41 - 7.31 8.76 9.7 11.30, 12.37, 13.1
# i = 8 not many lines: 7.46 7.503 11.30 12.37
# i = 9 weak lines: 7.46 7.503 11.30 12.37
# i = 10 almost no lines: 7.46 7.503
# i = 11 weak lines: 7.46 7.503 11.30 12.37
# i = 12 almost no lines: 12.30
#sel =~ (stellar_line1|stellar_line2|feature_sil|stellar_line3|stellar_line4) # can be expanded if more features will be added
sel =~ (stellar_line1|stellar_line2|feature_sil|stellar_line3|stellar_line4) # can be expanded if more features will be added
sel_lines =~ (stellar_line1|stellar_line2|stellar_line3|stellar_line4) # can be expanded if more features will be added
sel_lines_extreme =~ (stellar_line1|stellar_line2|stellar_line3|stellar_line4|stellar_line5|stellar_line6|stellar_line7)
#sel_lines = sel
sel_extreme =~ (stellar_line1|stellar_line2|feature_sil|stellar_line3|stellar_line4|stellar_line5|stellar_line6|stellar_line7)
if i==7:
tau_data = optical_depthfit(continuum_wav, continuum_flux, sel_extreme)
else: tau_data = optical_depthfit(continuum_wav, continuum_flux, sel)
if i==7:
modified_drude_model = fit_carbonyl_gauss_skewed(settings_array_skewed_gauss[i], continuum_wav, continuum_flux, sel_extreme, sel_lines_extreme)
else: modified_drude_model = fit_carbonyl_gauss_skewed(settings_array_skewed_gauss[i], continuum_wav, continuum_flux, sel, sel_lines)
#modified_drude_model = modified_drude(scale=settings_array_drude_modified[i,0], x_o=settings_array_drude_modified[i,1], gamma_o=settings_array_drude_modified[i,2], asym=settings_array_drude_modified[i,3])
#modified_drude_model = fit_carbonyl_gauss(settings_array_gauss[i], continuum_wav, continuum_flux, sel, sel_lines)
#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)
tau_model = modified_drude_model(continuum_wav)
#print(gl_fit_drude)
print(modified_drude_model)
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()
if i==7:
fit_result = fit_carbonyl_errors_gauss_skewed(names_miri[i], settings_array_skewed_gauss[i], continuum_wav, continuum_flux, sel_extreme, sel_lines_extreme)
else: fit_result = fit_carbonyl_errors_gauss_skewed(names_miri[i], settings_array_skewed_gauss[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 = np.std(tau_model2)
fig,axs=plt.subplots(2,1)
axs[0].errorbar(continuum_wav, tau_data, yerr=stddev, marker='',linestyle="-", color="grey", errorevery=60)
#axs[0].plot(continuum_wav, tau_data,".")
axs[0].plot(continuum_wav, tau_model2, color="orange", lw=2, zorder=10)
axs[0].axhline(y=0, color='k', ls=":", c="k")
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_model2)/stddev**2, color="grey")
axs[1].set_ylabel('residuals', fontsize=14)
axs[1].set_xlabel(r'Wavelength ($\mu$m)', fontsize=14)
axs[1].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_silicates.pdf",dpi=200, format="pdf")
plt.close()
amp_unc_ave = fit_mcmc_result.amplitude.unc
stddev_unc_ave = fit_mcmc_result.stddev.unc
mean_unc_ave = fit_mcmc_result.mean.unc
amplitude_array[i]=fit_mcmc_result.amplitude.value
amplitude_array_error[i]=amp_unc_ave
mean_array[i]=fit_mcmc_result.mean.value
mean_array_error[i]=mean_unc_ave
stddev_array[i] = fit_mcmc_result.stddev.value
stddev_array_error[i]=stddev_unc_ave
### calculate uncertainties using chains
chains = fit_result[1]
### integrated area and uncertainties
# We integrated the skewed gaussian function, for more information on the function and CDF see:
# https://en.wikipedia.org/wiki/Skew_normal_distribution
#x = 8. - 13. micron
amplitudes = chains[:,0]
means = chains[:,1]
stddeviation = chains[:,2]
alphas = chains[:,3]
h_min = (8.-means)/stddeviation
h_max = (13.-means)/stddeviation
xmax = (13.-means)/(stddeviation)
xmin = (8.-means)/(stddeviation)
Owens_max = 2.*(special.owens_t(h_max, alphas, out=None))
psi_max = 0.5 * (1.+special.erf(xmax/(np.sqrt(2))))
Owens_min = 2.*(special.owens_t(h_min, alphas, out=None))
psi_min = 0.5 * (1.+special.erf(xmin/(np.sqrt(2))))
surface_area = amplitudes*((psi_max - Owens_max)-(psi_min - Owens_min))
surface_area_unc_array = np.percentile(surface_area, [16, 50, 84])
surface_area_cm[i] = surface_area_unc_array[1]
surface_area_error[i] = 0.5 * (surface_area_unc_array[2] - surface_area_unc_array[0])
surface_area_unc_plus[i] = surface_area_unc_array[2] - surface_area_unc_array[1]
surface_area_unc_minus[i] = surface_area_unc_array[1] - surface_area_unc_array[0]
# Calculating the mean of a skewed gaussian
delta = alphas/(np.sqrt(1+alphas**2))
mean_skewed = means + stddeviation*delta*np.sqrt(2./np.pi)
means_skewed = np.percentile(mean_skewed, [16, 50, 84])
mean_array[i] = means_skewed[1]
mean_array_error[i] = 0.5 * (means_skewed[2] - means_skewed[0])
mean_array_unc_plus[i] = means_skewed[2] - means_skewed[1]
mean_array_unc_minus[i] = means_skewed[2] - means_skewed[0]
# Calculating the variance?
# Does the FWHM even mean anything here?
# Calculating the mode (maximum value)
m_0 = np.sqrt(2./np.pi)*delta - (1-np.pi/4.)*((np.sqrt(2./np.pi)*delta)**3./(1.-2./np.pi*delta**2.))-np.sign(alphas)/2.*np.exp(-2.*np.pi/np.absolute(alphas))
modes = means + stddeviation*m_0
mode = np.percentile(modes, [16, 50, 84])
mode_array[i] = mode[1]
mode_array_error[i] = 0.5 * (mode[2] - mode[0])
mode_array_unc_plus[i] = mode[2] - mode[1]
mode_array_unc_minus[i] = mode[1] - mode[0]
# half widht at half maximum
# find x at the mode and at half of the mode
y_modes = amplitudes * 2/stddeviation/np.sqrt(2 * np.pi) * np.exp(-0.5 * ((modes-means)/stddeviation)**2)*0.5 * (1+special.erf(alphas * (modes-means)/stddeviation/np.sqrt(2)))
y_mode = np.percentile(y_modes, [16, 50, 84])
optical_depths[i]=y_mode[1]
optical_depths_error[i] = 0.5 * (y_mode[2] - y_mode[0])
optical_depths_unc_plus[i] = y_mode[2] - y_mode[1]
optical_depths_unc_minus[i] = y_mode[1] - y_mode[0]
# looping over the chains
fwhms = np.empty(len(y_modes))
for m in range(0,len(y_modes)):
wave_min = continuum_wav[np.argsort((np.abs(tau_model2-y_modes[m]/2.)))]
fwhms[m] = np.max(wave_min[0:4]) - np.min(wave_min[0:4])
fwhm = np.percentile(fwhms, [16, 50, 84])
fwhm_array[i] = fwhm[1]
fwhm_array_error[i] = 0.5 * (fwhm[2] - fwhm[0])
fwhm_array_unc_plus[i] = fwhm[2] - fwhm[1]
fwhm_array_unc_minus[i] = fwhm[1] - fwhm[0]
# surface area in cm^-1
#function_cm = amplitude * np.exp(-0.5 * (x - mean) ** 2 / stddev**2)*(1-special.erf((gamma*(x-mean))/(stddev*math.sqrt(2))))*-1.*x^-2.
#function_cm_int = integrate.quad(lambda x: amplitudes * np.exp(-0.5 * (1/x - 1/means) ** 2 / 1/stddeviation**2)*(1-special.erf((alphas*(1/x-1/means))/(1/stddeviation*math.sqrt(2))))*-10**4.*x^-2. , 1./8., 1./13.)
#stddev_wavenm = 1./(chains[:,1]**2)*stddeviation*10000.
#surface_wavenm = chains[:,0]*stddev_wavenm * np.sqrt(2. * np.pi)
#surface_wavenm_unc_array = np.percentile(surface_wavenm, [16, 50, 84])
#surface_area_cm[i] = surface_wavenm_unc_array[1]
#surface_area_error[i] = 0.5 * (surface_wavenm_unc_array[2] - surface_wavenm_unc_array[0])
#surface_area_unc_plus[i] = surface_wavenm_unc_array[2] - surface_wavenm_unc_array[1]
#surface_area_unc_minus[i] = surface_wavenm_unc_array[1] - surface_wavenm_unc_array[0]
#break
plt.plot(Av_array,mode_array,linestyle="",marker="o")
plt.errorbar(Av_array,mode_array, yerr=mode_array_error, marker='o',linestyle="")
plt.ylabel(r'Mode (micron)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"modeAv.pdf",dpi=200)
plt.show()
plt.plot(Av_array,mean_array,linestyle="",marker="o")
plt.errorbar(Av_array,mean_array, yerr=mean_array_error, marker='o',linestyle="")
plt.ylabel(r'Mean (micron)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"meanAv.pdf",dpi=200)
plt.show()
plt.plot(Av_array,stddev_array,linestyle="",marker="o")
plt.errorbar(Av_array,stddev_array, yerr=stddev_array_error, marker='o',linestyle="")
plt.ylabel(r'stddev (micron)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"stddevAv.pdf",dpi=200)
plt.show()
#plt.plot(Av_array,fwhm_array,linestyle="",marker="o")
#plt.errorbar(Av_array,fwhm_array, yerr=fwhm_array_error, marker='o',linestyle="")
#plt.ylabel(r'FWHM (micron)', fontsize=14)
#plt.xlabel(r'Av', fontsize=14)
#plt.savefig(output_directory_plots+"fwhmAv.pdf",dpi=200)
#plt.show()
plt.plot(Av_array,surface_area_cm,linestyle="",marker="o")
plt.errorbar(Av_array,surface_area_cm, yerr=surface_area_error, marker='o',linestyle="")
plt.ylabel(r'integrated area ($micron$)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaAv.pdf",dpi=200)
plt.show()
surface_asymmetric_error = [surface_area_unc_minus, surface_area_unc_plus]
plt.plot(Av_array_new,surface_area_cm,linestyle="",marker="o")
plt.errorbar(Av_array_new,surface_area_cm, yerr=surface_asymmetric_error, marker='o',linestyle="")
plt.ylabel(r'ingrated area ($micron$)', fontsize=14)
plt.xlabel(r'Av', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaAvnew.pdf",dpi=200)
plt.show()
plt.plot(Rv_array_new,surface_area_cm,linestyle="",marker="o")
plt.errorbar(Rv_array_new,surface_area_cm, yerr=surface_asymmetric_error, marker='o',linestyle="")
plt.ylabel(r'integrated area ($cm^{-1}$)', fontsize=14)
plt.xlabel(r'Rv', fontsize=14)
plt.savefig(output_directory_plots+"surfaceareaRv.pdf",dpi=200)
plt.show()
print_results(names_miri, names_miri_2mass, Av_array, mode_array, mode_array_error, optical_depths, optical_depths_error, fwhm_array, fwhm_array_error, surface_area_cm, surface_area_error)