src/idl/dustem_activate_plugins.pro

Jean-Philippe Bernard
1 parent 755a4236
Showing 12 changed files with 152 additions and 149 deletions Show diff stats
src/idl/dustem_brute_force_fit.pro
src/idl/dustem_define_grid.pro
src/idl/dustem_init_params.pro
src/idl/dustem_init_plugins.pro
src/idl/dustem_make_sed_table.pro
src/idl/dustem_param_range2param_values.pro
src/idl/dustem_plot_dataset.pro
src/idl/dustem_plugin_phangs_stellar_continuum.pro
src/idl/dustem_read_cb19_stellar_templates.pro
src/idl/dustem_read_emiles_stellar_templates.pro
src/idl/dustem_set_params.pro
src/idl/dustem_write_isrf_release.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)
-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:
@@ -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
@@ -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)
@@ -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]
@@ -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
@@ -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
@@ -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:
@@ -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)
@@ -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
@@ -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
@@ -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: