PRO dustem_init,model=model $ ,wraptest=wraptest $ ,plot_it=plot_it $ ,noerase=noerase $ ,noobj=noobj $ ,st_model=st_model $ ,polarization=polarization $ ,verbose=verbose $ ,help=help $ ,grain_keywords=grain_keywords ;+ ; NAME: ; dustem_init ; ; PURPOSE: ; Defines system variables to run the DUSTEM code ; ; CATEGORY: ; DustEMWrap, Distributed, HighLevel, Initialization ; ; CALLING SEQUENCE: ; dustem_init,model=model, ; ; INPUTS: ; None ; ; OPTIONAL INPUT PARAMETERS: ; None ; ; OUTPUTS: ; None ; ; OPTIONAL OUTPUT PARAMETERS: ; ; st_model ; ; 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) ; 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 = set this variables for runs taking polarization into account ; noerase : set this keyword to not erase the temporary DAT files used by dustemwrap ; noobj : set this keyword to turn off object-oriented ; graphical output when using IDL. Object-oriented ; graphics are OFF by default for GDL and Fawlty. ; grain_keywords : Fortran grain keywords. Must be an array of Ngrains strings. Use '?' to keep default values. ; ; 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 DustEMWrap idl code must be installed ; ; PROCEDURES AND SUBROUTINES USED: ; ; 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 ;define which code is running, idl, gdl or fawlty defsysv,'!fl',exist=fawlty defsysv,'!gdl',exist=gdl IF fawlty THEN defsysv,'!dustemwrap_which_language','fawlty' IF gdl THEN defsysv,'!dustemwrap_which_language','gdl' IF not fawlty and not gdl THEN defsysv,'!dustemwrap_which_language','idl' ;Control of object-oriented graphics IF fawlty OR gdl THEN defsysv,'!dustem_noobj',1 ELSE defsysv,'!dustem_noobj',0 IF keyword_set(noobj) then defsysv,'!dustem_noobj',1 ;see if the user has set their dustem_which defsysv,'!dustem_which',exist=dwhich if not dwhich then begin defsysv,'!dustem_which','RELEASE' message,'!dustem_which not known from idl_startup, using RELEASE version',/info end un_mouchard={iteration:0L, $ chi2:0.d0,rchi2:0.d0, $ qchi2:0.d0,qrchi2:0.d0, $ uchi2:0.d0,urchi2:0.d0, $ parameters_descriptions:ptr_new(),parameters_values:ptr_new()} mouchard=ptr_new(replicate(un_mouchard,1)) defsysv,'!iteration_mouchard',mouchard 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 IF not keyword_set(polarization) THEN !run_pol=0 ;=== 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. defsysv,'!indef', exists=exists_indef IF NOT exists_indef THEN defsysv,'!indef',32768. ;=== read info about instrument file=!dustem_wrap_soft_dir+'instrument_description.xcat' st=read_xcat(file,/silent) defsysv,'!dustem_instrument_description',st defsysv, '!dustem_isrf_file', ptr_new() ;but there is only one system variable for this... Is this ok? defsysv, '!dustem_iter', {prv:1,act:1} ;define !dustem_fit dustem_fit_st={data:ptr_new(), $ ;,wavelength:ptr_new(), ;because wavelength arrays in the different datasets (data) can be different 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, $ qchi2:0.D0,qrchi2:0.D0, $ uchi2:0.D0,urchi2:0.D0, $ current_param_values:ptr_new(), current_param_errors:ptr_new() $ } defsysv, '!dustem_fit', ptr_new(dustem_fit_st) ;Data to fit dustem_data_st = dustem_define_dustem_data(polarization=polarization,rchi2_weight =rchi2_weight) defsysv, '!dustem_data', ptr_new(dustem_data_st) defsysv,'!dustem_show', ptr_new(dustem_data_st) 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 'FPOLEXT' 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', {" instr="fit_rchi2_weight_str={" 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) defsysv, '!fit_rchi2_weight',ptr_new(fit_rchi2_weight_str) ;#########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, '!dustem_dim', 0 ;after the new modifications I think that some of these pointers are not needed 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}, $ ; fpolext: {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_m:0.,txtwdth_x:0.}, $ ; plot title prms : {pl:0}, $ ; plgns : {pl:0}, $ ; runs : {pl:0, txtwdth:0.} $ ; } version={version:'V2.0',version_number:2.0} defsysv, '!dustem_version',version ;This is the current dustemwrap version ;NB: #This is definitely not a robust way to handle the _extra tags but it's a primariy workaround for the plotting tags of the data sets. defsysv, '!dustem_plot_range', { $ sed: {xr:[1.,5.0e5], yr:[5e-8,1.00e6]}, $ ; qsed: {xr:[1.,5.0e5], yr:[5e-8,1.00e6]}, $ ; used: {xr:[1.,5.0e5], yr:[5e-8,1.00e6]}, $ ; polsed: {xr:[1.,5.0e5], yr:[5e-8,1.00e6]}, $ ; polfrac: {xr:[1.,5.0e5], yr:[1.00E-04,50.00]}, $ ; psi_em: {xr:[1.,5.0e5], yr:[-90.,90.]}, $ ; ext: {xr:[0.01,30], yr:[5.00E-8,10]}, $ ; qext: {xr:[0.01,30], yr:[5.00E-8,10]}, $ ; uext: {xr:[0.01,30], yr:[5.00E-8,10]}, $ ; polext: {xr:[0.01,30], yr:[5.00E-8,10]}, $ ; fpolext: {xr:[0.01,30], yr:[1.00E-04,10.00]}, $ ; psi_ext: {xr:[0.01,30], yr:[-90.,90.]}, $ ; title_m : 'DUSTEMWRAP '+'('+!dustem_version.version+') '+'Spectral Energy Distribution', $ ; default title if no title is specified title_x : 'DUSTEMWRAP '+'('+!dustem_version.version+') '+'Dust Optical Depth' $ } defsysv, '!dustem_mlog', 0 ;necessary for the plotting of negative values on the axis in dstmwrp_exp.pro ;################################################################### ;These two sysvars allow for the polarization angle (whether in emission or extinction) to ;be updated in procedures (non-plugins) that create polarization data such as ;dustem_compute_stokes and dustem_compute_stokext defsysv, '!dustem_psi', 0. defsysv, '!dustem_psi_ext', 0. ;useful pointers (current dustem output structure and current mpfit params values) defsysv,'!dustem_current',ptr_new() ; defsysv,'!dustem_current_params',ptr_new() ;Is replaced by (*!dustem_fit).CURRENT_PARAM_VALUES ; defsysv, '!dustem_current_errors', ptr_new() ;Is replaced by (*!dustem_fit).current_param_errors=ptr_new(p_er_dim) defsysv, '!dustem_end', 0 ;Easiest and fastest way I could find to signal the end of the fit. If there is a better way please remove this. ;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 -> Array[3] ; GEMISSIV STRUCT -> Array[300] ; ISRF STRUCT -> Array[64] ; QABS STRUCT -> Array[64] ; CALOR STRUCT -> Array[5] ; SPECEM STRUCT -> Array[1] ; IONFRAC STRUCT -> 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 -> Array[1] ;; MIPS STRUCT -> Array[1] ;; ... etc 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) use_model='DEFAULT' 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 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 ;if the user modified the fortran GRAIN.DAT AND ALIGN.DAT files. 'USER_MODEL':BEGIN (*!dustem_inputs).grain='GRAIN.DAT' (*!dustem_inputs).align='ALIGN.DAT' END ELSE:BEGIN message,'model '+model+' unknown',/continue message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue stop END ENDCASE use_model=model ENDIF defsysv, '!dustem_model', use_model ;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 ;stop ;reading fortran files to set up !dustem_params defsysv, '!dustem_kwords', ptr_new() st_model=dustem_read_all(dir_in,/silent) !dustem_params=ptr_new(st_model) Ngrains=(*!dustem_params).Ngrains needed_keywords=strarr(Ngrains) ;replace dustem keywords requested by the user IF keyword_set(grain_keywords) THEN BEGIN IF n_elements(grain_keywords) NE Ngrains THEN BEGIN message,'Grain_keywords provided do not match the number of grains ('+strtrim(Ngrains,2)+') in the model',/continue message,'Please provide '+strtrim(Ngrains,2)+' keywords, using ? for those left to their default value',/continue stop ENDIF ELSE BEGIN FOR i=0L,Ngrains-1 DO BEGIN IF grain_keywords[i] NE '?' THEN BEGIN needed_keywords[i]=grain_keywords[i] ENDIF ELSE BEGIN needed_keywords[i]=(*!dustem_params).grains(i).TYPE_KEYWORDS ENDELSE ENDFOR ENDELSE ENDIF ELSE BEGIN FOR i=0L,Ngrains-1 DO BEGIN needed_keywords[i]=(*!dustem_params).grains(i).TYPE_KEYWORDS ENDFOR ENDELSE ;reading fortran inputs again, but now with replaced keywords defsysv,'!dustem_kwords', ptr_new(needed_keywords) st_model=dustem_read_all(dir_in,/silent) !dustem_params=ptr_new(st_model) ;stop ;=== create dynamical storage if needed IF not file_test(!dustem_dat,/dir) THEN BEGIN ;spawn,'mkdir '+!dustem_dat force_mkdir,!dustem_dat ENDIF IF not file_test(!dustem_res,/dir) THEN BEGIN ;spawn,'mkdir '+!dustem_dat force_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 ;window,0 & loadct,13 yr=[1e-14,1.e-6] xr=[1,5e3] tit='dustem_init wrapper test '+!dustem_model xtit=textoidl('\lambda (\mum)') ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)') yr=[1e-11,6.e-7] xr=[1,5e3] dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit ENDIF ENDIF ;Initialize filters dustem_filter_init the_end: END