dustem_init.pro 17.1 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464
PRO dustem_init,model=model $
;                ,dir=dir $
                ,wraptest=wraptest $
                ,plot_it=plot_it $
                ,noerase=noerase $
                ,st_model=st_model $
                ,polarization=polarization $
                ,verbose=verbose $
                ,help=help $
                ,kwords=kwords

;+
; NAME:
;   dustem_init
; PURPOSE:
;   Defines system variables to run the DUSTEM code
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, Initialization
; CALLING SEQUENCE:
;   dustem_init,model=model,
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    ;   st_model *** COMMENT AH -- what is this? END COMMENT AH ***
; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM
;           'MC10' model from Compiegne et al 2010 (default)
;           'DBP90' model from Desert et al 1990
;           'DL01' model from Draine & Li 2001
;           'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5
;           'DL07' model from Draine & Li 2007
;           'J13' model from Jones et al 2013, as updated in
;                 Koehler et al 2014
;           'G17_ModelA' model A from Guillet et al (2018). Includes
;                 polarisation. See Tables 2 and 3 of that paper for details.
;           'G17_ModelB' model B from Guillet et al (2018)
;           'G17_ModelC' model C from Guillet et al (2018)
;           'G17_ModelD' model A from Guillet et al (2018)
;   dir      = directory where to read data (overrides the default)
;   wraptest = test running dustem f90 through the wrapper
;   plot_it  = plots result of wraptest (ignored if wraptest not set)
;   help     = if set, print this help
;   verbose     = if set, output code diagnostic info (default OFF)
;   polarization *** COMMENT AH -- what is this? END COMMENT AH ***
;   noerase *** COMMENT AH -- what is this? END COMMENT AH ***
;   kwords *** COMMENT AH -- what are they? END COMMENT AH ***
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    1/ The following system variables are initialized to null values
;	or empty pointers:
;       !dustem_fit
;       !dustem_data
;       !dustem_inputs
;       !dustem_params
;       !run_ionfrac
;       !dustem_idl_continuum
;       !dustem_filters
;       !dustem_verbose
;       !dustem_show_plot
;    2/ optionally runs the fortran code (if wraptest)
;       This step calls the Fortran code to run a standard model
;    3/ !dustem_filters is initialized calling dustem_filter_init.pro.
; RESTRICTIONS:
;    The dustem fortran code must be installed
;    The dustem idl wrapper must be installed
; PROCEDURES AND SUBROUTINES USED:
;    *** COMMENT AH --> is this really NONE? ****
; EXAMPLES
;    dustem_init,model="MC10",/wrap,/plot_it
; MODIFICATION HISTORY:
;    Written by JPB,NF,DP Jan-2007
;    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_init'
  goto,the_end
ENDIF

dustem_define_la_common 

defsysv,'!dustem_verbose',0             ;1=verbose
IF keyword_set(verbose) THEN !dustem_verbose=1

defsysv,'!run_pol',0             ;0=no polar
IF keyword_set(polarization) THEN !run_pol=1

;=== color correction system variables
defsysv,'!dustem_previous_cc',ptr_new()   ;store previous color correction values
defsysv,'!dustem_do_cc',1   ;indicates to dustem if color corrections are needed. Will be changed when fitting
defsysv,'!dustem_never_do_cc',0   ;set to 1 to never do color corrections (not recommanded, but faster).

;=== polarization system variables
defsysv, '!run_circ', 0.  
defsysv, '!run_anis', 0.  
defsysv, '!run_rrf', 0.    
defsysv, '!run_univ', 0.    
defsysv, '!run_lin', 0.

defsysv, '!dustem_redshift', 0.    ;in case the SED needs to be shifted for a non zero redshift

;This is just to make sure that plugin system variables are created even if no plugins are used
pd_bidon=['NONE']
dustem_init_plugins,pd_bidon

defsysv,'!dustem_nocatch',0   ;set to catch errors in SED fitting.


;=== read info about instrument
file=!dustem_wrap_soft_dir+'instrument_description.xcat'
st=read_xcat(file,/silent)
defsysv,'!dustem_instrument_description',st

defsysv,'!indef', exists=exists_indef
IF NOT exists_indef THEN defsysv,'!indef',-32768.

;define !dustem_fit
dustem_fit_st={data:ptr_new(),wavelength:ptr_new(), $
               param_descs:ptr_new(),param_init_values:ptr_new(),param_func:ptr_new(), $
               fixed_param_descs:ptr_new(),fixed_param_init_values:ptr_new(), $
               chi2:0.D0,rchi2:0.D0, $
               current_param_values:ptr_new(), current_param_errors:ptr_new() $
               }
defsysv, '!dustem_fit', ptr_new(dustem_fit_st)    ;Data to fit

