improved

Jean-Philippe Bernard
1 parent d5838901
Showing 3 changed files with 93 additions and 44 deletions Show diff stats
LabTools/IRAP/JPB/fit_phangs_ngc0628_nir_continuum.pro
src/idl/dustem_brute_force_fit.pro
src/idl/dustem_define_grid.pro
@@ -859,13 +859,13 @@ grid_brute_force:
 ;stop
  
 Nvor=max(voronoi_id)
-GOs=fltarr(Nvor)+la_undef()
+G0s=fltarr(Nvor)+la_undef()
 Ypahs=fltarr(Nvor)+la_undef()
 YVSGs=fltarr(Nvor)+la_undef()
 facts=fltarr(Nvor)+la_undef()
 chi2s=fltarr(Nvor)+la_undef()
 rchi2s=fltarr(Nvor)+la_undef()
-dGOs=fltarr(Nvor)+la_undef()
+dG0s=fltarr(Nvor)+la_undef()
 dYpahs=fltarr(Nvor)+la_undef()
 dYVSGs=fltarr(Nvor)+la_undef() 
 use_NHmap=NHCO       ;used NH map in 1.e21 (from CO)
@@ -897,11 +897,26 @@ restore,file,/verb
 ind=where(vor_class NE 0)
 class_min=min(vor_class[ind])
 class_max=max(vor_class[ind])
+restore,data_dir+'ngc0628_isrf_min_prediction.sav',/verb
+;% RESTORE: Restored variable: ISRFS.
+;% RESTORE: Restored variable: OBJECT_DISTANCE.
+;% RESTORE: Restored variable: OBJECT_THICKNESS.
+;% RESTORE: Restored variable: SOURCE_NAME.
  
 ;For test
-class_min=14
+class_min=15
 class_max=15
  
+bidon=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
+
+fixed_parameters_description=['dustem_plugin_phangs_class_isrf_2']   ;This is the ISRF G0 factor
+
+file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT'
+;file='/Users/jpb/Soft_Libraries/dustem_fortran/data/ISRF.DAT'
+readcol,file,Mathis_wavs,Mathis_ISRF
+Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths)
+Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0)
+
 FOR isrf_class=class_min,class_max DO BEGIN
 	;==== This is to use the grid for this ISRF class
     isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
@@ -913,7 +928,10 @@ FOR isrf_class=class_min,class_max DO BEGIN
 			IF vvid mod 100 EQ 0 THEN BEGIN
 		   		message,'isrf#'+strtrim(isrf_class,2)+' Doing sed '+strtrim(vvid,2)+' '+strtrim(1.*vvid/Nvor_class*100.,2)+' %',/continue
 			ENDIF
+			;=== Get the G0 value for that vid
 			vid=ind_class[vvid]
+			G0=interpol(ISRFs[*,vid],isrf_wavelengths,1.0)/Mathis_1mic
+			fixed_parameters_values=[G0]
 	   		sed=*seds_ptr[vid]
 	   		;=== restrict sed to requested filters
 			sed=sed[ind_filters]
@@ -924,16 +942,25 @@ FOR isrf_class=class_min,class_max DO BEGIN
 	   		ind=where(sed.stokesI EQ la_undef(),count)
 	   		IF count EQ 0 THEN BEGIN
 			   	params=dustem_brute_force_fit(sed,table_name,use_filters,fact=fact,chi2=chi2,/normalize,rchi2=rchi2,show_sed=show_sed $
-		   								,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max)
-		   		GOs[vid]=params[0]
-		   		Ypahs[vid]=params[1]/fact ;extenssive quantities must be divided by normalization factor
-				Yvsgs[vid]=params[2]/fact
+		   								,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max $
+		   								,fixed_parameters_description=fixed_parameters_description,fixed_parameters_values=fixed_parameters_values)
+			   	;This would be if there was no fixed_parameters
+		   		;G0s[vid]=params[0]
+		   		;Ypahs[vid]=params[1]/fact ;extenssive quantities must be divided by normalization factor
+				;Yvsgs[vid]=params[2]/fact
+		   		G0s[vid]=G0
+		   		Ypahs[vid]=params[0]/fact ;extenssive quantities must be divided by normalization factor
+				Yvsgs[vid]=params[1]/fact
 			   	facts[vid]=fact
 			   	chi2s[vid]=chi2
 			   	rchi2s[vid]=rchi2
