improved

Jean-Philippe Bernard
1 parent eb263d61
Showing 3 changed files with 74 additions and 146 deletions Show diff stats
src/idl/dustem_fit_intensity_phangs_example.pro
src/idl/dustem_plugin_phangs_stellar_continuum.pro
src/idl/dustem_read_emiles_stellar_templates.pro
@@ -156,118 +156,12 @@ pd = [ $
 ; initial values vector for run without synchrotron plugin
 iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3,100.] 
  
-; initial values vector for run including the synchrotron plugin
-; [G0,PAH0,VSG,BG,alpha_CR,Amp_syn]
-;iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]
-
-; initial values vector for run including the synchrotron plugin and
-; fixing the G0 and BB continuum parameters
-;iv =   [2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]
-
 Npar=n_elements(pd)
 ulimed=replicate(0,Npar)
 llimed=replicate(1,Npar) 
 llims=replicate(1.e-15,Npar)
 fpd=[] & fiv=[]
  
-; example using fixed parameters (ISRF and a NIR continuum)
-
-; Uncomment the following lines to include these fixed parameters
-;fpd = ['(*!dustem_params).G0', $    ;G0
-;     'dustem_plugin_continuum_1', $                    ;Temperature of a BB
-;     'dustem_plugin_continuum_2']                      ;Peak amplitude of a BB
-
-; fixed values vector 
-; [G0,T_BB,Amp_BB]
-;fiv=[0.5,750,1.e-2]
-
-;; ;==================================
-;; ;=== EXAMPLES FOR OTHER DUST MODELS
-;; ;=== START BELOW HERE
-;; ;==================================
-
-
-;; ;=== AN EXAMPLE FOR DL07
-;; ;=== Here we fit the dust abundances of the model and the
-;; ;=== intensity of the ISRF (via gas.G0 since the model
-;; ;=== includes spinning dust)
-;; ;=== The free parameters are all lower-bounded at zero.
-;; use_model='DL07' ; you should specify this above, or in the command line
-;; pd = [  '(*!dustem_params).gas.G0'], $             
-;;         '(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH0 mass fraction
-;;         '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
-;;         '(*!dustem_params).grains(2).mdust_o_mh', $    ;Gra
-;;         '(*!dustem_params).grains(3).mdust_o_mh', $    ;Gra
-;;         '(*!dustem_params).grains(4).mdust_o_mh']      ;aSil
-;;    iv =   [1.5 ,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3]
-;;    Npar=n_elements(pd)
-;;    ulimed=replicate(0,Npar)
-;;    llimed=replicate(1,Npar)
-;;    llims=replicate(1.e-15,Npar)
-;; ;   fpd=['(*!dustem_params).gas.G0']
-;; ;   fiv=[2.]
-   
-
-;; ;; ;=== AN EXAMPLE FOR MC10
-;; ;; ;=== Here we fit the dust abundances of the MC10 model, the
-;; ;; ;=== intensity of the dust-heating radiation field as well as a plug-in:
-;; ;; ;===  (i) continuum due to a blackbody 
-;; ;; ;=== The intensity of the dust-heating radiation field is fixed to
-;; ;; ;=== 1.5*G0 and the tmperature of the blackbody is fixed to 1200K
-;; ;; ;=== The free parameters in the fit are lower-bounded at zero.
-;; use_model='MC10' ; you should specify this above, or in the command line
-;; pd = [ $
-;;      '(*!dustem_params).grains(0).mdust_o_mh'$      ;PAH0 mass fraction
-;;      ,'(*!dustem_params).grains(1).mdust_o_mh' $    ;PAH1 mass fraction
-;;      ,'(*!dustem_params).grains(2).mdust_o_mh' $    ;amCBEx
-;;      ,'(*!dustem_params).grains(3).mdust_o_mh' $    ;amCBEx
-;;      ,'(*!dustem_params).grains(4).mdust_o_mh' $    ;aSilx
-;;      ,'dustem_plugin_continuum_2' $                 ;Intensity at peak of the continuum
-;;      ]
-;; iv =   [ 7.8e-4, 7.8e-4, 1.65e-4, 1.45e-3, 6.7e-3, 0.003]
-;; Npar=n_elements(pd)
-;; ulimed=replicate(0,Npar)
-;; llimed=replicate(1,Npar)
-;; llims=replicate(1.e-15,Npar)
-;; fpd=[ $
-;;     '(*!dustem_params).G0'   $    ; ISRF intensity 
-;;     ,'dustem_plugin_continuum_1' $    ;temperature of blackbody the produces the continuum 
-;;     ]
-;; fiv=[1.5, 1200.]
-
-
-
-;; ;=== AN EXAMPLE FOR J13
-;; ;=== Here we fit the dust abundances of the J13 model, the
-;; ;=== intensity of the dust-heating radiation field as well as two plug-ins:
-;; ;=== (i) free-free emission
-;; ;=== (ii)continuum due to a blackbody 
-;; ;=== The temperature of the blackbody is fixed to 1000K
-;; ;=== The free parameters in the fit are all lower-bounded at zero.
-;; use_model='J13' ; you should specify this above, or in the command line
-;; pd = [ $
-;;      '(*!dustem_params).G0' $                                   ;G0
-;;      ,'(*!dustem_params).grains(0).mdust_o_mh'$                  ;CM20 -- power law size distribution
-;;      ,'(*!dustem_params).grains(1).mdust_o_mh'$                 ;CM20 -- logN size distribution
-;;      ,'(*!dustem_params).grains(2).mdust_o_mh' $                ;aPyM5
-;;      ,'(*!dustem_params).grains(3).mdust_o_mh' $                ;aOlM5
-;;      ,'dustem_plugin_freefree_1' $                              ;ionized gas temperature
-;;      ,'dustem_plugin_freefree_2' $                              ;free-free amplitude
-;;      ,'dustem_plugin_continuum_2']                              ;Intensity at peak of the BB continuum
-;; iv =   [1.2,1.7e-3, 6.3e-4, 2.55e-3, 2.55e-3, 7500., 0.4, 0.001]         
-;; Npar=n_elements(pd)
-;; ulimed=replicate(0,Npar)
-;; llimed=replicate(1,Npar)
-;; llims=replicate(1.e-15,Npar)
-;; fpd=[ 'dustem_plugin_continuum_1'] ; Temperature of the BB
-;; fiv=[1000.]                       
-
-
-;;Nfix=n_elements(fpd)
-;;if n_elements(fiv) ne Nfix then begin
-;;   message,'Number of fixed parameters (fpd) does not equal number of initial values of fixed parameters (fiv)',/info
-;;   stop
-;;end
  
 if keyword_set(wait) then begin
    message,'Finished setting dust model and plug-in parameters: '+use_model,/info
