dustem_set_data.pro 8.61 KB
FUNCTION dustem_set_data,help=help,ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI,qext=qext,uext=uext
;  qsed=qsed,used=used

;+
; NAME:
;    dustem_set_data
; PURPOSE:
;    Initializes the dustem data structure (!dustem_data)
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_set_data,spec[,help=]
; INPUTS:
;    spec: array of structures containing a SED
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    initializes !dustem_data
;    wavelengths are forced to be the central filter wavelength for
;    photometric data
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_init
;    dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
;    file=dir+'Gal_composite_spectrum.xcat'
;    spec=read_xcat(file,/silent)
;    dustem_set_data,spec
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
; 
;    V. Guillet (2012) : dustem_set_data is turned into a function which can be used 
;                        indifferently for sed, ext, polext
;-

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

;=== clean the data structure
!dustem_data.sed = ptr_new()
!dustem_data.ext = ptr_new()

if keyword_set(polext) then !dustem_data.polext = ptr_new()
if keyword_set(polsed) then !dustem_data.polsed = ptr_new()
if keyword_set(polfrac) then !dustem_data.polfrac = ptr_new()
indQ=where(sed.StokesQ NE la_undef(),NQ)
if NQ NE 0 then !dustem_data.qsed = ptr_new()
;if keyword_set(qsed) then !dustem_data.qsed = ptr_new()
indU=where(sed.StokesU NE la_undef(),NU)
if NU NE 0 then !dustem_data.used = ptr_new()
;if keyword_set(used) then !dustem_data.used = ptr_new()
if keyword_set(qext) then !dustem_data.qext = ptr_new()
if keyword_set(uext) then !dustem_data.uext = ptr_new()