-			   	dGOs[vid]=params_uncertainties[0]
-			   	dYpahs[vid]=params_uncertainties[1]/fact ;extenssive quantities must be devided by normalization factor
-			   	dYvsgs[vid]=params_uncertainties[2]/fact
+			   	;This would be if there was no fixed_parameters
+			   	;dG0s[vid]=params_uncertainties[0]
+			   	;dYpahs[vid]=params_uncertainties[1]/fact ;extenssive quantities must be devided by normalization factor
+			   	;dYvsgs[vid]=params_uncertainties[2]/fact
+			   	dG0s[vid]=0.   ;should use difference between fixed param and nearest table neighbor
+			   	dYpahs[vid]=params_uncertainties[0]/fact ;extenssive quantities must be devided by normalization factor
+			   	dYvsgs[vid]=params_uncertainties[1]/fact
 			   	;print,params[0],params[1],params[2],fact,chi2
 		   ENDIF
 		ENDFOR
@@ -945,7 +972,7 @@ ENDFOR
 ;stop
  
 ;save,GOs,Ypahs,Yvsgs,file=dir+source_name+'DBP90_JWST_G0_YPAH_YVSG.fits'
-save,GOs,Ypahs,Yvsgs,file=data_dir+source_name+'_DBP90_JWST_G0_YPAH_YVSG_ISRFclasses.fits'
+save,G0s,Ypahs,Yvsgs,facts,chi2s,rchi2s,file=data_dir+source_name+'_DBP90_JWST_G0_YPAH_YVSG_ISRFclasses.fits'
  
 ;Nvor=max(voronoi_id)
  
@@ -953,20 +980,20 @@ save,GOs,Ypahs,Yvsgs,file=data_dir+source_name+&#39;_DBP90_JWST_G0_YPAH_YVSG_ISRFcla
 ;=== make maps of parameters
 chi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 rchi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
-GO_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
+G0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 Ypah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 YVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
-dGO_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
+dG0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 dYpah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 dYVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 fact_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
 FOR vid=0LL,Nvor-1 DO BEGIN
 	IF vid mod 100 EQ 0 THEN print,(1.*vid)/Nvor*100.
 	;ind=where(voronoi_id EQ vid)
-    GO_map[*all_seds_indices[vid]]=GOs[vid]
+    G0_map[*all_seds_indices[vid]]=G0s[vid]
     Ypah_map[*all_seds_indices[vid]]=Ypahs[vid]
     YVSG_map[*all_seds_indices[vid]]=Yvsgs[vid]
-    dGO_map[*all_seds_indices[vid]]=dGOs[vid]
+    dG0_map[*all_seds_indices[vid]]=dG0s[vid]
     dYpah_map[*all_seds_indices[vid]]=dYpahs[vid]
     dYVSG_map[*all_seds_indices[vid]]=dYvsgs[vid]
     chi2_map[*all_seds_indices[vid]]=chi2s[vid]
@@ -980,13 +1007,27 @@ ENDFOR
 ;FOR vid=0LL,Nvor-1 DO chi2_map[*all_seds_indices[vid]]=chi2s[vid]
 stop
  
+;G0 histogram
 win=0L
+window,win,xsize=800,ysize=800 & win=win+1
+
+ind=where(G0s NE la_undef())
+res=histogram(G0s[ind],locations=xv)
+cgplot,xv,res,psym=10,title='G0s histogram',xtit='G0',ytit='Number',/ylog,yrange=[1,max(res)],/ysty
+
+ind=where(Ypahs NE la_undef())
+res=histogram(Ypahs[ind],locations=xv)
+cgplot,xv,res,psym=10,title='Ypahs histogram',xtit='Ypah',ytit='Number',/ylog,yrange=[1,max(res)],/ysty
+
 window,win,xsize=800,ysize=900 & win=win+1
  
 obp=[1.1,0.,1.15,1]
-image_cont20,GO_map,Href,/square,imrange=[0.01,100],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='G0'
-image_cont20,Ypah_map,Href,/square,imrange=[1,100]*1.e-3,image_color_table='jpbloadct',/silent
-image_cont20,YVSG_map,Href,/square,imrange=[1,100]*1.e-3,image_color_table='jpbloadct',/silent
+image_cont20,G0_map,Href,/square,imrange=[-2,40],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='G0'
+image_cont20,la_log10(G0_map),,Href,/square,imrange=[-0.2,2],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='G0'
+
+image_cont20,la_log10(Ypah_map),Href,/square,imrange=[-2.,6],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='log10(Ypah)'
+
+image_cont20,la_log10(YVSG_map),Href,/square,imrange=[-2.,6],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='log10(Yvsg)'
  
 image_cont20,chi2_map,Href,/square,imrange=[0.01,100]*1.e4,image_color_table='jpbloadct',/silent
  
@@ -1004,7 +1045,7 @@ image_cont20,jwst_4_coutours[*,*,0],href,/square,imrange=[-1.,10],image_color_ta
  
 window,0,xsize=800,ysize=1000
 aa=alog10(GO_map)
-wset,0 & image_cont20,GO_map,href,/square,imrange=[-1.,100],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1
+wset,0 & image_cont20,GO_map,href,/square,imrange=[-2.,100],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1
 wset,0 & image_cont20,aa,href,/square,imrange=[0.5,2.5],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1,levels=levs
 window,1,xsize=800,ysize=1000
 wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1,levels=levs
 FUNCTION dustem_brute_force_fit,sed $
