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

Jean-Philippe Bernard
2 parents 21d1a616 6cde48ea
Showing 1 changed file with 92 additions and 0 deletions Show diff stats
src/idl/dustem_plugin_modify_dust_pol.pro
@@ -0,0 +1,92 @@
+Function dustem_plugin_modify_dust_pol, key=key, val=val, scope=scope,paramtag=paramtag,help=help
+
+;+
+; NAME:
+;    dustem_plugin_modify_dust_pol
+; PURPOSE:
+;    Produces stokes emission parameters for the dust and synchrotron components
+; CATEGORY:
+;    DUSTEM Wrapper
+; CALLING SEQUENCE:
+;    dustem_create_stokes(st,key=key,val=val)
+; INPUTS:
+;
+; OPTIONAL INPUT PARAMETERS:
+;    key  = input parameter number
+;    val  = input parameter value
+; OUTPUTS:
+;    out = array containing the stokes emission parameters associated to the dust/synchrotron component or a concatenated version for both
+; OPTIONAL OUTPUT PARAMETERS:
+;    None
+; ACCEPTED KEY-WORDS:
+;    help                  = if set, print this help
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    The dustem fortran code must be installed
+;    The dustem idl wrapper must be installed
+; PROCEDURE:
+;    This is a dustem plugin
+;-
+
+IF keyword_set(help) THEN BEGIN
+  doc_library,'dustem_create_stokes'
+  out=0.
+  goto,the_end
+ENDIF
+
+
+paramtag=['p','Psi (deg)']
+
+IF  ptr_valid((*!dustem_fit).CURRENT_PARAM_VALUES) THEN BEGIN
+    p_di=(*(*!dustem_fit).CURRENT_PARAM_VALUES)
+    
+ENDIF ELSE BEGIN
+          p_di=(*(*!dustem_fit).PARAM_INIT_VALUES)    
+      ENDELSE 
+          
+st=dustem_run(p_di)      
+;stop
+
+psi=0.
+smallp=1.
+IF keyword_set(key) THEN BEGIN 
+  ind1=where(key EQ 1,count1)
+  ind2=where(key EQ 2,count2)
+  IF count1 NE 0 THEN smallp=val[ind1[0]] ; setting smallp from pd
+  IF count2 NE 0 THEN psi=val[ind2[0]] ; setting psi from pd   
+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
+I=((st.sed).em_tot)*fact
+
+
+
+Nwaves=(size(I))[1]
+
+frac=P/I*smallp
+tes=where(finite(frac) eq 0)
+frac(tes)=0.
+
+
+polar_ippsi2iqu,I,Q,U,frac,replicate(psi,Nwaves) 
+
+out=fltarr(Nwaves,3)
+
+out[*,0]=I
+out[*,1]=Q
+out[*,2]=U
+
+scope='REPLACE_QSED+REPLACE_USED'
+
+return, out
+
+the_end:
+
+
+end
+