PRO dustem_fit,help=help,noplot=noplot,mode=mode

;+
; NAME:
;    dustem_fit
; PURPOSE:
;    Runs conjointly a fit on extinction, emission and polarization by extinction
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    sed=dustem_fit([/noplot,/help])
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    noplot : do not display figures
; 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:
;    TLS added by D. Paradis.
;    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

; Verbose
defsysv, '!dustem_verbose', 1

data_files = { $
	;	sed: 	!dustem_wrap_soft_dir + 'Data/SEDs/Gal_composite_spectrum.xcat', $
	sed: 	!dustem_wrap_soft_dir + 'Data/SEDs/sedsk-70116_2.xcat', $
        ;        sed: 	!dustem_wrap_soft_dir + 'Data/SEDs/spec_firas_chi2.xcat', $
		ext: 	!dustem_wrap_soft_dir + 'Data/EXTs/spec_sk-70116-2.xcat' $
	;	polext:	!dustem_wrap_soft_dir + 'Data/POLEXTs/Serkowski_DISM.xcat', $
	;	polsed:	'' $
             }

; Local variable, to define global !fit_rchi2_weight
rchi2_weight = { $
		sed:	1.,	$
 		ext:    1.	$
 	;	polext:	1.	$
 	;	polsed:	1. 	$
	       }

; Flags to determine what to fit
defsysv, '!dustem_data', {    $	;Data to fit
 				sed:	ptr_new(),	$
  				ext:	ptr_new()	$
  			;	polext:	ptr_new() 	$
;				polsed:	ptr_new() 	$
			 }

; Define global !fit_rchi2_weight with the same structure and field order like !dustem_data
tagnames = tag_names(!dustem_data)
instr="defsysv, '!fit_rchi2_weight', {"
for i = 0, n_elements(tagnames)-1 do instr+=tagnames(i)+': rchi2_weight.'+tagnames(i)+', '
instr+='empty: 0. }
toto=execute(instr)

;=== Set weight for reduced chi2 ponderation between emission, extinction, polarized extinction and emission ===
; Calculate total according to fileds in !dustem_data structure 
total = 0
for i = 0, n_tags(!dustem_data)-1 do total = total + !fit_rchi2_weight.(i)
; values are relative to total
for i = 0, n_tags(!dustem_data)-1 do !fit_rchi2_weight.(i) /= total

dustem_init,mode=use_mode

if keyword_set(noplot) then defsysv,'!dustem_show_plot', 0

; Check that polarisation has been run (-pol in GRAIN.DAT) if fit on polarization is asked for
if (tag_exist(!dustem_data,'polext') or tag_exist(!dustem_data,'polsed')) and !run_pol eq 0. then begin
	message,"pol option must be set in GRAIN.DAT to fit polarization",/info
	stop
endif

; Check that DTLS has been run (-dtls in GRAIN.DAT) 
if (tag_exist(!dustem_data,'dtls') ) and !run_tls eq 0. then begin
	message,"dtls option must be set in GRAIN.DAT to use the TLS model",/info
	stop
endif

;=== Read observational data ===
for i = 0, n_elements(tagnames)-1 do begin
	instr='!dustem_data.('+strtrim(i,2)+') = dustem_set_data(read_xcat(data_files.'+tagnames(i)+',/silent))'
	toto=execute(instr)
endfor

;=== Parameters for fit ===
pd_mass = [ $
       '(*!dustem_params).grains(0).mdust_o_mh', $    ;PAH0 mass fraction
       '(*!dustem_params).grains(1).mdust_o_mh', $    ;PAH1 mass fraction
       '(*!dustem_params).grains(2).mdust_o_mh', $    ;Small amCBEx
       '(*!dustem_params).grains(3).mdust_o_mh' $     ;Large amCBEx
  ;     '(*!dustem_params).grains(4).mdust_o_mh'  $     ;aSil
     ]

