diff --git a/src/idl/dustem_brute_force_fit.pro b/src/idl/dustem_brute_force_fit.pro index fc9c50c..0675d34 100644 --- a/src/idl/dustem_brute_force_fit.pro +++ b/src/idl/dustem_brute_force_fit.pro @@ -1,7 +1,7 @@ FUNCTION dustem_brute_force_fit,sed,table_name,filters,normalize=normalize,fact=fact,chi2=chi2,rchi2=rchi2,show_sed=show_sed ;stop -;sed must be normalized to NH=1e21 H/cm2 +;sed must be normalized to NH=1e20 H/cm2 defsysv,'!dustem_grid',0,exist=exist IF exist EQ 0 THEN BEGIN diff --git a/src/idl/dustem_compute_sed.pro b/src/idl/dustem_compute_sed.pro index e03d2dd..38a1de6 100755 --- a/src/idl/dustem_compute_sed.pro +++ b/src/idl/dustem_compute_sed.pro @@ -1,6 +1,7 @@ FUNCTION dustem_compute_sed,p_dim,$ st=st,$ sed_spec=sed_spec, $ + use_previous_fortran=use_previous_fortran, $ help=help ;+ @@ -20,7 +21,7 @@ FUNCTION dustem_compute_sed,p_dim,$ ; p_dim = parameter values ; ; OPTIONAL INPUT PARAMETERS: -; st = Dustem output structure +; st = Dustem Fortran output structure ; ; OUTPUTS: ; sed = computed SED for filters in !dustem_data @@ -30,6 +31,7 @@ FUNCTION dustem_compute_sed,p_dim,$ ; sed_spec = dustem emission spectrum ; ; ACCEPTED KEY-WORDS: +; use_previous_fortran = if set, uses the output of the previous fortran run, and does not run the fortran. ; help = If set, print this help ; ; COMMON BLOCKS: @@ -59,7 +61,7 @@ IF keyword_set(help) THEN BEGIN ENDIF IF not keyword_set(st) THEN BEGIN - dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st,use_previous_fortran=use_previous_fortran ENDIF ; Convert into MJy/sr/1d20 @@ -67,20 +69,22 @@ fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is SED_spec = st.sed.em_tot * fact ;ADDING PLUGIN TO SPECTRUM---------------- -scopes=tag_names((*!dustem_plugin)) +Nplugins=n_elements(*!dustem_plugin) +scopes=(*!dustem_plugin).scope IF scopes[0] NE 'NONE' THEN BEGIN -;IF ptr_valid(!dustem_plugin) THEN BEGIN - for i=0L,n_tags(*!dustem_plugin)-1 do begin - if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_SED') then SED_spec+=(*(*!dustem_plugin).(i).spec)[*,0] - endfor + FOR i=0L,Nplugins-1 DO BEGIN + IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) eq 'ADD_SED') THEN BEGIN + SED_spec+=(*(*!dustem_plugin)[i].spec)[*,0] + ENDIF + ENDFOR ENDIF ;------------------------------------------ -if not isarray(SED_spec) THEN stop +IF not isarray(SED_spec) THEN stop ;COMPUTE THE MODEL SED - +;JPB: problem there: we should be able to compute an SED even when there is no data in !dustem_data dustem_sed = (*(*!dustem_data).sed).values * 0. IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN diff --git a/src/idl/dustem_fit_intensity_example.pro b/src/idl/dustem_fit_intensity_example.pro index 1be4180..59ee653 100644 --- a/src/idl/dustem_fit_intensity_example.pro +++ b/src/idl/dustem_fit_intensity_example.pro @@ -257,7 +257,8 @@ if keyword_set(wait) then begin end ;== INITIALISE DUSTEM -dustem_init,model=use_model,polarization=use_polarization +dustem_init,model=use_model,polarization=use_polarization,/show_plot +;stop !dustem_nocatch=1 !dustem_verbose=use_verbose IF keyword_set(noobj) THEN !dustem_noobj=1 @@ -303,7 +304,7 @@ tit='FIT INTENSITY EXAMPLE' ; plot title ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title xtit=textoidl('\lambda (\mum)') ; x-axis title - +;stop ;=== RUN THE FIT t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $ diff --git a/src/idl/dustem_fit_intensity_phangs_example.pro b/src/idl/dustem_fit_intensity_phangs_example.pro index 90a5f22..8bc7d2d 100644 --- a/src/idl/dustem_fit_intensity_phangs_example.pro +++ b/src/idl/dustem_fit_intensity_phangs_example.pro @@ -10,7 +10,7 @@ PRO dustem_fit_intensity_phangs_example,model=model $ ;+ ; NAME: -; dustem_fit_intensity_example +; dustem_fit_intensity_phangs_example ; ; PURPOSE: ; This routine is an example of how to fit an observational SED @@ -169,7 +169,9 @@ if keyword_set(wait) then begin end ;== INITIALISE DUSTEM -dustem_init,model=use_model,polarization=use_polarization +;!dustem_show_plot=1 +dustem_init,model=use_model,polarization=use_polarization,/show_plot +;stop !dustem_nocatch=1 !dustem_verbose=use_verbose IF keyword_set(noobj) THEN !dustem_noobj=1 diff --git a/src/idl/dustem_plugin_phangs_stellar_continuum.pro b/src/idl/dustem_plugin_phangs_stellar_continuum.pro index d2c41c2..2d62321 100644 --- a/src/idl/dustem_plugin_phangs_stellar_continuum.pro +++ b/src/idl/dustem_plugin_phangs_stellar_continuum.pro @@ -31,6 +31,7 @@ FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $ ; paramtag = plugin parameter names as strings ; ACCEPTED KEY-WORDS: ; help = if set, print this help +; method = SSP used. can be EMILES or CB19. default='EMILES' ; COMMON BLOCKS: ; None ; SIDE EFFECTS: @@ -60,31 +61,25 @@ IF keyword_set(scope) THEN BEGIN goto, the_scope ENDIF -;IF keyword_set(paramtag) THEN BEGIN -; out=0 -; goto, the_paramtag -;ENDIF - -;stop - -;default values of input parameters +;use_method='EMILES' +use_method='CB19' age_values=[0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5] metalicity_values=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397100] ;read template info -info_st=dustem_read_emiles_stellar_templates(age_values,metalicity_values,/info_only) -;stop -age_values_emiles=info_st.age -order=sort(age_values_emiles) -age_values_emiles=age_values_emiles[order] -un=uniq(age_values_emiles) -age_values_emiles=age_values_emiles[un] -met_values_emiles=info_st.z -order=sort(met_values_emiles) -met_values_emiles=met_values_emiles[order] -un=uniq(met_values_emiles) -met_values_emiles=met_values_emiles[un] +;==== This is actually not used +;info_st=dustem_read_emiles_stellar_templates(age_values,metalicity_values,/info_only) +;age_values_emiles=info_st.age +;order=sort(age_values_emiles) +;age_values_emiles=age_values_emiles[order] +;un=uniq(age_values_emiles) +;age_values_emiles=age_values_emiles[un] +;met_values_emiles=info_st.z +;order=sort(met_values_emiles) +;met_values_emiles=met_values_emiles[order] +;un=uniq(met_values_emiles) +;met_values_emiles=met_values_emiles[un] ;==== This does not work because the info files has many more entries than there are available ssps. ;age_values=age_values_emiles @@ -132,8 +127,18 @@ IF exist EQ 0 THEN BEGIN defsysv,'!dustem_plugin_phangs_stellar_continuum',st ENDIF IF not ptr_valid(!dustem_plugin_phangs_stellar_continuum.ssps) THEN BEGIN - message,'Initializing SSP templates',/continue - template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars) + message,'Initializing SSP templates',/info + CASE use_method of + 'EMILES': BEGIN + template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars) + END + 'CB19': BEGIN + ;== this is just to get the Mstars + bidon=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars) + ;== This is to read in the CB19 templates + template_flux=dustem_read_cb19_stellar_templates(age_values=age_values,metalicity_values=metalicity_values,template_wav=template_wav,Mstars=Mstars) + END + ENDCASE !dustem_plugin_phangs_stellar_continuum.ssps=ptr_new(template_flux) !dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav) !dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars) @@ -157,7 +162,7 @@ FOR i=2L,Nparam-1 DO BEGIN IF count NE 0 THEN BEGIN ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-2] output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind]) - message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/continue + message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/info ENDIF ENDIF ENDFOR @@ -178,12 +183,12 @@ fact3=1./NHref ;[cm2/H] to go from ergs/s/cm2 to ergs/s/H fact4=1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/lambir)*1.e20/1.e7 ;This goes from ergs/s/H to MJy/sr (see dustem_plot_dataset) fact=fact2*fact3*fact4 -print,'flux(1mic)=',interpol(flux,lambir,1.),' ergs/s/cm2/AA' -print,'fact2=',interpol(fact2,lambir,1.) -print,'fact3=',fact3 -print,'fact4=',interpol(fact4,lambir,1.) -print,'fact=',interpol(fact,lambir,1.) -print,'flux(1mic)=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' MJy/sr' +;print,'flux(1mic)=',interpol(flux,lambir,1.),' ergs/s/cm2/AA' +;print,'fact2=',interpol(fact2,lambir,1.) +;print,'fact3=',fact3 +;print,'fact4=',interpol(fact4,lambir,1.) +;print,'fact=',interpol(fact,lambir,1.) +;print,'flux(1mic)=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' MJy/sr' ;goto,no_unred @@ -213,7 +218,7 @@ output[*,0]=output[*,0]*fact output[*,1]=output[*,1]*fact output[*,2]=output[*,2]*fact -print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr' +;print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr' ;stop the_scope: diff --git a/src/idl/dustem_read_emiles_stellar_templates.pro b/src/idl/dustem_read_emiles_stellar_templates.pro index c3f5fc5..af16ee5 100644 --- a/src/idl/dustem_read_emiles_stellar_templates.pro +++ b/src/idl/dustem_read_emiles_stellar_templates.pro @@ -39,7 +39,7 @@ FUNCTION dustem_read_emiles_stellar_templates,age_values, $ ; PROCEDURES AND SUBROUTINES USED: ; ; EXAMPLES -; +; st=dustem_read_emiles_stellar_templates(age_values,Z_values) ; MODIFICATION HISTORY: ; Written by JPB June-2023 ; Evolution details on the DustEMWrap gitlab. @@ -91,8 +91,8 @@ metalicity_values_4str=string(metalicity_values_4str,format='(F5.2)') FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],' ','p') FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],'-','m') -print,age_values -print,metalicity_values +;print,age_values +;print,metalicity_values ;stop ;Ntemplates=n_elements(paramvalues)-1 @@ -119,11 +119,11 @@ FOR i=1L,Ntemplates DO BEGIN met=metalicity_values[imet] met_str=metalicity_values_4str[imet] age_str=age_values_4str[iage] - message,'reading template for Age '+age_str+' Met '+met_str,/continue + message,'reading template for Age '+age_str+' Met '+met_str,/info file=dir_templates+file_start+'Z'+met_str+'T'+age_str+file_end fi=file_info(file) IF fi.exists THEN BEGIN - message,'Reading SSP '+file,/continue + message,'Reading SSP '+file,/info spec=readfits(file,h,2) templates[ii]=ptr_new(spec) make_axis,h,wavs diff --git a/src/idl/dustem_sed_extractor.pro b/src/idl/dustem_sed_extractor.pro index 568b63c..31e53ee 100644 --- a/src/idl/dustem_sed_extractor.pro +++ b/src/idl/dustem_sed_extractor.pro @@ -1,6 +1,6 @@ -FUNCTION dustem_sed_extractor,maps $ - ,index_roi $ - ,filters, $ +FUNCTION dustem_sed_extractor,maps, $ + index_roi, $ + filters, $ comments=comments, $ maps_order=maps_order, $ method=method, $ @@ -25,7 +25,7 @@ FUNCTION dustem_sed_extractor,maps $ ; sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only]) ; ; INPUTS: -; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets] +; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]. Caution, maps are changed if set_undefined is set ! ; index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest) ; filters : DustEMWrap filters corresponding to each dataset [Ndata_sets] ; @@ -117,9 +117,10 @@ ENDIF ELSE BEGIN ENDIF ENDELSE -used_maps=maps +;used_maps=maps IF keyword_set(set_undefined) THEN BEGIN + used_maps=maps ;===== do a map of saturated stokesI pixels mask_I_undef=maps[*,*,0,0] mask_I_undef[*]=0 @@ -132,8 +133,6 @@ IF keyword_set(set_undefined) THEN BEGIN mask_I_undef[ind]=1 ENDIF ENDFOR - - ;used_maps=maps ;set saturated stokesI pixels to la_undef in indbad=where(mask_I_undef EQ 1,count) IF count NE 0 THEN BEGIN @@ -146,34 +145,36 @@ IF keyword_set(set_undefined) THEN BEGIN ENDFOR ENDFOR ENDIF + maps=used_maps + used_maps=0 ENDIF CASE use_method OF 'MEAN':BEGIN - FOR j=0L,Nfilt-1 DO BEGIN - this_map=used_maps[*,*,indI,j] + FOR j=0L,Nfilt-1 DO BEGIN + this_map=maps[*,*,indI,j] sed[j].stokesI=la_mean(this_map[index_roi]) - this_map=used_maps[*,*,indII,j] + this_map=maps[*,*,indII,j] sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi) IF not keyword_set(total_intensity_only) THEN BEGIN - this_map=used_maps[*,*,indQ,j] + this_map=maps[*,*,indQ,j] sed[j].stokesQ=la_mean(this_map[index_roi]) - this_map=used_maps[*,*,indU,j] + this_map=maps[*,*,indU,j] sed[j].stokesU=la_mean(this_map[index_roi]) - this_map=used_maps[*,*,indQQ,j] + this_map=maps[*,*,indQQ,j] sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi) - this_map=used_maps[*,*,indUU,j] + this_map=maps[*,*,indUU,j] sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi) - this_map=used_maps[*,*,indIQ,j] + this_map=maps[*,*,indIQ,j] sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi) - this_map=used_maps[*,*,indIU,j] + this_map=maps[*,*,indIU,j] sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi) - this_map=used_maps[*,*,indQU,j] + this_map=maps[*,*,indQU,j] sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi) ENDIF - ENDFOR - END + ENDFOR + END 'MEAN_ONOFF': BEGIN ;This is mean ON - mean MEAN_ONOFF IF not keyword_set(index_off) THEN BEGIN message,'index_off keyword must be set with method '+use_method,/continue @@ -182,36 +183,36 @@ CASE use_method OF count_off=n_elements(index_off) comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels'] FOR j=0L,Nfilt-1 DO BEGIN - this_map=used_maps[*,*,indI,j] + this_map=maps[*,*,indI,j] sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) - this_map=used_maps[*,*,indII,j] + this_map=maps[*,*,indII,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaII=la_add(sig_on,sig_off) IF not keyword_set(total_intensity_only) THEN BEGIN - this_map=used_maps[*,*,indQ,j] + this_map=maps[*,*,indQ,j] sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) - this_map=used_maps[*,*,indU,j] + this_map=maps[*,*,indU,j] sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) ;variances - this_map=used_maps[*,*,indQQ,j] + this_map=maps[*,*,indQQ,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaQQ=la_add(sig_on,sig_off) - this_map=used_maps[*,*,indUU,j] + this_map=maps[*,*,indUU,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaUU=la_add(sig_on,sig_off) ;co-variances - this_map=used_maps[*,*,indIQ,j] + this_map=maps[*,*,indIQ,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaIQ=la_add(sig_on,sig_off) - this_map=used_maps[*,*,indIU,j] + this_map=maps[*,*,indIU,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaIU=la_add(sig_on,sig_off) - this_map=used_maps[*,*,indQU,j] + this_map=maps[*,*,indQU,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaQU=la_add(sig_on,sig_off) @@ -232,21 +233,21 @@ CASE use_method OF message,'reference filter '+reference_filter+' not found in provided filter list',/continue stop ENDIF - mapI_ref=used_maps[*,*,indI,indf] - mapQ_ref=used_maps[*,*,indQ,indf] - mapU_ref=used_maps[*,*,indU,indf] - mapII_ref=used_maps[*,*,indII,indf] - mapQQ_ref=used_maps[*,*,indQQ,indf] - mapUU_ref=used_maps[*,*,indUU,indf] + mapI_ref=maps[*,*,indI,indf] + mapQ_ref=maps[*,*,indQ,indf] + mapU_ref=maps[*,*,indU,indf] + mapII_ref=maps[*,*,indII,indf] + mapQQ_ref=maps[*,*,indQQ,indf] + mapUU_ref=maps[*,*,indUU,indf] Iref=la_mean(mapI_ref[index_roi]) Qref=la_mean(mapQ_ref[index_roi]) Uref=la_mean(mapU_ref[index_roi]) FOR j=0L,Nfilt-1 DO BEGIN ;==== do the correlation in Stokes I x=mapI_ref[index_roi] - y=(used_maps[*,*,indI,j])[index_roi] + y=(maps[*,*,indI,j])[index_roi] sx=la_power(mapII_ref[index_roi],0.5) - sy=la_power((used_maps[*,*,indI,j])[index_roi],0.5) + sy=la_power((maps[*,*,indI,j])[index_roi],0.5) ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good) comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels'] IF count_good NE 0 THEN BEGIN @@ -268,9 +269,9 @@ CASE use_method OF ENDIF ;Do the correlation between Q maps x=[mapQ_ref[index_roi],mapU_ref[index_roi]] - y=[(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]] + y=[(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]] sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5) - sy=la_power([(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]],0.5) + sy=la_power([(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]],0.5) ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good) comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels'] IF count_good NE 0 THEN BEGIN diff --git a/src/idl/dustem_write_chrg.pro b/src/idl/dustem_write_chrg.pro index 8ccf24a..be14472 100755 --- a/src/idl/dustem_write_chrg.pro +++ b/src/idl/dustem_write_chrg.pro @@ -94,7 +94,7 @@ FOR i=0L,Nst-1 DO BEGIN fv=str_sep(ffile,'/') file=dir+fv(n_elements(fv)-1) ;Open file - message,'Will write CHRG file '+file,/continue + message,'Will write CHRG file '+file,/info openw,unit,file,/get_lun ;Print comments IF Ncomments NE 0 THEN BEGIN -- libgit2 0.21.2