if keyword_set(polarization) then begin
	defsysv, '!dustem_data', {    $ ;Data to fit
                                sed:    ptr_new(),      $ ; 
                                ext:    ptr_new(),      $ ;
                                qsed: ptr_new(),        $ ; 
                                used: ptr_new(),        $ ; 
                                polsed: ptr_new(),      $ ; 
                                polfrac: ptr_new(),     $ ;   
                                psi_em: ptr_new(),      $ ; 
                                qext: ptr_new(),        $ ; 
                                uext: ptr_new(),        $ ;
                                polext: ptr_new(),      $ ;  
                                psi_ext: ptr_new()      $ ; 
                                 
                         }
    
         
         rchi2_weight = {sed: 0.,ext: 0.,qsed:0.,used:0.,polsed:0.,polfrac:0.,psi_em:0.,qext:0.,uext:0.,polext:0.,psi_ext:0.}
         
    defsysv,'!dustem_show', {    $ ;Data to fit
                                sed:    ptr_new(),      $ ; 
                                ext:    ptr_new(),      $ ;
                                qsed: ptr_new(),        $ ; 
                                used: ptr_new(),        $ ; 
                                polsed: ptr_new(),      $ ; 
                                polfrac: ptr_new(),     $ ;   
                                psi_em: ptr_new(),      $ ; 
                                qext: ptr_new(),        $ ; 
                                uext: ptr_new(),        $ ;
                                polext: ptr_new(),      $ ;  
                                psi_ext: ptr_new()      $ ; 
                                 
                         }     
         
endif else begin
	
	defsysv, '!dustem_data', {    $ ;Data to fit
                                sed:    ptr_new(),      $
                                ext:    ptr_new()       $
                         }
	rchi2_weight = {sed: 0.,ext: 0.}
	
	defsysv, '!dustem_show', {    $ ;Data to fit
                                sed:    ptr_new(),      $
                                ext:    ptr_new()       $
                         }
endelse

if !run_pol then begin
    tagnames=tag_names(!dustem_data)
    testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'POLSED' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT'
    ind_data=where(~testpol,ctestpol)
    if ctestpol ne 0 then tagnames =(tag_names(!dustem_data))[ind_data] 

endif else tagnames=tag_names(!dustem_data)

instr="defsysv, '!fit_rchi2_weight', {"
FOR i = 0, n_elements(tagnames)-1 DO BEGIN
  instr+=tagnames(i)+': rchi2_weight.'+tagnames(i)
  IF i NE n_elements(tagnames)-1 THEN instr+=','
ENDFOR
instr+='}'
toto=execute(instr)

;#########system variables necessary for cgwindow plotting##########     
defsysv, '!dustemcgwin_id', {    $ ;IDs of windows to plot                  
                            sed:    la_undef(),      $                          
                            ext:    la_undef(),      $
                            prms:   la_undef(),      $
                            plgns:  la_undef()       $
                     }                                                      

defsysv, '!dustemcgwin_ncmds', {    $ ;Data to fit
                            sed:     {pl:0, nrm:0},      $ ; 
                            ext:     {pl:0, nrm:0},      $ ;
                            qsed:    {pl:0, nrm:0},      $ ; 
                            used:    {pl:0, nrm:0},      $ ; 
                            polsed:  {pl:0, nrm:0},      $ ; 
                            polfrac: {pl:0},             $ ;   
                            psi_em:  {pl:0},             $ ; 
                            qext:    {pl:0, nrm:0},      $ ; 
                            uext:    {pl:0, nrm:0},      $ ;
                            polext:  {pl:0, nrm:0},      $ ;  
                            psi_ext: {pl:0},             $ ;
                            chi2 :   {pl:0, txtwdth:0.}, $ ; ;beacause chi2 needs to be updated (needing the command index)
                            rchi2 :  {pl:0, txtwdth:0.}, $ ;
                            pltit :  {pl:0, txtwdth:0.}, $ ; plot title
                            prms :   {pl:0},             $ ;
                            plgns :  {pl:0}             $ ;
                     }

defsysv, '!dustem_mlog', 0      ;necessary for the plotting of negative values on the axis in dstmwrp_exp.pro

;###################################################################


defsysv, '!dustem_psi', 0.


                     
;IDL> help,*!dustem_data,/str
;** Structure <2173b20>, 4 tags, length=448, data length=448, refs=1:
;   INSTRU_NAMES    STRING    Array[14]
;   FILT_NAMES      STRING    Array[14]
;   VALUES          FLOAT     Array[14]
;   SIGMA           FLOAT     Array[14]
defsysv,'!dustem_inputs',ptr_new()    ;dustem input files

