diff --git a/src/idl/dustem_fillup_systvar_from_fits.pro b/src/idl/dustem_fillup_systvar_from_fits.pro index 72865ae..b6d4097 100644 --- a/src/idl/dustem_fillup_systvar_from_fits.pro +++ b/src/idl/dustem_fillup_systvar_from_fits.pro @@ -1,75 +1,151 @@ PRO dustem_fillup_systvar_from_fits,syst_var,str_input_SED,used_pol - - Nsed=n_elements(str_input_SED) + + ;IC: some way for the procedure to be aware of emission and extinction + test_m = fix(total(tag_names (str_input_sed) EQ 'STOKESI')) ;a proper boolean + + Nsed=n_elements(str_input_SED) sed=dustem_initialize_internal_sed(Nsed) - - (syst_var).sed=ptr_new(sed) - - (*syst_var.sed).instru_names=dustem_filter2instru(str_input_SED.filter) - (*syst_var.sed).filt_names=strtrim(str_input_SED.filter,2) - (*syst_var.sed).wav=str_input_SED.wavelength - (*syst_var.sed).values=str_input_SED.STOKESI - (*syst_var.sed).sigma=la_power(str_input_SED.varianceII,0.5) - + + IF test_m THEN BEGIN + (syst_var).sed=ptr_new(sed) + (*syst_var.sed).instru_names=dustem_filter2instru(str_input_SED.filter) + (*syst_var.sed).filt_names=strtrim(str_input_SED.filter,2) + (*syst_var.sed).wav=str_input_SED.wavelength + (*syst_var.sed).values=str_input_SED.STOKESI + (*syst_var.sed).sigma=la_power(str_input_SED.varianceII,0.5) + ENDIF ELSE BEGIN + (syst_var).ext=ptr_new(sed) + (*syst_var.ext).instru_names=dustem_filter2instru(str_input_SED.filter) + (*syst_var.ext).filt_names=strtrim(str_input_SED.filter,2) + (*syst_var.ext).wav=str_input_SED.wavelength + (*syst_var.ext).values=str_input_SED.STOKESI_EXT + (*syst_var.ext).sigma=la_power(str_input_SED.varianceII_EXT,0.5) + ENDELSE + IF used_pol THEN BEGIN ;stop - ;===== note: those should include only points with polarization data, not all sed points - indq=where(str_input_SED.STOKESQ NE la_undef(),Nqsed) - indu=where(str_input_SED.STOKESU NE la_undef(),Nused) - qsed=dustem_initialize_internal_sed(Nqsed) - used=dustem_initialize_internal_sed(Nused) - indpolar=where(str_input_SED.STOKESU NE la_undef() AND str_input_SED.STOKESQ NE la_undef(),Npolar) - - psi_sed=dustem_initialize_internal_sed(Npolar) - smallp_sed=dustem_initialize_internal_sed(Npolar) - largep_sed=dustem_initialize_internal_sed(Npolar) - - syst_var.qsed=ptr_new(qsed) - syst_var.used=ptr_new(used) - (*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED[indq].filter) - (*syst_var.qsed).filt_names=strtrim(str_input_SED[indq].filter,2) - (*syst_var.qsed).wav=str_input_SED[indq].wavelength - (*syst_var.qsed).values=str_input_SED[indq].STOKESQ - (*syst_var.qsed).sigma=la_power(str_input_SED[indq].varianceQQ,0.5) - ;now U - (*syst_var.used).instru_names=dustem_filter2instru(str_input_SED[indu].filter) - (*syst_var.used).filt_names=strtrim(str_input_SED[indu].filter,2) - (*syst_var.used).wav=str_input_SED[indu].wavelength - (*syst_var.used).values=str_input_SED[indu].STOKESU - (*syst_var.used).sigma=la_power(str_input_SED[indu].varianceUU,0.5) - ;=== add up dependent fields - syst_var.psi_em=ptr_new(psi_sed) - syst_var.polfrac=ptr_new(smallp_sed) - syst_var.polsed=ptr_new(largep_sed) - polar_iqu2ippsi,str_input_SED[indpolar].STOKESI,str_input_SED[indpolar].STOKESQ,str_input_SED[indpolar].STOKESU,p,psi,largep=largep - (*syst_var.psi_em).values=psi - (*syst_var.polfrac).values=p - (*syst_var.polsed).values=largep - CIQ=str_input_SED[indpolar].VARIANCEII & CIQ[*]=0. ;because input SEDs don't include covariances, yet - CIU=str_input_SED[indpolar].VARIANCEII & CIU[*]=0. ;because input SEDs don't include covariances, yet - CQU=str_input_SED[indpolar].VARIANCEII & CQU[*]=0. ;because input SEDs don't include covariances, yet - polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI $ - ,str_input_SED[indpolar].STOKESQ $ - ,str_input_SED[indpolar].STOKESU $ - ,str_input_SED[indpolar].VARIANCEII $ - ,str_input_SED[indpolar].VARIANCEQQ $ - ,str_input_SED[indpolar].VARIANCEUU $ - ,CIQ,CIU,CQU,variance_smallp,variance_psi,variance_largep,/lac,snr_largeP=snr_largeP - - (*syst_var.psi_em).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter) - (*syst_var.psi_em).filt_names=strtrim(str_input_SED[indpolar].filter,2) - (*syst_var.psi_em).wav=str_input_SED[indpolar].wavelength - - (*syst_var.polfrac).instru_names = (*syst_var.psi_em).instru_names - (*syst_var.polfrac).filt_names = (*syst_var.psi_em).filt_names - (*syst_var.polfrac).wav = (*syst_var.psi_em).wav - - (*syst_var.polsed).instru_names = (*syst_var.psi_em).instru_names - (*syst_var.polsed).filt_names = (*syst_var.psi_em).filt_names - (*syst_var.polsed).wav = (*syst_var.psi_em).wav + IF test_m THEN BEGIN;IC: quickest way I found to make this work + ;===== note: those should include only points with polarization data, not all sed points + indq=where(str_input_SED.STOKESQ NE la_undef(),Nqsed) + indu=where(str_input_SED.STOKESU NE la_undef(),Nused) + qsed=dustem_initialize_internal_sed(Nqsed) + used=dustem_initialize_internal_sed(Nused) + indpolar=where(str_input_SED.STOKESU NE la_undef() AND str_input_SED.STOKESQ NE la_undef(),Npolar) + + psi_sed=dustem_initialize_internal_sed(Npolar) + smallp_sed=dustem_initialize_internal_sed(Npolar) + largep_sed=dustem_initialize_internal_sed(Npolar) + + syst_var.qsed=ptr_new(qsed) + syst_var.used=ptr_new(used) + (*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED[indq].filter) + (*syst_var.qsed).filt_names=strtrim(str_input_SED[indq].filter,2) + (*syst_var.qsed).wav=str_input_SED[indq].wavelength + (*syst_var.qsed).values=str_input_SED[indq].STOKESQ + (*syst_var.qsed).sigma=la_power(str_input_SED[indq].varianceQQ,0.5) + ;now U + (*syst_var.used).instru_names=dustem_filter2instru(str_input_SED[indu].filter) + (*syst_var.used).filt_names=strtrim(str_input_SED[indu].filter,2) + (*syst_var.used).wav=str_input_SED[indu].wavelength + (*syst_var.used).values=str_input_SED[indu].STOKESU + (*syst_var.used).sigma=la_power(str_input_SED[indu].varianceUU,0.5) + ;=== add up dependent fields + syst_var.psi_em=ptr_new(psi_sed) + syst_var.polfrac=ptr_new(smallp_sed) + syst_var.polsed=ptr_new(largep_sed) + polar_iqu2ippsi,str_input_SED[indpolar].STOKESI,str_input_SED[indpolar].STOKESQ,str_input_SED[indpolar].STOKESU,p,psi,largep=largep + (*syst_var.psi_em).values=psi + (*syst_var.polfrac).values=p + (*syst_var.polsed).values=largep + CIQ=str_input_SED[indpolar].VARIANCEII & CIQ[*]=0. ;because input SEDs don't include covariances, yet + CIU=str_input_SED[indpolar].VARIANCEII & CIU[*]=0. ;because input SEDs don't include covariances, yet + CQU=str_input_SED[indpolar].VARIANCEII & CQU[*]=0. ;because input SEDs don't include covariances, yet + polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI $ + ,str_input_SED[indpolar].STOKESQ $ + ,str_input_SED[indpolar].STOKESU $ + ,str_input_SED[indpolar].VARIANCEII $ + ,str_input_SED[indpolar].VARIANCEQQ $ + ,str_input_SED[indpolar].VARIANCEUU $ + ,CIQ,CIU,CQU,variance_smallp,variance_psi,variance_largep,/lac,snr_largeP=snr_largeP + + (*syst_var.psi_em).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter) + (*syst_var.psi_em).filt_names=strtrim(str_input_SED[indpolar].filter,2) + (*syst_var.psi_em).wav=str_input_SED[indpolar].wavelength + + (*syst_var.polfrac).instru_names = (*syst_var.psi_em).instru_names + (*syst_var.polfrac).filt_names = (*syst_var.psi_em).filt_names + (*syst_var.polfrac).wav = (*syst_var.psi_em).wav + + (*syst_var.polsed).instru_names = (*syst_var.psi_em).instru_names + (*syst_var.polsed).filt_names = (*syst_var.psi_em).filt_names + (*syst_var.polsed).wav = (*syst_var.psi_em).wav + + (*syst_var.psi_em).sigma=la_power(variance_psi,0.5) + (*syst_var.polfrac).sigma=la_power(variance_smallp,0.5) + (*syst_var.polsed).sigma=la_power(variance_largep,0.5) + ENDIF ELSE BEGIN - (*syst_var.psi_em).sigma=la_power(variance_psi,0.5) - (*syst_var.polfrac).sigma=la_power(variance_smallp,0.5) - (*syst_var.polsed).sigma=la_power(variance_largep,0.5) + ;===== note: those should include only points with polarization data, not all sed points + indq=where(str_input_SED.STOKESQ_EXT NE la_undef(),Nqsed) + indu=where(str_input_SED.STOKESU_EXT NE la_undef(),Nused) + qsed=dustem_initialize_internal_sed(Nqsed) + used=dustem_initialize_internal_sed(Nused) + indpolar=where(str_input_SED.STOKESU_EXT NE la_undef() AND str_input_SED.STOKESQ_EXT NE la_undef(),Npolar) + + psi_sed=dustem_initialize_internal_sed(Npolar) + smallp_sed=dustem_initialize_internal_sed(Npolar) + largep_sed=dustem_initialize_internal_sed(Npolar) + + syst_var.qext=ptr_new(qsed) + syst_var.uext=ptr_new(used) + (*syst_var.qext).instru_names=dustem_filter2instru(str_input_SED[indq].filter) + (*syst_var.qext).filt_names=strtrim(str_input_SED[indq].filter,2) + (*syst_var.qext).wav=str_input_SED[indq].wavelength + (*syst_var.qext).values=str_input_SED[indq].STOKESQ_EXT + (*syst_var.qext).sigma=la_power(str_input_SED[indq].varianceQQ_EXT,0.5) + ;now U + (*syst_var.uext).instru_names=dustem_filter2instru(str_input_SED[indu].filter) + (*syst_var.uext).filt_names=strtrim(str_input_SED[indu].filter,2) + (*syst_var.uext).wav=str_input_SED[indu].wavelength + (*syst_var.uext).values=str_input_SED[indu].STOKESU_EXT + (*syst_var.uext).sigma=la_power(str_input_SED[indu].varianceUU_EXT,0.5) + ;=== add up dependent fields + syst_var.psi_ext=ptr_new(psi_sed) + syst_var.fpolext=ptr_new(smallp_sed) + syst_var.polext=ptr_new(largep_sed) + polar_iqu2ippsi,str_input_SED[indpolar].STOKESI_EXT,str_input_SED[indpolar].STOKESQ_EXT,str_input_SED[indpolar].STOKESU_EXT,p,psi,largep=largep + (*syst_var.psi_ext).values=psi + (*syst_var.fpolext).values=p + (*syst_var.polext).values=largep + CIQ_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CIQ_EXT[*]=0. ;because input SEDs don't include covariances, yet + CIU_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CIU_EXT[*]=0. ;because input SEDs don't include covariances, yet + CQU_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CQU_EXT[*]=0. ;because input SEDs don't include covariances, yet + polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI_EXT $ + ,str_input_SED[indpolar].STOKESQ_EXT $ + ,str_input_SED[indpolar].STOKESU_EXT $ + ,str_input_SED[indpolar].VARIANCEII_EXT $ + ,str_input_SED[indpolar].VARIANCEQQ_EXT $ + ,str_input_SED[indpolar].VARIANCEUU_EXT $ + ,CIQ_EXT,CIU_EXT,CQU_EXT,variance_smallp_EXT,variance_psi_EXT,variance_largep_EXT,/lac,snr_largeP=snr_largeP_EXT + + (*syst_var.psi_ext).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter) + (*syst_var.psi_ext).filt_names=strtrim(str_input_SED[indpolar].filter,2) + (*syst_var.psi_ext).wav=str_input_SED[indpolar].wavelength + + (*syst_var.fpolext).instru_names = (*syst_var.psi_ext).instru_names + (*syst_var.fpolext).filt_names = (*syst_var.psi_ext).filt_names + (*syst_var.fpolext).wav = (*syst_var.psi_ext).wav + + (*syst_var.polext).instru_names = (*syst_var.psi_ext).instru_names + (*syst_var.polext).filt_names = (*syst_var.psi_ext).filt_names + (*syst_var.polext).wav = (*syst_var.psi_ext).wav + + (*syst_var.psi_ext).sigma=la_power(variance_psi_EXT,0.5) + (*syst_var.fpolext).sigma=la_power(variance_smallp_EXT,0.5) + (*syst_var.polext).sigma=la_power(variance_largep_EXT,0.5) + + ENDELSE + + ENDIF END \ No newline at end of file diff --git a/src/idl/dustem_make_fits_predicted_ext.pro b/src/idl/dustem_make_fits_predicted_ext.pro index 07a202d..29c4a1f 100644 --- a/src/idl/dustem_make_fits_predicted_ext.pro +++ b/src/idl/dustem_make_fits_predicted_ext.pro @@ -51,16 +51,36 @@ ENDIF ;===== define structures containing results -one_str_input_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)} -one_str_predicted_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)} +one_str_input_EXT={FILTER:'',Wavelength:la_undef(4),StokesI_EXT:la_undef(5),StokesQ_EXT:la_undef(5),StokesU_EXT:la_undef(5),varianceII_EXT:la_undef(5),varianceQQ_EXT:la_undef(5),varianceUU_EXT:la_undef(5)} +one_str_predicted_EXT={FILTER:'',Wavelength:la_undef(4),StokesI_EXT:la_undef(5),StokesQ_EXT:la_undef(5),StokesU_EXT:la_undef(5)} Next=n_elements((*syst_var.ext).filt_names) str_input_EXT=replicate(one_str_input_EXT,Next) str_input_EXT.FILTER=(*syst_var.ext).filt_names str_input_EXT.Wavelength=(*syst_var.ext).wav -str_input_EXT.I=(*syst_var.ext).values -str_input_EXT.varianceII=la_power((*syst_var.ext).sigma,2.) - +str_input_EXT.StokesI_EXT=(*syst_var.ext).values +str_input_EXT.varianceII_EXT=la_power((*syst_var.ext).sigma,2.) +stop +;IC; if I didn't misunderstand the str_input_EXT output lacks this information +IF !run_pol THEN BEGIN + qext=aos2soa(*syst_var.qext) + FOR i=0L,n_elements(qext.wav)-1 DO BEGIN + ind=where(str_input_EXT.FILTER EQ qext[i].filt_names and str_input_EXT.wavelength EQ qext[i].wav,count) + IF count NE 0 THEN BEGIN + str_input_EXT[ind[0]].STOKESQ_EXT=qext[i].values + str_input_EXT[ind[0]].varianceQQ_EXT=la_power(qext[i].sigma,2.) + ENDIF + ENDFOR + uext=aos2soa(*syst_var.uext) + FOR i=0L,n_elements(uext.wav)-1 DO BEGIN + ind=where(str_input_EXT.FILTER EQ uext[i].filt_names and str_input_EXT.wavelength EQ uext[i].wav,count) + IF count NE 0 THEN BEGIN + str_input_EXT[ind[0]].STOKESU_EXT=uext[i].values + str_input_EXT[ind[0]].varianceUU_EXT=la_power(uext[i].sigma,2.) + ENDIF + ENDFOR +ENDIF +;IC: If we make an analogy with emission, why aren't there any predicted variances? ;N_predicted_ext=n_elements(dustem_predicted_polext[0]) N_predicted_ext=n_elements(dustem_predicted_ext) str_predicted_EXT=replicate(one_str_predicted_EXT,N_predicted_EXT) @@ -71,15 +91,15 @@ IF N_predicted_ext NE Next THEN BEGIN ;This is for cases when the extinction da ENDIF ELSE BEGIN str_predicted_EXT.FILTER=(*syst_var.ext).filt_names ;taken from input SED str_predicted_EXT.Wavelength=(*syst_var.ext).wav ;taken from input SED - str_predicted_EXT.I=dustem_predicted_ext + str_predicted_EXT.StokesI_EXT=dustem_predicted_ext ENDELSE -;stop + IF !run_pol THEN BEGIN - str_predicted_EXT.Q=interpol(dustem_predicted_Qext,(*syst_var.qext).wav,(*syst_var.ext).wav) - str_predicted_EXT.U=interpol(dustem_predicted_Uext,(*syst_var.uext).wav,(*syst_var.ext).wav) + str_predicted_EXT.StokesQ_EXT=interpol(dustem_predicted_Qext,(*syst_var.qext).wav,(*syst_var.ext).wav) + str_predicted_EXT.StokesU_EXT=interpol(dustem_predicted_Uext,(*syst_var.uext).wav,(*syst_var.ext).wav) ENDIF ELSE BEGIN - str_predicted_EXT.Q=fully_undefined_ext.Q - str_predicted_EXT.U=fully_undefined_ext.U + str_predicted_EXT.StokesQ_EXT=fully_undefined_ext.StokesQ_EXT + str_predicted_EXT.StokesU_EXT=fully_undefined_ext.StokesU_EXT ENDELSE RETURN,str_predicted_EXT diff --git a/src/idl/dustem_plot_dataset.pro b/src/idl/dustem_plot_dataset.pro index 2289e4a..3354b8e 100644 --- a/src/idl/dustem_plot_dataset.pro +++ b/src/idl/dustem_plot_dataset.pro @@ -363,7 +363,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang if keyword_set(refresh) then begin ;The data points in the plot that are being refreshed ;making sure we're in the normalized plot so that overplotting happens here - cgplot,1/vectw,vectx,/xlog,/nodata,/ys,xs=1,pos=position[1],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,charsize=1.0 + cgplot,1/vectw,vectx,/xlog,/nodata,ys=8,xs=8,pos=position[1],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,charsize=1.0 IF ct_spec NE 0 THEN BEGIN xx=((*(*!dustem_data).ext).wav)[idx_spec] @@ -410,7 +410,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang ;making sure we're in the normal plot so that overplotting happens here - cgplot,1/vectw,vectx,/xlog,/ylog,/nodata,/ys,xs=1,pos=position[0],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=yr,yticks=2,ymino=2,xticklen=0.1,charsize=1.0 + cgplot,1/vectw,vectx,/xlog,/ylog,/nodata,ys=8,xs=8,pos=position[0],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=yr,yticks=2,ymino=2,xticklen=0.1,charsize=1.0 ;Plotting of the spectra of the dust species @@ -463,7 +463,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang if ct_spec ne 0 then begin ;Plotting of spectrum data points (to be fitted) - cgplot,1/((*(*!dustem_data).ext).wav)(idx_spec),((*(*!dustem_data).ext).values)(idx_spec),/ylog,/xlog,/ys,xs=9,pos=position[0],noerase=1,xtit=' ',charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,yr=yr,psym=16,syms=0.8,ytickformat='dstmwrp_exp' + cgplot,1/((*(*!dustem_data).ext).wav)(idx_spec),((*(*!dustem_data).ext).values)(idx_spec),/ylog,/xlog,ys=8,xs=8,pos=position[0],noerase=1,xtit=' ',charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,yr=yr,psym=16,syms=0.8,ytickformat='dstmwrp_exp' ;Plotting of the spectrum error points rms=2.*((*(*!dustem_data).ext).sigma)(idx_spec)/2. @@ -475,7 +475,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang ;Plotting of filter data points (to be fitted) + testing for prior plotting if ct_spec ne 0 then cgoplot,1/((*(*!dustem_data).ext).wav)(idx_filt),((*(*!dustem_data).ext).values)(idx_filt),charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position[0],/ys,xs=9,noerase=1,xtickformat='(A1)',xr=xr,yr=yr,ytickformat='dstmwrp_exp',/ylog,/xlog else $ - cgplot,1/((*(*!dustem_data).ext).wav)(idx_filt),((*(*!dustem_data).ext).values)(idx_filt),charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position[0],/ys,xs=9,noerase=1,xtickformat='(A1)',xtit=' ',xr=xr,yr=yr,ytickformat='dstmwrp_exp',/ylog,/xlog + cgplot,1/((*(*!dustem_data).ext).wav)(idx_filt),((*(*!dustem_data).ext).values)(idx_filt),charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position[0],ys=8,xs=8,noerase=1,xtickformat='(A1)',xtit=' ',xr=xr,yr=yr,ytickformat='dstmwrp_exp',/ylog,/xlog ;Plotting of the filter error points @@ -502,7 +502,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang ;stop if ct_spec_hdn ne 0 then begin ;Plotting of hidden spectrum data points - cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_spec_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_spec_hdn),pos=position[0],/ylog,/xlog,/ys,xs=9,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8 + cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_spec_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_spec_hdn),pos=position[0],/ylog,/xlog,ys=8,xs=8,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8 ;Plotting of hidden spectrum error points rms=2.*(((*(*!dustem_show).ext).sigma)[idx_rmv_sed])(idx_spec_hdn)/2. @@ -511,7 +511,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang if ct_filt_hdn then begin ;Plotting of hidden filter data points - cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_filt_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_filt_hdn),pos=position[0],/ylog,/xlog,/ys,xs=9,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8 + cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_filt_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_filt_hdn),pos=position[0],/ylog,/xlog,ys=8,xs=8,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8 ;Plotting of hidden filter error point rms=2.*(((*(*!dustem_show).ext).sigma)[idx_rmv_sed])(idx_filt_hdn)/2. @@ -1635,7 +1635,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang rms=2.*((*(*!dustem_data).qsed).sigma)(idx_filt)/2. - if ct_spec ne 0 then dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,xs=1,xr=xr,noerase=1,/positive_only,xtickformat='(A1)',xtit='',/xlog,/ylog,yr=yr,ytickformat='(A1)';,ytickformat='(A1)' + if ct_spec ne 0 then dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,xs=1,xr=xr,noerase=1,/positive_only,xtickformat='(A1)',xtit='',/xlog,/ylog,yr=yr,/overplot,ytickformat='(A1)';,ytickformat='(A1)' if ct_spec eq 0 then dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,xs=1,xr=xr,noerase=1,/positive_only,xtickformat='(A1)',xtit='',/xlog,/ylog,yr=yr;,ytickformat='(A1)' dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],color='Dodger Blue',/positive_only,noerase=1,/overplot,rms=rms,xtit='',charsize=1.15,xtickformat='(A1)',xs=1,psym=16,syms=0.8,xr=xr,/xlog,/ylog,yr=yr,ytickformat='(A1)' @@ -2218,7 +2218,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang if keyword_set(refresh) then begin ;The data points in the plot that are being refreshed ;Making sure we're placed in the right plot to refresh. Here normalized plot for QEXT. - cgplot,1/vectw,vectx,/nodata,/xlog,/ys,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8 + cgplot,1/vectw,vectx,/nodata,/xlog,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8 IF ct_spec NE 0 THEN BEGIN @@ -2255,7 +2255,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang xtit=textoidl('1/\lambda (\mum^{-1})') - cgplot,1/vectw,vectx,/xlog,/ys,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 + cgplot,1/vectw,vectx,/xlog,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 xyouts,pospltxt_n[0],pospltxt_n[1],textoidl('norm'),color=0,/normal,charsize=1.1 endelse @@ -2574,7 +2574,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang if keyword_set(refresh) then begin ;The data points in the plot that are being refreshed ;Making sure we're placed in the right plot to refresh. Here normalized plot for QEXT. - cgplot,1/vectw,vectx,/nodata,/xlog,/ys,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8 + cgplot,1/vectw,vectx,/nodata,/xlog,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8 IF ct_spec NE 0 THEN BEGIN @@ -2607,7 +2607,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang xtit=textoidl('1/\lambda (\mum^{-1})') ;if !run_pol then xtit = '' - cgplot,1/vectw,vectx,/xlog,/ys,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 + cgplot,1/vectw,vectx,/xlog,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 xyouts,pospltxt_n[0],pospltxt_n[1],textoidl('norm'),color=0,/normal,charsize=1.1 endelse diff --git a/src/idl/dustem_read_fits_table.pro b/src/idl/dustem_read_fits_table.pro index f5fcfdf..1c6406b 100644 --- a/src/idl/dustem_read_fits_table.pro +++ b/src/idl/dustem_read_fits_table.pro @@ -224,7 +224,7 @@ IF unit_dustem_predicted_extinction_sed NE -1 THEN BEGIN ENDIF ELSE BEGIN ENDELSE IF unit_dustem_shown_extinction_sed NE -1 THEN BEGIN - str_shown_EXT=mrdfits(file,unit_shown_extinction_sed,header_shown_ext) + str_shown_EXT=mrdfits(file,unit_dustem_shown_extinction_sed,header_shown_ext) ENDIF ELSE BEGIN ENDELSE ;=== read plugin informations and plugins headers @@ -252,7 +252,8 @@ ENDIF ELSE BEGIN dustem_init,model=0 ENDELSE -defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH'))) +IF unit_input_emission_sed NE -1 THEN defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH'))) +IF unit_input_extinction_sed NE -1 THEN defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_ext,'NH'))) !dustem_model=used_model IF !dustem_version.version NE used_version THEN BEGIN @@ -308,7 +309,6 @@ IF Nplugin NE 0 THEN BEGIN (*!dustem_plugin).(i).PARAMTAG=ptr_new(paramtags) ENDFOR ENDIF -;stop ;==== This is to create (*!dustem_data).sed in case it does not exist already ; it doesn't exist at all times because dustem_init has been run @@ -319,7 +319,6 @@ IF unit_input_emission_sed NE -1 THEN BEGIN ENDIF ;help,str_shown_SED -;stop ;==== This is to create (*!dustem_show).sed in case it does not exist already IF unit_dustem_shown_sed NE -1 THEN BEGIN diff --git a/src/idl/dustem_write_fits_table.pro b/src/idl/dustem_write_fits_table.pro index 26c39f4..729ec44 100755 --- a/src/idl/dustem_write_fits_table.pro +++ b/src/idl/dustem_write_fits_table.pro @@ -187,7 +187,7 @@ IF Nplugins NE 0 THEN BEGIN str_plugins[i]=plugin_tags[i] unit=unit+1 ENDFOR -ENDIF ELSE BEGIN +ENDIF ELSE BEGIN ;IC: This should be commented right? These quantities do not exist. unit_plugins[i]=-1 str_plugins[i]='' ENDELSE @@ -290,29 +290,51 @@ FOR i=0L,Nparams-1 DO BEGIN ENDFOR sxaddpar,hdr,'COMMENT','===================================================',after='PARF'+strtrim(i,2) -;stop -;===== complement extension header for input SED -sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT' -sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name' -sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]' -sxaddpar,header_input_sed,'TTYPE3','I','Observed Total Stokes I Intensity [??]' -sxaddpar,header_input_sed,'TTYPE4','Q','Observed Stokes Q Intensity [??]' -sxaddpar,header_input_sed,'TTYPE5','U','Observed Stokes U Intensity [??]' -sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]' -sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]' -sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]' -sxaddpar,header_input_sed,'NH',*!dustem_hcd,'Column density [H/cm2]' +IF unit_input_sed NE -1 THEN BEGIN -;stop + ;===== complement extension header for input SED + sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT' + sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name' + sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]' + sxaddpar,header_input_sed,'TTYPE3','STOKESI','Observed Total Stokes I Intensity [??]' + sxaddpar,header_input_sed,'TTYPE4','STOKESQ','Observed Stokes Q Intensity [??]' + sxaddpar,header_input_sed,'TTYPE5','STOKESU','Observed Stokes U Intensity [??]' + sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]' + sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]' + sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]' + sxaddpar,header_input_sed,'NH',*!dustem_hcd,'Column density [H/cm2]' + + ;===== complement extension header for predicted SED + sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP' + sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name' + sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]' + sxaddpar,header_predicted_SED,'TTYPE3','STOKESI','Observed Total Stokes I Intensity [??]' + sxaddpar,header_predicted_SED,'TTYPE4','STOKESQ','Observed Stokes Q Intensity [??]' + sxaddpar,header_predicted_SED,'TTYPE5','STOKESU','Observed Stokes U Intensity [??]' +ENDIF ELSE BEGIN;IC: hesitated between this and 'IF unit_input_ext NE -1 THEN BEGIN'. This is under the assumtion that we always fit something. + + ;===== complement extension header for input extinction SED + sxaddpar,header_input_ext,'TITLE','OBSERVED INPUT EXTINCTION SED TO DUSTEMWRAP',before='COMMENT' + sxaddpar,header_input_ext,'TTYPE1','FILTER','Dustemwrap Filter name' + sxaddpar,header_input_ext,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]' + sxaddpar,header_input_ext,'TTYPE3','STOKESI_EXT','Observed Total Stokes I optical depth []' + sxaddpar,header_input_ext,'TTYPE4','STOKESQ_EXT','Observed Stokes Q optical depth []' + sxaddpar,header_input_ext,'TTYPE5','STOKESU_EXT','Observed Stokes U optical depth []' + sxaddpar,header_input_ext,'TTYPE6','VARIANCEII_EXT','Observed Total Stokes I optical depth Variance [??]' + sxaddpar,header_input_ext,'TTYPE7','VARIANCEQQ_EXT','Observed Stokes Q optical depth Variance [??]' + sxaddpar,header_input_ext,'TTYPE8','VARIANCEUU_EXT','Observed Stokes U optical depth Variance [??]' + sxaddpar,header_input_ext,'NH',*!dustem_hcd,'Column density [H/cm2]' + + ;===== complement extension header for predicted extinction SED + sxaddpar,header_predicted_EXT,'TITLE','OUTPUT BEST FIT EXTINCTION SED BY DUSTEMWRAP' + sxaddpar,header_predicted_EXT,'TTYPE1','FILTER','Dustemwrap Filter name' + sxaddpar,header_predicted_EXT,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]' + sxaddpar,header_predicted_EXT,'TTYPE3','STOKESI_EXT','Observed Total Stokes I optical depth []' + sxaddpar,header_predicted_EXT,'TTYPE4','STOKESQ_EXT','Observed Stokes Q optical depth []' + sxaddpar,header_predicted_EXT,'TTYPE5','STOKESU_EXT','Observed Stokes U optical depth []' -;===== complement extension header for predicted SED -sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP' -sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name' -sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]' -sxaddpar,header_predicted_SED,'TTYPE3','I','Observed Total Stokes I Intensity [??]' -sxaddpar,header_predicted_SED,'TTYPE4','Q','Observed Stokes Q Intensity [??]' -sxaddpar,header_predicted_SED,'TTYPE5','U','Observed Stokes U Intensity [??]' +ENDELSE ;===== complement extension header for predicted spectrum total sxaddpar,header_predicted_spectra,'TITLE','OUTPUT BEST FIT PER-GRAIN AND TOTAL EMISSION SPECTRA BY DUSTEMWRAP' -- libgit2 0.21.2