tls=['(*!dustem_params).tls.a_dtls', $
     '(*!dustem_params).tls.lc' $
   ;  '(*!dustem_params).tls.vm' $
]

iv_tls=[10., 8.]

;iv_mass = [ $
;	7.8e-4, $
;	7.8e-4, $
;	1.65e-3, $
;	1.45e-3, $
;	7.8e-2 $         ;7.8e-3
;	]  	; DUST MASSES
;iv_mass = [ $
;	9.3E-4, $
;	0.000449, $
;	2.2e-3 $
;;	0.0008 $
;;	7.8e-4 $         ;7.8e-3
;	]  	; DUST MASSES
iv_mass = [ $
	1.3e-3, $
	1.8e-9, $
	6.4e-3, $
	0.0015 $
	]  	; DUST MASSES

pd_sdist = [ $
;;       '(*!dustem_params).grains(0).amin', $         ;PAH0 minmal size 
;;       '(*!dustem_params).grains(0).amax' $         ;PAH0 maximal size 
       '(*!dustem_params).grains(0).alpha_o_a0', $   ;PAH0 mean size
       '(*!dustem_params).grains(2).alpha_o_a0' $   ;PAH0 mean size

;;       '(*!dustem_params).grains(0).at', $           ;PAH0 sigma size
;       '(*!dustem_params).grains(0).amin', $         ;PAH1 minmal size 
;       '(*!dustem_params).grains(0).amax' $         ;PAH1 maximal size 
;       '(*!dustem_params).grains(2).alpha_o_a0' $   ;PAH1 mean size
;;       '(*!dustem_params).grains(1).at', $           ;PAH1 sigma size
;;       '(*!dustem_params).grains(2).amin', $         ;VSG minimal size
;;       '(*!dustem_params).grains(2).amax', $         ;VSG maximal size
;        '(*!dustem_params).grains(3).amin', $        ;Carbon BG minimal size 
;;       '(*!dustem_params).grains(3).amax', $         ;Carbon BG maximal size
;       '(*!dustem_params).grains(2).alpha_o_a0', $   ;Carbon BG power law slope
;       '(*!dustem_params).grains(3).alpha_o_a0' $   ;Carbon BG power law slope

;       '(*!dustem_params).grains(3).at', $    	     ;Carbon BG threshold size (exp. decay)
;       '(*!dustem_params).grains(3).ac', $           ;Carbon BG cut-off size (exp. decay)
;       '(*!dustem_params).grains(3).gamma', $        ;Carbon BG exp. decay slope
;       '(*!dustem_params).grains(4).amin', $         ;Silicate BG minimal size 
;;       '(*!dustem_params).grains(4).amax', $         ;Silicate BG maximal size
;       '(*!dustem_params).grains(4).alpha_o_a0', $   ;Silicate BG power law slope
;       '(*!dustem_params).grains(4).at', $    	     ;Silicate BG threshold size (exp. decay)
;       '(*!dustem_params).grains(4).ac', $           ;Silicate BG cut-off size (exp. decay)
;       '(*!dustem_params).grains(4).gamma' $        ;Silicate BG exp. decay slope
     ]
iv_sdist = [	$
;3.5000E-08, 1.2000E-07, 6.4000E-08, 1.0000E-01 $	; PAH0
;;3.5000E-08, 1.2000E-07, 6.4000E-08, 1.0000E-01, $	; PAH1
-4.98, 3.12e-7  $
;-4.16 $
;6.4e-8,-3.4	$		; VSG
;;6.0000E-07, 5.0000E-04,-2.6000E+00, 0.7000E-05,  0.7000E-05, 0.9000E+00, $ Carbon BG
;;6.0000E-07, 5.0000E-04,-3.4000E+00, 3.0000E-05,  1.0000E-05, 0.4200E+00  $ Silicate BG
;6.0000E-07, -2.6000E+00, 0.7000E-05,  0.7000E-05, 0.9000E+00, $ Carbon BG
;6.0000E-07, -3.4000E+00, 3.0000E-05,  1.0000E-05, 0.4200E+0   $ Silicate BG
	   ]

