FUNCTION dustem_mpfit_run,x,p_min,niter=niter_completed,_EXTRA=_extra

;+
; NAME:
;    dustem_mpfit_run
;
; 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
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

;note that x is not used.

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_mpfit_run'
  sed=0.
  goto,the_end
ENDIF

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

IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()

;RUN THE MODEL AND GET THE SPECTRUM

dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st ;st is an output...
st_model=dustem_read_all(!dustem_soft_dir)


chi2_sed      = 0  
chi2_ext      = 0  
chi2_qsed      =0
chi2_used      =0 
chi2_qext      =0
chi2_uext      =0 

rchi2_sed     = 0
rchi2_ext     = 0
rchi2_qsed     =0
rchi2_used     =0
rchi2_qext     =0
rchi2_uext     =0

wrchi2_sed    = 0
wrchi2_ext    = 0
wrchi2_qsed    =0
wrich2_used    =0
wrchi2_qext    =0
wrich2_uext    =0

n_sed         = 0
n_ext         = 0
n_qsed         =0
n_used         =0
n_qext         =0
n_uext         =0

win = 0

tagnames=tag_names((*!dustem_data))
dustem_out = [0]
ind_data_tag=dustem_data_tag(count=count_data_tag)


FOR i_tag=0,count_data_tag-1 DO BEGIN

    CASE tagnames(ind_data_tag(i_tag)) OF
        
		'SED'	: BEGIN 
	        dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
          ;help,*(*!dustem_fit).CURRENT_PARAM_VALUES
          ;stop
	        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
             
            ;These two lines were put here so that the chi2 is updated everytime we plot
            ;Moving it to the bottom so that all chi2 are taken into account  
;  	        (*!dustem_fit).chi2=chi2_sed
;  	        (*!dustem_fit).rchi2=rchi2_sed    
	        
	        dustem_out = [ dustem_out, dustem_sed ]
            
		END

		'EXT' 	: BEGIN     		
            dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec)   
	        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

		 	dustem_out = [ dustem_out, dustem_ext ]

		 	END

		'POLEXT': BEGIN 

			
			;Haven't investigated !run_circ=1 yet... 
			;---------------the procedure within this block has not been updated for a vert long time
			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 = dustem_compute_polext(p_dim,st=st,POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext,dustem_ext=dustem_ext)

			ENDIF
			END

        'QSED': BEGIN 
        		    
		        stokes = 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 = stokes[0]  
		        ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
		        ;dustem_qsed(-1)=dustem_qsed(-2)
