Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Annie Hughes
2 parents ab7bdcb3 6ad05c38
Showing 12 changed files with 119 additions and 106 deletions Show diff stats
LabTools/IRAP/JPB/srun/make_phangs_grids.pro
plugin_scopes_definition.xcat
src/idl/dustem_activate_plugins.pro
src/idl/dustem_compute_stokes.pro
src/idl/dustem_fit_intensity_mbb_example.pro
src/idl/dustem_init_plugins.pro
src/idl/dustem_plugin_freefree.pro
src/idl/dustem_plugin_mbbdy.pro
src/idl/dustem_plugin_modify_dust_pol.pro
src/idl/dustem_plugin_synchrotron.pro
src/idl/dustem_read_dustem_lv.pro
src/idl/dustem_run_plugins.pro
@@ -13,6 +13,9 @@ PRO make_phangs_grids,bidon
 ;more {$phangs_srundir}/make_phangs_grids.log | grep "iteration "
 ;more {$phangs_srundir}/make_phangs_grids.log | grep "Will produce"
+;This won't work, because it needs the Fortran to set files to the same area ...
+;defsysv,'!DUSTEM_RES','/tmp/jbernard/Dustem_nohup'
+
 t0=systime(/sec)
 ;make grid 
@@ -13,13 +13,13 @@
 \------------------------------------------------
 \This scope is for plugins adding to the output SEDs 
 ADD           SED        1
-ADD           POLSED     1
+\ADD           POLSED     1
 \This scope is for plugins replacing the ISRF before fortran is run 
 REPLACE       ISRF       2
 ADD           ISRF       2
 \This scope is for plugins replacing the fortran SEDs 
 REPLACE       SED        3
-REPLACE       POLSED     3
+\REPLACE       POLSED     3
 REPLACE       POLEXT     3
 \This scope is for plugins which fully replace the Fortran, including by calling Fortran several times ...
 REPLACE       FORTRAN    4
@@ -78,36 +78,27 @@ ENDELSE
 ;stop
-;Run the plugins, according to their scopes 
+;===== Run the plugins, according to their scopes (and the Fortran)
-;normal (additive) plugins come first as they can be used by other plugins 
-;The ISRF plugin is next - we can always change its scope however we please to. DustEm runs here. This is with the /force_dustem_run keyword. 
-;The plugins that use the DustEM output come next - we can always change their 'common' scope 
+;The 'ADD_SED' plugins are run first. Note that their output could be used by plugins run later. 
+;The plugins related to 'ISRF' are run second, and the Fortrn is run at that stage, just after the corresponding plugins.
+;The plugins which Replace the output of the fortran are run last.
 ;IF /avoid is used the specified scopes are avoided. They're considered when /avoid is missing.
 ;using '**' so not to write the entire scope
-;=== This runs plugins with scopes including STELLAR, REPLACE and ISRF
-;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['ADD','REPLACE'], avoid=1,use_previous_fortran=use_previous_fortran
-dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['ADD_SED'], avoid=1,use_previous_fortran=use_previous_fortran
+;=== This runs the plugins
+;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['ADD_SED'], avoid=1,use_previous_fortran=use_previous_fortran
+dustem_run_plugins, p_dim ,param_descs, param_values, param_func, 1, avoid=1,use_previous_fortran=use_previous_fortran
-;=== The line below makes fawlty go Segmentation fault. st is undefined at that stage
-;=== it would run fine if avoid was set to 1
-;IC: I think the segmentation fault comes from the fact that the user wants to use models that have the /pol keyword in the GRAIN.DAT file without being in polarization mode.
+;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['ADD_ISRF','REPLACE_ISRF'] ,force_dustem_run=1,st=st,use_previous_fortran=use_previous_fortran ;dustem output is available at this step.
+dustem_run_plugins, p_dim ,param_descs, param_values, param_func, 2,force_dustem_run=1,st=st,use_previous_fortran=use_previous_fortran ;dustem output is available at this step.
-;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['STELLAR','ISRF'] ,force_dustem_run=1,st=st,use_previous_fortran=use_previous_fortran ;dustem output is available at this step.
-dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['ADD_ISRF','REPLACE_ISRF'] ,force_dustem_run=1,st=st,use_previous_fortran=use_previous_fortran ;dustem output is available at this step.
+;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['REPLACE_SED'],use_previous_fortran=use_previous_fortran
+dustem_run_plugins, p_dim ,param_descs, param_values, param_func, 3,use_previous_fortran=use_previous_fortran
-;added this small condition to ensure that dustem is run
-;(if no ISRF plugin is used)  
-; if ~isa(st) then begin
-;     st=dustem_run(p_dim)
-;     !dustem_current=ptr_new(st)
-; endif
-;THE ABOVE CONDITION IS UNNECESSARY SINCE DUSTEM IS CALLED ABOVE WHETHER PLUGINS ARE RUN OR NOT.
-
-;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['REPLACE'],use_previous_fortran=use_previous_fortran
-dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['REPLACE_SED'],use_previous_fortran=use_previous_fortran
+;dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['REPLACE_SED'],use_previous_fortran=use_previous_fortran
+dustem_run_plugins, p_dim ,param_descs, param_values, param_func, 4,use_previous_fortran=use_previous_fortran
 the_end:
@@ -75,18 +75,25 @@ P = st.polsed.em_tot * fact   ;This is Polarized intensity P
 Nwaves=(size(stI))[1]
-frac=P/stI
+;stop
+
+;This some times leads to division by 0
+frac=P/stI       
 tes=where(finite(frac) eq 0)
 frac(tes)=0.
-scopes=tag_names((*!dustem_plugin))
+scopes=(*!dustem_plugin).scope
+n_plugins=n_elements(scopes)
+
 IF scopes[0] NE 'NONE' THEN BEGIN
-   FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
-      IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN
-          Q_spec = (*(*!dustem_plugin).(i).spec)[*,1]
-          U_spec = (*(*!dustem_plugin).(i).spec)[*,2]
-      ENDIF 
-  ENDFOR    
+   FOR i=0L,n_plugins-1 DO BEGIN
+      ;IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN
+      ;IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN
+      IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) EQ 'REPLACE_SED') THEN BEGIN
+          Q_spec = (*(*!dustem_plugin)[i].spec)[*,1]
+          U_spec = (*(*!dustem_plugin)[i].spec)[*,2]
+      ENDIF
+  ENDFOR
 ENDIF
   IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN
@@ -94,18 +101,14 @@ ENDIF
   ENDIF    
 IF scopes[0] NE 'NONE' THEN BEGIN          
-  FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
-      
-      IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'ADD_POLSED') THEN BEGIN
-          
-          Q_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
-          U_spec+=(*(*!dustem_plugin).(i).spec)[*,2]
-          
-      
+  FOR i=0L,n_plugins-1 DO BEGIN
+      ;IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'ADD_POLSED') THEN BEGIN
+      ;IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) EQ 'ADD_POLSED') THEN BEGIN
+      IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) EQ 'ADD_SED') THEN BEGIN
+          Q_spec+=(*(*!dustem_plugin)[i].spec)[*,1]
+          U_spec+=(*(*!dustem_plugin)[i].spec)[*,2]
       ENDIF
-
   ENDFOR
-
 ENDIF
 ;-----------------------------------------
 ind=where(finite(U_spec) NE 1 OR finite(Q_spec) NE 1,count)
@@ -162,7 +162,7 @@ if keyword_set(wait) then begin
 end
 ;== INITIALISE DUSTEM
-dustem_init,model=use_model,polarization=use_polarization
+dustem_init,model=use_model,polarization=use_polarization,/show_plots
 !dustem_nocatch=1
 !dustem_verbose=use_verbose
 IF keyword_set(noobj) THEN !dustem_noobj=1
@@ -88,7 +88,13 @@ plugin_names=plugin_names[un]
 parameter_types=parameter_types[un]
 ind=where(parameter_types EQ 'PLUGIN',Nplugins)
-one_plugin_st={name:'',spec:ptr_new(),scope:'',paramtag:ptr_new()}
+;===== This is the dustem structure describing the plugin scopes
+one_plugin_st={name:'', $           ;name of the plugin
+               spec:ptr_new(), $    ;output of the plugin (a spectrum)
+               scope:'', $          ;scope of the plugin
+               run_order:0L,  $      ;plugin run order
+               paramtag:ptr_new() $ ;parameter names
+               }
 IF Nplugins EQ 0 THEN BEGIN
     plugin_st=replicate(one_plugin_st,1)
