initial procedures to plot polsed and polext seperately

Ilyes Choubani
1 parent b7f936fc
Showing 5 changed files with 100 additions and 49 deletions Show diff stats
src/idl/dustem_compute_polext.pro
src/idl/dustem_plot_fit_sed.pro
src/idl/dustem_plot_polext.pro
src/idl/dustem_plot_polsed.pro
src/idl/dustem_sed_plot.pro
@@ -6,8 +6,8 @@ FUNCTION dustem_compute_polext ,p_dim,sti,_extra=extra,out_st=out_st,dustem_qext
 ;written to homogonize DustEM
  
 IF not keyword_set(sti) THEN BEGIN
+  dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)
   sti=dustem_run(p_dim)
-  dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,sti
 ENDIF
  
  
@@ -20,31 +20,39 @@ If ptr_valid(*!dustem_plugin) THEN BEGIN
  ;endif
  ;------------------------------------------
  
- ;MODIFYING SPECTRUM VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE 
-
- for i=0L,n_tags(*!dustem_scope)-1 do begin
-
- if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'MODIFY_POLEXT') then begin
-
- ;making sure it's the stokes plugin that gets used here
-
- iidx=where(tag_names(*!dustem_scope) eq 'STEXT')
+ 
+; ;=============================THIS BLOCK IS WRONG YOU HAVE TO CHANGE IT=============================
+;  
+;  
+;  ;MODIFYING SPECTRUM VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE 
+; 
+;  for i=0L,n_tags(*!dustem_scope)-1 do begin
+; 
+;  if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'MODIFY_POLEXT') then begin
+; 
+;  ;making sure it's the stokes plugin that gets used here
+; 
+;  iidx=where(tag_names(*!dustem_scope) eq 'STEXT')
+; 
+; 
+;  IF iidx eq i THEN BEGIN ;not sure about this condition but you get the spirit
+; 
+;  Q_ext=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
+;  U_ext=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
+;  dustem_qext=interpol(Q_ext,sti.polext.wav,(*!dustem_data.polext).wav)
+;  dustem_uext=interpol(U_ext,sti.polext.wav,(*!dustem_data.polext).wav)
+;  
+;  endif 
+;  endif
+;  endfor
+;  ;endif
+;  ;-------------------------------------------
+; ;=======================================================================================================
  
+ENDIF
  
- IF iidx eq i THEN BEGIN ;not sure about this condition but you get the spirit
  
- Q_ext=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
- U_ext=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
- dustem_qext=interpol(Q_ext,sti.polext.wav,(*!dustem_data.polext).wav)
- dustem_uext=interpol(U_ext,sti.polext.wav,(*!dustem_data.polext).wav)
- 
- endif 
- endif
- endfor
- ;endif
- ;-------------------------------------------
  
-ENDIF
 dustem_polext = interpol(sti.polext.ext_tot,sti.polext.wav,(*!dustem_data.polext).wav) 
  
  
@@ -167,13 +167,13 @@ IF count_filt NE 0 THEN BEGIN
   plotsym,8
   xx=((*!dustem_data.sed).wav)[ind_filt]
   yy=sed[ind_filt]/norm[ind_filt]
-  cgoplot,xx,yy,color=use_col_sed_filt,psym=8,syms=2
+  cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
 ENDIF
 IF count_spec NE 0 THEN BEGIN
   plotsym,0
   xx=((*!dustem_data.sed).wav)[ind_spec]
   yy=sed[ind_spec]/norm[ind_spec]
-  cgoplot,xx,yy,color=use_col_sed_spec,psym=8,syms=0.5
+  cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
 ENDIF
  
  
@@ -247,11 +247,15 @@ IF keyword_set(res) THEN BEGIN
       str=str+textoidl(' \pm ')+string(errors(i),format=frmt1)
     ENDIF
  
-    IF STRUPCASE(strmid(strtrim(parameter_description,2),0,6)) eq 'DUSTEM' THEN BEGIN
-       
     xxpos=use_legend_xpos*0.07
     yypos=use_legend_ypos*0.03
     xyouts,xxpos,yypos,'Model: '+!dustem_model,color=0,/normal,charsize=legend_charsize
