def plot_3orders_ms(fig,axes,solarx, solary,x_grid,y_grid,x_grid_bin, y_grid_bin,
file_grid, path, file_df,wvl_slices,flatfields,
fit_matrices, fit_matrices_err, fit_bin_matrices, fit_bin_matrices_err,
cont_slice_1 = slice(0,10), cont_slice_2 = slice(40,50),
params_prec = {"int":2,"wvl":1,"fwhm":1},title=None,title_ycenter=None):
y_select_pixel,x_select_pixel = find_nearest_pixel(solarx, solary,
x_grid,y_grid)
y_select_pixel_bin,x_select_pixel_bin = find_nearest_pixel(solarx, solary,
x_grid_bin,y_grid_bin)
solarx_select, solary_select = x_grid[y_select_pixel,x_select_pixel], y_grid[y_select_pixel,x_select_pixel]
solarx_select_bin, solary_select_bin = (x_grid_bin[y_select_pixel_bin,x_select_pixel_bin],
y_grid_bin[y_select_pixel_bin,x_select_pixel_bin])
file_index = file_grid[x_select_pixel]
print(file_df.iloc[file_index,0], 69 - y_select_pixel_bin)
frame = CCDData.read(os.path.join(path,file_df.iloc[file_index,0]),hdu=0,unit="adu")
wvl = CCDData.read(os.path.join(path,file_df.iloc[file_index,0]),hdu=2,unit="adu").data
if title is None:
title_text = r"$({:.0f}'',{:.0f}'')$".format(solarx_select_bin, solary_select_bin)
else:
title_text = r"{:} $({:.0f}'',{:.0f}'')$".format(title,solarx_select_bin, solary_select_bin)
fig.text(1,title_ycenter,title_text,fontsize=20,rotation=270)
show_single_fit_ms(axes[0,0],frame=frame,wvl=wvl,order=61,flatfield_1d=flatfields[0],y_select=y_select_pixel,
x_select=x_select_pixel,solarx=solarx_select,solary=solary_select,fit_matrix=fit_matrices[0],
fit_matrix_err=fit_matrices_err[0],wvl_slice=wvl_slices[0],cont_slice_1=cont_slice_1,
cont_slice_2=cont_slice_2,params_prec=params_prec)
show_single_fit_ms(axes[0,1],frame=frame,wvl=wvl,order=62,flatfield_1d=flatfields[1],y_select=y_select_pixel,
x_select=x_select_pixel,solarx=solarx_select,solary=solary_select,fit_matrix=fit_matrices[1],
fit_matrix_err=fit_matrices_err[1],wvl_slice=wvl_slices[1],cont_slice_1=cont_slice_1,
cont_slice_2=cont_slice_2,params_prec=params_prec)
show_single_fit_ms(axes[0,2],frame=frame,wvl=wvl,order=63,flatfield_1d=flatfields[2],y_select=y_select_pixel,
x_select=x_select_pixel,solarx=solarx_select,solary=solary_select,fit_matrix=fit_matrices[2],
fit_matrix_err=fit_matrices_err[2],wvl_slice=wvl_slices[2],cont_slice_1=cont_slice_1,
cont_slice_2=cont_slice_2,params_prec=params_prec)
show_single_fit_mix_ms(axes[0,3],frame=frame,wvl=wvl,orders=[61,62,63],flatfields=flatfields,y_select=y_select_pixel,
x_select=x_select_pixel,solarx=solarx_select,solary=solary_select,fit_matrices=fit_matrices,
fit_matrices_err=fit_matrices_err,wvl_slices=wvl_slices,cont_slice_1=slice(0,10),
cont_slice_2=slice(30,40),params_prec=params_prec)
show_single_fit_ms(axes[1,0],frame=frame,wvl=wvl,order=61,flatfield_1d=flatfields[0],y_select=y_select_pixel_bin,
x_select=x_select_pixel_bin,solarx=solarx_select_bin,solary=solary_select_bin,fit_matrix=fit_bin_matrices[0],
fit_matrix_err=fit_bin_matrices_err[0],wvl_slice=wvl_slices[0],cont_slice_1=cont_slice_1,
cont_slice_2=cont_slice_2,params_prec=params_prec,bin=True)
show_single_fit_ms(axes[1,1],frame=frame,wvl=wvl,order=62,flatfield_1d=flatfields[1],y_select=y_select_pixel_bin,
x_select=x_select_pixel_bin,solarx=solarx_select_bin,solary=solary_select_bin,fit_matrix=fit_bin_matrices[1],
fit_matrix_err=fit_bin_matrices_err[1],wvl_slice=wvl_slices[1],cont_slice_1=cont_slice_1,
cont_slice_2=cont_slice_2,params_prec=params_prec,bin=True)
show_single_fit_ms(axes[1,2],frame=frame,wvl=wvl,order=63,flatfield_1d=flatfields[2],y_select=y_select_pixel_bin,
x_select=x_select_pixel_bin,solarx=solarx_select_bin,solary=solary_select_bin,fit_matrix=fit_bin_matrices[2],
fit_matrix_err=fit_bin_matrices_err[2],wvl_slice=wvl_slices[2],cont_slice_1=cont_slice_1,
cont_slice_2=cont_slice_2,params_prec=params_prec,bin=True)
show_single_fit_mix_ms(axes[1,3],frame=frame,wvl=wvl,orders=[61,62,63],flatfields=flatfields,y_select=y_select_pixel_bin,
x_select=x_select_pixel_bin,solarx=solarx_select_bin,solary=solary_select_bin,fit_matrices=fit_bin_matrices,
fit_matrices_err=fit_bin_matrices_err,wvl_slices=wvl_slices,cont_slice_1=slice(0,10),
cont_slice_2=slice(30,40),params_prec=params_prec,bin=True)
def show_single_fit_ms(ax_spec, frame, wvl, order, flatfield_1d, y_select, x_select, solarx, solary,
fit_matrix,fit_matrix_err,wvl_slice,cont_slice_1,cont_slice_2,params_prec,bin=False,
):
frame_wavelength = wvl/order/10.
profile_to_fit = frame.data/flatfield_1d[:,np.newaxis]/frame.header["exptime"]
profile_to_fit = np.flip(profile_to_fit,axis=0)
readout_noise = np.float64(frame.header["ronoise"])
profile_err = np.sqrt(frame.data*frame.header["exptime"] + readout_noise**2)/flatfield_1d[:,np.newaxis]/frame.header["exptime"]
profile_err = np.flip(profile_err,axis=0)
if bin == True:
profile_to_fit = np.mean(profile_to_fit[y_select*5:(y_select+1)*5,wvl_slice],axis=0)
profile_err = np.sqrt(np.nanmean((profile_err**2)[y_select*5:(y_select+1)*5,wvl_slice],axis=0))/np.sqrt(5)
else:
profile_to_fit = profile_to_fit[y_select,wvl_slice]
profile_err = profile_err[y_select,wvl_slice]
frame_wavelength_sliced = frame_wavelength[wvl_slice]
cont_wvl = frame_wavelength_sliced[np.r_[cont_slice_1, cont_slice_2]]
cont_int = profile_to_fit[np.r_[cont_slice_1, cont_slice_2]]
cont_fit_param = np.polyfit(cont_wvl, cont_int, deg = 1)
cont_fit_poly = np.poly1d(cont_fit_param)
profile_res = profile_to_fit - cont_fit_poly(frame_wavelength_sliced)
profile_fit = gaussian(frame_wavelength_sliced,*fit_matrix[:3,y_select,x_select]) + \
fit_matrix[3,y_select,x_select] - np.mean(cont_fit_poly(frame_wavelength_sliced))
wvl_to_plot = np.linspace(frame_wavelength_sliced[0],frame_wavelength_sliced[-1],301)
profile_fit_to_plot = gaussian(wvl_to_plot,*fit_matrix[:3,y_select,x_select]) + \
fit_matrix[3,y_select,x_select]
fit_res = profile_res - profile_fit
plot_xlim = fit_matrix[0,y_select,x_select] + \
3*np.array([-fit_matrix[2,y_select,x_select],
fit_matrix[2,y_select,x_select]])
ax_spec.errorbar(frame_wavelength_sliced,profile_res,profile_err,ds="steps-mid",capsize=2,
color="#E87A90",label = r"$I_{\rm obs}$",lw=2,zorder=15)
ax_spec.fill_between(frame_wavelength_sliced,
np.ones_like(frame_wavelength_sliced)*np.min(profile_res),profile_res,
step='mid',color="#FEDFE1",alpha=0.6)
ax_spec.plot(wvl_to_plot,profile_fit_to_plot,color="black",ls="-",label = r"$I_{\rm fit}$",lw=2,
zorder=16,alpha=0.7)
ax_spec.set_xlim(plot_xlim)
ax_spec.tick_params(labelsize=16)
order_colors = {61:"#C1272D",62:"#008176",63:"#0000A7"}
ax_spec.text(0.97,0.87,r"\textbf{{{:2d}}}".format(order),fontsize=20,
transform=ax_spec.transAxes,ha="right",va="center",color=order_colors[order])
# ax_spec.text(0.97,0.85,r"{:.2f}\,$R_\odot$".format(np.sqrt(solarx**2 + solary**2)/950.),fontsize=20,
# transform=ax_spec.transAxes,ha="right",va="center")
pixel_size_green = np.abs(np.mean(np.diff(frame_wavelength_sliced)))
inst_width_pix_green = 1.86
inst_width_nm_green = pixel_size_green*inst_width_pix_green
vlos_plot = -(fit_matrix[0,y_select,x_select] - 530.29)/530.29*const.c.cgs.value*1e-5
vlos_err_plot = -(fit_matrix_err[0,y_select,x_select] - 530.29)/530.29*const.c.cgs.value*1e-5
veff_plot = np.sqrt(fit_matrix[2,y_select,x_select]**2 - inst_width_nm_green**2)/530.29*const.c.cgs.value*1e-5/np.sqrt(4*np.log(2))
veff_err_plot = fit_matrix[2,y_select,x_select]/np.sqrt(fit_matrix[2,y_select,x_select]**2 - inst_width_nm_green**2)* \
fit_matrix_err[2,y_select,x_select]/530.29*const.c.cgs.value*1e-5/np.sqrt(4*np.log(2))
vlos_text_fmt = r'$v : {:#.{vlos_data_prec}g}$' #\pm{:#.{vlos_err_prec}g}$'
veff_text_fmt = r'$v_{{\rm eff}} : {:#.{veff_data_prec}g}$' #\pm{:#.{veff_err_prec}g}$'
vlos_data_prec = np.ceil(np.log10(np.abs(vlos_plot))).astype("int") - \
np.ceil(np.log10(vlos_err_plot)).astype("int") + params_prec["wvl"]
ax_spec.text(0.03,0.87,vlos_text_fmt.format(num2tex(vlos_plot),
vlos_data_prec = 2),ha = 'left',va = 'center',
color = 'black',fontsize = 20,linespacing=1.5,transform=ax_spec.transAxes)
veff_data_prec = np.ceil(np.log10(np.abs(veff_plot))).astype("int") - \
np.ceil(np.log10(veff_err_plot)).astype("int") + params_prec["fwhm"]
ax_spec.text(0.03,0.67,veff_text_fmt.format(num2tex(veff_plot),
veff_data_prec = 3),ha = 'left',va = 'center',
color = 'black',fontsize = 20,linespacing=1.5,transform=ax_spec.transAxes)
# ax_spec.set_title(r"Fe \textsc{{xiv}} {:d} order".format(order),fontsize=20)
def show_single_fit_mix_ms(ax_spec, frame, wvl, orders, flatfields, y_select, x_select, solarx, solary,
fit_matrices,fit_matrices_err,wvl_slices,cont_slice_1,cont_slice_2,params_prec,bin=False):
readout_noise = frame.header["ronoise"]
frame_noise = np.sqrt(frame.data*frame.header["exptime"] + readout_noise**2)/frame.header["exptime"]
if bin == False:
images_sliced = [np.flip(frame.data[:,wvl_slice_]/flatfield_[:,np.newaxis]/frame.header["exptime"],axis=0)[y_select,:] \
for (wvl_slice_,flatfield_) in zip(wvl_slices,flatfields)]
errs_sliced = [np.flip(frame_noise[:,wvl_slice_]/flatfield_[:,np.newaxis]/frame.header["exptime"],axis=0)[y_select,:] \
for (wvl_slice_,flatfield_) in zip(wvl_slices,flatfields)]
else:
images_sliced = [np.mean(np.flip(frame.data[:, wvl_slice_]/flatfield_[:,np.newaxis]/frame.header["exptime"],axis=0)[y_select*5:(y_select+1)*5,:],axis=0) \
for (wvl_slice_,flatfield_) in zip(wvl_slices,flatfields)]
errs_sliced = [np.sqrt(np.mean(np.flip(frame_noise[:, wvl_slice_]/flatfield_[:,np.newaxis]/frame.header["exptime"],axis=0)[y_select*5:(y_select+1)*5,:]**2,axis=0))/np.sqrt(5) \
for (wvl_slice_,flatfield_) in zip(wvl_slices,flatfields)]
images_sliced_interp = np.array(images_sliced)
errs_sliced_interp = np.array(errs_sliced)
wvls_sliced = [wvl[wvl_slice_]/order_/10. for \
(wvl_slice_ , order_)in zip(wvl_slices, orders)]
interp_specs = []
for ii in range(3):
wvl_grid = wvls_sliced[ii]*u.nm
wvl_grid_interp = wvls_sliced[2]*u.nm
interp_spline = SplineInterpolatedResampler()
spec_grid = images_sliced[ii][:]*u.adu
spec_err = StdDevUncertainty(errs_sliced[ii][:],unit=u.adu)
if (wvl_grid[1] - wvl_grid[0]) < 0:
input_spec = Spectrum1D(flux=np.flip(spec_grid),spectral_axis=np.flip(wvl_grid),uncertainty=np.flip(spec_err))
else:
input_spec = Spectrum1D(flux=spec_grid,spectral_axis=wvl_grid,uncertainty=spec_err)
interp_spec = interp_spline(input_spec,wvl_grid_interp)
images_sliced_interp[ii,:] = interp_spec.flux.value
errs_sliced_interp[ii,:] = interp_spec.uncertainty.array
interp_specs.append(interp_spec)
# interp_func = interp1d(wvls_sliced[ii], images_sliced[ii],
# kind="cubic",fill_value="extrapolate")
# images_sliced_interp[ii,:] = interp_func(wvls_sliced[2])
# interp_func_err = interp1d(wvls_sliced[ii], errs_sliced[ii],
# kind="linear",fill_value="extrapolate")
# errs_sliced_interp[ii,:] = interp_func_err(wvls_sliced[2])
#now don't normalize the profiles before summation
# image_combined = np.nanmean(images_sliced_interp/np.array(fit_matrices)[:3,1,y_select,x_select,np.newaxis]*fit_matrices[2][1,y_select,x_select],axis=0)[5:-5]
# err_f1 = fit_matrices[2][1,y_select,x_select]/(3*fit_matrices[0][1,y_select,x_select])
# err_f2 = fit_matrices[2][1,y_select,x_select]/(3*fit_matrices[1][1,y_select,x_select])
# err_f3 = 1./3
# err_f4 = - images_sliced_interp[0,:]*fit_matrices[2][1,y_select,x_select,np.newaxis]/(3*fit_matrices[0][1,y_select,x_select,np.newaxis]**2)
# err_f5 = - images_sliced_interp[1,:]*fit_matrices[2][1,y_select,x_select,np.newaxis]/(3*fit_matrices[1][1,y_select,x_select,np.newaxis]**2)
# err_f6 = 1./3*(images_sliced_interp[0,:]/fit_matrices[0][1,y_select,x_select,np.newaxis]+ \
# images_sliced_interp[1,:]/fit_matrices[1][1,y_select,x_select,np.newaxis])
# err_combined = np.sqrt((err_f1*errs_sliced_interp[0,:])**2 + \
# (err_f2*errs_sliced_interp[1,:])**2 + \
# (err_f3*errs_sliced_interp[2,:])**2 + \
# (err_f4*fit_matrices_err[0][1,y_select,x_select,np.newaxis])**2 + \
# (err_f5*fit_matrices_err[1][1,y_select,x_select,np.newaxis])**2 + \
# (err_f6*fit_matrices_err[2][1,y_select,x_select,np.newaxis])**2)
# peak_index = np.argmax(image_combined) + 5
# err_each_term = np.array([(err_f1*errs_sliced_interp[0,peak_index])**2,
# (err_f2*errs_sliced_interp[1,peak_index])**2,
# (err_f3*errs_sliced_interp[2,peak_index])**2,
# (err_f4[peak_index]*fit_matrices_err[0][1,y_select,x_select])**2,
# (err_f5[peak_index]*fit_matrices_err[1][1,y_select,x_select])**2,
# (err_f6[peak_index]*fit_matrices_err[2][1,y_select,x_select])**2])
# print(err_each_term)
image_combined = np.nansum(images_sliced_interp,axis=0)[5:-5]
err_combined = np.sqrt(errs_sliced_interp[0,:]**2 + \
errs_sliced_interp[1,:]**2 + \
errs_sliced_interp[2,:]**2)
err_combined = err_combined[5:-5]
frame_wavelength_sliced = wvls_sliced[2][5:-5]
cont_wvl = frame_wavelength_sliced[np.r_[cont_slice_1, cont_slice_2]]
cont_int = image_combined[np.r_[cont_slice_1, cont_slice_2]]
cont_fit_param = np.polyfit(cont_wvl, cont_int, deg = 1)
cont_fit_poly = np.poly1d(cont_fit_param)
profile_res = image_combined - cont_fit_poly(frame_wavelength_sliced)
wvl_to_plot = np.linspace(frame_wavelength_sliced[0],frame_wavelength_sliced[-1],301)
profile_fit_to_plot = gaussian(wvl_to_plot,*fit_matrices[3][:3,y_select,x_select]) + \
fit_matrices[3][3,y_select,x_select]
plot_xlim = fit_matrices[3][0,y_select,x_select] + \
3*np.array([-fit_matrices[3][2,y_select,x_select],
fit_matrices[3][2,y_select,x_select]])
# fit_model = SpectrumFitSingle(data=profile_res,err=err_combined, wvl=frame_wavelength_sliced,line_number=1,
# line_wvl_init=frame_wavelength_sliced[np.argmax(profile_res)],int_max_init=profile_res.max(),fwhm_init=0.1)
# fit_model.run_lse(absolute_sigma=False)
# fit_model.plot()
# ax_spec.step(frame_wavelength_sliced,image_combined,where="mid",
# color="#E87A90",label = r"$I_{\rm obs}$",lw=2,zorder=15)
# ax_spec.fill_between(frame_wavelength_sliced,
# np.ones_like(frame_wavelength_sliced)*np.min(image_combined),image_combined,
# step='mid',color="#FEDFE1",alpha=0.6)
ax_spec.errorbar(frame_wavelength_sliced,profile_res,err_combined,ds="steps-mid",capsize=2,
color="#E87A90",label = r"$I_{\rm obs}$",lw=2,zorder=15)
ax_spec.fill_between(frame_wavelength_sliced,
np.ones_like(frame_wavelength_sliced)*np.min(profile_res),profile_res,
step='mid',color="#FEDFE1",alpha=0.6)
ax_spec.plot(wvl_to_plot,profile_fit_to_plot,color="black",ls="-",label = r"$I_{\rm fit}$",lw=2,
zorder=16,alpha=0.7)
ax_spec.set_xlim(plot_xlim)
ax_spec.tick_params(labelsize=16)
ax_spec.text(0.97,0.87,r"\textbf{add up}",fontsize=20,
transform=ax_spec.transAxes,ha="right",va="center",)
ax_spec.text(0.97,0.67,r"{:.2f}\,$R_\odot$".format(np.sqrt(solarx**2 + solary**2)/950.),fontsize=20,
transform=ax_spec.transAxes,ha="right",va="center")
vlos_plot = -(fit_matrices[3][0,y_select,x_select] - 530.29)/530.29*const.c.cgs.value*1e-5
vlos_err_plot = -(fit_matrices_err[3][0,y_select,x_select] - 530.29)/530.29*const.c.cgs.value*1e-5
veff_plot = np.sqrt(fit_matrices[3][2,y_select,x_select]**2 - inst_width_nm_green**2)/530.29*const.c.cgs.value*1e-5/np.sqrt(4*np.log(2))
veff_err_plot = fit_matrices[3][2,y_select,x_select]/np.sqrt(fit_matrices[3][2,y_select,x_select]**2 - inst_width_nm_green**2)* \
fit_matrices_err[3][2,y_select,x_select]/530.29*const.c.cgs.value*1e-5/np.sqrt(4*np.log(2))
vlos_text_fmt = r'$v : {:#.{vlos_data_prec}g}$' #\pm{:#.{vlos_err_prec}g}$'
veff_text_fmt = r'$v_{{\rm eff}} : {:#.{veff_data_prec}g}$' #\pm{:#.{veff_err_prec}g}$'
vlos_data_prec = np.ceil(np.log10(np.abs(vlos_plot))).astype("int") - \
np.ceil(np.log10(vlos_err_plot)).astype("int") + params_prec["wvl"]
ax_spec.text(0.03,0.87,vlos_text_fmt.format(num2tex(vlos_plot),
vlos_data_prec = 2),ha = 'left',va = 'center',
color = 'black',fontsize = 20,linespacing=1.5,transform=ax_spec.transAxes)
veff_data_prec = np.ceil(np.log10(np.abs(veff_plot))).astype("int") - \
np.ceil(np.log10(veff_err_plot)).astype("int") + params_prec["fwhm"]
ax_spec.text(0.03,0.67,veff_text_fmt.format(num2tex(veff_plot),
veff_data_prec = 3),ha = 'left',va = 'center',
color = 'black',fontsize = 20,linespacing=1.5,transform=ax_spec.transAxes)
# ax_spec.set_title(r"Fe \textsc{{xiv}} added-up".format(),fontsize=20)