; 		        testspass1 = ((*(*!dustem_data).qsed).filt_names)(-1) eq 'SPASS1' 
;  			    if testspass1 then dustem_qsed(-1)=dustem_qsed(-2)
;  			    
			    chi2_qsed = total(((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2) 
			    n_qsed = n_elements((*(*!dustem_data).qsed).values)
			    rchi2_qsed  = chi2_qsed / n_qsed
	        	wrchi2_qsed  = rchi2_qsed * (*!fit_rchi2_weight).qsed 
	            dustem_out = [ dustem_out, dustem_qsed ] ; classic command
	            message,"chi2 QSED = "+strtrim(chi2_qsed,2)+"  rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
	        	print,"chi2 per wl: ",((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2
		        
		    END   

        'USED': BEGIN 
    		    dustem_used = stokes[1] 
 			    chi2_used  = total(((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2) 
 			    n_used = n_elements((*(*!dustem_data).used).values)
 			    rchi2_used  = chi2_used / n_used
  	        	wrchi2_used  = rchi2_used * (*!fit_rchi2_weight).used
 	            dustem_out = [ dustem_out, dustem_used] ; classic command
 	            message,"chi2 USED = "+strtrim(chi2_used,2)+"  rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  	        	print,"chi2 per wl: ",((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2
                
  		END
  		
        'POLSED': BEGIN
 			
             	  IF !run_pol and !run_lin THEN BEGIN
                        ;the keyword for dustem_sed here is used so that the function uses it to compute the polarization fraction instead of
                        ;running dustem_compute_sed again in the procedure. Which it used to do.  
            	       	dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac,dustem_sed=dustem_sed)
            	      	
         		  ENDIF
		END
        	
 		
  		'QEXT': BEGIN ;UPDATED  
                
                Stokext = dustem_compute_stokext(p_dim,st=st,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spec=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) 
                dustem_qext = stokext[0]
                chi2_qext   =total(((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2) 
   			    n_qext = n_elements((*(*!dustem_data).qext).values)
   			    rchi2_qext  = chi2_qext / n_qext
  	        	wrchi2_qext  = rchi2_qext * (*!fit_rchi2_weight).qext
   	            dustem_out = [ dustem_out, dustem_qext] ; classic command
   	            message,"chi2 QEXT = "+strtrim(chi2_qext,2)+"  rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  	        	print,"chi2 per wl: ",((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2 
     	           
  		END
;  		
  		'UEXT': BEGIN ;UPDATED
          		dustem_uext = stokext[1]
 			    chi2_uext   =total(((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2) 
   			    n_uext = n_elements((*(*!dustem_data).uext).values)
   			    rchi2_uext  = chi2_uext / n_uext
  	        	wrchi2_uext  = rchi2_uext * (*!fit_rchi2_weight).uext
   	            dustem_out = [ dustem_out, dustem_uext] ; classic command
   	            message,"chi2 UEXT = "+strtrim(chi2_uext,2)+"  rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  	        	print,"chi2 per wl: ",((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2
            
 	           
  		END
      	
	    ELSE : ;Do nothing
	    
    ENDCASE

ENDFOR

;=== computing chi^2
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).qsed ne 0 then begin
    	total_n += n_qsed
    	total_chi2 += chi2_qsed
    	total_rchi2 += rchi2_qsed
    endif
    if (*!fit_rchi2_weight).used ne 0 then begin
    	total_n += n_used
    	total_chi2 += chi2_used
    	total_rchi2 += rchi2_used
    endif
    if (*!fit_rchi2_weight).qext ne 0 then begin
    	total_n += n_qext
    	total_chi2 += chi2_qext
    	total_rchi2 += rchi2_qext
    endif
    if (*!fit_rchi2_weight).uext ne 0 then begin
    	total_n += n_uext
    	total_chi2 += chi2_uext
    	total_rchi2 += rchi2_uext
    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 $
             		 + chi2_qsed * (*!fit_rchi2_weight).qsed $
             		 + chi2_used * (*!fit_rchi2_weight).used $
             		 + chi2_qext * (*!fit_rchi2_weight).qext $
             		 + chi2_uext * (*!fit_rchi2_weight).uext  
endif

;Using this expression for now for the weighting. 
(*!dustem_fit).chi2 = total_wchi2
(*!dustem_fit).rchi2 = total_wchi2/(total_n-DOF)


message,'======================================================================',/continue
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);total_rchi2/4*total_n/(total_n-DOF) original ; I do not understand this line. Maybe redo some computation.
message,'======================================================================',/continue
;the following is needed under fawlty (not IDL) which message,/reset does not work like in IDL
!Error_State.msg=''

If !dustem_show_plot THEN BEGIN
    ;=== plotting the fit
    if !dustem_iter.act NE !dustem_iter.prv or !dustem_iter.prv EQ 1 then begin
        IF !dustem_noobj EQ 0 THEN BEGIN
          dustemwrap_plot,p_dim,st,dustem_sed,SED_spec,dustem_qsed,Q_spec,dustem_used,U_spec,dustem_polsed,P_spec,dustem_polfrac,SP_spec, dustem_psi_em, PSI_spec,dustem_ext,EXT_spec,dustem_qext,QEXT_spec,dustem_uext,UEXT_spec,dustem_polext,POLEXT_spec,dustem_fpolext,SPEXT_spec,dustem_psi_ext,PSIEXT_spec,_extra=_extra
        ENDIF ELSE BEGIN
          dustemwrap_plot_noobj,p_dim,st,dustem_sed,SED_spec,dustem_qsed,Q_spec,dustem_used,U_spec,dustem_polsed,P_spec,dustem_polfrac,SP_spec, dustem_psi_em, PSI_spec,dustem_ext,EXT_spec,dustem_qext,QEXT_spec,dustem_uext,UEXT_spec,dustem_polext,POLEXT_spec,dustem_fpolext,SPEXT_spec,dustem_psi_ext,PSIEXT_spec,_extra=_extra
        ENDELSE
        if !dustem_iter.act EQ 1 then begin
            !dustem_iter.prv=0
            !dustem_iter.act=0 
        endif else !dustem_iter.prv=!dustem_iter.act
    endif
ENDIF

;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