dustem_read_all.pro 10 KB
FUNCTION dustem_read_all,dir_in,silent=silent,help=help

;+
; NAME:
;    dustem_read_all
; PURPOSE:
;    Reads all Dustem input data
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    st=dustem_read_all(dir_in,silent=silent,help=help)
; INPUTS:
;    dir_in: input directory
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    st: input structures
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    silent    = If set, function is silent
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_init
;    dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/'
;    st=dustem_read_all(dir)
; 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.
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_read_all'
  st=0.
  goto,the_end
ENDIF

;IF getenv('DUSTEM_WHICH') EQ 'WEB3p8' THEN BEGIN
IF !dustem_which EQ 'WEB3p8' THEN BEGIN     
  st=dustem_read_all_web3p8(dir_in,silent=silent)
  goto,ionfrac
ENDIF

;IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
  st=dustem_read_all_lv(dir_in,silent=silent)
  goto,ionfrac
ENDIF

;== note: the output keywords are set only for LV version
file=dir_in+'GRAIN.DAT'
st_grains=dustem_read_grain(file,silent=silent)

;stop

file=dir_in+'TAILLE.DAT' ;I put this back to its original form 
st_taille=dustem_read_taille(file,silent=silent)

file=dir_in+'PROPMASS.DAT'
st_propmass=dustem_read_propmass(file,silent=silent)

;=== This does not exist in MC's new version
;IF getenv('DUSTEM_WHICH') EQ 'DESERT' THEN BEGIN
IF !dustem_which EQ 'DESERT' THEN BEGIN
  file=dir_in+'GEMISSIV.DAT'
  st_gemissiv=dustem_read_gemissiv(file,silent=silent)
ENDIF

file=dir_in+'ISRF.DAT'
st_isrf=dustem_read_isrf(file,silent=silent,Nisrf=Nisrf)

;stop

;=== This does not exist in XD's new version
;IF getenv('DUSTEM_WHICH') EQ 'COMPIEGNE' THEN BEGIN
IF !dustem_which EQ 'COMPIEGNE' THEN BEGIN
  Qabs_files=findfile(dir_in+'QABS_*.DAT',count=Nfiles)
  st_qabs=ptrarr(Nfiles)
  FOR i=0,Nfiles-1 DO BEGIN
    st=dustem_read_qabs(Qabs_files(i),silent=silent)
    st_qabs(i)=ptr_new(st)
  ENDFOR
ENDIF
;IF getenv('DUSTEM_WHICH') EQ 'DESERT' THEN BEGIN
IF !dustem_which EQ 'DESERT' THEN BEGIN
  file=dir_in+'QABS.DAT'
  st_qabs=dustem_read_qabs_desert(file,silent=silent)
ENDIF


file=dir_in+'SPECEM.DAT'
st_specem=dustem_read_specem(file,silent=silent)

file=dir_in+'CALOR.DAT'
st_calor=dustem_read_calor(file,silent=silent)

;=== This does not exist in MC's new version
;IF getenv('DUSTEM_WHICH') EQ 'DESERT' THEN BEGIN
IF !dustem_which EQ 'DESERT' THEN BEGIN
  file=dir_in+'IONFRAC.DAT'
  st_ionfrac=dustem_read_ionfrac(file,silent=silent)
ENDIF

;=== This does not exist in XD's new version
;IF getenv('DUSTEM_WHICH') EQ 'COMPIEGNE' THEN BEGIN
IF !dustem_which EQ 'COMPIEGNE' THEN BEGIN
  file=dir_in+'LAMBDA.DAT'
  st_lambda=dustem_read_lambda(file,silent=silent)
ENDIF


Ngrains=n_elements(st_grains)

one_stt_grains={ $
type:st_grains(0).type, $
albmax:st_grains(0).albmax, $
densgr:st_grains(0).densgr, $
ngem:st_grains(0).ngem, $
hydro:st_grains(0).hydro, $
ioni:st_grains(0).ioni, $
alpha:st_taille(0).alpha, $
tmin:st_taille(0).tmin, $
tmax:st_taille(0).tmax, $
ndiscr:st_taille(0).ndiscr, $
propmass:st_propmass(0).propmass}

