rest of plotting changes

Ilyes Choubani
1 parent 9388fd62
Showing 4 changed files with 95 additions and 141 deletions Show diff stats
src/idl/dustem_fit_sed_ext_pol_example.pro
src/idl/dustem_plugin_stellar_population.pro
src/idl/dustem_write_isrf_lv.pro
src/idl/dustemwrap_plot_noobj.pro
@@ -322,6 +322,7 @@ dustem_init,model=use_model,/pol;,kwords=[&#39;logn&#39;,&#39;plaw&#39;,&#39;plaw&#39;];,polarization=us
 !dustem_nocatch=1
 !dustem_verbose=use_verbose
 ;!dustem_show_plot=1
+;!dustem_noobj=1
 !EXCEPT=2 ;so I can locate the plotting error...
  
  
@@ -441,16 +442,23 @@ yr_x = [5.00E-8,10]
 xr_m = [1.,5e5]
 yr_m = [5e-8,1.00e6]
  
+;Because we're in emission+extinction 
+;two seperate titles have to be provided
+;if only one is present, the default titles are kept.
+tit_m='Spectral Energy Distribution' ; combinations that work: anthing that starts with 'tit' and ends with '_m' or '_emission'
+tit_x='Dust Optical Depth' ;same for emission but ends with '_x' or '_extinction'. Not case sensitive.
  
-tit='Spectral Energy Distribution'
+
+;These two are not needed anymore for any of the plottings.... 
 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
 xtit=textoidl('\lambda (\mum)')
+
 ;Set show_plot to 0 to hide plot
 ;Commented or set to 1 is the same since !dustem_show_plot (existing sysvar) is initialized to 1 in dustem_init
 ;show_plot = 0 
 t1=systime(0,/sec)
 res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
-                      ,/xlog,/ylog,xr=xr,yr=yr,xr_m=xr_m,yr_m=yr_m,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,title=tit $
+                      ,/xlog,/ylog,xr_m=xr_m,yr_m=yr_m,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,tit_m=tit_m,tit_x=tit_x $
                       ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
                       ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
 t2=systime(0,/sec)
@@ -465,9 +473,9 @@ 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,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (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)'
+    dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
   ENDELSE
   IF keyword_set(wait) THEN BEGIN
      message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info
@@ -3590,32 +3590,36 @@ FOR i=0L,n_elements(comp_pop.popid)-1 DO BEGIN ; Looping over all the stellar po
  
     Inu = planck(wave_angstrom,(comp_pop.temperature)[i])*(wave_angstrom)^2/c2a ; ergs/cm2/s/Hz/sr
  
-    ;THIS PART IS ONLY FOR ORION. 
-    ;IT SHOULD NOT BE PROVIDED IN THE RELEASE 
-    if isa(!dustem_current) then begin
+;     ;THIS PART IS ONLY FOR ORION. 
+;     ;IT SHOULD NOT BE PROVIDED IN THE RELEASE 
+;     if isa(!dustem_current) then begin
+;     
+;         idff = where(tag_names(*!dustem_plugin) EQ 'FREEFREE',ct_ff)
+;     
+;         if ct_ff ne 0 then begin
+;     
+;         ;look for freefree temperature in the pd vector (sysvar)
+;         ;we need an index
+;     
+;         id_fftmp = where(strupcase(strmid((*(*!dustem_fit).param_descs),14)) EQ 'FREEFREE_1',ct_fftmp)
+;         if ct_fftmp ne 0 then Te = ((*(*!dustem_fit).current_param_values)[id_fftmp])[0]*(((*(*!dustem_fit).param_init_values))[id_fftmp])[0]
+;            ; test if Te is correct
+;     
+;         Beta2 = 2.e-10*(Te)^(-3/4) ;cm^3/s-1
+;         Rs_p3 = 1.1355d+57 ;Rs^3 for orion nebula
+;         nH_stromgren = sqrt(Rs_p3*4*!pi*Beta2/(3*(!const.sigma)*1.0d-04*((comp_pop.temperature)[i])^4*!pi))
+;     
+;     
+;         Inu *= exp(((*!dustem_current).ext.ext_tot)*(-nh_stromgren)*4/3*Rs_p3^(1/3)/1.0d21)
+;     
+;            
+;         endif
+;     
+;     ENDIF
+    
+    
+    If !dustem_dim and isa(!dustem_current) then Inu *= exp(((*!dustem_current).ext.ext_tot)*(-(*!dustem_HCD))/1.0d21)
  
-        idff = where(tag_names(*!dustem_plugin) EQ 'FREEFREE',ct_ff)
-    
-        if ct_ff ne 0 then begin
-    
-        ;look for freefree temperature in the pd vector (sysvar)
-        ;we need an index
-    
-        id_fftmp = where(strupcase(strmid((*(*!dustem_fit).param_descs),14)) EQ 'FREEFREE_1',ct_fftmp)
-        if ct_fftmp ne 0 then Te = ((*(*!dustem_fit).current_param_values)[id_fftmp])[0]*(((*(*!dustem_fit).param_init_values))[id_fftmp])[0]
-           ; test if Te is correct
-    
-        Beta2 = 2.e-10*(Te)^(-3/4) ;cm^3/s-1
-        Rs_p3 = 1.1355d+57 ;Rs^3 for orion nebula
-        nH_stromgren = sqrt(Rs_p3*4*!pi*Beta2/(3*(!const.sigma)*1.0d-04*((comp_pop.temperature)[i])^4*!pi))
-    
-    
-        Inu *= exp(((*!dustem_current).ext.ext_tot)*(-nh_stromgren)*4/3*Rs_p3^(1/3)/1.0d21)
-    
-           
-        endif
-    
-    ENDIF
  
     stellar_component=stellar_component+(comp_pop.nstars)[i]*omega*Inu ;no *4 because of Lambert's cosine law.
  
@@ -104,6 +104,7 @@ IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN
         ;Let's hope they keep releasing it otherwise we'll have to include it ourselves
         ma_isrf_dir=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
         ma_isrf=dustem_read_isrf(ma_isrf_dir)
+        stop
         final_isrf=final_isrf/G0_mathis[0]+ma_isrf.isrf
         Ncomments+=1
         c=strarr(Ncomments)
@@ -91,71 +91,40 @@ 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. 
+;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=st  ; 0/0 division case?
+
+ENDIF 
  
-  ;# 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?
+;# Emission
+
+if not keyword_set(dustem_sed) and isa((*!dustem_data).sed) then begin ;better test?
  
-      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=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
+  dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
  
-  ;# 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=st,EXT_spec=EXT_spec)
-      
-      if !run_pol && !run_lin then begin
-          
-; 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
+  if !run_pol && !run_lin then begin
+     dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac)
+     toto = dustem_compute_stokes(p_dim,st=st,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em) 
+     dustem_qsed = toto[0]
+     dustem_used = toto[1]   
+  endif         
+endif
+
+;# Extinction
+
+if not keyword_set(dustem_ext) and isa((*!dustem_data).ext) then begin
+    dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec)
+  
+    if !run_pol && !run_lin then begin
  
