diff --git a/src/idl/dustem_brute_force_fit.pro b/src/idl/dustem_brute_force_fit.pro index 0b7eb77..0e8f894 100644 --- a/src/idl/dustem_brute_force_fit.pro +++ b/src/idl/dustem_brute_force_fit.pro @@ -25,7 +25,7 @@ FUNCTION dustem_brute_force_fit,sed $ ; CALLING SEQUENCE: ; res=dustem_brute_force_fit(sed,table_name,filters[,/normalize][,fact=][,chi2=][,rchi2=][rchi2][best_sed=][,/show_sed][,/nostop][,/reset]) ; INPUTS: -; sed = sed to be fit. Should contain all filters initialy provided to the routine +; sed = sed to be fit. Should contain all filters provided to the routine trhough input variable filter ; table_name = name of the fits file to be used as a grid (see dustem_make_sed_table.pro for creation) ; filters = filters to be used in the grid table (only used if !dustem_grid does not exist or if /reset) ; OPTIONAL INPUT PARAMETERS: @@ -73,6 +73,7 @@ ENDIF ;sed must be normalized to NH=1e20 H/cm2 ;==== test for the existence of a loaded grid and initialize if needed. +;==== only filters contained in the filters variable are kept in the grid, and the sum of grid seds are recomputed accordingly defsysv,'!dustem_grid',0,exist=exist IF exist EQ 0 or keyword_set(reset)THEN BEGIN dustem_define_grid,table_name,filters @@ -82,6 +83,11 @@ seds=!dustem_grid.seds Nseds=!dustem_grid.Nsed Nfilters=!dustem_grid.Nfilters Nparams=!dustem_grid.Nparams +table_params=!dustem_grid.params +table_pmins=!dustem_grid.pmin_values +table_pmaxs=!dustem_grid.pmax_values +;dustem_grid={Nsed:Nsed,Nfilters:Nfilters,Nparams:Ntparams,filters:filters,params:table_params,seds:st,pmin_values:table_pmins,pmax_values:table_pmaxs} + chi2s=dblarr(Nseds) grid_seds=dblarr(Nseds,Nfilters) total_seds=dblarr(Nseds) @@ -128,21 +134,30 @@ ind=where(chi2s EQ chi2,count) params=grid_param_values[ind[0],*] ;these are the best fit parameters fact=facts[ind[0]] ; This is the corresponding scaling factor -Nfreedom=Nseds-Nparams-1 ;This is the degree of freedom +Nfreedom=Nfilters-Nparams-1 ;This is the degree of freedom rchi2=chi2/Nfreedom ;This is the reduced chi2 best_grid_sed=reform(grid_seds[ind[0],*]) ;This is the best SED fround in the grid. best_sed=best_grid_sed*fact ;This is the best scaled grid SED matching the input SED. ;==== should check bumping to table edges there -;params_hit=intarr(Nparams) -;FOR i=0L,Nparams-1 DO BEGIN -; IF params[i] EQ minvalues[i] THEN params_hit[i]=-1 -; IF params[i] EQ maxvalues[i] THEN params_hit[i]=1 -;ENDFOR +params_hit=intarr(Nparams) +FOR i=0L,Nparams-1 DO BEGIN + IF params[i] EQ table_pmins[i] THEN params_hit[i]=-1 + IF params[i] EQ table_pmaxs[i] THEN params_hit[i]=1 +ENDFOR ;==== should compute parameter uncertainties there +;stop conf_level=90./100. -dchi2=delta_chi2(Nfreedom,conf_level) +;print,Nfreedom,conf_level +IF Nfreedom NE 0. THEN BEGIN + dchi2=delta_chi2(Nfreedom,conf_level) +ENDIF ELSE BEGIN + message,'Nfreedom = 0. Beware of meaningless fit',/continue + stop + message,'seting Nfreedom = 1. Beware of meaningless fit',/continue + dchi2=delta_chi2(1.,conf_level) +ENDELSE ind=where(chi2s LE chi2+dchi2,count) params_min=fltarr(Nparams) params_max=fltarr(Nparams) diff --git a/src/idl/dustem_define_grid.pro b/src/idl/dustem_define_grid.pro index d5b9fb0..79958de 100644 --- a/src/idl/dustem_define_grid.pro +++ b/src/idl/dustem_define_grid.pro @@ -1,4 +1,5 @@ -PRO dustem_define_grid,table_name,filters,help=help +PRO dustem_define_grid,table_name,filters,help=help, $ + table_params=table_params,table_pmins=table_pmins,table_pmaxs=table_pmaxs,table_pnvals=table_pnvals,table_plogs=table_plogs ;+ ; NAME: @@ -8,7 +9,7 @@ PRO dustem_define_grid,table_name,filters,help=help ; CATEGORY: ; Dustem ; CALLING SEQUENCE: -; dustem_define_grid,table_name,filters +; dustem_define_grid,table_name,filters[,table_params=][,table_pmins=][,table_pmaxs=][,table_pnvals=][,table_plogs=] ; INPUTS: ; table_name : fits file name containing the grid ; filters : name of dustemwrap filters to be included in the GRID @@ -17,7 +18,11 @@ PRO dustem_define_grid,table_name,filters,help=help ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: -; None +; table_params= parameter names in the grid table +; table_pmins= parameter min values in the grid table +; table_pmaxs= parameter max values in the grid table +; table_pnvals= parameter number of values in the grid table +; table_plogs= 1=parameter values in log in the grid table, 0=linear ; ACCEPTED KEY-WORDS: ; help = If set, print this help ; COMMON BLOCKS: diff --git a/src/idl/dustem_init_params.pro b/src/idl/dustem_init_params.pro index f74b20e..4b29224 100644 --- a/src/idl/dustem_init_params.pro +++ b/src/idl/dustem_init_params.pro @@ -194,6 +194,7 @@ if keyword_set(polarization) and polexists eq 0 then $ ; basic sanity checking of plugins ; to be written +;stop ;== INITIALIZE DUST MODEL PARAMETERS and PLUGIN FREE PARAMETERS dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims,fixed=fixed,tied=tied @@ -203,7 +204,9 @@ IF not keyword_set(no_reset_plugin_structure) THEN BEGIN ENDIF ;== INITIALIZE ANY FIXED PARAMETERS FOR DUST MODEL AND PLUGINS -if keyword_set(fpd) then dustem_init_fixed_params,fpd,fiv +IF keyword_set(fpd) THEN BEGIN + dustem_init_fixed_params,fpd,fiv +ENDIF message,'Finished initializing free and fixed parameters of dust model and plugins',/info diff --git a/src/idl/dustem_init_plugins.pro b/src/idl/dustem_init_plugins.pro index febbada..b27f71e 100644 --- a/src/idl/dustem_init_plugins.pro +++ b/src/idl/dustem_init_plugins.pro @@ -103,11 +103,12 @@ ENDELSE ;==== invok each plugin to get their scopes and parameter tag names FOR i=0L,Nplugins-1 DO BEGIN + scope=1 str='toto='+plugin_st[i].name+'(scope=scope)' str=str[0] toto=execute(str) plugin_st[i].scope=scope - + paramtag=1 ;otherwise, paramtag may nt be returned by plugin str='toto='+plugin_st[i].name+'(paramtag=paramtag)' str=str[0] toto=execute(str) diff --git a/src/idl/dustem_make_sed_table.pro b/src/idl/dustem_make_sed_table.pro index 9e3cbb3..1e8c57e 100644 --- a/src/idl/dustem_make_sed_table.pro +++ b/src/idl/dustem_make_sed_table.pro @@ -9,6 +9,8 @@ PRO dustem_make_sed_table,model, $ filters=filters, $ log=log, $ show_seds=show_seds, $ + use_isrf_class=use_isrf_class, $ + isrf_class=isrf_class, $ help=help ;+ @@ -31,6 +33,7 @@ PRO dustem_make_sed_table,model, $ ; filters : name of dustemwrap filters to be included in the SED calculations (default = IRAS filters) ; fpd : fixed parameter description ; fiv : fixed parameter values +; isrf_class ; If set use_isrf_class is set, uses ISRF from this PHANGs ISRF class ; filename : output fits file name for the table (default ='/tmp/dustem_seds.fits') ; log : array indicating if a given parameter is to be sampled in linear (0) or log scale (1) ; show_seds : if set, plots SEDs as they are computed. @@ -39,6 +42,7 @@ PRO dustem_make_sed_table,model, $ ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: +; use_isrf_class = If set, uses ISRF from a given PHANGs ISRF class ; help = If set, print this help ; COMMON BLOCKS: ; None @@ -132,14 +136,15 @@ dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext ;=== initialize at least one parameter as 'free' ;=== this is only to keep dustemwrap happy in its management of variables ;=== no fitting is going to be done -pd = ['(*!dustem_params).gas.G0'] ; G0 -pval = [1.] +;pd = ['(*!dustem_params).gas.G0'] ; G0 +;pval = [-1.e-20] +pd=parameters_description ;=== declare some other parameters as fixed ;=== this shows how to add any non-dust plugins to the observed spectra/SED ;=== for the dust-model parameters, it is not necessary, but allows the ;=== user to (i) change default values and (2) see the adopted values in the graphical output windows -fpd=parameters_description +;fpd=parameters_description seds=ptrarr(Nc) yrange=[1.e-5,1.e5] colors=long(range_gen(Nc,[0,255])) @@ -149,13 +154,16 @@ FOR i=0L,Nc-1 DO BEGIN ;dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext dummy=0 ;otherwise model is not recomputed !! - fpval = *parameter_values[i] - + ;fpval = *parameter_values[i] + pval = *parameter_values[i] message,'================= computing model '+strtrim(i,2)+'/'+strtrim(Nc,2)+' ('+strtrim(1.*i/Nc*100.)+'%)',/continue - print,'parameter_values=',fpval + ;print,'parameter_values=',fpval + print,'parameter_values=',pval ;=== initialise all this information into dustemwrap's brain - dustem_init_params,model,pd,pval,fpd=fpd,fiv=fpval,pol=use_polarisation + ;stop + ;dustem_init_params,model,pd,pval,fpd=fpd,fiv=fpval,pol=use_polarisation + dustem_init_params,model,pd,pval,fpd=fpd,fiv=fiv,pol=use_polarisation ;=== compute the predicted SED and extinction curve for the dust-model ;=== and any plugins @@ -196,11 +204,11 @@ FOR i=0L,Nc-1 DO BEGIN ;==== fill in dependent columns of the SED. sed=dustem_fill_sed_dependent_columns(sed) seds[i]=ptr_new(sed) - tit='Grid SEDs model '+model - xtit='Wavelength [mic]' - ytit='SED [MJy/sr for NH=1.e20 H/cm^2]' - IF keyword_set(show_seds) THEN BEGIN + IF keyword_set(show_seds) THEN BEGIN IF i EQ 0 THEN BEGIN + tit='Grid SEDs model '+model + xtit='Wavelength [mic]' + ytit='SED [MJy/sr for NH=1.e20 H/cm^2]' cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog ENDIF ELSE BEGIN cgoplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,color=colors[i] diff --git a/src/idl/dustem_param_range2param_values.pro b/src/idl/dustem_param_range2param_values.pro index f457ced..8bf9612 100644 --- a/src/idl/dustem_param_range2param_values.pro +++ b/src/idl/dustem_param_range2param_values.pro @@ -22,7 +22,7 @@ FOR i=0L,Nc-1 DO BEGIN vec[j]=(*par_values[j])[ij[0,j]] ENDFOR ivs[i]=ptr_new(vec) - print,vec + ;print,vec ;stop ENDFOR diff --git a/src/idl/dustem_plot_dataset.pro b/src/idl/dustem_plot_dataset.pro index 26373e2..170a9cc 100644 --- a/src/idl/dustem_plot_dataset.pro +++ b/src/idl/dustem_plot_dataset.pro @@ -3183,6 +3183,7 @@ if keyword_set(dataset) then begin ; testpop = STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' IF parameter_type EQ 'PLUGIN' THEN BEGIN + ;stop mm = where((*!dustem_plugin).name EQ plugin_name,count) ; ; Selecting a plugin through matching the string name of the plugin form the scope system variable with the one read from the parameter description vector prmtg = (*(*!dustem_plugin)[mm].paramtag) ;these are all plugin parameter names indtg=key diff --git a/src/idl/dustem_plugin_phangs_stellar_continuum.pro b/src/idl/dustem_plugin_phangs_stellar_continuum.pro index 44de1fa..3dbc84c 100644 --- a/src/idl/dustem_plugin_phangs_stellar_continuum.pro +++ b/src/idl/dustem_plugin_phangs_stellar_continuum.pro @@ -15,20 +15,21 @@ FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $ ; CATEGORY: ; DustEM, Distributed, Mid-Level, Plugin ; CALLING SEQUENCE: -; stellar_brightness=dustem_plugin_phangs_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help]) +; stellar_ISRF=dustem_plugin_phangs_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help]) ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; key = input parameter number -; First parameter: E(B-V) for extinction applied to the template, as used in Phangs +; First parameter: Amplitude factor for the emission +; Second parameter: E(B-V) for extinction to stars applied to the template, as used in Phangs ; Following 13*6 parameters are weights by which to multiply the MILES templates [Msun/pc^2] ; val = corresponding input parameter value ; OUTPUTS: -; stellar_brightness = stellar spectrum spectrum resampled on dustem wavelengths [MJy/sr] +; stellar_isrf = stellar ISRF spectrum resampled on dustem wavelengths [MJy/sr] ; OPTIONAL OUTPUT PARAMETERS: -; scope = scope of the plugin +; scope = if set, return the scope of the plugin ; paramdefault = default values of parameters -; paramtag = plugin parameter names as strings +; paramtag = if set, return the return plugin parameter names as strings ; ACCEPTED KEY-WORDS: ; help = if set, print this help ; method = SSP used. can be EMILES or CB19. default='EMILES' @@ -56,7 +57,8 @@ IF keyword_set(help) THEN BEGIN goto,the_end ENDIF -IF keyword_set(scope) THEN BEGIN +IF keyword_set(scope) THEN BEGIN + scope='ADD_SED' out=0 goto, the_scope ENDIF @@ -95,7 +97,7 @@ metalicity_values=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397 Nage=n_elements(age_values) Nmetalicity=n_elements(metalicity_values) -Nparam=Nage*Nmetalicity+2 ;+2 is for Amplitude and opacity +Nparam=Nage*Nmetalicity+2 ;+2 is for Amplitude and E(B-V) ;==== define parameter tags IF keyword_set(paramtag) THEN BEGIN @@ -170,6 +172,7 @@ FOR i=2L,Nparam-1 DO BEGIN output[*,0]=output[*,0]+paramvalues[i]*Mstar*ssp message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/info ;ENDIF + ;stop ENDIF ENDFOR output[*,0]=output[*,0]*paramvalues[0] @@ -201,7 +204,11 @@ fact=fact2*fact3*fact4 ;==== reden spectrum using Calzetti's extinction law (which is what Phangs does) ;cgplot,lambir,flux,/xlog,/ylog ;factor -1 is to reden, instead of unreden. factor 1e4 is to go from microns to Angstroem +;This uses the prescription of Calzetti et al 2000. It is applied from 912AA to 2.2 mic, no correction applied before and beyound that. +;Note that the supplied color excess should be that derived for the stellar continuum, EBV(stars), which (in Calzetti 2000) is related to the +;reddening derived from the gas, EBV(gas), via the Balmer decrement by EBV(stars) = 0.44*EBV(gas) calz_unred, lambir*1.e4, output[*,0], -1.*paramvalues[1], flux_red, R_V = R_V +;Note: we could actually use our own dustemwrap extinction curve to do that. output[*,0]=flux_red no_unred: @@ -214,7 +221,6 @@ output[*,0]=output[*,0]*fact ;stop the_scope: -scope='ADD_SED' the_paramtag: diff --git a/src/idl/dustem_read_cb19_stellar_templates.pro b/src/idl/dustem_read_cb19_stellar_templates.pro index bf12fca..c0b9f4c 100644 --- a/src/idl/dustem_read_cb19_stellar_templates.pro +++ b/src/idl/dustem_read_cb19_stellar_templates.pro @@ -67,7 +67,8 @@ Z=alog10(Z/Z0) ; -2.18299 -1.88196 -1.48401 -1.18298 -0.881955 -0.580925 -0.404834 -0.279895 -0.182985 -0.0368569 0.0474640 0.118045 0.294136 0.419075 0.595166 NZ=n_elements(Z) -dir_templates=!dustem_wrap_soft_dir+'/Data/SSPs/CB2019/CB19_Chabrier/Mu100/' +dir_templates=!phangs_data_dir+'/SSPs/CB2019/CB19_Chabrier/Mu100/' +;dir_templates=!dustem_wrap_soft_dir+'/Data/SSPs/CB2019/CB19_Chabrier/Mu100/' files=file_search(dir_templates+'*.fits') Nfiles=n_elements(files) ;metalicities=fltarr(Nfiles) diff --git a/src/idl/dustem_read_emiles_stellar_templates.pro b/src/idl/dustem_read_emiles_stellar_templates.pro index 9fd16fd..81b81d4 100644 --- a/src/idl/dustem_read_emiles_stellar_templates.pro +++ b/src/idl/dustem_read_emiles_stellar_templates.pro @@ -53,7 +53,8 @@ IF keyword_set(help) THEN BEGIN END ;dir_templates=!dustem_wrap_soft_dir+'/Data/SSPs/ssp_emiles-ch/' -dir_templates=!dustem_wrap_soft_dir+'/Data/SSPs/ssp_emiles-ch_extended/' +dir_templates=!phangs_data_dir+'/SSPs/ssp_emiles-ch_extended/' +;dir_templates=!dustem_wrap_soft_dir+'/Data/SSPs/ssp_emiles-ch_extended/' ;read templates info file ;IMF slope Z Age Mtotal M(*+remn) M* Mremn Mgas M(*+remn)/Lv M*/Lv Mv unknown diff --git a/src/idl/dustem_set_params.pro b/src/idl/dustem_set_params.pro index ee0533b..640c629 100755 --- a/src/idl/dustem_set_params.pro +++ b/src/idl/dustem_set_params.pro @@ -86,13 +86,8 @@ FOR i=0L,Nparams-1 DO BEGIN ;Not to be used for the moment. - - END ENDCASE - - - ENDFOR dustem_write_all,*!dustem_params,!dustem_dat diff --git a/src/idl/dustem_write_isrf_release.pro b/src/idl/dustem_write_isrf_release.pro index 9ca1dcb..7f57b43 100755 --- a/src/idl/dustem_write_isrf_release.pro +++ b/src/idl/dustem_write_isrf_release.pro @@ -59,7 +59,7 @@ ENDIF ;IF PLUGINS AREN'T SET. WRITE THE DEFAULT DUSTEM ISRF.DAT FILE Nplugins=n_elements(*!dustem_plugin) scopes = strarr(Nplugins) -IF (*!dustem_plugin)[0].name EQ 'NONE' THEN GOTO, the_after +IF (*!dustem_plugin)[0].name EQ 'NONE' THEN GOTO,write_isrf_file ;IF PLUGINS ARE SET, LOOP OVER THEIR SCOPES. ;Getting the actual scopes @@ -69,141 +69,108 @@ FOR i=0L, Nplugins-1 DO BEGIN ;ENDIF ELSE goto, the_after ENDFOR -;TESTING FOR THE TWO CURRENT ISRF PLUGINS -tstpop = where(scopes EQ 'STELLAR_POPULATION',ctpop) -;tstusrisrf = where(scopes EQ 'USER_ISRF',ctusrisrf) -tstusrisrf = where(scopes EQ 'REPLACE_ISRF',ctusrisrf) +Ncomments=3 ;default number of comments in the .DAT file (minus the ISRF used) +c=strarr(Ncomments) +c[*]='#' +c[0]='# DUSTEM_WRAP: Interstellar radiation field ' + +;Testing for the two current ISRF plugin scopes +ind_replace_isrf = where(scopes EQ 'REPLACE_ISRF',count_replace_isrf) +ind_add_isrf = where(scopes EQ 'ADD_ISRF',count_add_isrf) -IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN - ;stop +;stop +IF count_add_isrf NE 0 OR count_replace_isrf NE 0 THEN BEGIN + lambir=dustem_get_wavelengths() + ;IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN ;NB: Since we're using G0_mathis=1 by default now. We will test on the !dustem_params sysvar instead of the parameter description vectors G0_mathis = 1. ;Default value - IF ((*!dustem_params).g0) NE 1. then G0_mathis = ((*!dustem_params).g0) -; Fortran now uses negative numbers to specify Gas.dat G0 -; If positive G0 present in Gas.dat, it takes precedence over grain.dat G0 -; IF ((*!dustem_params).gas.g0) NE 1. then G0_mathis = float((*!dustem_params).gas.g0) - IF ((*!dustem_params).gas.g0) gt 0. then G0_mathis = float((*!dustem_params).gas.g0) - + IF ((*!dustem_params).g0) NE 1. THEN BEGIN + G0_mathis = ((*!dustem_params).g0) + ENDIF + ; Fortran now uses negative numbers to specify Gas.dat G0 + ; If positive G0 present in Gas.dat, it takes precedence over the G0 value in GRAIN.DAT + ; IF ((*!dustem_params).gas.g0) NE 1. then G0_mathis = float((*!dustem_params).gas.g0) + IF ((*!dustem_params).gas.g0) gt 0. THEN BEGIN + G0_mathis = float((*!dustem_params).gas.g0) + ENDIF ;Getting the ISRF default wavelengths ;using this structure - st=((*!dustem_params).isrf) + ;st=((*!dustem_params).isrf) final_isrf = fltarr(n_elements(st)) ;Initializing the ISRF vector - Ncomments=3 ;default number of comments in the .DAT file (minus the ISRF used) - c=strarr(Ncomments) - c[0]='# DUSTEM_WRAP: exciting radiation field featuring' - ;=== ADD to the ISRF - IF ctpop NE 0. THEN final_isrf=final_isrf+(*(*!dustem_plugin)[tstpop].spec) - IF ctusrisrf NE 0. THEN final_isrf=final_isrf+(*(*!dustem_plugin)[tstusrisrf].spec) - - IF ctpop NE 0. and ctusrisrf EQ 0. THEN c[1] = '# Stellar population ISRF' - IF ctpop EQ 0. and ctusrisrf NE 0. THEN c[1] = '# User-defined ISRF' - - IF ctpop NE 0. and ctusrisrf NE 0. THEN BEGIN - Ncomments+=1 - c=strarr(Ncomments) - c[0]='# DUSTEM_WRAP: exciting radiation field featuring' - c[1] = '# Stellar population ISRF' - c[2] = '# User-defined ISRF' + ;IF ctpop NE 0. THEN final_isrf=final_isrf+(*(*!dustem_plugin)[tstpop].spec) + ;JPB: This line below can have an undefined plugin ISRF if the REPLACE_ISRF plugin has not been called earlier except with /scope or paramtags keywords + IF count_replace_isrf NE 0. THEN BEGIN + IF ptr_valid((*!dustem_plugin)[ind_replace_isrf].spec) THEN BEGIN + c[1]='# DUSTEM_WRAP: replaced by radiation field from plugin '+(*!dustem_plugin)[ind_replace_isrf].name + ;final_isrf=final_isrf+(*(*!dustem_plugin)[tstusrisrf].spec) + final_isrf=interpol(*(*!dustem_plugin)[ind_replace_isrf].spec,lambir,st.lambisrf) + ENDIF ELSE BEGIN + message,'Replacement ISRF not found. Using Mathis',/continue + ;stop + mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT' + sst=dustem_read_isrf(mathis_isrf_datfile) + final_isrf=sst.isrf + ;stop + ;LEFT blank intensionally + ENDELSE ENDIF - - ;WHEN G0 IS USED (MATHIS FIELD ON) - IF G0_mathis GT 1E-10 THEN BEGIN ;discussed with JP - ;This is reading the mathis isrf from the fortran release - mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT' - mathis_isrf=dustem_read_isrf(mathis_isrf_datfile) - final_isrf=final_isrf/G0_mathis+mathis_isrf.isrf ;[0] - Ncomments=4 - c=strarr(Ncomments) - c[0]='# DUSTEM_WRAP: exciting radiation field featuring' - - IF ctpop NE 0 and ctusrisrf EQ 0 THEN c[1] = '# Stellar population ISRF' - IF ctpop EQ 0 and ctusrisrf NE 0 THEN c[1] = '# User-Defined ISRF' - - IF ctpop NE 0 and ctusrisrf NE 0 THEN BEGIN - Ncomments+=2 - c=strarr(Ncomments) - c[0]='# DUSTEM_WRAP: exciting radiation field featuring' - c[1] = '# Stellar population ISRF' - c[2] ='# User-defined ISRF' - c[3] = '# Mathis ISRF' - ENDIF ELSE c[Ncomments-3] = '# Mathis ISRF' - - ENDIF ELSE final_isrf = final_isrf/G0_mathis ;WHEN G0 ISN'T USED (MATHIS FIELD OFF) ie: lower values than 1.0E-10. - + IF count_add_isrf NE 0 THEN BEGIN + IF ptr_valid((*!dustem_plugin)[ind_add_isrf].spec) THEN BEGIN + c[2]='# DUSTEM_WRAP: added radiation field from plugin '+(*!dustem_plugin)[ind_replace_isrf].name + ;final_isrf=final_isrf+(*(*!dustem_plugin)[ind_add_isrf].spec) + final_isrf=final_isrf+interpol(*(*!dustem_plugin)[ind_add_isrf].spec,lambir,st.lambisrf) + ENDIF ELSE BEGIN + ;stop + ;LEFT blank intensionally + ENDELSE + ENDIF + ;JPB: I don't think this loop is in fact needed, as the default value (without plugins) is the Mathis ISRF and ADD_ISRF will just add to this + ;IF G0_mathis GT 1E-10 THEN BEGIN ;discussed with JP + ; ;This is reading the mathis isrf from the fortran release + ; mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT' + ; mathis_isrf=dustem_read_isrf(mathis_isrf_datfile) + ; final_isrf=final_isrf/G0_mathis+mathis_isrf.isrf ;This is because Fortran uses ISRF*G0 which is now final_isrf+G0*mathis_isrf + ; c=[c,'# Mathis ISRF'] + ; Ncomments=Ncomments+1 + ;ENDIF ELSE BEGIN + ; ;final_isrf = final_isrf/G0_mathis ;This is because Fortran uses ISRF*G0 + ;ENDELSE + ;=== set the ISRF in the st structure st.isrf = final_isrf +ENDIF ELSE BEGIN ;This is Mathis field only + ;==== default comment lines for ISRF.DAT + Ncomments=4 + c=strarr(Ncomments) + c[1]='# Mathis ISRF' +ENDELSE - ;Last .DAT header comments lines - c[Ncomments-2]='# Nbr of points' - c[Ncomments-1]='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)' - - ;Computing the 'G0' of the new ISRF (integral(ISRF)/integral(MATHIS_ISRF)) -; un=uniq(ma_isrf.isrf) ;because of the zeros -; -; int_mathis=INT_TABULATED((ma_isrf.lambisrf)[un],(ma_isrf.isrf)[un]) -; -; -; ;integrating the composite isrf -; int_isrf=INT_TABULATED((ma_isrf.lambisrf)[un],(st.isrf)[un]) -; print, 'G0_stellar_population=' -; print, int_isrf/(int_mathis) - - -;WRITING THE NEW ISRF - - file=!dustem_dat+'data/ISRF.DAT' - openw,unit,file,/get_lun - - FOR i=0,Ncomments-1 DO BEGIN - printf,unit,c[i] - ENDFOR - - n_waves=n_elements(st) - printf,unit,n_waves - - FOR i=0L,n_waves-1 DO BEGIN - printf,unit,st(i).lambisrf,st(i).isrf - ENDFOR - - close,unit - free_lun,unit - - goto, the_end -ENDIF - - - -the_after: ;Defaut ISRF writing (MATHIS) - -;This was changed because 6 comments aren't needed -Ncomments=4 -c=strarr(Ncomments) +c[Ncomments-2]='# Nbr of points' +c[Ncomments-1]='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)' -;First lines of the ISRF.DAT file -c[0]='# DUSTEM_WRAP: exciting radiation field featuring' -c[1]='# Mathis ISRF' -c[2]='# Nbr of points' -c[3]='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)' +;==== Write the ISRF file +write_isrf_file: openw,unit,file,/get_lun - FOR i=0,Ncomments-1 DO BEGIN printf,unit,c[i] ENDFOR - n_waves=n_elements(st) printf,unit,n_waves - FOR i=0L,n_waves-1 DO BEGIN printf,unit,st[i].lambisrf,st[i].isrf ENDFOR - close,unit free_lun,unit +;print,(*!dustem_params).g0 +;print,(*!dustem_params).gas.g0 +;cgplot,st.lambisrf,st.isrf,/xlog,/ylog +;stop + the_end: -- libgit2 0.21.2