@@ -2,6 +2,8 @@ FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $
                                 			   ,val=val $
                                 			   ,scope=scope $
                                 			   ,paramtag=paramtag $
+                                			   ,age_values=age_values $
+                                			   ,metalicity_values=metalicity_values $
                                 			   ,paramdefault=paramdefault $
                                 			   ,help=help
  
@@ -87,28 +89,29 @@ 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
 ;metalicity_values=met_values_emiles
+;This does not work either
+;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=Nmetalicity,bins=bins)
+;print,age_values
+;     0.030000000     0.050000000     0.080000000      0.15000000      0.25000000      0.40000000      0.60000000       1.0000000       1.7500000       3.0000000       5.0000000       8.5000000       13.500000
+;print,metalicity_values
+;     -1.49000    -0.960000    -0.350000    0.0600000     0.260000     0.400000
+
+;stop
  
 Nage=n_elements(age_values)
 Nmetalicity=n_elements(metalicity_values)
-Nparam=Nage*Nmetalicity+1    ;+1 is for opacity
+Nparam=Nage*Nmetalicity+2    ;+2 is for Amplitude and opacity
 paramdefault=dblarr(Nparam)
 paramdefault[*]=0.D0         ;default parameter values is 0.
  
 ;==== define parameter tags
-
-;If the test below is activated (ie paramtags computed only when requested) the pluggin does not work.
-;IF keyword_set(paramtag) THEN BEGIN
-	paramtag=strarr(Nparam)
-	paramtag[0]='E(B-V)'
-	ii=1L
-	FOR i=1L,Nparam-1 DO BEGIN
-		ij=index2ij([i-1],[Nage,Nmetalicity])
-		paramtag[ii]='Age'+strtrim(ij[0,0],2)+'Met'+strtrim(ij[0,1],2)
-		ii=ii+1
-	ENDFOR
-    ;out=0.
-    ;GOTO,the_paramtag
-;ENDIF
+paramtag=strarr(Nparam)
+paramtag[0]='Amplitude'
+paramtag[1]='E(B-V)'
+FOR i=2L,Nparam-1 DO BEGIN
+	ij=index2ij([i-2],[Nage,Nmetalicity])
+	paramtag[i]='Age'+strtrim(ij[0,0],2)+'Met'+strtrim(ij[0,1],2)
+ENDFOR
  
 ;==== decode key keyword. 
 ;stop