;=== This is not a system variable anymore
;defsysv, '!dustem_ext', ptr_new()    ;Contains the values of extinction
;=== This is not a system variable anymore
;defsysv, '!dustem_specem', ptr_new() ;Contains the values of emission

defsysv, '!dustem_params', ptr_new() ;Contains the values of all Desert Model parameters (except extinction !)
;IDL> help,*!dustem_params,/str
;** Structure <2926e10>, 7 tags, length=13312, data length=13282, refs=1:
;   GRAINS          STRUCT    -> <Anonymous> Array[3]
;   GEMISSIV        STRUCT    -> <Anonymous> Array[300]
;   ISRF            STRUCT    -> <Anonymous> Array[64]
;   QABS            STRUCT    -> <Anonymous> Array[64]
;   CALOR           STRUCT    -> <Anonymous> Array[5]
;   SPECEM          STRUCT    -> <Anonymous> Array[1]
;   IONFRAC         STRUCT    -> <Anonymous> Array[10]

;defsysv, '!dustem_fit_sed',ptr_new()          ;Best fit SED

;The following is not used in the RELEASE version of the code, only in VERSTRAETE
;IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
  defsysv, '!run_ionfrac', 0.     ; 0:no call to dustem_create_ionfrac UNLESS it is called as a param
                                  ; 1:call to dustem_create_ionfrac with the MIX_xxx.DAT files
                                  ; 2:call to dustem_create_ionfrac with the SIZE_xxx.DAT files
;ENDIF

defsysv,'!dustem_filters',ptr_new()     ;filter information (filled by dustem_filter_init.pro)
;; IDL> help,*!dustem_filters,/str
;; ** Structure <1b95608>, 17 tags, length=8232, data length=8118, refs=1:
;;    IRAC            STRUCT    -> <Anonymous> Array[1]
;;    MIPS            STRUCT    -> <Anonymous> Array[1]
;;    MSX             STRUCT    -> <Anonymous> Array[1]
;;    IRAS            STRUCT    -> <Anonymous> Array[1]
;;    DIRBE           STRUCT    -> <Anonymous> Array[1]
;;    SPM             STRUCT    -> <Anonymous> Array[1]
;;    ISOCAM          STRUCT    -> <Anonymous> Array[1]
;;    ISOPHOTP        STRUCT    -> <Anonymous> Array[1]
;;    ISOPHOTC        STRUCT    -> <Anonymous> Array[1]
;;    EFIRAS          STRUCT    -> <Anonymous> Array[1]
;;    ARCHEOPS        STRUCT    -> <Anonymous> Array[1]
;;    HFI             STRUCT    -> <Anonymous> Array[1]
;;    LFI             STRUCT    -> <Anonymous> Array[1]
;;    WMAP            STRUCT    -> <Anonymous> Array[1]
;;    SPIRE           STRUCT    -> <Anonymous> Array[1]
;;    PACS            STRUCT    -> <Anonymous> Array[1]
;;    AKARI           STRUCT    -> <Anonymous> Array[1]

defsysv,'!dustem_verbose',0                                  ;1=verbose
defsysv,'!dustem_show_plot',1           ;0=show no plot


;apparently not needed anymore. Was removed.
;defsysv,'!ki',0.  ;apparently needed by mpfit, but what is it ?

;== This is for the default mode
dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'}
!dustem_inputs=ptr_new(dustem_inputs)


defsysv, '!dustem_HCD', ptr_new(1.00E20) ;I do not see the reason behind using a pointer in this case.


;===This is the Desert option

IF keyword_set(model) THEN BEGIN

defsysv, '!dustem_model', model


  CASE model OF
    'MC10':BEGIN
      (*!dustem_inputs).grain='GRAIN_MC10.DAT'
    END
    'DBP90':BEGIN
      (*!dustem_inputs).grain='GRAIN_DBP90.DAT'
    END
    'DL01':BEGIN
      (*!dustem_inputs).grain='GRAIN_DL01.DAT'
    END
    'WD01_RV5p5B':BEGIN
      (*!dustem_inputs).grain='GRAIN_WD01_Rv5.5B.DAT'
    END
    'DL07':BEGIN
      (*!dustem_inputs).grain='GRAIN_DL07.DAT'
    END
    'J13':BEGIN
      (*!dustem_inputs).grain='GRAIN_J13.DAT'
    END
    'G17_MODELA':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELA.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT'
    END
    'G17_MODELB':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELB.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELB.DAT'
    END
    'G17_MODELC':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELC.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELC.DAT'
    END
    'G17_MODELD':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELD.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELD.DAT'
    END
    'THEMIS':BEGIN
      (*!dustem_inputs).grain='GRAIN_THEMIS.DAT'
      (*!dustem_inputs).align='ALIGN_THEMIS.DAT'
    END
    'NY_MODELA':BEGIN ;TB removed?
      (*!dustem_inputs).grain='GRAIN_NY_MODELA.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT' 
    END
    ;This is a user-defined model ;This will probably need to be modified.    
    'USER_MODEL':BEGIN ;TB removed? 
      (*!dustem_inputs).grain='GRAIN_USER.DAT' 
    END
    ELSE:BEGIN
      message,'model '+model+' unknown',/continue
      message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
      stop
    END
  ENDCASE
