dustem_define_grid.pro 5.96 KB
PRO dustem_define_grid,table_name $
										  ,filters=filters $
										  ,help=help $
					   					,table_params=table_params $
					   					,table_pmins=table_pmins $
					   					,table_pmaxs=table_pmaxs $
					   					,table_pnvals=table_pnvals $
					   					,table_plogs=table_plogs $
					   					,status=status

;+
; NAME:
;    dustem_define_grid
; PURPOSE:
;    loads an SED Grid table into !dustem_grid
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    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
; OPTIONAL INPUT PARAMETERS:
;    NONE
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    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:
;    None
; SIDE EFFECTS:
;    !dustem_grid is updated
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_make_sed_table,'DBP90',['(*!dustem_params).G0','(*!dustem_params).grains(0).mdust_o_mh'], $
;								[0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',plog=[1,0],/show_seds, $
;                              	filters=[dustem_instru2filters('IRAC'),dustem_instru2filters('SPIRE')]
;    dustem_define_grid,'/tmp/IRAS_GO.fits',['IRAC3','SPIRE2']
;    help,!dustem_grid
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (2023)
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_define_grid'
  goto,the_end
ENDIF

;=== read the grid table
message,'reading grid data from '+table_name,/continue
;stop
info=file_info(table_name)
status=1
IF info.exists THEN BEGIN
	st0=mrdfits(table_name,0,h0)   ;This gets the main header
	st_params=mrdfits(table_name,1,h_params)
	st_seds=mrdfits(table_name,2,h_seds)
ENDIF ELSE BEGIN
	message,'grid data '+table_name+' not found',/continue
	status=-1
	goto,the_end
	;stop
ENDELSE

;stop
;=== get filters and params names from grid table header
i=0L
table_filters=['']
fname='bidon'
count=1
WHILE count NE 0 DO BEGIN
	fname=sxpar(h0,'FNAME'+strtrim(i+1,2),count=count)
  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)

;stop

table_params=['']
table_pmins=[0.d0]
table_pmaxs=[0.d0]
table_pnvals=[0L]
table_plogs=[0]
count=1
i=0L
WHILE count NE 0 DO BEGIN
	pname=sxpar(h0,'PNAME'+strtrim(i+1,2),count=count)
    IF count NE 0 THEN BEGIN
    	table_params=[table_params,strtrim(pname,2)]
    	pmin=sxpar(h0,'PMIN'+strtrim(i+1,2),count=count)
    	pmax=sxpar(h0,'PMAX'+strtrim(i+1,2),count=count)
    	pnval=sxpar(h0,'PNVAL'+strtrim(i+1,2),count=count)
    	plog=sxpar(h0,'PLOG'+strtrim(i+1,2),count=count)
    	table_pmins=[table_pmins,pmin]
    	table_pmaxs=[table_pmaxs,pmax]
    	table_pnvals=[table_pnvals,pnval]
    	table_plogs=[table_plogs,plog]
    ENDIF
    i=i+1
ENDWHILE
table_params=table_params[1:*]
table_pmins=table_pmins[1:*]
table_pmaxs=table_pmaxs[1:*]
table_pnvals=table_pnvals[1:*]
table_plogs=table_plogs[1:*]
Ntparams=n_elements(table_params)
;stop

;=== select only requested filters
;stop
use_filters=table_filters   ;default is to keep all filters
IF keyword_set(filters) THEN use_filters=filters

Nfilters=n_elements(use_filters)
keep_tfilters=intarr(Ntfilters)   ;filters to be kept in SED
add_tfilters=intarr(Nfilters)   ;filters to be added to the table SEDs
tags=tag_names(st_seds[0])
FOR i=0L,Nfilters-1 DO BEGIN
    ind=where(table_filters EQ use_filters[i],count)
    IF count NE 0 THEN BEGIN
    	keep_tfilters[ind]=1
    ENDIF ELSE BEGIN
    	add_tfilters[i]=1
    ENDELSE
ENDFOR
ind_add=where(add_tfilters EQ 1,count_add)
;stop

;==== remove un-needed filters from the table SEDs
FOR i=0L,Ntfilters-1 DO BEGIN
	IF keep_tfilters[i] EQ 0 THEN BEGIN
		col_name=tags[i]
		st_seds=structure_remove_column(st_seds,col_name)
	ENDIF
ENDFOR

;stop
previous_st_seds=st_seds

;==== add filters not present in the table SEDs (setting corresponding intensities to zero)
;==== while keeping column TOTAL at the end.
IF count_add NE 0 THEN BEGIN
	totals=st_seds.total
	st_seds=structure_remove_column(st_seds,'TOTAL')
	FOR i=0L,count_add-1 DO BEGIN
    col_name='I'+use_filters[i]
		st_seds=structure_add_column(st_seds,col_name,0.)
		table_filters=[table_filters,use_filters[i]]   ;== modify table_filters accordingly
  ENDFOR
  new_filters=strmid(tag_names(st_seds[0]),1,1000)
	Ntfilters=n_elements(table_filters)
  wavs=dustem_filter2wav(new_filters)
  order=sort(wavs)
  ;stop
  st_seds=structure_order_column(st_seds,order)
  table_filters=table_filters[order]
  ;stop
	st_seds=structure_add_column(st_seds,'TOTAL',0.)
	st_seds.total=totals
ENDIF

sed_totals=st_seds.(Nfilters-1)

;stop

;==== recompute totals
Nsed=n_elements(st_seds)
;==== recompute totals (needed if some columns have been removed)
FOR i=0L,Nsed-1 DO BEGIN
	this_total=0.d0
	FOR j=0L,Nfilters-1 DO BEGIN
		;this_total=this_total+(st[i]).(j+Ntparams)
		this_total=this_total+st_seds[i].(j)
	ENDFOR
	sed_totals[i]=this_total
ENDFOR

;stop

dustem_grid={Nsed:Nsed $ ;
					  ,Nfilters:Nfilters $ ;
					  ,Nparams:Ntparams $  ;
					  ,filters:use_filters $   ;
					  ,params:table_params $;
					  ,st_seds:st_seds $          ;
					  ,st_params:st_params $      ;
					  ,sed_totals:sed_totals $     ;
					  ,pmin_values:table_pmins $  ;
					  ,pmax_values:table_pmaxs $  ;
					  }

;defsysv,'!dustem_grid',dustem_grid
defsysv,'!dustem_grid',ptr_new(dustem_grid)
message,'Defined !dustem_grid',/continue

the_end:

END