-							   ,table_name $
-							   ,filters $
-							   ,normalize=normalize $
-							   ,params_hit=params_hit $
-							   ,params_uncertainties=params_uncertainties $
-							   ,params_min=params_min $
-							   ,params_max=params_max $
-							   ,fact=fact $
-							   ,chi2=chi2 $
-							   ,rchi2=rchi2 $
-							   ,best_sed=best_sed $
-							   ,best_grid_sed=best_grid_sed $
-							   ,show_sed=show_sed $
-							   ,nostop=nostop $
-							   ,reset=reset
+							   							 ,table_name $
+							   							 ,filters $
+							   							 ,fixed_parameters_description=fixed_parameters_description $
+							   							 ,fixed_parameters_values=fixed_parameters_values $
+							   							 ,normalize=normalize $
+							   							 ,params_hit=params_hit $
+							   							 ,params_uncertainties=params_uncertainties $
+							   							 ,params_min=params_min $
+							   							 ,params_max=params_max $
+							   							 ,fact=fact $
+							   							 ,chi2=chi2 $
+							   							 ,rchi2=rchi2 $
+							   							 ,best_sed=best_sed $
+							   							 ,best_grid_sed=best_grid_sed $
+							   							 ,show_sed=show_sed $
+							   							 ,nostop=nostop $
+							   							 ,reset=reset
  
 ;+
 ; NAME:
@@ -29,7 +31,8 @@ FUNCTION dustem_brute_force_fit,sed $
 ;    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:
-;    None
+;    fixed_parameters_description = if set, describes parameters of the table that should be considered fixed.
+;    fixed_parameters_values      = values of parameters of the table that should be considered fixed.
 ; OUTPUTS:
 ;    res = best values of the free parameters. Not that linear parameters to the sed should be multiplied by fact.
 ; OPTIONAL OUTPUT PARAMETERS:
@@ -79,6 +82,7 @@ IF exist EQ 0 or keyword_set(reset)THEN BEGIN
     dustem_define_grid,table_name,filters
 ENDIF
  
+;=== Should marginalize table for fixed parameters, if any
 seds=!dustem_grid.seds
 Nseds=!dustem_grid.Nsed
 Nfilters=!dustem_grid.Nfilters
@@ -86,8 +90,12 @@ 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}
  
+IF keyword_set(fixed_parameters_description) THEN BEGIN
+	dustem_marginalize_grid,fixed_parameters_description,fixed_parameters_values,seds,Nseds,Nfilters,Nparams,table_params,table_pmins,table_pmaxs
+ENDIF
+
+;=== These are the full arrays needed for calculations
 chi2s=dblarr(Nseds)
 grid_seds=dblarr(Nseds,Nfilters)
 total_seds=dblarr(Nseds)
@@ -109,7 +117,7 @@ ENDFOR
 ;stop
  
 ;=== loop over grid lines to compute chi2s
-;we should actually store the total of grid sed for each grid sed into the grid fits file (!)
+;we do actually store the total of grid sed for each grid sed into the grid fits file
 ;we could also store the total of observed seds in the observed sed save fits file (!)
 facts=dblarr(Nseds)
 facts[*]=1.d0
@@ -119,12 +127,12 @@ IF keyword_set(normalize) THEN BEGIN
 		facts[i]=total(sed.stokesI)/total_seds[i]
 		this_sed=sed.stokesI/facts[i]
 		this_var=sed.sigmaII/facts[i]^2
-	    ;chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
-	    chi2s[i]=total((this_sed-grid_seds[i,*])^2/this_var)
+	  ;chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
+	  chi2s[i]=total((this_sed-grid_seds[i,*])^2/this_var)
 	ENDFOR
 ENDIF ELSE BEGIN
 	FOR i=0L,Nseds-1 DO BEGIN
-	    chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
+	   chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
 	ENDFOR
 ENDELSE
 ;stop
@@ -139,14 +147,14 @@ 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
+;==== Check bumping to table edges.
 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
+;==== Compute parameter uncertainties
 ;stop
 conf_level=90./100.
 ;print,Nfreedom,conf_level
@@ -69,8 +69,8 @@ fname=&#39;bidon&#39;
 count=1
 WHILE count NE 0 DO BEGIN
 	fname=sxpar(h,'FNAME'+strtrim(i+1,2),count=count)
-    IF count NE 0 THEN table_filters=[table_filters,strtrim(fname,2)]
-    i=i+1
+  IF count NE 0 THEN table_filters=[table_filters,strtrim(fname,2)]
+  i=i+1
 ENDWHILE
 table_filters=table_filters[1:*]
 Ntfilters=n_elements(table_filters)