ENDIF

;Initialize !dustem_params structure
IF not keyword_set(dir) THEN BEGIN
;  CASE getenv('DUSTEM_WHICH') OF
  CASE !dustem_which OF
    'DESERT':dir_in=!dustem_wrap_soft_dir+'/Data/DESERT_POST_ISO_MORE/'
    'COMPIEGNE': dir_in=!dustem_wrap_soft_dir+'/Data/MC_DAT/'
    'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/LV_DAT/'
    'RELEASE': dir_in=!dustem_soft_dir
    ELSE: dir_in=!dustem_wrap_soft_dir+'/Data/les_DAT/';!dustem_wrap_soft_dir+'/Data/les_DAT/'
  ENDCASE
ENDIF ELSE BEGIN
  dir_in=dir
ENDELSE
;stop
;dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/'+dir
;dir_in=!dustem_soft_dir+'/src/'+dir

if keyword_set(kwords) then defsysv, '!dustem_kwords', ptr_new(kwords) else defsysv, '!dustem_kwords', ptr_new()

st_model=dustem_read_all(dir_in,/silent)
!dustem_params=ptr_new(st_model)

;=== create dynamical storage if needed
IF not file_test(!dustem_dat,/dir) THEN BEGIN
  spawn,'mkdir '+!dustem_dat
ENDIF
IF not file_test(!dustem_res,/dir) THEN BEGIN
  spawn,'mkdir '+!dustem_dat
ENDIF

;=== remove dynamical storage directories
;=== This is needed for iterative use of DUSTEM from IDL
if not keyword_set(noerase) then begin
CASE !dustem_which OF
  'VERSTRAETE':BEGIN
    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'rm -rf '+!dustem_dat+'/les_DAT'
      spawn,'rm -rf '+!dustem_dat+'/les_QABS'
      spawn,'rm -rf '+!dustem_dat+'/les_CAPA'
    ENDIF
    IF file_test(!dustem_res,/dir) THEN BEGIN
      spawn,'rm -rf '+!dustem_res+'/les_RES'
    ENDIF
  END
  'RELEASE':BEGIN 
    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'rm -rf '+!dustem_dat+'/data'
      spawn,'rm -rf '+!dustem_dat+'/hcap'
      spawn,'rm -rf '+!dustem_dat+'/oprop'
      spawn,'rm -rf '+!dustem_dat+'/out'
      ;spawn,'rm -rf '+!dustem_dat+'/plots' TO BE TESTED. NOT SURE ABOUT WHAT I AM DOING HERE. ASK JP ABOUT THIS
    ENDIF
  END
  ELSE: BEGIN
  END
ENDCASE

;=== create storage directories
;=== This is needed for iterative use of DUSTEM from IDL
CASE !dustem_which OF
  'VERSTRAETE':BEGIN
    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'mkdir '+!dustem_dat+'/les_DAT'
      spawn,'mkdir '+!dustem_dat+'/les_QABS'
      spawn,'mkdir '+!dustem_dat+'/les_CAPA'
    ENDIF
    IF file_test(!dustem_res,/dir) THEN BEGIN
      spawn,'mkdir '+!dustem_res+'/les_RES'
    ENDIF
  END
  'RELEASE':BEGIN  ; IN CHANGED THIS AS WELL
    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'mkdir '+!dustem_dat+'/data'
      spawn,'mkdir '+!dustem_dat+'/hcap'
      spawn,'mkdir '+!dustem_dat+'/oprop'
      spawn,'mkdir '+!dustem_dat+'/out'
      ;spawn,'mkdir '+!dustem_dat+'/plots' TO BE TESTED
    ENDIF
  END
  ELSE: BEGIN
  END
ENDCASE
endif ; /noerase

;=== If needed, test runing the model through the wrapper
IF keyword_set(wraptest) THEN BEGIN
  dir_out=!dustem_dat
  dustem_write_all,st_model,dir_out
  st=dustem_run()
  IF keyword_set(plot_it) THEN BEGIN
    loadct,13
    dustem_plot_nuinu_em,st,yr=[1e-14,1.e-6],/ysty,xr=[1,5e3],/xsty
  ENDIF
ENDIF

;Initialize filters
dustem_filter_init

the_end:

END