IF keyword_set(sed) THEN BEGIN
  wavs=sed.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(sed.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(sed(ind).filter)
  ENDIF
  ;=== define observations structure
  ;obs_st={instru_names:sed.instru,filt_names:sed.filter,wav:wavs, values:sed.spec,sigma:sed.error} 
  obs_st={instru_names:sed.instru,filt_names:sed.filter,wav:wavs, values:sed.StokesI,sigma:sqrt(sed.sigmaII)} 
   ;=== fill !dustem_data
  !dustem_data.sed = ptr_new(obs_st)
  ;VGb
  ; If f_HI is specified, multiply the data by f_HI
  if keyword_set(f_HI) then begin
	if f_HI gt 0 then (*!dustem_data.sed).values *= f_HI else f_HI = 1.
  endif else f_HI = 1.
  defsysv,'!dustem_f_HI',f_HI
  ;VGe
ENDIF

;===== not sure if this below is vere happening anymore
IF keyword_set(polsed) THEN BEGIN
  wavs=polsed.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(polsed.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(polsed(ind).filter)
  ENDIF
  ;=== define observations structure
  ;obs_st={instru_names:polsed.instru,filt_names:polsed.filter,wav:wavs, values:polsed.spec,sigma:polsed.error} 
  obs_st={instru_names:polsed.instru,filt_names:polsed.filter,wav:wavs, values:polsed.largeP,sigma:sqrt(sed.sigma_largeP)} 
  ;=== fill !dustem_data
  !dustem_data.polsed = ptr_new(obs_st)
  ;VGe
  ; If f_HI is specified, multiply the data by f_HI
  if keyword_set(f_HI) then if f_HI gt 0 then $
	(*!dustem_data.polsed).values *= f_HI
  ;VGe
ENDIF


;IF keyword_set(qsed) THEN BEGIN
IF NQ NE 0 THEN BEGIN
  qsed=sed[indQ]
  wavs=qsed.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(qsed.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs[ind]=dustem_filter2wav(qsed[ind].filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:qsed.instru,filt_names:qsed.filter,wav:wavs, values:qsed.StokesQ,sigma:sqrt(qsed.SIGMAQQ)} 
  ;=== fill !dustem_data
  !dustem_data.qsed = ptr_new(obs_st)
  ;VGe
  ; If f_HI is specified, multiply the data by f_HI
  if keyword_set(f_HI) then if f_HI gt 0 then $
	(*!dustem_data.qsed).values *= f_HI
  ;VGe
ENDIF

IF NU NE 0 THEN BEGIN
;IF keyword_set(used) THEN BEGIN
  used=sed[indU]
  wavs=used.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(used.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(used(ind).filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:used.instru,filt_names:used.filter,wav:wavs, values:used.StokesU,sigma:sqrt(qsed.SIGMAUU)} 
  ;=== fill !dustem_data
  !dustem_data.used = ptr_new(obs_st)
  ;VGe
  ; If f_HI is specified, multiply the data by f_HI
  if keyword_set(f_HI) then if f_HI gt 0 then $
	(*!dustem_data.used).values *= f_HI
  ;VGe
ENDIF

IF keyword_set(qext) THEN BEGIN
  wavs=qext.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(qext.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(qext(ind).filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:qext.instru,filt_names:qext.filter,wav:wavs, values:qext.spec,sigma:qext.error} 
  ;=== fill !dustem_data
  !dustem_data.qext = ptr_new(obs_st)
  ;VGe
  ; If f_HI is specified, multiply the data by f_HI
  if keyword_set(f_HI) then if f_HI gt 0 then $
 	(*!dustem_data.qext).values *= f_HI
  ;VGe
ENDIF

IF keyword_set(uext) THEN BEGIN
  wavs=uext.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(uext.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(uext(ind).filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:uext.instru,filt_names:uext.filter,wav:wavs, values:uext.spec,sigma:uext.error} 
  ;=== fill !dustem_data
  !dustem_data.uext = ptr_new(obs_st)
  ;VGe
  ; If f_HI is specified, multiply the data by f_HI
  if keyword_set(f_HI) then if f_HI gt 0 then $
 	(*!dustem_data.uext).values *= f_HI
  ;VGe
ENDIF

;===== not sure if this below is vere happening anymore
IF keyword_set(polfrac) THEN BEGIN
  wavs=polfrac.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(polfrac.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(polfrac(ind).filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:polfrac.instru,filt_names:polfrac.filter,wav:wavs, values:polfrac.spec,sigma:polfrac.error} 
  ;=== fill !dustem_data
  !dustem_data.polfrac = ptr_new(obs_st)
ENDIF

IF keyword_set(polext) THEN BEGIN
   wavs=polext.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(polext.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(polext(ind).filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:polext.instru,filt_names:polext.filter,wav:wavs, values:polext.spec,sigma:polext.error} 
                                ;=== fill !dustem_data
  !dustem_data.polext = ptr_new(obs_st)
ENDIF

IF keyword_set(ext) THEN BEGIN
  wavs=ext.wave
  ;=== Impose central wavelengths for photometric channels
  ind=where(ext.filter NE 'SPECTRUM',count)
  IF count NE 0 THEN BEGIN
    wavs(ind)=dustem_filter2wav(ext(ind).filter)
  ENDIF
  ;=== define observations structure
  obs_st={instru_names:ext.instru,filt_names:ext.filter,wav:wavs, values:ext.spec,sigma:ext.error} 
  ;=== fill !dustem_data
  !dustem_data.ext = ptr_new(obs_st)
ENDIF


!fit_rchi2_weight.sed=1.
!fit_rchi2_weight.ext=1.
if keyword_set(polext) then !fit_rchi2_weight.polext=1.
if keyword_set(polsed) then !fit_rchi2_weight.polsed=1.
if keyword_set(polfrac) then !fit_rchi2_weight.polfrac=1.
if keyword_set(qsed) then !fit_rchi2_weight.qsed=1.
if keyword_set(used) then !fit_rchi2_weight.used=1.
if keyword_set(qext) then !fit_rchi2_weight.qext=1.
if keyword_set(uext) then !fit_rchi2_weight.uext=1.

IF keyword_set(rchi2_weight) THEN begin
	!fit_rchi2_weight.sed     = rchi2_weight.sed
	!fit_rchi2_weight.ext     = rchi2_weight.ext
	if keyword_set(polext) then !fit_rchi2_weight.polext  = rchi2_weight.polext
	if keyword_set(polsed) then !fit_rchi2_weight.polsed  = rchi2_weight.polsed
	if keyword_set(polfrac) then !fit_rchi2_weight.polfrac = rchi2_weight.polfrac
    if keyword_set(qsed) then !fit_rchi2_weight.qsed = rchi2_weight.qsed
    if keyword_set(used) then !fit_rchi2_weight.used = rchi2_weight.used
    if keyword_set(qext) then !fit_rchi2_weight.qext = rchi2_weight.qext
    if keyword_set(uext) then !fit_rchi2_weight.uext = rchi2_weight.uext
endif

return,ptr_new(obs_st)


the_end:

END