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

Ilyes Choubani
2 parents 25140f99 474c4b97
Showing 13 changed files with 52 additions and 60 deletions Show diff stats
instrument_description.xcat
src/idl/dustem_compute_polsed.pro
src/idl/dustem_compute_stokes.pro
src/idl/dustem_filter2instru.pro
src/idl/dustem_fit_intensity_example.pro
src/idl/dustem_fit_sed_ext_pol_example.pro
src/idl/dustem_get_band_flux.pro
src/idl/dustem_plugin_mbbdy.pro
src/idl/dustem_plugin_modify_dust_pol.pro
src/idl/dustem_plugin_synchrotron.pro
src/idl/dustem_write_fits_table.pro
src/idl/dustemwrap_plot.pro
src/idl/dustemwrap_plot_noobj.pro
@@ -395,4 +395,5 @@ GISMO      GISMO1     2000.00    MJy/sr      nuInu=cste   0.0046389
 \Everything should be checked.
 BOLOCAM    BOLOCAM1   1105.4     MJy/sr      nuInu=cste   0
 \================ ================================================================================
-
+\This last dummy filter needs to be included to cope for spectrum data
+UNKNOWN    SPECTRUM   -32768.    MJy/sr      NONE         0
 FUNCTION dustem_compute_polsed,p_dim,$
-st=st,$
-P_spec=P_spec,$
-SP_spec=SP_spec,$
-dustem_polfrac=dustem_polfrac,$
-out_st=out_st,$
-dustem_sed=dustem_sed,$
-_extra=extra
-
+                               st=st,$
+                               P_spec=P_spec,$
+                               SP_spec=SP_spec,$
+                               dustem_polfrac=dustem_polfrac,$
+                               dustem_sed=dustem_sed
 ;+
 ; NAME:
@@ -16,7 +13,7 @@ _extra=extra
 ; CATEGORY:
 ;    Dustem
 ; CALLING SEQUENCE:
-;    sed=dustem_compute_polsed(p_dim[,st=][,P_spec=],[SP_spec=][dustem_polfrac=][out_st=][dustem_sed=][,_extra=][,/help])
+;    sed=dustem_compute_polsed(p_dim[,st=][,P_spec=],[SP_spec=][dustem_polfrac=][dustem_sed=][,/help])
 ; INPUTS:
 ;    p_dim      = parameter values
 ; OPTIONAL INPUT PARAMETERS:
@@ -29,7 +26,8 @@ _extra=extra
 ;    SP_spec       = Dustem polarized emission fraction
 ;    dustem_polfrac = Computed polarization emission fraction for spectrum points in !dustem_data
 ; OPTIONAL OUTPUT PARAMETERS:
-;    out_st        = Dustem output structure
+;    st        = Dustem structure
+;    dustem_sed = Computed emission for spectrum points in !dustem_data
 ; ACCEPTED KEY-WORDS:
 ;    help      = If set, print this help
 ; COMMON BLOCKS:
@@ -162,7 +160,7 @@ out_st=st
 ;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac)
 if !run_lin then begin
