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

; 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

;stop

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

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,''

IF !dustem_noobj THEN BEGIN

	dustemwrap_plot_noobj,p_dim,stp,_extra=_extra

     ENDIF ELSE BEGIN

	dustemwrap_plot,p_dim,stp,_extra=_extra

ENDELSE

the_end:

RETURN, p_dim

END