FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra

;+
; NAME:

;    dustem_mpfit_run_vg

; PURPOSE:
;    function used to fit dustem model with dustem_mpfit_sed
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:

;    res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])

; INPUTS:
;    x         = necessary for interface but unused
;    params    = Dustem Model parameters normalized to initial values
; OPTIONAL INPUT PARAMETERS:
;    cc_sed    = SED color correction factors
;    _extra    = extra parameters for the plot routine
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    SED and model are plotted if !dustem_show_plot
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES

;    res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])

; MODIFICATION HISTORY:
;    Adapted to extinction and polarization by Vincent Guillet (IAS)
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;note that x is not used.

;COMMON MPFIT_ERROR, ERROR_CODE

IF keyword_set(help) THEN BEGIN

  doc_library,'dustem_mpfit_run_vg'

  sed=0.
  goto,the_end
ENDIF

p_dim=p_min*(*(*!dustem_fit).param_init_values)

;Run dustem with as an interface to mpfitfun

IF !dustem_idl_continuum THEN cont=dustem_create_continuum()

;CALL THE DUSTEM_CREATE FUNCTIONS
;Actually, res is never used
;THIS line below was added by NF, but removed by JPB on Dec 5 2009
;dustem_set_params,p_dim         ;NF on October 2009 (params have to be updated before calling some functions !)
;dustem_create_functions,p_dim/(*!IDO),cont=cont,res=res
dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,res=res

;if ISRF has been modified and IONFRAC has not been modified, then call dustem_create_ionfrac
;NF on October 2009 : handling of DUSTEM_CREATE_IONFRAC has been changed
;this call only if IONFRACPAH master-keyword in GRAIN.DAT
;Removed for now (JPB). Use of ionfrac should be fully revisited.
IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()

;RUN THE MODEL AND GET THE SPECTRUM
st = dustem_run(p_dim)

st_model=dustem_read_all(!dustem_soft_dir)

chi2_sed      = 0  
chi2_ext      = 0  
chi2_polext   = 0 

chi2_polsed   = 0 
chi2_polfrac   = 0 

rchi2_sed     = 0
rchi2_ext     = 0
rchi2_polext  = 0

rchi2_polsed  = 0
rchi2_polfrac  = 0

wrchi2_sed    = 0
wrchi2_ext    = 0
wrchi2_polext = 0

wrchi2_polsed = 0
wrchi2_polfrac = 0

n_sed         = 0
n_ext         = 0
n_polext      = 0

n_polsed      = 0
n_polfrac      = 0

win = 0

tagnames=tag_names(!dustem_data)
dustem_out = [0]

ind_data_tag=dustem_data_tag(count=count_data_tag)

;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1

for i_tag=0,count_data_tag-1 do begin

;stop