-    if not keyword_set(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
+    if not keyword_set(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
     If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin                            
         if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin 
@@ -4,8 +4,6 @@ FUNCTION dustem_compute_stokes,p_dim,$
                                U_spec=U_spec,$
                                PSI_spec=PSI_spec,$
                                dustem_psi_em=dustem_psi_em,$
-                               out_st=out_st,$
-                               _extra=extra,$
                                help=help
 ;+
@@ -16,7 +14,7 @@ FUNCTION dustem_compute_stokes,p_dim,$
 ; CATEGORY:
 ;    Dustem
 ; CALLING SEQUENCE:
-;    out=dustem_compute_stokes(p_dim[,st][,Q_spec][,U_spec][,PSI_spec][,dustem_psi_em][,out_st=][,_extra=]
+;    out=dustem_compute_stokes(p_dim[,st][,Q_spec][,U_spec][,PSI_spec][,dustem_psi_em]
 ; INPUTS:
 ;    p_dim      = parameter values
 ; OPTIONAL INPUT PARAMETERS:
@@ -26,7 +24,7 @@ FUNCTION dustem_compute_stokes,p_dim,$
 ;    out[0]     = Computed Stokes Q emission for spectrum points in !dustem_data
 ;    out[1]     = Computed Stokes U emission for spectrum points in !dustem_data 
 ; OPTIONAL OUTPUT PARAMETERS:
-;    out_st        = Dustem output structure
+;    st          = Dustem output structure
 ;    Q_spec        = Dustem Stokes Q emission 
 ;    U_spec        = Dustem Stokes U emission
 ;    PSI_spec      = Dustem Polarization angle output in emission
@@ -55,8 +53,7 @@ IF keyword_set(help) THEN BEGIN
 ENDIF
 IF not keyword_set(st) THEN BEGIN
-    dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st
-    ;st=dustem_run(p_dim)
+    dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st
 ENDIF
@@ -174,8 +171,6 @@ IF count_uspec NE 0 THEN BEGIN
   if isa(dustem_used) then dustem_used[ind_uspec]=interpol(U_spec,st.polsed.wav,(((*(*!dustem_data).used).wav)[ind_uspec]))
 ENDIF
-out_st=st
-
 ;generating interpolates for Psi_em
 dustem_psi_em = 0.5*atan(dustem_used,dustem_qsed)/!dtor
@@ -44,11 +44,11 @@ Nfilt=n_elements(filter)
 ind1=lonarr(Nfilt)
 FOR i=0L,Nfilt-1 DO BEGIN
-  ind=where(!dustem_instrument_description.filter EQ strtrim(filter(i),2),count)
+  ind=where(!dustem_instrument_description.filter EQ strtrim(filter[i],2),count)
   IF count EQ 0 THEN BEGIN
     message,'Filter: '+strtrim(filter(i),2)+' not found in instrument_description.xcat'
   ENDIF
-  ind1(i)=ind(0)
+  ind1[i]=ind[0]
 ENDFOR
 instru=!dustem_instrument_description(ind1).instru
@@ -387,7 +387,7 @@ IF keyword_set(fits_save_and_restore) THEN BEGIN
   ;==== plot result taken from the saved fits table
   res=*(*!dustem_fit).CURRENT_PARAM_VALUES
   IF !dustem_noobj THEN BEGIN
-     dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+    dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
   ENDIF ELSE BEGIN
     dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
   ENDELSE
@@ -418,7 +418,7 @@ ext.EXT_Q = toto[0]
 ext.EXT_U = toto[1]
 sed.StokesI = dustem_compute_sed(rv,st=st)
-toto = dustem_compute_stokes(rv,st=st,Qspec=dustem_qsed,Uspec=dustem_used)
+toto = dustem_compute_stokes(rv,st=st,Q_spec=dustem_qsed,U_spec=dustem_used)
 sed.StokesQ = toto[0]
 sed.StokesU = toto[1]
@@ -238,7 +238,8 @@ FUNCTION DUSTEM_GET_BAND_FLUX, dpath, inst, xs=xs, ys=ys, smi=smi, unit=unit, VE
            for itp=0,ntype-1 do begin
               inu = 1d17 * yw(ixw,itp)/nu/1d9 ; MJy/sr 
               y0 = INTERPOL(yw(ixw,itp),nu,wbd(k))
-              cc = HFI_COLOR_CORRECTION( bp,bnm(k),bolo_cc, nu,inu )
+;              cc = HFI_COLOR_CORRECTION( bp,bnm(k),bolo_cc, nu,inu )
+              cc = HFI_COLOR_CORRECTION_DUSTEM( bp.name,bp.freq/1e9,bp.trans,bnm(k),bolo_cc, nu,inu )
               smdat.(itx).ym(k,itp)  = cc.cc * y0                  ; SED (erg/cm2/s/sr) @ band center
               smdat.(itx).flx(k,itp) = 1d17*cc.cc*y0/wbd(k)/1d9    ; in MJy/sr
               smdat.(itx).cc(k,itp)  = cc.cc
@@ -42,7 +42,7 @@ ENDIF
 ;default values of input parameters
 amp = 0.5 ; Amplitude of the modified black-body emission  
-temp = 19 ; Temperature of the dust species
+temp = 20.0 ; Temperature of the dust
 beta = 1.8 ; Emissivity index
 w0=100. ; Theoretical wavelength (3THz) where optical depth equals unity
@@ -53,7 +53,7 @@ smallp=0.
 psi=0.
 scope='ADD_SED+ADD_POLSED'
-paramtag=['Amp','T_MBB (K)','beta (Plaw_index)','p','Psi (deg)']
+paramtag=['Amp',textoidl('T_{MBB}')+' [K]',textoidl('\beta'),'p','Psi [deg]']
 IF keyword_set(key) THEN BEGIN 
@@ -8,7 +8,7 @@ FUNCTION dustem_plugin_modify_dust_pol, key=key, val=val, scope=scope, paramtag=
 ; CATEGORY:
 ;    DUSTEM Wrapper
 ; CALLING SEQUENCE:
-;    a=dustem_plugin_modify_dust_pol(st,[key=][val=][,scope=][,paramtag=][,/help])
+;    a=dustem_plugin_modify_dust_pol(st,[key=][val=][,scope=][,paramtag=][,paramdefault=][,/help])
 ; INPUTS:
 ;    st = dustem structure
 ; OPTIONAL INPUT PARAMETERS:
@@ -38,8 +38,6 @@ IF keyword_set(help) THEN BEGIN
   GOTO,the_end
 ENDIF
-
-
 IF keyword_set(scope) THEN BEGIN
     out=0
     GOTO, the_end
@@ -57,7 +55,7 @@ if ~isa(!dustem_current) then begin
 ENDIF else st=(*!dustem_current)
 ;below are the default values for the plugin parameters
-smallp_fact=1.     ;This is the default multiplicative factore to the dust polarization
+smallp_fact=1.     ;This is the default dust polarization factor
 psi=0.             ;This is the default polarization angle
 IF keyword_set(key) THEN BEGIN 
@@ -102,7 +100,7 @@ RETURN, out
 the_end:
 scope='REPLACE_POLSED'
-paramtag=['p','Psi (deg)']
+paramtag=['p',textoidl('\psi')+' [deg]']
 END
@@ -52,8 +52,8 @@ psi=0.        ;default polarization angle
 smallp=0.3;0.3    ;default polarization fraction
 scope='ADD_SED+ADD_POLSED'
 ;paramtag=['s (plaw_index)','Amp','p','Psi (deg)']
-paramtag=[textoidl('\alpha_{CR}'),'Amp','p','Psi [deg]']
-paramdefault=[s,A,psi,smallp]
+paramtag=[textoidl('\alpha_{CR}'),'Amp','p',textoidl('\psi')+' [deg]']
+paramdefault=[s,A,smallp,psi]
 IF keyword_set(key) THEN BEGIN 
   ind1=where(key EQ 1,count1)
@@ -55,19 +55,16 @@ IF keyword_set(help) THEN BEGIN
 ENDIF
 ;===== run dustem to get predicted spectrum and sed
-IF ptr_valid((*!dustem_data).sed) THEN BEGIN
-	stop
+IF ptr_valid((*!dustem_data).sed) THEN $
 	dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)
-ENDIF ELSE BEGIN
-ENDELSE
 IF !run_pol THEN BEGIN
 ;This is needed to compute Q,U of both SED and spectra
-   dustem_predicted_polsed=dustem_compute_stokes(*(*!dustem_fit).current_param_values,bidon, $
-                                                 dustem_predicted_Qsed,dustem_predicted_Used, $
-                                                 dustem_predicted_Q_spec,dustem_predicted_U_spec,dustem_predicted_PSI_spec, $
-                                                 dustem_psi_em,out_st=out_st)
-
+   dustem_predicted_polsed=dustem_compute_stokes(*(*!dustem_fit).current_param_values, $
+                                                 Q_spec=dustem_predicted_Q_spec,U_spec=dustem_predicted_U_spec,PSI_spec=dustem_predicted_PSI_spec, $
+                                                 dustem_psi_em=dustem_psi_em)
+   dustem_predicted_Qsed=dustem_predicted_polsed[0]
+   dustem_predicted_Used=dustem_predicted_polsed[1]
 ;Same in extinction, but does not work in polarization extinction was not fitted.
 ;   dustem_predicted_polext=dustem_compute_stokext(*(*!dustem_fit).current_param_values,bidon, $
 ;											  dustem_predicted_Qext,dustem_predicted_Uext, $
@@ -188,7 +188,7 @@ if isa((*!dustem_data).sed) then begin
 ;Leaving as is.
 ;JPB: need here to store values of current parameter values (including plugins) into !iteration_mouchard
-    help,(*!dustem_fit).CURRENT_PARAM_VALUES   ;This is not currently filled up ...
+;    help,(*!dustem_fit).CURRENT_PARAM_VALUES   ;This is not currently filled up ...
 ;Below will work only for non plugin parameters
     goto,skip1
@@ -92,32 +92,33 @@ END
 ;This is also necessary for the plotting of the results of the fit (Last iteration) 
 IF not keyword_set(st) THEN BEGIN
   ;Activation of the plugins is needed because of the plotting of the total emission model that includes them. 
-  dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st  ; 0/0 division case?
+  dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st  ; 0/0 division case?
   ;# Emission
   ;if isa(!dustem_data.sed) then begin ;using !dustem_data instead of !dustem_show because in this release the user is compelled to fit IQU whether in M or X
   if not keyword_set(dustem_sed) and isa((*!dustem_data).sed) then begin ;better test?
-      dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
+      dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
       if !run_pol && !run_lin then begin
-          
-          dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec,dustem_polfrac)
-          toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em) 
+         dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac)
+; the following call need to change to include out instead of dustem_qsed/dustem_used?
+         toto = dustem_compute_stokes(p_dim,st=st,dustem_qsed,dustem_used,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em) 
-          
       endif         
   endif
   ;# Extinction
 ;   if isa(!dustem_data.ext) then begin
   if not keyword_set(dustem_ext) and isa((*!dustem_data).ext) then begin
-      dustem_ext = dustem_compute_ext(p_dim,st,EXT_spec)
+      dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec)
       if !run_pol && !run_lin then begin