sst_grains=replicate(one_stt_grains,Ngrains)

sst_grains.type=st_grains.type
sst_grains.albmax=st_grains.albmax
sst_grains.densgr=st_grains.densgr
sst_grains.ngem=st_grains.ngem
sst_grains.hydro=st_grains.hydro
sst_grains.ioni=st_grains.ioni
sst_grains.alpha=st_taille.alpha
sst_grains.tmin=st_taille.tmin
sst_grains.tmax=st_taille.tmax
sst_grains.ndiscr=st_taille.ndiscr
sst_grains.propmass=st_propmass.propmass


;IF getenv('DUSTEM_WHICH') EQ 'DESERT' THEN BEGIN
IF !dustem_which EQ 'DESERT' THEN BEGIN
  st={Ngrains:Ngrains,grains:sst_grains, $
      gemissiv:st_gemissiv,isrf:st_isrf,qabs:st_qabs,calor:st_calor,specem:st_specem,ionfrac:st_ionfrac}
ENDIF
;IF getenv('DUSTEM_WHICH') EQ 'COMPIEGNE' THEN BEGIN
IF !dustem_which EQ 'COMPIEGNE' THEN BEGIN
  st={Ngrains:Ngrains,grains:sst_grains, $
      isrf:st_isrf,qabs:st_qabs,calor:st_calor,specem:st_specem,lambda:st_lambda}
ENDIF
;IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
  st={Ngrains:Ngrains,grains:sst_grains, $
      isrf:st_isrf,qabs:st_qabs,calor:st_calor,specem:st_specem,lambda:st_lambda}
ENDIF

ionfrac:

;=== JPB added tests on !dustem_verbose for messages
;stop
;NF on October 2009 : test IONFRACPAH keyword (only if 'VERSTRAETE' model)
;IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
   IF !dustem_verbose NE 0 THEN message, "============= IONization FRACtion of the PAHs =============", /info
;Get master-keywords
   master_keywords = strsplit(st.keywords,' ',/regex,/extr)
;check presence of IONFRACPAH in master-keywords
   tmp = where(master_keywords eq 'IONFRACPAH',test_ionfracpah)
;check presence of dustem_create_ionfrac in the *!param_desc
   ;JPB: added test on validity of param_desc
   IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN
     tmp2 = total(stregex(*(*!dustem_fit).param_descs,'ionfrac',/bool,/fold_case))
   ENDIF ELSE BEGIN
     tmp2=0
   ENDELSE
   test_dustemcreateionfrac = (tmp2 GT 0) ? 1 : 0
;if IONFRACPAH or "dustem_create_ionfrac" is present ...
   IF test_ionfracpah or test_dustemcreateionfrac THEN BEGIN
      IF !dustem_verbose NE 0 THEN message, "You want the IONization FRACtion of the PAHs to be computed by the model", /info
;  Check if there are both PAH1 and PAH0 provided
      pah0 = (where(strmid(st.grains.type,0,4) eq 'PAH0', count0))(0)
      pah1 = (where(strmid(st.grains.type,0,4) eq 'PAH1', count1))(0)
      IF count0 eq count1 and count0 ne 0 THEN BEGIN
         IF !dustem_verbose NE 0 THEN BEGIN
           message, "and you have provided me with both PAH0 and PAH1 components.", /info
           message, "I will calculate the fractions.", /info
           message, "I will assume you have forced them to keep the same abundance (initial values + tied).", /info
         ENDIF
      ENDIF ELSE BEGIN
         IF !dustem_verbose NE 0 THEN BEGIN
           message, "but you have NOT provided me with both PAH0 and PAH1 components.", /info
           message, "Please correct and rerun the code."
         ENDIF
      ENDELSE
;  Check presence of SIZE in master-keywords
      tmp = where(master_keywords eq 'SIZE',test_size)
