FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI

;+
; 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
;-



;fully deprecate the use of polfrac?

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

;=== clean the data structure
;This step detects polarization without the polarization keyword of dustem_init

!dustem_data.sed = ptr_new()
!dustem_data.ext = ptr_new()

NP=0
NQ=0
NU=0
IF !run_pol THEN BEGIN
  ind=where(sed.largeP NE la_undef(),NP)
  if NP NE 0 then !dustem_data.polsed = ptr_new()

  ind=where(sed.StokesQ NE la_undef(),NQ)
  if NQ NE 0 then !dustem_data.qsed = ptr_new()

  ind=where(sed.StokesU NE la_undef(),NU)
  if NU NE 0 then !dustem_data.used = ptr_new()
ENDIF

;NOTA BENE: if we want to make a similar change to the ext structure we have to find a name 
;that suits polext and the cos and sin components

; ind=where(ext.largeP NE la_undef(),NPext)
; if NPext NE 0 then !dustem_data.polext = ptr_new()
; ind=where(sed.StokesQext NE la_undef(),NQext)
; if NQext NE 0 then !dustem_data.qext = ptr_new()
; ind=where(sed.StokesUext NE la_undef(),NUext)
; if NUext NE 0 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.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

    if NP NE 0 then BEGIN
        ;=== define observations structure
        obs_st.values=sed.largeP
        obs_st.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 NQ NE 0 then BEGIN
        ;=== define observations structure
        ;obs_st={instru_names:polsed.instru,filt_names:polsed.filter,wav:wavs, values:polsed.spec,sigma:polsed.error} 
        obs_st.values=sed.StokesQ
        obs_st.sigma=sqrt(sed.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
        ;=== define observations structure
        ;obs_st={instru_names:polsed.instru,filt_names:polsed.filter,wav:wavs, values:polsed.spec,sigma:polsed.error} 
        obs_st.values=sed.StokesU
        obs_st.sigma=sqrt(sed.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


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

;I'm not so sure about the extinction


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 NP NE 0 THEN !fit_rchi2_weight.polsed=1.
IF NQ NE 0 THEN !fit_rchi2_weight.qsed=1.
IF NU NE 0 THEN !fit_rchi2_weight.used=1.

;needs the rest

; 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 NP NE 0 THEN !fit_rchi2_weight.polsed  = rchi2_weight.polsed
	IF NQ NE 0 THEN !fit_rchi2_weight.qsed  = rchi2_weight.qsed
	IF NU NE 0 THEN !fit_rchi2_weight.used  = rchi2_weight.used
	
	;you need to add the rest too
	
; 	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