-          dustem_polext = dustem_compute_polext(p_dim,st,POLEXT_spec,SPEXT_spec,dustem_fpolext)
-          toto = dustem_compute_stokext(p_dim,st,dustem_qext,dustem_uext,QEXT_spec,UEXT_spec,PSIEXT_spec,dustem_psi_ext) 
+; the following call need to change to include out instead of dustem_fpolext?
+          dustem_polext = dustem_compute_polext(p_dim,st=st,polext_spec=POLEXT_spec,spext_spec=SPEXT_spec,dustem_fpolext)
+; the following call need to change to include out instead of QEXT/UEXT-spec?
+          toto = dustem_compute_stokext(p_dim,st=st,dustem_qext,dustem_uext,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spe=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) 
       endif
   endif
@@ -128,12 +129,13 @@ ENDIF ELSE BEGIN ;st is provided
     ;# Emission
     ;if isa(!dustem_data.sed) then begin
     if not keyword_set(dustem_sed) and isa((*!dustem_data).sed) then begin
-        dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
+        dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
         if !run_pol && !run_lin then begin
-            dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec,dustem_polfrac)
-            toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em) 
+            dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac)
+; the following call need to change to include out instead of dustem_q/used?
+            toto = dustem_compute_stokes(p_dim,st=st,dustem_qsed,dustem_used,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em) 
             ;stop    
         endif
     endif