;  Check if MIX is a keyword provided in GRAIN.DAT for PAH0 and PAH1
      pah0_keyword = st.grains(pah0).flag
      pah1_keyword = st.grains(pah1).flag
      tmp = where(pah0_keyword eq 'MIX',test_mix_pah0)
      tmp = where(pah1_keyword eq 'MIX',test_mix_pah1)
; if SIZE is not present ...
      IF not test_size THEN BEGIN
;   if MIX is not present for both PAH components ...
         IF not test_mix_pah0 or not test_mix_pah1 THEN BEGIN
            IF !dustem_verbose NE 0 THEN BEGIN
              message, "You did not provide me with neither SIZE master-keyword nor MIX keyword for both PAHs components.", /info
              message, "I will add the MIX keyword for both PAHs components and generate 50%-50% MIX_xxx.DAT files.", /info
              message, "I will assume you have forced them to have the same size distributions.", /info
            ENDIF
;      Add MIX keyword to both PAH components (had to concatenate them
;      to order the F90 code to read them as 1 column)
            if not test_mix_pah0 then $
               st.grains(pah0).flag = st.grains(pah0).flag+'MIX'
            if not test_mix_pah1 then $
               st.grains(pah1).flag = st.grains(pah1).flag+'MIX'
;      Update structure and generate MIX_xxx.DAT files
;            full_st0 = ptr_new({file:getenv('DUSTEM_DAT')+'/MIX_'+st.grains(pah0).type+'.DAT',fmix:fltarr(st.grains(pah0).ndiscr)+0.5})
;            full_st1 = ptr_new({file:getenv('DUSTEM_DAT')+'/MIX_'+st.grains(pah1).type+'.DAT',fmix:fltarr(st.grains(pah1).ndiscr)+0.5})
            full_st0 = ptr_new({file:!dustem_dat+'/MIX_'+st.grains(pah0).type+'.DAT',fmix:fltarr(st.grains(pah0).ndiscr)+0.5})
            full_st1 = ptr_new({file:!dustem_dat+'/MIX_'+st.grains(pah1).type+'.DAT',fmix:fltarr(st.grains(pah1).ndiscr)+0.5})
            full_st = ptrarr(2) & full_st(0)=full_st0 & full_st(1)=full_st1
;            dustem_write_mix, getenv('DUSTEM_DAT'), full_st
            dustem_write_mix,!dustem_dat, full_st
            st.mix(pah0) = ptr_new(dustem_read_mix((*full_st0).file))
            st.mix(pah1) = ptr_new(dustem_read_mix((*full_st1).file))
            !run_ionfrac = 1.
         ENDIF ELSE BEGIN
            IF !dustem_verbose NE 0 THEN BEGIN
              message, "You provided me with the MIX keyword for both PAHs components.", /info
              message, "I will use your MIX_xxx.DAT files to update the fractions.", /info
            ENDIF
            !run_ionfrac = 1.
         ENDELSE
; if SIZE is present
      ENDIF ELSE BEGIN
         IF !dustem_verbose NE 0 THEN BEGIN
           message, "You provided me with the SIZE master-keyword.", /info
           message, "I will use your SIZE_xxx.DAT files to update the fractions.", /info
           message, "I will assume you have forced them to have the same size distributions.", /info
         ENDIF
         !run_ionfrac = 2.
      ENDELSE
;if IONFRACPAH is NOT present ...
   ENDIF ELSE BEGIN
      IF !dustem_verbose NE 0 THEN BEGIN
        message, "You do not want the IONization FRACtion of the PAHs to be computed by the model", /info
      ENDIF
      !run_ionfrac = 0.
   ENDELSE
   
   IF test_dustemcreateionfrac THEN BEGIN
      !run_ionfrac = -1. * !run_ionfrac
   ENDIF
   IF !dustem_verbose NE 0 THEN BEGIN
     message, "============= IONization FRACtion of the PAHs =============", /info
     wait, 5                      ;might be a bit long, but I want to be sure the user sees it
   ENDIF
  ;JPB: The above is stupid, specially before I added the !dustem_verbose test
ENDIF
;NF on October 2009 : test IONFRACPAH keyword

the_end:

RETURN,st

END