@@ -140,51 +143,78 @@ ENDIF
  
 lambir=dustem_get_wavelengths()
 Nwavs=n_elements(lambir)
-output=fltarr(Nwavs,3)
+output=dblarr(Nwavs,3)
 wavs=*!dustem_plugin_phangs_stellar_continuum.wavs
 Mstars=*!dustem_plugin_phangs_stellar_continuum.mstars
 ;==== sum up the stellar spectrum
 ;stop
-FOR i=1L,Nparam-1 DO BEGIN
+FOR i=2L,Nparam-1 DO BEGIN
 	IF paramvalues[i] NE 0. THEN BEGIN
 		;stop
 		;only interpolate within the template wavelengths
 		ind=where(lambir LE max(wavs) AND lambir GE min(wavs),count)
-		Mstar=Mstars[i-1]
+		Mstar=Mstars[i-2]
 		IF count NE 0 THEN BEGIN
-			ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-1]
+			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
 		ENDIF
 	ENDIF
 ENDFOR
-;==== Caution: at this point, output is unredenned Flambda in Lsun/pc^2/Angstroem
+output[*,0]=output[*,0]*paramvalues[0]
+flux=output[*,0]
+;no polarization for now
+output[*,1]=0.
+output[*,2]=0. 
+
+
+;==== Caution:
+;output is Flambda in ergs/s/cm2/AA
+;at this point, output is unredenned 
+;==== must be transformed into MJy/sr
+NHref=*!DUSTEM_HCD   ;NH dustemwrap reference [H/cm^2]
+fact2=lambir*1.e4      ;[AA] to go from F_lambda in ergs/s/cm2/AA to Fnu ergs/s/cm2
+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'
+
+;goto,no_unred
  
-stop
 ;==== reden spectrum using Calzetti's extinction law (which is what Phangs does)
-flux=output[*,0]
+;cgplot,lambir,flux,/xlog,/ylog
 ;factor -1 is to reden, instead of unreden. factor 1e4 is to go from microns to Angstroem
-calz_unred, lambir*1.e4, flux, -1.*paramvalues[0], flux_red, R_V = R_V
+calz_unred, lambir*1.e4, flux, -1.*paramvalues[1], flux_red, R_V = R_V
 output[*,0]=flux_red
  
-;no polarization for now
-output[*,1]=0.
-output[*,2]=0. 
+no_unred:
+
+;fact=Lsun*watts2ergs/Mpc2cm*1.e-4*4.*!pi/Mpc2cm/NHref   ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr to MJy/sr
+
+;fact=Lsun/pc*(lambir/1.e6)^2/c/4./!pi/pc*1.e20   ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr to MJy/sr
+;print,'flux=',interpol(flux,lambir,1.),' Lsun/Mpc^2/AA'
+;print,'fact1=',fact1
+;print,'fact2=',interpol(fact2,lambir,1.)
+;print,'fact3=',fact3
+;print,'fact=',interpol(fact,lambir,1.)
+;print,'flux=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' ergs/s/H'
+;print,'fact4=',interpol(fact4,lambir,1.)
+;print,'fact5=',fact5
+;print,'flux=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.)*interpol(fact4,lambir,1.),' MJy/sr'
  
-;at that stage, output is F_lambda in Lsun/A/pc^2
-;going to W/m2/sr/hz
-pc=3.085677581e16 ;m
-Lsun=3.828e26  ;Watts
-c=3.e5*1e3     ;m/s
-;Fact=Lsun       ;Watts/A
-;fact=fact*c     ;Watts/hz
-;fact=fact/4./!pi  ;Watts/Hz/sr
-;fact=fact/pc^2   ;Watts/m2/Hz/sr
-
-fact=Lsun/pc*c/4./!pi/pc   ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr
-output=output*fact
-;print,fact
-;      9.59803
+;output=output*fact
+output[*,0]=output[*,0]*fact
+output[*,1]=output[*,1]*fact
+output[*,2]=output[*,2]*fact
+
+print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr'
+;stop
  
 the_scope:
 scope='ADD_SED'
@@ -91,6 +91,10 @@ metalicity_values_4str=string(metalicity_values_4str,format=&#39;(F5.2)&#39;)
 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
+;stop
+
 ;Ntemplates=n_elements(paramvalues)-1
 Ntemplates=n_elements(age_values)*n_elements(metalicity_values)
 Nage=n_elements(age_values)