@@ -114,6 +120,18 @@ FOR i=0L,Nplugins-1 DO BEGIN
     toto=execute(str)
     plugin_st[i].paramtag=ptr_new(paramtag)
 ENDFOR
+;derive the run order from the scopes
+file_scopes=!dustem_wrap_soft_dir+'plugin_scopes_definition.xcat'
+scope_st=read_xcat(file_scopes,/silent)
+possible_scopes=scope_st.objective+'_'+scope_st.target
+FOR i=0L,Nplugins-1 DO BEGIN
+    ind=where(possible_scopes EQ plugin_st[i].scope,count)
+    IF count EQ 0 THEN BEGIN
+        message,'scope '+plugin_st[i].scope+' not recognized',/continue
+        stop
+    ENDIF
+    plugin_st[i].run_order=scope_st[ind[0]].run_order
+ENDFOR
 ;creates the !dustem_plugin structure
 defsysv, '!dustem_plugin', ptr_new(plugin_st)
@@ -23,8 +23,11 @@ FUNCTION dustem_plugin_freefree,key=key $
 ;
 ; OPTIONAL INPUT PARAMETERS:
 ;    key  = input parameter number
+;         1 : Tgas [K]
+;         2 : free-free Amplitude []
+;         3 : polarization fraction [%] (default=0)
+;         4 : polarization angle [deg]
 ;    val  = input parameter value
-;
 ; OUTPUTS:
 ;    freefree = free-free spectrum (on dustem wavelengths)
 ;
@@ -77,29 +80,21 @@ IF keyword_set(key) THEN BEGIN
   b=where(key EQ 2,count2)
   c=where(key EQ 3,count3) ;default polarization fraction -newly added
   d=where(key EQ 4,count4) ;default polarization angle -newly added
-  IF count1 NE 0 then Tgas=(val(a))(0)
-  IF count2 NE 0 then Amplitude=(val(b))(0)
-  IF count3 NE 0 then smallp=val[c[0]]   ;-newly added
-  IF count4 NE 0 then psi=val[d[0]]    ;-newly added
+  IF count1 NE 0 then Tgas=val[a[0]]
+  IF count2 NE 0 then Amplitude=val[b[0]]
+  IF count3 NE 0 then smallp=val[c[0]]
+  IF count4 NE 0 then psi=val[d[0]]
 ENDIF
-;IF !dustem_which EQ 'DESERT' THEN BEGIN
-;  lambir=((*!dustem_params).gemissiv.lambir)
-;ENDIF ELSE BEGIN
-;  lambir=((*!dustem_params).lambda.lambda)
-;ENDELSE
-;replaced by the following to ease the life of plugins writters
 lambir=dustem_get_wavelengths()
 Nwavs=n_elements(lambir)
 cmic=3.e14
 nu=cmic/lambir   ;Hz
 mjy=1           ;output is in MJy/sr
-lambir_ref=10000.
-
+lambir_ref=10000.  ;in mic (this is 1cm)
-output=fltarr(Nwavs,3) ;newly added
+output=fltarr(Nwavs,3)
-;stop
 ;use_method='Deschenes2008'
 ;use_method='WallsGabaud1998'
 use_method='Dickinson2003_norm'
@@ -140,10 +135,7 @@ ENDCASE
 polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
 output[*,1]=Q
-output[*,2]=U
-;output[*,3]=replicate(smallp,Nwavs)
-;output[*,4]=replicate(psi,Nwavs)
- 
+output[*,2]=U 
 the_end:
 RETURN,output
@@ -52,7 +52,8 @@ w0=100. ; Theoretical wavelength (3THz) where optical depth equals unity
 smallp=0.
 psi=0.
-scope='ADD_SED+ADD_POLSED'
+;scope='ADD_SED+ADD_POLSED'
+scope='ADD_SED'
 paramtag=['Amp',textoidl('T_{MBB}')+' [K]',textoidl('\beta'),'p','Psi [deg]']
 IF keyword_set(key) THEN BEGIN 
@@ -5,31 +5,33 @@ FUNCTION dustem_plugin_modify_dust_pol, key=key, val=val, scope=scope, paramtag=
 ;    dustem_plugin_modify_dust_pol
 ;  
 ; PURPOSE:
-;    modifies Emission Stokes Q,U values for dust polarization in st according to the given keywords and values
-;
+;    Replaces the emission in polarization (Stokes Q,U) with a fixed polarization angle and fraction
+;    Q,U to be replaced are computed from the total intensity as:
+;    Q=I*val[0]*p*cos(2*val[1])
+;    U=I*val[0]*p*sin(2*val[1])
+;    where I is the total intensity in the model and p is the polarization fraction in the model 
 ; CATEGORY:
 ;    DustEMWrap, Plugin, Mid-level, Distributed
 ;
 ; CALLING SEQUENCE:
-;    a=dustem_plugin_modify_dust_pol(st,[key=][val=][,scope=][,paramtag=][,paramdefault=][,/help])
+;    a=dustem_plugin_modify_dust_pol([key=][val=][,scope=][,paramtag=][,paramdefault=][,/help])
 ;
 ; INPUTS:
 ;    st = dustem structure
 ;
 ; OPTIONAL INPUT PARAMETERS:
-;    key  = input parameter numbers (first = polarization fraction in %, default=1.%, second=polarization angle, default=0.)
+;    key  = input parameter numbers
+;           first  = multiplication factor for the polarization fraction
+;           second = polarization angle [deg], default=0.
 ;    val  = input parameter values
-;
 ; OUTPUTS:
-;    out = array containing the stokes emission parameters associated to the dust/synchrotron component or a concatenated version for both
+;    out  = array containing the stokes emission parameters associated to the dust/synchrotron component or a concatenated version for both
 ;
 ; OPTIONAL OUTPUT PARAMETERS:
 ;    scope = if set, returns only the scope of the pluggin
 ;    paramtag = if set, returns only the parameter tags
-;
 ; ACCEPTED KEY-WORDS:
 ;    help                  = if set, print this help
-;
 ; COMMON BLOCKS:
 ;    None
 ;
@@ -53,20 +55,25 @@ IF keyword_set(help) THEN BEGIN
 ENDIF
 IF keyword_set(scope) THEN BEGIN
+    ;scope='REPLACE_POLSED'
+    scope='REPLACE_SED'
     out=0
     GOTO, the_end
 ENDIF 
 IF keyword_set(paramtag) THEN BEGIN
+    paramtag=['p',textoidl('\psi')+' [deg]']
     out=0.
     GOTO, the_end
 ENDIF 
 ;Retrieving this sytem variable  (dustem output)  
-if ~isa(!dustem_current) then begin
-    out=0 
+IF ~isa(!dustem_current) then begin
+    out=0
     goto, the_end 
-ENDIF else st=(*!dustem_current)
+ENDIF ELSE BEGIN
+    st=(*!dustem_current)
+ENDELSE
 ;below are the default values for the plugin parameters
 smallp_fact=1.     ;This is the default dust polarization factor
@@ -76,12 +83,13 @@ IF keyword_set(key) THEN BEGIN
   ind1=where(key EQ 1,count1)
   ind2=where(key EQ 2,count2)
   IF count1 NE 0 THEN smallp_fact=val[ind1[0]] ; setting smallp from pd - this is another polarization fraction (constant) that is applied to the total dust emission
-  IF count2 NE 0 THEN psi=val[ind2[0]] & !dustem_psi = psi ; setting psi from pd.  !dustem_psi here helps for the plotting.   
+  IF count2 NE 0 THEN psi=val[ind2[0]]
+  !dustem_psi = psi ; setting psi from pd.  !dustem_psi here helps for the plotting.   
 ENDIF
 fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/((st.polsed).wav))*1.e20/1.e7
-P=((st.polsed).em_tot)*fact   ; This is the polarized emission P
-I=((st.sed).em_tot)*fact      ; This is the total intensity emission I
+P=((st.polsed).em_tot)*fact   ; This is the polarized emission P currently in the model
+I=((st.sed).em_tot)*fact      ; This is the total intensity emission I currently in the model
 Nwaves=(size(I))[1]
@@ -90,7 +98,7 @@ Nwaves=(size(I))[1]
 frac_model=P*0.
 indI = where(I ne 0, countI)
-if countI ne 0 then frac_model[indI] = P[indI]/I[indI] ;This is the polarization fraction in the model
+if countI ne 0 then frac_model[indI] = P[indI]/I[indI] ;This is the polarization fraction currently in the model
 frac_used=frac_model*smallp_fact
@@ -98,22 +106,15 @@ psi_used = replicate(psi,Nwaves)
 polar_ippsi2iqu,I,Q,U,frac_used,psi_used
