Added Modified-Black-body plugin

Ilyes Choubani
1 parent 53711375
Showing 6 changed files with 130 additions and 27 deletions Show diff stats
src/idl/dustem_activate_plugins.pro
src/idl/dustem_cc.pro
src/idl/dustem_compute_sed.pro
src/idl/dustem_init_plugins.pro
src/idl/dustem_plot_fit_sed.pro
src/idl/dustem_plugin_mbbdy.pro
@@ -102,8 +102,6 @@ FOR i=0L,n_elements(p_min)-1 DO BEGIN
  
             (*!dustem_scope).(k)=ptr_new(scope)
  
-            tgnms=tag_names((*!dustem_scope))
-            
             str='((*!dustem_plugin).('+strtrim(k,2)+'))=ptr_new('+ftn+'(key=index,val=value)'+')' & str=str(0)
             toto=execute(str) & IF !dustem_verbose NE 0 THEN message,strupcase(strmid(ftn,7)),/info
  
@@ -132,7 +132,7 @@ FOR i=0L,Nbands-1 DO BEGIN
       'NUINU=CSTE': BEGIN
         spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
         num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
-        den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)            
+        den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         cc(i)=num/den/spec0*(st.central_wavelengths(ii))
       END
       'MIPS': BEGIN
@@ -59,6 +59,9 @@ spec = st.sed.em_tot * fact
  
 ;ADDING PLUGIN TO SPECTRUM----------------
 ;if n_tags(!dustem_data.sed) gt 1 then begin ; executing only when !dustem_data.sed is set meaning that this allows for the creation of fake seds
+
+
+
 scopes=tag_names((*!dustem_scope))
 IF scopes[0] NE 'NONE' THEN BEGIN
   for i=0L,n_tags(*!dustem_scope)-1 do begin
@@ -23,30 +23,30 @@ for i=0L,n_elements(pd)-1 do begin
 ENDFOR
  
  
-
-IF KEYWORD_SET(fpd) THEN BEGIN
-    
-    for i=0L,n_elements(fpd)-1 do begin
-        fi=strtrim(strmid(fpd(i),0,strlen(plugind_detect_string)),2)
-        if fi eq plugind_detect_string then begin
-            if isa((*!dustem_fit).fixed_param_descs) then begin
-            
-                ftn = strmid((*(*!dustem_fit).fixed_param_descs)(i),0) ; String containing the name of the plugin and the keyword used (ie: dustem_create_continuum_2)    
-            
-            endif else  ftn = strmid((*(*!dustem_fit).param_descs)(i),0) 
-            
-                
-            ii = strsplit(ftn,'_',count=countx)
-            ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
-            ftn = strmid(ftn,14,ii-14) ; String containing the name of the plugin without the associated keyword
-            plugin_names=[plugin_names,ftn]
-            Nplugins=Nplugins+1
-        ENDIF
-    ENDFOR
-
-    
-ENDIF
-
+;Dustem_init_fixed_params works so we do not need this block 
+; 
+; IF KEYWORD_SET(fpd) THEN BEGIN
+;     
+;     for i=0L,n_elements(fpd)-1 do begin
+;         fi=strtrim(strmid(fpd(i),0,strlen(plugind_detect_string)),2)
+;         if fi eq plugind_detect_string then begin
+;             if isa((*!dustem_fit).fixed_param_descs) then begin
+;             
+;                 ftn = strmid((*(*!dustem_fit).fixed_param_descs)(i),0) ; String containing the name of the plugin and the keyword used (ie: dustem_create_continuum_2)    
+;             
+;             endif else  ftn = strmid((*(*!dustem_fit).param_descs)(i),0) 
+;             
+;                 
+;             ii = strsplit(ftn,'_',count=countx)
+;             ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
+;             ftn = strmid(ftn,14,ii-14) ; String containing the name of the plugin without the associated keyword
+;             plugin_names=[plugin_names,ftn]
+;             Nplugins=Nplugins+1
+;         ENDIF
+;     ENDFOR
+; 
+;     
+; ENDIF
  
 IF Nplugins EQ 0 THEN BEGIN
     plugin_names=['NONE']
