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 fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;
; MODIFICATION HISTORY:
;    Written by V. Guillet (2012)
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_mpfit_data'
  goto,the_end
ENDIF

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
previous_show_plot = !dustem_show_plot

; Build array concatenating observational data

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' or tagnames EQ 'FPOLEXT'
    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.
        ;we also only need the points that match.
        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

IF tag_exist(_extra,'show_plot') THEN !dustem_show_plot = _extra.show_plot
IF ~tag_exist(_extra,'show_plot') THEN !dustem_show_plot = 0

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)

;stop                 
;!dustem_iter.act=niter_completed

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

;PRINT RESULTS
print,''
print,''
print,'<o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><'
print,'<o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><'
print,''
print,'FINISHED DUSTEMWRAP FITTING...'
print,''
print,'The final values of the free parameters (pd) are : '
print, p_dim
print,''
print,'<o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><'
print,'<o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><  <o)))><'
print,''
print,''


; we assume that the user knows if they want plots or not
IF !dustem_show_plot gt 0 then begin
   IF !dustem_noobj THEN BEGIN
      dustemwrap_plot_noobj,p_dim,st=stp,_extra=_extra
   ENDIF ELSE BEGIN
      dustemwrap_plot,p_dim,st=stp,_extra=_extra
   ENDELSE
end

!dustem_show_plot = previous_show_plot
the_end:

RETURN, p_dim

END