-out=fltarr(Nwaves,3) ; modified this. This is the only plugin that has this number of outputs.
-
-;degtorad = !pi/180
+out=fltarr(Nwaves,3)
 out[*,0]=I
 out[*,1]=Q
 out[*,2]=U
-; out[*,3]=frac_used
-; ;out[*,4]=0.5*atan(U,Q)/degtorad
-; out[*,4]=psi_used
 RETURN, out
 the_end:
-scope='REPLACE_POLSED'
-paramtag=['p',textoidl('\psi')+' [deg]']
 END
@@ -64,7 +64,8 @@ s=3           ;power spectrum of CR E distribution
 A=1           ;Synchrotron radiation amplitude at 10 mm
 psi=0.        ;default polarization angle
 smallp=0.3;0.3    ;default polarization fraction
-scope='ADD_SED+ADD_POLSED'
+;scope='ADD_SED+ADD_POLSED'
+scope='ADD_SED'
 ;paramtag=['s (plaw_index)','Amp','p','Psi (deg)']
 paramtag=[textoidl('\alpha_{CR}'),'Amp','p',textoidl('\psi')+' [deg]']
 paramdefault=[s,A,smallp,psi]
@@ -79,6 +79,7 @@ ENDWHILE
 ;readf,unit,str
 CLOSE,unit
 free_lun,unit
+
 str=strcompress(str)
 str=strtrim(str,2)
 vv=str_sep(str,' ')
@@ -100,12 +101,11 @@ Nlines=n_elements(wav)
 instruc='one_st={wav:0.,'
 FOR i=0L,Ngrains-1 DO BEGIN
-  instruc=instruc+'em_grain_'+strtrim(i+1,2)+':0.,'
+  instruc=instruc+'em_grain_'+strtrim(i+1,2)+':0.D0,'
 ENDFOR
-instruc=instruc+'em_tot:0.}'
+instruc=instruc+'em_tot:0.D0}'
 message,'Executing '+instruc,/info
 toto=execute(instruc)
-
 st=replicate(one_st,Nlines)
 st.wav=wav
@@ -113,6 +113,7 @@ FOR i=0L,Ngrains-1 DO BEGIN
   instruc='st.(i+1)='+'em_grain_'+strtrim(i+1,2)
   message,'Executing '+instruc,/info
   toto=execute(instruc)
+  ;stop
 ENDFOR
 st.em_tot=em_tot
-PRO dustem_run_plugins, p_dim ,$
-                        param_descs,$
-                        param_values,$
-                        param_func,$
-                        scopes,$
-                        st=st,$
-                        avoid=avoid,$
-                        force_dustem_run=force_dustem_run,$
-                        use_previous_fortran=use_previous_fortran, $
-                        help=help
+PRO dustem_run_plugins, p_dim $
+                        ,param_descs $
+                        ,param_values $
+                        ,param_func $
+                        ,run_order $
+                        ,st=st $
+                        ,avoid=avoid $
+                        ,force_dustem_run=force_dustem_run $
+                        ,use_previous_fortran=use_previous_fortran  $
+                        ,help=help
 ;+
 ; NAME:
@@ -28,7 +28,7 @@ PRO dustem_run_plugins, p_dim ,$
 ;    param_descs = parameter description vector
 ;    param_values = current parameter values
 ;    param-func   = plugin indices array 
-;    scopes       = scopes of plugins 
+;    run_order    = run order of the plugins to be run
 ;
 ; OPTIONAL INPUT PARAMETERS:
 ;    avoid      = scopes should be avoided
@@ -85,12 +85,14 @@ Nplugins=n_elements(*!dustem_plugin)
 must_be_run=intarr(Nplugins)
 FOR i=0L,Nplugins-1 DO BEGIN
-	FOR j=0L,n_elements(scopes)-1 DO BEGIN
-		pos=strposmulti((*!dustem_plugin)[i].scope,scopes[j],count)
-		IF count NE 0 THEN must_be_run[i]=1
-  	ENDFOR
+		IF (*!dustem_plugin)[i].run_order EQ run_order THEN BEGIN
+			must_be_run[i]=1
+		ENDIF
 ENDFOR
+;print,must_be_run
+;stop
+
 ;=== run plugins that need be run
 FOR i=0L,Nplugins-1 DO BEGIN
 	IF must_be_run[i] THEN BEGIN