-ENDIF ELSE BEGIN ;st is provided 
-    ;stop
+      dustem_polext = dustem_compute_polext(p_dim,st=st,polext_spec=POLEXT_spec,spext_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext)
  
-    ;# 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=st,SED_spec=SED_spec)
-        
-        if !run_pol && !run_lin then begin
-            
-            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
+      toto = dustem_compute_stokext(p_dim,st=st,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spe=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) 
+      dustem_qext = toto[0]
+      dustem_uext = toto[1]   
     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=st,EXT_spec=EXT_spec)
-        
-        if !run_pol && !run_lin then begin
-            
-            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
-    
-ENDELSE 
+endif
+  
  
 ;IC: Will come back for this to include all datasets and their associated data.
 ;ADDING THIS BLOCK FOR COMPLETENESS
@@ -231,35 +200,8 @@ frmt2=&#39;(10F10.2)&#39;
 ;Instead of setting the extra tags for plotting in dustemcgwin, we'll place them in this (parent) procedure. 
 tgnms_extra = tag_names(_extra)
  
+dustem_set_plot_range,test_sed, test_ext,_extra= _extra
  
-
-;I have decided to keep track of the command index manually as I did not find any command online that does this.
-cmdind_m = 0 ;for emission
-cmdind_x = 0 ;for extinction
-;trying to do the same thing for parameters and plugins
-cmdind_prms = 0 ;for parameters
-cmdind_plgns = 0 ;for plugins
-
-;ADDING PLUGIN(S) TO SPECTRUM----------------
-
-;these following indices will be used to change the keep the same information in the _extra structure
-;this is an initial solution until I finish coding the _extra-filtering procedure. 
-
-iswinsed = !dustemcgwin_id.sed EQ la_undef()
-iswinext = !dustemcgwin_id.ext EQ la_undef()
-
-;trying to do the same thing for parameters and plugins
-;But we need to know if the extinction or the emission mode is on.
-;so the use of ths pointer will be different than the others
-iswinprms = !dustemcgwin_id.prms EQ la_undef()
-iswinplgns = !dustemcgwin_id.plgns EQ la_undef()
-
-
-;plotstrct = !p
-
-;I think I will test over the presence of the errors because they are not computed at that point.
-
-;errors = ((*!dustem_fit).current_param_errors)
 if isa((*!dustem_fit).current_param_errors) then begin
     errors = (*(*!dustem_fit).current_param_errors)
 endif else begin