@@ -193,6 +193,11 @@ IF tag_exist(*!dustem_scope,&#39;SYNCHROTRON&#39;) THEN BEGIN
   cgoplot,st.sed.wav,(*(*!dustem_plugin).synchrotron)[*,0],color='Crimson',linestyle=3
 ENDIF
  
+IF tag_exist(*!dustem_scope,'MBBDY') THEN BEGIN
+  cgoplot,st.sed.wav,(*(*!dustem_plugin).mbbdy)[*,0],color='Magenta',linestyle=3
+ENDIF
+
+
 cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
 ;plot the normalized data as well. 
 ;----------------------------------------
@@ -0,0 +1,97 @@
+FUNCTION dustem_plugin_mbbdy ,key=key,val=val,scope=scope,help=help
+
+;+
+; NAME:
+;    dustem_plugin_mbbdy
+; PURPOSE:
+;    DUSTEM plugin to compute modified black-body emission
+; CATEGORY:
+;    DUSTEM Wrapper
+; CALLING SEQUENCE:
+;    output=dustem_plugin_mbbdy([,key=][,val=])
+; INPUTS:
+;    None
+; OPTIONAL INPUT PARAMETERS:
+;    key  = input parameter number
+;    val  = input parameter value
+; OUTPUTS:
+;    output = MBB spectrum (on dustem wavelengths) (output[*,0] is Stokes I, output[*,1] Stokes Q and output[*,2] Stokes U)
+; OPTIONAL OUTPUT PARAMETERS:
+;    None
+; ACCEPTED KEY-WORDS:
+;    help                  = if set, print this help
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    None
+; PROCEDURE:
+;    This is a dustem plugin
+; EXAMPLES
+;
+; MODIFICATION HISTORY:
+;    Written by JPB 
+;-
+
+IF keyword_set(help) THEN BEGIN
+  doc_library,'dustem_plugin_mbbdy'
+  output=0.
+  goto,the_end
+ENDIF
+
+;default values of input parameters
+amp = 0.5 ; Amplitude of the modified black-body emission  
+temp = 19 ; Temperature of the dust species
+beta = 1.8 ; Emissivity index
+
+w0=100. ; Theoretical wavelength (3THz) where optical depth equals unity
+;check formula of modified black body in the optically thin case
+
+;default polarization values 
+smallp=0.
+psi=0.
+
+
+IF keyword_set(key) THEN BEGIN 
+  
+  ind1=where(key EQ 1,count1) ; reading amp
+  ind2=where(key EQ 2,count2) ; reading temp
+  ind3=where(key EQ 3,count3) ; reading beta
+  
+  ;In the case this plugin is used in polarization 
+  ind4=where(key EQ 4,count4) ; reading smallp
+  ind5=where(key EQ 5,count5) ; reading psi  
+  
+      
+  IF count1 NE 0 then amp=val[ind1[0]] ;else default value
+  IF count2 NE 0 then temp=val[ind2[0]] ; // //
+  IF count3 NE 0 then beta=val[ind3[0]] ; // //
+  IF count4 NE 0 then smallp=val[ind4[0]] ; // //
+  IF count5 NE 0 then psi=val[ind5[0]] ; // //
+   
+ENDIF
+
+lambir=dustem_get_wavelengths()
+Nwavs=n_elements(lambir)
+spec=amp*(lambir/w0)^(-1.*beta)*dustem_planck_function(temp,lambir)
+
+
+;see Deschenes et al 2008, eq 6
+output=fltarr(Nwavs,3)
+output[*,0]=spec
+
+;polarization this actually need to 
+polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
+
+output[*,1]=Q
+output[*,2]=U
+
+
+scope='ADD_SED+ADD_QSED+ADD_USED'
+
+
+the_end:
+RETURN,output
+  
+END