+   
+    IF STRUPCASE(strmid(strtrim(parameter_description,2),0,6)) eq 'DUSTEM' THEN BEGIN
+       
+;     xxpos=use_legend_xpos*0.07
+;     yypos=use_legend_ypos*0.03
+;     xyouts,xxpos,yypos,'Model: '+!dustem_model,color=0,/normal,charsize=legend_charsize
  
     xxpos=use_legend_xpos*0.75
     yypos=use_legend_ypos*1.03
-PRO dustem_plot_polext, st, p_dim, aligned=aligned
+PRO dustem_plot_polext,st,p_dim,dustem_polext,aligned=aligned,win=win
  
 IF keyword_set(ps) THEN BEGIN
     set_plot, 'PS'
@@ -9,15 +9,35 @@ ENDIF ELSE BEGIN
 ENDELSE
  
  
-fact = (3.e10/(1.e-4*(st.polsed.wav)))/1.d40
-fact1 = (3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40
+titstq='Polarization Cross-section in extinction'
+ytitstq=textoidl('\sigma_{pol} (cm^2/H) for N_H=10^{20} H/cm^2')
+xtit=textoidl('\lambda^{-1} (\mum^{-1})')
+xr=[1e-4,5.0]
+yr=[1e-5,2.]
+Ngrains=(*!dustem_params).Ngrains
+Nwaves=(size(st.polext.ext_tot))[1]
  
+fact=1E+01
+Pext=st.polext.ext_tot*fact
+normpol1=dustem_polext*fact
  
-for i = 0, Ngrains-1 do begin
- if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1) ;investigate this line of the plotting of the different components in polarization
-endfor
-
+if not keyword_set(st) then begin
+    ;activate the plugins here 
+    st=dustem_run(p_dim)
+    dustem_polext = dustem_compute_polext(p_dim,st,dustem_qext,dustem_uext,Q_ext,U_ext)
+ENDIF
  
+ylog=1.
+cgplot,1/st.polext.wav,Pext,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',charsize=1.3
+cgoplot,1/st.polext.wav,Pext,color='black'
+cgoplot,1/(*!dustem_data.polext).wav,(*!dustem_data.polext).values*fact,color='Tomato',psym=16,symsize=1,thick=2
+err_bar,1/(*!dustem_data.polext).wav,(*!dustem_data.polext).values*fact,yrms=3.*((*!dustem_data.polext).sigma)/2.*fact,color=cgColor('Tomato')
+cgoplot,1/(*!dustem_data.polext).wav,dustem_polext*fact,psym=6,color='red',symsize=0;2
+
+cgplot,1/st.polext.wav,Pext/Pext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3    
+cgoplot,1/st.polext.wav,Pext/Pext,color='black'
+cgoplot,1/(*!dustem_data.polext).wav,(*!dustem_data.polext).values*fact/normpol1,color='Tomato',psym=16,symsize=1,thick=2
+err_bar,1/(*!dustem_data.polext).wav,(*!dustem_data.polext).values*fact/normpol1,yrms=3.*((*!dustem_data.polext).sigma)/2./normpol1*fact,color=cgColor('Tomato')
  
  
 IF keyword_set(ps) THEN BEGIN
@@ -14,12 +14,19 @@ ENDELSE
 ;This is because the default procedures for dustem
  
 titstq='Polarization Intensity'
-ytitstq=textoidl('P_\nu (W/m^2/sr) for N_H=10^20 H/cm^2')
+ytitstq=textoidl('P_\nu (W/m^2/sr) for N_H=10^{20} H/cm^2')
 xtit=textoidl('\lambda (\mum)')
 xr=[10,2e4]
 yr=[1e-30,1e-19]
 Ngrains=(*!dustem_params).Ngrains
 Nwaves=(size(st.polsed.em_tot))[1]
+
+
+if not keyword_set(st) then begin
+    ;activate the plugins here 
+    st=dustem_run(p_dim)
+    dustem_polsed = dustem_compute_polsed(p_dim,st)
+ENDIF
  
  
 fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17