;print,i_tag,ind_data_tag(i_tag)
;print,tagnames(ind_data_tag(i_tag))
    CASE tagnames(ind_data_tag(i_tag)) OF

	'SED'	: BEGIN

        dustem_sed = dustem_compute_sed(p_dim,st,cont=cont) 

		; We normalize to 353
		;x=min(((*!dustem_data.sed).wav-850)^2,j353)
		;dustem_sed = dustem_sed/dustem_sed(j353)*(*!dustem_data.sed).values(j353)

        chi2_sed   = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)

		n_sed = n_elements((*!dustem_data.sed).values)

       	rchi2_sed  = chi2_sed / n_sed
       	wrchi2_sed  = rchi2_sed * !fit_rchi2_weight.sed
       	message,"chi2 SED = "+strtrim(chi2_sed,2)+"  rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed

		;print,'wl ',(*!dustem_data.sed).wav
		;print,"Delta SED per wl: ",(*!dustem_data.sed).values-dustem_sed
		;print,"error per wl: ",(*!dustem_data.sed).sigma
		;print,"chi2 per wl: ",((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2

        alpha = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $

              		/ total(((*!dustem_data.sed).values)^2/(*!dustem_data.sed).sigma^2)

        beta = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $

             		/ total(dustem_sed^2/(*!dustem_data.sed).sigma^2)
	

        print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta

        	;dustem_plot_nuinu_em,st,dustem_sed,cont,/print_radiance,p_dim=p_dim
               
		; Plot if needed

        IF !dustem_show_plot NE 0 THEN BEGIN
            window,0
            dustem_plot_fit_sed,st,dustem_sed,cont,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
    		;ytit='Total intensity'
            ;xtit=textoidl('\lambda (\mum)')

			;win+=1
			;xr=[10,5e3]
			;yr=[5e-33,5.e-28]
			;nodiff = 0
			;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1,nodiff=nodiff
			;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2,nodiff=nodiff

	 	ENDIF

        dustem_out = [ dustem_out, dustem_sed ]

	END

	'EXT' 	: BEGIN 

		dustem_ext = interpol(st.ext.ext_tot,st.ext.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;," weighted rchi2 EXT=",wrchi2_ext
        	;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2

		; Plot if needed
		IF !dustem_show_plot NE 0 THEN BEGIN

			win+=1
                	;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
			xr=[0.3,40]
			yr=[1.e-26,1.e-21]
			dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
			dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
	
			win+=1
                	;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
			xr=[0.2,10]
			yr=[0,3e-21]
			dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
			dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2

		ENDIF

		dustem_plot_ext,st,/print_Rv

	 	dustem_out = [ dustem_out, dustem_ext ]

	 	END

	'POLEXT': BEGIN

		IF !dustem_show_plot NE 0 and !run_circ then begin
                   win+=1
                 
			dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03]
		endif

		IF !run_lin then begin

			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_polext  = rchi2_polext * !fit_rchi2_weight.polext
        		print,"chi2 POLEXT = ",chi2_polext,"  rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
        		;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2

        		; Fit in (pmax,lmax,K)
        		x=st.polext.wav
        		logx=alog(x)
        		logy=alog(st.polext.ext_tot)
        		;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)

        		; Whittet+1992 (Table 1)
        		Uband = 0.36
        		Bband = 0.43
        		Vband = 0.55
        		Rband = 0.63
        		Iband = 0.78
        		Jband = 1.21
        		Hband = 1.64
        		Kband = 2.04
        		bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
        		nbands = n_elements(bands)
        		jband = indgen(nbands)*0
        		for k=0,nbands-1 do begin
                		xx=min(abs(st.polext.wav-bands(k)),j)
                		jband(k) = j
        		endfor
        		;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
        		res = poly_fit(logx(jband),logy(jband),2)
		
        		print
        		print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
        		print,"-------------------------"
        		K    = -res(2)
        		lmax = exp(res(1)/(2*K))
        		pmax = exp(res(0)+K*alog(lmax)^2)
        		print,"* lmax = ",lmax," micron"
        		print,"* K    = ",K
        		print,"* pmax = ",pmax," for NH=1e21"
        		print
		
			; Plot if needed
			IF !dustem_show_plot NE 0 THEN BEGIN

       				; Polarization fraction in the UV-IR
                           win+=1
                           
       				dustem_plot_polar,st,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,4],/xsty,win=win
				
   				; Polarisation curve (Serkowski)
				win+=1
				xr=[0.1,50]
				yr=[5e-25,5e-23]
				dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,/donotclose,multi=1
				dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,multi=2
	
				win+=1
				dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
				
                                
                             ENDIF

			;dustem_plot_polar,st,/print_ratio,cont=cont
                     
			dustem_out = [ dustem_out, dustem_polext ]

		ENDIF

		END

	'POLFRAC': BEGIN

		IF !run_lin then begin

                dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont) 
		dustem_sed_tmp = dustem_compute_sed(p_dim,st,cont=cont) 
		dustem_sed = dustem_polsed
		ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
		for i = 0, count_filt-1 do begin
			j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
			if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
		endfor 
		dustem_polfrac = dustem_polsed/dustem_sed

		; We normalize at 353GHz
		;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
		;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)

        	chi2_polfrac   = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
		n_polfrac = n_elements((*!dustem_data.polfrac).values)
        	rchi2_polfrac  = chi2_polfrac / n_polfrac
        	wrchi2_polfrac  = rchi2_polfrac * !fit_rchi2_weight.polfrac
        	print,"chi2 POLFRAC = ",chi2_polfrac,"  rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
        	;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2

                dustem_out = [ dustem_out, dustem_polfrac ]

		; Plot if needed
		IF !dustem_show_plot NE 0 THEN BEGIN

			; Also show prediction for polarized P/I
			win+=1
       			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,yr=[0,0.25],/ysty,xr=[10,5d3]

	 	ENDIF

	 	ENDIF

		END

	'POLSED': BEGIN

		IF !run_lin then begin

                dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont) 

		; We normalize to 353
		;x=min(((*!dustem_data.polsed).wav-850)^2,j353)
		;dustem_polsed = dustem_polsed/dustem_polsed(j353)*(*!dustem_data.polsed).values(j353)

        	chi2_polsed   = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
		n_polsed = n_elements((*!dustem_data.polsed).values)
        	rchi2_polsed  = chi2_polsed / n_polsed
        	wrchi2_polsed  = rchi2_polsed * !fit_rchi2_weight.polsed
        	print,"chi2 POLSED = ",chi2_polsed,"  rchi2 POLSED = ",rchi2_polsed;," weighted rchi2 POLSED=",wrchi2_polsed
        	print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
        	;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma

                dustem_out = [ dustem_out, dustem_polsed ]

		; Plot if needed
		IF !dustem_show_plot NE 0 THEN BEGIN

			win+=1
			xr=[10,5e3]
			yr=[1e-34,1e-28]
			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,xr=xr,yr=yr,/xsty,/donotclose,multi=1
			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,xr=xr,yr=yr,/xsty,multi=2

		ENDIF

		ENDIF

		END

        ELSE :	stop

    ENDCASE