@@ -408,18 +350,21 @@ if test_m then begin
  
     ;NB issue here because the user can't modify the title after the dustem_mpfit_data call at the main procedure level (user proc)
     ;Getting the width of the title 
-    tit=!dustem_plot_range.title_m;
-
-    cgtext, 0.3, 0.95,tit+' (RUNNING)',color=0,/normal,charsize=-1,width=thiswidth 
-    widthtot = thiswidth
-    xxpos = (1 - thiswidth)/2
-    yypos = 0.96
  
-    cgtext, 0.3, 0.95,tit,color=0,/normal,charsize=-1,width=thiswidth 
+    tit=!dustem_plot_range.title_m
+    
+    xxpos = 0.5
+    yypos =0.96
+    
+    ;the runs inbetween
+    if !dustem_end EQ 0 then cgtext, xxpos , yypos,tit+' (RUNNING)',color=0,/normal,charsize = 1.6,alignment=0.5
+    
+    ;subsequent iterations to the end of the fit
+    if !dustem_end EQ -1 then cgtext, xxpos , yypos,tit,color=0,/normal,charsize = 1.6,alignment=0.5
+    
+    ;last iteration of the fit
+    if !dustem_end EQ 1 then cgtext, xxpos , yypos,tit+' (FINAL RUN)',color=0,/normal,charsize = 1.6,alignment=0.5
  
-    cgtext, xxpos, yypos,tit,color=0,/normal,charsize = 1.6 & cmdind_m+=1
-    xxpos+=thiswidth*1.6
-    cgtext, xxpos, yypos,' (RUNNING)',color=0,/normal,charsize = 1.6 
  
    if !run_pol then begin
         ;,wback='grey'
@@ -670,28 +615,22 @@ if test_x then begin
     cgtext, xxpos + widthtext , yypos,strtrim(string(rchi2,format=frmt2),2),color=0,/normal,charsize = 1.0
     ;saving the command id to replace it
  
+   ;plotting of title 
+    tit=!dustem_plot_range.title_x
  
-    ;Plotting of the title of the dashboard (plot(s))
-    
-    ;NB: JP is right, the user should be able to set the entirety of the title.
-    ;Make this block aware of the saving of the data  (from the fits table)
+    xxpos = 0.5
+    yypos =0.96
  
+    ;the runs inbetween
+    if !dustem_end EQ 0 then cgtext, xxpos , yypos,tit+' (RUNNING)',color=0,/normal,charsize = 1.6,alignment=0.5
  
-    ;Getting the width of the title 
-    tit=!dustem_plot_range.title_x;
+    ;subsequent iterations to the end of the fit
+    if !dustem_end EQ -1 then cgtext, xxpos , yypos,tit,color=0,/normal,charsize = 1.6,alignment=0.5
  
-    cgtext, 0.3, 0.95,tit+' (RUNNING)',color=0,/normal,charsize=-1,width=thiswidth
-    widthtot = thiswidth
-    xxpos = (1 - thiswidth)/2
-    yypos = 0.96
+    ;last iteration of the fit
+    if !dustem_end EQ 1 then cgtext, xxpos , yypos,tit+' (FINAL RUN)',color=0,/normal,charsize = 1.6,alignment=0.5
  
-    cgtext, 0.3, 0.95,tit,color=0,/normal,charsize=-1,width=thiswidth
-      
-    cgtext, xxpos, yypos,tit,color=0,/normal,charsize = 1.6
-    xxpos+=thiswidth*1.6
-    cgtext, xxpos, yypos,' (RUNNING)',color=0,/normal,charsize = 1.6
  
-   
     if !run_pol then begin
  
         ;position arrays for plot and normalized graph
@@ -853,6 +792,8 @@ if test_x then begin
 ENDIF 
  
  
+if !dustem_end then !dustem_end = -1 ;so that the last iteration is only plotted once.
+
 ;PLOTTING OF THE PARAMETERS
 window,2,xpos=param_wxpos,ypos=param_wypos,xsize=param_wxsize,ysize=param_wysize,title='DUSTEMWRAP '+!dustem_version.version+' (PARAMETERS)' ; keeping the same y dimension of the emission/extinction plot? 
 dustem_plot_dataset, st, p_dim, errors,dataset='PARAMETERS',_extra=_extra ;& cmdind_prms+=1