@@ -74,6 +81,7 @@ ENDIF
  
 normpol1=dustem_polsed*fact1
 normpol=P
+;normpol=P*0.+1. ;and then use this to plot the normalized data
  
 ; for i = 0, Ngrains-1 do begin
 ;  if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1) ;investigate this line of the plotting of the different components in polarization
@@ -83,14 +91,14 @@ normpol=P
 ylog=1.
 cgplot,st.polsed.wav,P,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',charsize=1.3
 cgoplot,st.polsed.wav,P,color='black'
-cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*1E-20,color='Tomato',psym=16,symsize=1,thick=2
-err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*1E-20,yrms=3.*((*!dustem_data.polsed).sigma)/2.*1E-20,color=cgColor('Tomato')
+cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
+err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1,yrms=3.*((*!dustem_data.polsed).sigma)/2.*fact1,color=cgColor('Tomato')
 cgoplot,(*!dustem_data.polsed).wav,dustem_polsed*fact1,psym=6,color='red',symsize=2
  
 cgplot,st.polsed.wav,P/P,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3    
 cgoplot,st.polsed.wav,P/P,color='black'
-cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*1E-20/normpol1,color='Tomato',psym=16,symsize=1,thick=2
-err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*1E-20/normpol1,yrms=3.*((*!dustem_data.polsed).sigma)/2./normpol1*1E-20,color=cgColor('Tomato')
+cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
+err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.polsed).sigma)/2./normpol1*1E-20,color=cgColor('Tomato')
  
 ;stop
  
 PRO dustem_sed_plot,p_min,_extra=_extra,function_name=function_name,pol=pol,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png
  
 IF not keyword_set(function_name) THEN BEGIN
-  ;Run dustem with as an interface to mpfitfun
+  ;Run dustem with as an interface to mpfitfun; what does this mean? This whole procedure will have to be changed
 ;stop
   ;== JPB: st used to come out from here
   ;even when not set. Had to replavce by out_st ...
-  
-  ;dustem_sed=dustem_compute_sed(p_min,st=st,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
-  dustem_sed=dustem_compute_sed(p_min,st=st,out_st=out_st)
-  
+ 
+;   act_wins=0   
+;   FOR i=0L,n_tags(!dustem_data)-1 DO BEGIN
+;       IF isa(!dustem_data.(i)) THEN act_wins+=1   
+;   ENDFOR   
+ 
+  dustem_sed=dustem_compute_sed(p_min,st=st,out_st=out_st) 
   st=out_st
   ;dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos
-  
-  IF keyword_set(ps) THEN dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,use_model=use_model ELSE dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos,use_model=use_model
+  win=1;win-act_wins+1
+  IF keyword_set(ps) THEN dustem_plot_fit_sed,st,dustem_sed,win=win,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,use_model=use_model ELSE dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,win=win,legend_xpos=legend_xpos,legend_ypos=legend_ypos,use_model=use_model
  
   ;I do not understand the reason begind the stop below - either way polarization is handled by other procedures
   IF keyword_set(pol) THEN BEGIN
     ;stop
-    ;dustem_polsed=dustem_compute_polsed(p_min,st=pst,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
-    ;dustem_plot_polsed,st,
-    dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win;,_Extra=extra
-    ;dustem_plot_polar,st,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,cont=cont,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth
+    IF isa(*!dustem_data.polsed) THEN BEGIN
+        win+=1
+        dustem_polsed=dustem_compute_polsed(p_min,st=st)
+        dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win,_Extra=extra
+    ENDIF
+    IF isa(*!dustem_data.polext) THEN BEGIN
+        win+=1
+        dustem_polext=dustem_compute_polext(p_min,st=st)
+        dustem_plot_polext, st, p_dim, dustem_polext, aligned=aligned, win=win,_Extra=extra
+    ENDIF
+    
+     ;dustem_plot_polar,st,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,cont=cont,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth
   ENDIF
 ENDIF ELSE BEGIN
   CASE function_name OF