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 * 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()
tuple_return2 = [tau_data, power_cont(continuum_wav)]
return(tuple_return2)
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=5000, 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*5000))
#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,stddevs]
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/rebin/'
output_directory_plots_emcee = '/home/zeegers/wisci_first_shot/carbonyl_github/carbonyl_24_3/carbonyl/silicon_feature/plots_fshot/rebin/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()
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.59) & (wavelength_merged_new <= 13.1))
new_waves= wavelength_merged_new[continuum_sel]
new_fluxs = flux_merged_new[continuum_sel]
continuum_error = uncs_merged_new[continuum_sel]
# Let's rebin the spectrum here
rebin = [6.6,13.0,1000]
continuum_wav = np.arange(rebin[0],rebin[1],
(rebin[1]-0.5*(rebin[1]-rebin[0]))/rebin[2])
continuum_flux = spectres(continuum_wav,new_waves, new_fluxs)
# We need to get the signal to noise from the ETC calculation, probably the best we've got so far?
# Line list for CPD: 66371.2
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))
feature_sil = ((continuum_wav > 8.0) & (continuum_wav <= 12.3))
stellar_line1 = ((continuum_wav > 6.92) & (continuum_wav <= 6.953))
stellar_line2 = ((continuum_wav > 7.45) & (continuum_wav <= 7.466))
stellar_line3 = ((continuum_wav > 7.424) & (continuum_wav <=7.53))
stellar_line4 = ((continuum_wav > 7.76) & (continuum_wav <=7.785))
stellar_line5 = ((continuum_wav > 8.124) & (continuum_wav <=8.179))
stellar_line6 = ((continuum_wav > 11.25) & (continuum_wav <= 11.35))
stellar_line7 = ((continuum_wav > 12.30) & (continuum_wav <= 12.40))
stellar_line8 = ((continuum_wav > 8.7) & (continuum_wav <= 8.8))
stellar_line9 = ((continuum_wav > 9.6) & (continuum_wav <= 9.8))
stellar_line10 = ((continuum_wav > 12.565) & (continuum_wav <= 12.631))
stellar_line11 = ((continuum_wav > 13.0) & (continuum_wav <= 13.2))
#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 =~ (stellar_line2|stellar_line3|feature_sil|stellar_line6|stellar_line7) # can be expanded if more features will be added
sel_lines =~ (stellar_line2|stellar_line3|stellar_line6|stellar_line7) # 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|stellar_line8|stellar_line9|stellar_line10|stellar_line11)
#sel_lines = sel
sel_extreme =~ (stellar_line1|stellar_line2|stellar_line3|stellar_line4|stellar_line5|feature_sil|stellar_line6|stellar_line7|stellar_line8|stellar_line9|stellar_line10|stellar_line11)
if i==7:
optical_depth_return = optical_depthfit(continuum_wav, continuum_flux, sel_extreme)
else: optical_depth_return = optical_depthfit(continuum_wav, continuum_flux, sel)
tau_data = optical_depth_return[0]
power_cont_array = optical_depth_return[1]
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].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 = fit_result[2]
fig,axs=plt.subplots(3,1, height_ratios=[2, 3, 1], figsize=(3.6, 5.5))
axs[0].plot(continuum_wav, continuum_flux*continuum_wav**2., color = "black")
axs[0].plot(continuum_wav[sel_extreme], continuum_flux[sel_extreme]*continuum_wav[sel_extreme]**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\cdot$ F($\lambda$) [$\mu\mathrm{m}^2$ Jy]', fontsize=10)
axs[0].yaxis.set_tick_params(labelsize=9)
axs[0].set_xlim([6.6, 13.0])
axs[0].set_xticks([])
axs[1].errorbar(continuum_wav, tau_data, yerr=stddev, marker='',linestyle="-", color="black", ecolor="grey")
axs[1].plot(continuum_wav, tau_model2, color="red", lw=2, zorder=10, label = r'10 $\mu$m silicate')
axs[1].axhline(y=0, color='k', ls=":", c="k")
axs[1].set_ylabel('optical depth', fontsize=10)
axs[1].set_ylim([-0.06, 0.26])
axs[1].yaxis.set_tick_params(labelsize=9)
axs[1].set_xlim([6.6, 13.0])
axs[1].set_xticks([])
axs[2].plot(continuum_wav, (tau_data-tau_model2)/stddev, color="grey")
axs[2].axhline(y=0, color='k', ls=":", c="k")
axs[2].set_ylim([-6.4, 6.4])
axs[2].set_xlim([6.6, 13.0])
axs[2].xaxis.set_tick_params(labelsize=9)
axs[2].yaxis.set_tick_params(labelsize=9)
axs[2].set_ylabel('residuals', fontsize=9)
axs[2].set_xlabel(r'Wavelength [$\mu$m]',fontsize=9)
plt.subplots_adjust(hspace=0,left=0.22, right=0.96, top=0.98, bottom = 0.13)
plt.show()
fig.savefig(output_directory_plots+names_miri[i]+"plot_res_carbonyl.pdf",dpi=300)
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()