PRO phangs_plot_fitted_seds,vid, $
							all_wavs, $
							fit_wavs, $
							sed_norm, $
							total_predicted_sed, $
							total_sed_fit, $
							sed_fit, $
							weighted_sed, $
							best_grid_sed, $
							sed_stars, $
							sed_stars_reddened, $
							subtract_star_light=subtract_star_light, $
							plot_dust_only_best_fit=plot_dust_only_best_fit, $
							plot_dust_only_best_weighted_fit=plot_dust_only_best_weighted_fit, $
							show_dustem_stellar_spectrum=show_dustem_stellar_spectrum, $
							show_dustem_redenned_stellar_spectrum=show_dustem_redenned_stellar_spectrum, $
							show_best_total_emission_spectrum=show_best_total_emission_spectrum, $
							show_best_grain_emission_spectra=show_best_grain_emission_spectra

mmvec=[sed_norm.stokesI,total_sed_fit.stokesI,sed_fit.stokesI,best_grid_sed]
ind=where(mmvec NE 0)
yrange=[1.e-4,max(mmvec)]
xtit='Wavelength [mic]'
ytit='Inu for NH=1.e20 H/cm2'
tit='Voronoi bin '+strtrim(vid,2)

order_all=sort(all_wavs)
order_fit=sort(fit_wavs)

use_all_wavs=all_wavs[order_all]
use_all_data=sed_norm[order_all].stokesI
use_total_predicted_sed=total_predicted_sed[order_all]
use_fit_wavs=fit_wavs[order_fit]
use_total_sed_fit=total_sed_fit[order_fit]
use_sed_fit=sed_fit[order_fit]
use_weighted_sed=weighted_sed[order_fit]
use_sed_stars=sed_stars[order_all]
use_sed_stars_reddened=sed_stars_reddened[order_all]

ind_data_ok=where(finite(use_all_data) EQ 1,count_data_ok)

use_all_wavs=use_all_wavs[ind_data_ok]
use_all_data=use_all_data[ind_data_ok]
use_total_predicted_sed=use_total_predicted_sed[ind_data_ok]
;xr=minmax(use_all_wavs[ind_data_ok])
;yr=minmax(use_all_data[ind_data_ok])

;enlarge_range
xr=[min(use_all_wavs)*0.8,max(use_all_wavs)*1.1]
yr=[min(use_all_data)*0.8,max(use_all_data)*1.1]

;stop

;=== to do: colors per missions
color_stellar='blue'
color_reddened_stellar='orange'
color_dust='red'
color_dust_only='pink'
color_total_fit='red'
color_data='black'
color_spectrum_dust_only='red'
color_spectrum_total='black'

;=== This is the full data SED, normalized to NH 
cgplot,use_all_wavs,use_all_data,color=color_data,psym='Filled Circle',/xlog,/ylog,xtit=xtit,ytit=ytit,tit=tit,xr=xr,/xsty,yr=yr,/ysty
;cgoplot,use_all_wavs,use_all_data,color=color_data,linestyle=0

;=== This is the total best fit (star + dust)
cgoplot,use_all_wavs,use_total_predicted_sed.stokesI,color=color_total_fit,psym='Open Square'
;cgoplot,use_all_wavs,use_total_predicted_sed.stokesI,color=color_total_fit,linestyle=0

;== This is the Fitted SED without stellar subtraction, with uncertainties
;cgoplot,use_fit_wavs,use_total_sed_fit.stokesI,color='pink',psym='Filled Circle'
;err_bar,use_fit_wavs,use_total_sed_fit.stokesI,yrms=2.*sqrt(total_sed_fit.sigmaII),color=1

;This is the Fitted SED stellar subtracted, with uncertainties
;cgoplot,use_fit_wavs,use_sed_fit.stokesI,color='blue',psym='Filled Circle'
;cgoplot,use_fit_wavs,use_sed_fit.stokesI,color='blue',linestyle=0
;err_bar,use_fit_wavs,use_sed_fit.stokesI,yrms=2.*sqrt(sed_fit.sigmaII),color=1

IF keyword_set(plot_dust_only_best_fit) THEN BEGIN
	;=== This is the best dust SED fit (nearest in grid)
	cgoplot,fit_wavs,best_sed,color=color_dust_only,psym='Open Square'
	;cgoplot,fit_wavs,best_sed,color='blue',linestyle=0
ENDIF

IF keyword_set(plot_dust_only_best_weighted_fit) THEN BEGIN
	;=== This is the best dust SED fit (likelyhood weighted)
	cgoplot,use_fit_wavs,use_weighted_sed,color=color_dust_only,psym='Open Square'
	;cgoplot,use_fit_wavs,use_weighted_sed,color='blue',linestyle=0
ENDIF

;=== This is the stellar redenned and unredenned SEDs
IF keyword_set(subtract_star_light) THEN BEGIN
	cgoplot,use_all_wavs,use_sed_stars,color=color_stellar,psym='Open Square'
	;cgoplot,use_all_wavs,use_sed_stars,color='orange',linestyle=2
	cgoplot,use_all_wavs,use_sed_stars_reddened,color=color_reddened_stellar,psym='Open Square'
	;cgoplot,use_all_wavs,use_sed_stars_reddened,color='red',linestyle=1
ENDIF

IF keyword_set(show_dustem_stellar_spectrum) THEN BEGIN
	dustem_all_wavs=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
	cgoplot,isrf_wavelengths,show_dustem_stellar_spectrum,color=color_stellar,linestyle=2
	;stop
ENDIF
IF keyword_set(show_dustem_redenned_stellar_spectrum) THEN BEGIN
	dustem_all_wavs=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
	cgoplot,isrf_wavelengths,show_dustem_redenned_stellar_spectrum,color=color_reddened_stellar,linestyle=2
	;stop
ENDIF
IF keyword_set(show_best_total_emission_spectrum) THEN BEGIN
	dustem_all_wavs=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
	ntags=n_tags(show_best_total_emission_spectrum)
	spectrum=fltarr(ntags)
	fact=dustem_get_emission_conversion_factor()
	FOR i=0L,ntags-1 DO spectrum[i]=show_best_total_emission_spectrum.(i)
	;spectrum=spectrum/fact ;==== Why divide by Fact here ???
	cgoplot,dustem_all_wavs,spectrum,color=color_spectrum_dust_only
	;==== overplot total (dust+star) spectrum
	IF keyword_set(show_dustem_stellar_spectrum) THEN BEGIN
		;stop
		dust_and_stars_spectrum=spectrum+interpol(show_dustem_redenned_stellar_spectrum,isrf_wavelengths,dustem_all_wavs)
		cgoplot,dustem_all_wavs,dust_and_stars_spectrum,color=color_spectrum_total
	ENDIF
ENDIF

IF keyword_set(show_best_grain_emission_spectra) THEN BEGIN
	;stop
	Ngrains=n_elements(show_best_grain_emission_spectra)
	dustem_all_wavs=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)
	ntags=n_tags(*show_best_grain_emission_spectra[0])
	fact=dustem_get_emission_conversion_factor()
	colors=dustem_grains_colors(Ngrains,ls,/cgplot)
	FOR k=0L,Ngrains-1 DO BEGIN
		spectrum=fltarr(ntags)
		FOR i=0L,ntags-1 DO spectrum[i]=(*show_best_grain_emission_spectra[k]).(i)
		;spectrum=spectrum/fact ;==== Why divide by Fact here ???
		cgoplot,dustem_all_wavs,spectrum,color=colors[k],linestyle=ls[k]
		print,k,minmax(spectrum)
		;stop
	ENDFOR
	;stop
ENDIF

END