# 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() fit3 = fitting.SimplexLSQFitter() 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 = fit3(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_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_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)}) #, 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 = fit3(gl_init_gaus1+gl_init_gaus2+gl_init_gaus3+gl_init_gaus4+gl_init_gaus5+stellar_line_init, 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=20000) tau_data = np.log(power_cont(continuum_wav)/continuum_flux) stddevs = np.std(tau_data[sel]) weights = 1.0 / (stddevs) 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) # only the dust and feature gl_init_gaus_fit1=models.Gaussian1D(amplitude=fit_result[0].amplitude_0.value, stddev=fit_result[0].stddev_0.value, mean=fit_result[0].mean_0.value) gl_init_gaus_fit2=models.Gaussian1D(amplitude=fit_result[0].amplitude_1.value, stddev=fit_result[0].stddev_1.value, mean=fit_result[0].mean_1.value) gl_init_gaus_fit3=models.Gaussian1D(amplitude=fit_result[0].amplitude_2.value, stddev=fit_result[0].stddev_2.value, mean=fit_result[0].mean_2.value) gl_init_gaus_fit4=models.Gaussian1D(amplitude=fit_result[0].amplitude_3.value, stddev=fit_result[0].stddev_3.value, mean=fit_result[0].mean_3.value) gl_init_gaus_fit5=models.Gaussian1D(amplitude=fit_result[0].amplitude_4.value, stddev=fit_result[0].stddev_4.value, mean=fit_result[0].mean_4.value) model_dust = gl_init_gaus_fit1 + gl_init_gaus_fit2 + gl_init_gaus_fit3 + gl_init_gaus_fit4 + gl_init_gaus_fit5 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], 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\cdot$ F($\lambda$) [$\mu\mathrm{m}^2$ Jy]', fontsize=10) axs[0].yaxis.set_tick_params(labelsize=9) axs[0].set_xlim([2.6, 3.9]) 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_data, marker='',linestyle="-", color = "black", zorder=10) axs[1].plot(continuum_wav, tau_model2, color="cyan", lw=2, zorder=50, label = 'stellar line') axs[1].plot(continuum_wav, model_dust(continuum_wav), color="red", lw=2, zorder=100, label = r'carbon') axs[1].axhline(y=0, color='k', ls=":", c="k") axs[1].set_ylabel('optical depth', fontsize=10) axs[1].yaxis.set_tick_params(labelsize=9) axs[1].set_xlim([2.6, 3.9]) axs[1].set_xticks([]) axs[1].legend(fontsize=8, facecolor='white', framealpha=0.8) axs[2].plot(continuum_wav, (tau_data-tau_model2)/stddev, color="grey") axs[2].set_ylabel('residuals', fontsize=9) axs[2].set_xlabel(r'Wavelength [$\mu$m]',fontsize=10) axs[2].axhline(y=0, color='k', ls=":", c="k") axs[2].set_xlim([2.6, 3.9]) axs[2].xaxis.set_tick_params(labelsize=9) axs[2].yaxis.set_tick_params(labelsize=9) plt.subplots_adjust(hspace=0,left=0.22, right=0.98, top=0.98, bottom = 0.13) 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()