added all grain types and play nice with system variables

Annie Hughes
1 parent abb4c515
Showing 1 changed file with 58 additions and 13 deletions Show diff stats
src/idl/dustem_plugin_dl07_isrf_model.pro
@@ -61,21 +61,21 @@ IF keyword_set(help) THEN BEGIN
   goto,the_end
 ENDIF
  
+spec=0.
+
 IF keyword_set(scope) THEN BEGIN 
-    out=0
+    out=0 & st=0
     goto, the_end
 ENDIF 
  
 IF keyword_set(paramtag) THEN BEGIN
-    out=0
+    out=0 & st=0
     goto, the_end 
 ENDIF
  
 ;dustem_print_params
 ;stop
  
-spec=0.
-
 ;==== default parameter values
 Umin=1.
 Umax=1.e6
@@ -106,19 +106,47 @@ Umean=(1.-gamma)*Umin+gamma*Umin*alog(Umax/Umin)/(1.-Umin/Umax)
 lambir=dustem_get_wavelengths()
  
 Nwavs=n_elements(lambir)
-spec=fltarr(Nwavs,3)
-spec2=fltarr(Nwavs,3)
+;spec=fltarr(Nwavs,3) ; if ever one day we enforce IQU outputs
+;spec2=fltarr(Nwavs,3)
+spec=fltarr(Nwavs) 
+spec2=fltarr(Nwavs)
+
+Ngrains=(*!dustem_params).ngrains
+spec_grains=fltarr(Nwavs,ngrains) 
+spec2_grains=fltarr(Nwavs,ngrains)
+
+
+;==== save everything to put back later
+saved_pd=0 & saved_pinfo=0 & saved_fpd=0
+IF ptr_valid((*!dustem_fit).param_descs) THEN BEGIN
+  pd=*(*!dustem_fit).param_descs
+  iv=*(*!dustem_fit).param_init_values
+  cv=*(*!dustem_fit).current_param_values
+  saved_pd=1
+ENDIF
+IF ptr_valid(!dustem_parinfo) THEN BEGIN
+   parinfo=*(!dustem_parinfo)
+  saved_pinfo=1
+ENDIF
+IF ptr_valid((*!dustem_fit).fixed_param_descs) THEN BEGIN
+  fpd=*(*!dustem_fit).fixed_param_descs
+  fiv=*(*!dustem_fit).fixed_param_init_values
+  saved_fpd=1
+ENDIF
  
 ;==== Run dustem with G0(a.k.a. U)=1.
 param_desc=['(*!dustem_params).G0']
 param_values=[Umin]
-
 dustem_init_params,!dustem_model,param_desc,param_values
  
 ;stop
 st=dustem_run(param_values)
 spec1=st.sed.em_tot
 spec=(1.-gamma)*spec1
+for k=0,ngrains-1 do begin
+   this_spec=st.sed.(k+1)    ; position 0 is the wavelengths
+   spec_grains[*,k]=(1.-gamma)*this_spec
+endfor
  
 ;==== add contribution of dust seeing a power law of U
 spec2[*]=0.
@@ -128,7 +156,9 @@ NUs=10L   ;This seems sufficient
 ;U values over which integration will be performed
 Us=range_gen(NUs,[Umin,Umax],/log)
 specs=fltarr(NUs,Nwavs)
+specs_grains=fltarr(NUs,Nwavs,Ngrains)
 ;stop
+
 ffacts=fltarr(NUs)
 IF gamma NE 0. THEN BEGIN
     fact=gamma*(alpha-1.)/(Umin^(1.-alpha)-Umax^(1.-alpha))
@@ -138,29 +168,44 @@ IF gamma NE 0. THEN BEGIN
         st=dustem_run(param_values)
         ffacts[i]=Us[i]^(-1.*alpha)*fact
         specs[i,*]=ffacts[i]*st.sed.em_tot
+        for k=0,ngrains-1 do specs_grains[i,*,k]=ffacts[i]*st.sed.(k+1)
     ENDFOR
     ;integrate over U
     ;stop
     FOR j=0L,Nwavs-1 DO BEGIN
         integ=integral(Us,reform(specs[*,j]),Umin,Umax)
         spec[j]=spec[j]+integ[0]
+        for k=0,ngrains-1 do begin
+           integ_grains=integral(Us,reform(specs_grains[*,j,k]),Umin,Umax)
+           spec_grains[j,k]=spec_grains[j,k]+integ_grains[0]
+        endfor
     ENDFOR
 ENDIF
  
-;stop
-
-;This is to update keyword st with final spectrum. Useless, as there is no st keyword !
+;This is to update keyword st with the final spectra. 
 ;Not necessary here, but to be used in plugins with scope 'REPLACE_FORTRAN' which actually do not call the Fortran at all.
 ;Ngrains=3
 ;st=dustem_make_empty_fortran_st(Ngrains)
 st.sed.em_tot=spec
+for k=0,ngrains-1 do st.sed.(k+1)=spec_grains[*,k]
+
+;restore params to their original values
+if saved_pd eq 1 then begin
+  (*!dustem_fit).param_descs=ptr_new(pd)
+  (*!dustem_fit).param_init_values=ptr_new(iv)
+  (*!dustem_fit).param_init_values=ptr_new(cv)
+end
+if saved_fpd eq 1 then begin
+  (*!dustem_fit).fixed_param_descs=ptr_new(fpd)
+  (*!dustem_fit).fixed_param_init_values=ptr_new(fiv)
+end
+if saved_pinfo eq 1 then begin
+   !dustem_parinfo=ptr_new(parinfo)
+end
  
-;dustem_print_params
-;stop
  
 the_end:
 scope='REPLACE_FORTRAN'
-
 ;parameter tags. 
 paramtag=['gamma','alpha','Umin','Umax']