@@ -143,12 +145,12 @@ ENDIF ELSE BEGIN ;st is provided
     if not keyword_set(dustem_ext) and isa((*!dustem_data).ext) then begin
-        dustem_ext = dustem_compute_ext(p_dim,st,EXT_spec)
+        dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec)
         if !run_pol && !run_lin then begin
-            dustem_polext = dustem_compute_polext(p_dim,st,POLEXT_spec,SPEXT_spec,dustem_fpolext)
-            toto = dustem_compute_stokext(p_dim,st,dustem_qext,dustem_uext,QEXT_spec,UEXT_spec,PSIEXT_spec,dustem_psi_ext) 
+            dustem_polext = dustem_compute_polext(p_dim,st=st,POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext)
+            toto = dustem_compute_stokext(p_dim,st=st,dustem_qext,dustem_uext,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spec=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) 
         endif
     endif
@@ -176,7 +178,7 @@ new_mouchard.rchi2=(*!dustem_fit).rchi2
 ;Leaving as is.
 ;JPB: need here to store values of current parameter values (including plugins) into !iteration_mouchard
-help,(*!dustem_fit).CURRENT_PARAM_VALUES   ;This is not currently filled up ...
+;help,(*!dustem_fit).CURRENT_PARAM_VALUES   ;This is not currently filled up ...
 ;Below will work only for non plugin parameters
 goto,skip1