FUNCTION dustem_mpfit_data,tol=tol,Nitermax=Nitermax,gtol=gtol,function_name=function_name,xtol=xtol,_extra=_extra

;+
; NAME:
;    dustem_mpfit_data
; PURPOSE:
;    call MPFIT to fit observational datas using ponderated reduced chi2
;    between extinction, emission, and polarisation as defined in !fit_rchi2_weight
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol) 
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;
; MODIFICATION HISTORY:
;    Written by V. Guillet (2012)
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

IF keyword_set(mode) THEN BEGIN
  use_mode=strupcase(mode)
ENDIF ELSE BEGIN
  use_mode='COMPIEGNE_ETAL2010'    ;Default is last dustem model
ENDELSE


disp_str='====================================================='
IF not keyword_set(function_name) THEN BEGIN
  func_name='dustem_mpfit_run'
ENDIF ELSE BEGIN
  func_name=function_name
ENDELSE

;functargs={plot_it:1}

use_tol=1.e-14   ;Default tolerance value
use_Nitermax=3 ;Default maximum number of iterations

IF keyword_set(tol) THEN use_tol=tol
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax

; Build array concatenating observational datas
; Order : that of !dustem_data (sed, ext, polext, polsed)

if !run_pol then begin
    
    ind_data_tag=dustem_data_tag(count=count_data_tag)
    tagnames=tag_names(!dustem_data)

    testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT'

    data_ind=where(~testpol,ctestpol)
    if ctestpol ne 0 then begin
        match,ind_data_tag,data_ind,ind1,ind2 ;not using match2 because when matching indices, if the first is 0 it doesn't count which is wrong.
        ;I also only need the points that match.
        ;testind1 = where(ind1 ne -1,ctind1)
        ;if ctind1 ne 0 then ind_data_tag = ind_data_tag[testind1]
        ind_data_tag = ind_data_tag[ind1]
        count_data_tag = n_elements(ind_data_tag) 
    endif
endif else ind_data_tag=dustem_data_tag(count=count_data_tag)

fields=['wav','values','sigma']

FOR j=0,n_elements(fields)-1 DO BEGIN

	instruc = fields(j)+' = ([0'	; we add 0 first to allow for ',' after it

	for i=0, count_data_tag-1 do begin $	; Add each data set

		; Adapt sigma to the reduced chi2 ponderation defined in !fit_rchi2_weight between sed, ext, polext and polsed
		; Total_reduced_chi2_driving_mpfit  = sum_i ( chi2_i * (weight_i/n_i) )
		;	n 	: the total number of observational variables in the data
		;	n_i 	: the number of obs. variable in data i 
		;	w_i	: the weight for data i from !fit_rchi2_weight for the ponderation of reduced chi2
		; ex : 	w_sed = 0.8, w_ext = 0.2 , i.e. rchi2_sed counts 4 times as much as rchi2_ext
		;	sed : 314 pts,	ext : 52 pts, i.e. 6 times much more constraints in sed than in ext
		; sigma_sed *= sqrt(314./366/0.8) = 1.03
		; sigma_ext *= sqrt(52./366./0.2) = 0.84
		; By increasing sigma by a factor w, we diminish the weight of the corresponding point by a factor w^2
		if count_data_tag gt 1 and j eq 2 then begin ; sigma ponderation = sqrt(ni/wi)
			if !fit_rchi2_weight.(i) gt 0 then begin
				weight = sqrt(1./!fit_rchi2_weight.(i) ) 
			endif else begin
				weight = 1.d99
			endelse
		endif else begin
			weight = 1.
		endelse
		;print,"WEIGHT",i,j,weight,n_elements((*!dustem_data.(i)).sigma),"/",n_elements(wav)

		instruc += ', (*!dustem_data.(' + strtrim(ind_data_tag(i),2) + ')).'+fields(j)+'*'+strtrim(weight,2)

		;print,"INSTRUC",instruc

	endfor
	instruc += '])(1:*)'	; from 1 to exclude first 0
	;print,'-> INSTRUC: ',instruc

    toto=execute(instruc)	; Variables wav, values & sigma are created

ENDFOR

;stop

nocatch=!dustem_nocatch

p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
                 parinfo=(*!dustem_parinfo), perror=perror, yfit=yfit, $
                 bestnorm=bestnorm, dof=dof,quiet=quiet, status=status, xtol=xtol, niter=niter_completed, ftol=use_tol, $

                 maxiter=use_Nitermax,gtol=gtol,errmsg=errmsg,functargs=_extra,nocatch=nocatch)

message,errmsg,/info
message,'mpfitfun Status='+strtrim(status,2),/info
p_dim=p_min*(*(*!dustem_fit).param_init_values)

p_er_dim=perror*(*(*!dustem_fit).param_init_values)
;PRINT RESULTS
;print,'<°)))><  <°)))><  <°)))><  <°)))><  <°)))><  <°)))><  <°)))><  <°)))><'
print, 'Final parameters : ', p_dim
;=== store best fit values in syst variables
(*!dustem_fit).chi2=bestnorm
(*!dustem_fit).rchi2=bestnorm/dof

(*!dustem_fit).current_param_values=ptr_new(p_dim)
;The above was removed because it seems it was not used anywhere.
;!dustem_fit_sed=ptr_new(dustem_compute_sed(p_dim,st,cont=cont,_extra=extra))
(*!dustem_fit).current_param_errors=ptr_new(p_er_dim)
RETURN, p_dim

END