PRO plot_lastrun_vg,run=run

if keyword_set(run) then dustem_init,/wraptest

loadct,13

defsysv,'!dustem_inputs',ptr_new()
dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT'} 
!dustem_inputs=ptr_new(dustem_inputs)

;=== Read model ===
st_model=dustem_read_all(!dustem_soft_dir)
defsysv,'!dustem_params' , ptr_new(st_model)
st=dustem_read_all_res(!dustem_res,/silent)

;=== Read observational data ===
data_files = { $
		sed: 	!dustem_wrap_soft_dir + 'Data/SEDs/Gal_composite_spectrum.xcat', $
		ext: 	!dustem_wrap_soft_dir + 'Data/EXTs/Mathis1990_DISM.xcat', $
		polext:	!dustem_wrap_soft_dir + 'Data/POLEXTs/Serkowski_DISM.xcat', $
		polsed:	'' $
	     }

defsysv, '!dustem_data', {    $ ;Data to fit
                                sed:    ptr_new(),      $
                                ext:    ptr_new(),      $
                                polext: ptr_new()       $
;                                polsed: ptr_new()       $
                         }

tagnames = tag_names(!dustem_data)
for i = 0, n_elements(tagnames)-1 do begin
        instr='!dustem_data.('+strtrim(i,2)+') = dustem_set_data(read_xcat(data_files.'+tagnames(i)+',/silent))'
        toto=execute(instr)
endfor

if tag_exist(st,"sed") then begin

	window,1
	;data_sed = read_xcat(data_files.sed,/silent)
	fact = 1./4./!pi*1.e4/1.e7*1.e20 
	dustem_plot_nuinu_em,st,yr=[1e-11,1.e-7],/ysty,xr=[1,5e3],/xsty

	dustem_sed = interpol(st.sed.em_tot*fact,st.sed.wav,(*!dustem_data.sed).wav)
	chi2_sed   = total(((*!dustem_data.sed).values*(3.e8/1.e-6/(*!dustem_data.sed).wav)/1.d20-dustem_sed)^2/((*!dustem_data.sed).sigma*(3.e8/1.e-6/(*!dustem_data.sed).wav)/1.d20)^2)
	n_sed = n_elements((*!dustem_data.sed).values)
	rchi2_sed  = chi2_sed / n_sed
	;wrchi2_sed  = rchi2_sed * !fit_rchi2_weight.sed
	print,"chi2 SED = ",chi2_sed,"  rchi2 SED = ",rchi2_sed

endif

if tag_exist(st,"ext") then begin

	window,2
	dustem_plot_ext,st,st_model,yr=[1.e-6,1.e0],/ysty,xr=[0.3,400],/xsty

	window,3
	dustem_plot_ext,st,st_model,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty

	dustem_ext = interpol(st.ext.ext_tot,st.sed.wav,(*!dustem_data.ext).wav)
	chi2_ext   = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
	n_ext = n_elements((*!dustem_data.ext).values)
	rchi2_ext  = chi2_ext / n_ext
	;wrchi2_ext  = rchi2_ext * !fit_rchi2_weight.ext
	print,"chi2 EXT = ",chi2_ext,"  rchi2 EXT = ",rchi2_ext

endif

if tag_exist(st,"polext") and tag_exist(st,"polsed") then begin

	; Polarization fraction in the UV-IR
	window,4
	dustem_plot_polar,st,st_model,yr=[0,0.08],/ysty,xr=[0.1,4],/xsty
	
	; Polarisation curve (Serkowski)
	window,5
	dustem_plot_polar,st,st_model,/UV,yr=[1e-3,3e-2],/ysty,xr=[0.1,10],/xsty
	
	; Polarization fraction from the UV to the submm, extinction and emission
	window,6
	dustem_plot_polar,st,st_model,yr=[0,0.30],/ysty,xr=[0.1,5e3],/xsty

	dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
	chi2_polext   = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
	n_polext = n_elements((*!dustem_data.polext).values)
	rchi2_polext  = chi2_polext / n_polext
	;wrchi2_sed  = rchi2_sed * !fit_rchi2_weight.polext
	print,"chi2 POLEXT = ",chi2_polext,"  rchi2 POLEXT = ",rchi2_polext

endif

dustem_plot_sdist

stop

END