endfor

; 4 bands fit : IRAS 100 and 857 down to 353 GHz
;fit_sed,"~/tmp/out/SED.RES",['IRAS','HFI'],iband=[0,0,0,1,1,1,1,1,1,0],beta=fixbeta,/noplot
DOF = n_elements(p_min)
total_chi2 = 0
total_rchi2 = 0
total_n = 0
if !fit_rchi2_weight.ext ne 0 then begin
	total_n += n_ext
	total_chi2 +=  chi2_ext
	total_rchi2 += rchi2_ext
endif
if !fit_rchi2_weight.sed ne 0 then begin
	total_n += n_sed
	total_chi2 += chi2_sed
	total_rchi2 += rchi2_sed
endif
if !run_pol then begin
if !fit_rchi2_weight.polext ne 0 then begin
	total_n += n_polext
	total_chi2 += chi2_polext
	total_rchi2 += rchi2_polext
endif
if !fit_rchi2_weight.polsed ne 0 then begin
	total_n += n_polsed
	total_chi2 += chi2_polsed
	total_rchi2 += rchi2_polsed
endif
if !fit_rchi2_weight.polfrac ne 0 then begin
	total_n += n_polfrac
	total_chi2 += chi2_polfrac
	total_rchi2 += rchi2_polfrac
endif
endif

total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $
            + chi2_sed * !fit_rchi2_weight.sed 
if !run_pol then begin
	total_wchi2 	+= chi2_polext  * !fit_rchi2_weight.polext $
             		 + chi2_polsed  * !fit_rchi2_weight.polsed $
             		 + chi2_polfrac * !fit_rchi2_weight.polfrac 
endif

print
print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
print,"Nb of data points = ", total_n
print,"DOF = ", DOF
print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext,rchi2_polsed,rchi2_polfrac
print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext,wrchi2_polsed,wrchi2_polfrac
if !run_pol then 	print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ 
else 			print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed
print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF)

;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF)
; Reduced chi2 is like if there were the same nb of data points for each curve
;print,"Mean reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac) * total_n/(total_n-DOF)

; Plot size distribution
win+=1
xr=[0.3,5e3]
yr=[0.1,100]
IF !dustem_show_plot NE 0 THEN dustem_plot_sdist,xr=xr,yr=yr,ps=filename,win=win;,/donotclose

; Suppress the artificial [0] at the beginning 
dustem_out = dustem_out[1:*]

IF defined(extra) THEN BEGIN
  IF keyword_set(extra.plot_it) THEN BEGIN
;    message,'params='+strtrim(p_min,2),/info,/continue
;    stop

    print,'dustem_mpfit_run_vg: Current parameters values:',p_dim

  ENDIF
ENDIF

the_end:

RETURN,dustem_out

END