pd_pol_alignment = [ $
            '(*!dustem_params).pol(4).atresh', $     ; Threshold size for grain alignment
            '(*!dustem_params).pol(4).pstiff', $     ; Stiffness of transition between non-aligned and aligned grains sizes
       	    '(*!dustem_params).pol(4).plev'  $       ; Maximal alignment efficiency of grains
            ]
iv_pol_alignment = [ 	$
		0.1, 	$
		0.5, 	$
		1.0 	$
		]

pd_G0=['(*!dustem_params).G0'] ;Fixed parameter description
iv_G0=[17.]

pd_cont = ['dustem_create_continuum_2']
iv_cont = [12.e-4]
;stop
;pd = [ pd_pol_alignment, pd_mass ]
;iv = [ iv_pol_alignment, iv_mass ]

; Parameters
;pd = [ pd_G0, pd_mass, pd_pol_alignment, pd_sdist, pd_cont ]
;iv = [ iv_G0, iv_mass, iv_pol_alignment, iv_sdist, iv_cont ]
;pd = [ pd_mass, pd_sdist, pd_cont ]
;iv = [ iv_mass, iv_sdist, iv_cont ]
;pd = [ pd_mass, tls,pd_cont ]
;iv = [ iv_mass, iv_tls,iv_cont ]
pd = [pd_G0,  pd_mass, pd_cont,pd_sdist];, tls]
iv = [iv_G0, iv_mass, iv_cont,iv_sdist];,iv_tls]
;pd = [ pd_cont,tls,pd_G0 ]
;iv = [ iv_cont,iv_tls,iv_G0 ]

Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(0.,Npar)
ulims=replicate(0.,Npar)

; Upper limit on alignment efficiency
;ulimed(1+n_elements(pd_mass)+2) = 1
;ulims(1+n_elements(pd_mass)+2) = 1.
; Upper limit on grain size
; amCBEx
;ulimed(1+n_elements(pd_mass)+n_elements(pd_pol_alignment)+9) = 1
;ulims(1+n_elements(pd_mass)+n_elements(pd_pol_alignment)+9) = 9.9d-4

;Npar=n_elements(pd)
;ulimed=replicate(0,Npar)
;llimed=replicate(1,Npar)
;llims=replicate(0.,Npar)
;ulims=replicate(0.,Npar)
;stop
;     ulimed=[0   , 0.           ,0   ,0, 0, 0,0,1]  ; aucun des parametres n'a une valeur superieure limite
;ulims=[0   , 0.   ,0           ,0, 0, 0, 0,5.e-6]
     llimed=[1 ,         1   ,1, 1,  1., 1,0.,1];,1, 1]  ; tous les parametres ont une valeur inferieure limite
    llims =[0           ,0   ,0,0, 0,0., 0, 1.e-8];,0,3.]


;=== Fixed parameters
;fpd=['(*!dustem_params).G0'] ;Fixed parameter description
;fiv=[1000]
;dustem_init_fixed_params,fpd,fiv
;fpd=['(*!dustem_params).grains(1).mdust_o_mh'] ;Fixed parameter description
;fiv=[1.e-9]
;dustem_init_fixed_params,fpd,fiv

;== SET THE FITTED PARAMETERS
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims

;=== RUN fit
tol=1.e-13
Nitermax=20               ;maximum number of iteration. This is the criterium which will stop the fit procedure

loadct,13
!y.range=[1.e-3,100.]              ;This is to ajust plot range from outside the routine
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax) ;,func_name=func_name,cf_min=cf_min)
t2=systime(0,/sec)
print,'executed in ',t2-t1,' sec'

;st=dustem_read_all_res(!dustem_res,/silent)

END