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)