FUNCTION dustem_brute_force_fit,sed $
							   							 ,table_name $
							   							 ,filters $
							   							 ,fixed_parameters_description=fixed_parameters_description $
							   							 ,fixed_parameters_values=fixed_parameters_values $
							   							 ,normalize=normalize $
							   							 ,params_hit=params_hit $
							   							 ,params_uncertainties=params_uncertainties $
							   							 ,params_min=params_min $
							   							 ,params_max=params_max $
							   							 ,fact=fact $
							   							 ,chi2=chi2 $
							   							 ,rchi2=rchi2 $
							   							 ,best_sed=best_sed $
							   							 ,best_grid_sed=best_grid_sed $
							   							 ,show_sed=show_sed $
							   							 ,nostop=nostop $
							   							 ,reset=reset $
							   							 ,weighted_params=weighted_params $
							   							 ,all_chi2s=all_chi2s $
							   							 ,status=status

;+
; NAME:
;    dustem_brute_force_fit
; PURPOSE:
;    does brut force fit from a pre-computed model grid
; CATEGORY:
;    DustEM
; CALLING SEQUENCE:
;    res=dustem_brute_force_fit(sed,table_name,filters[,/normalize][,fact=][,chi2=][,rchi2=][rchi2][best_sed=][,/show_sed][,/nostop][,/reset])
; INPUTS:
;    sed  = sed to be fit. Should contain all filters provided to the routine trhough input variable filter
;    table_name = name of the fits file to be used as a grid (see dustem_make_sed_table.pro for creation)
;    filters    = filters to be used in the grid table (only used if !dustem_grid does not exist or if /reset)
; OPTIONAL INPUT PARAMETERS:
;    fixed_parameters_description = if set, describes parameters of the table that should be considered fixed.
;    fixed_parameters_values      = values of parameters of the table that should be considered fixed.
; OUTPUTS:
;    res = best values of the free parameters. Not that linear parameters to the sed should be multiplied by fact.
; OPTIONAL OUTPUT PARAMETERS:
;    fact = best scaling factor. Such that best_sed*fact is for NH=1.e20 H/cm2 (as in the grid). IF sed is already normalized to NH=1e20 H/cm2, no need to set /normalize
;    chi2 = chi2 of the best fit
;    rchi2 = reduced chi2 of the best fit
;    best_grid_sed = best SED found in the grid
;    best_sed = best fit SED matching input SED (such that best_sed=best_grid_sed*fact)
;    params_hit= 1 if parameter hits the top, -1 if hits the bottom, 0 otherwise
;    params_uncertainties = 90% confidence level uncertainty on parameters
;    params_min = 90% lower limit of parameters
;    params_max = 90% upper limit of parameters
;    weighted_params = likelyhood (1/chi2) weighted best parameter values
; ACCEPTED KEY-WORDS:
;    help                  = if set, print this help
;    normalize = if set, the best match is done based on the shape of the SEDs only, instead of their absolute values. fact will be computed.
;    show_sed = if set, plot the seds
;    nostop   = if set, do not stop
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
;    Scans the grid SEDs to find the best fit (lowest chi2)
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by JPB Jan 2024
;    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_brute_force_fit'
  params=0.
  goto,the_end
ENDIF

;stop
;sed must be normalized to NH=1e20 H/cm2

;==== test for the existence of a loaded grid and initialize if needed.
;==== only filters contained in the filters variable are kept in the grid, and the sum of grid seds are recomputed accordingly
defsysv,'!dustem_grid',0,exist=exist
IF exist EQ 0 or keyword_set(reset)THEN BEGIN
    dustem_define_grid,table_name,filters=filters,status=status
    IF status NE 1 THEN BEGIN
    	params=-1
    	goto,the_end
    ENDIF
ENDIF

;seds=!dustem_grid.seds
seds=(*!dustem_grid).st_seds
params=(*!dustem_grid).st_params
totals=(*!dustem_grid).sed_totals
Nseds=(*!dustem_grid).Nsed
Nfilters=(*!dustem_grid).Nfilters
Nparams=(*!dustem_grid).Nparams
table_params=(*!dustem_grid).params
table_pmins=(*!dustem_grid).pmin_values
table_pmaxs=(*!dustem_grid).pmax_values

;=== marginalize table for fixed parameters, if any
IF keyword_set(fixed_parameters_description) THEN BEGIN
	dustem_marginalize_grid,fixed_parameters_description,fixed_parameters_values,seds,Nseds,params,Nparams,table_params,table_pmins,table_pmaxs
ENDIF

;help,seds,Nseds,params,Nparams,table_params,table_pmins,table_pmaxs
;SEDS            STRUCT    = -> <Anonymous> Array[400]
;NSEDS           LONG      =          400
;PARAMS          STRUCT    = -> <Anonymous> Array[400]
;NPARAMS         LONG      =            3
;TABLE_PARAMS    STRING    = Array[3]
;TABLE_PMINS     DOUBLE    = Array[3]
;TABLE_PMAXS     DOUBLE    = Array[3]

;stop

;=== These are the full arrays needed for calculations
chi2s=dblarr(Nseds)
grid_seds=dblarr(Nseds,Nfilters)
total_seds=dblarr(Nseds)
grid_param_values=dblarr(Nseds,Nparams)

;=== compute sed grid Array
FOR i=0L,Nseds-1 DO BEGIN
	FOR j=0L,Nfilters-1 DO BEGIN
		;grid_seds[i,j]=seds[i].(j+Nparams)
		grid_seds[i,j]=seds[i].(j)
	ENDFOR
	;total_seds[i]=seds[i].total
	total_seds[i]=totals[i]
ENDFOR
;=== compute the parameter values Array
FOR i=0L,Nseds-1 DO BEGIN
	FOR j=0L,Nparams-1 DO BEGIN
		;grid_param_values[i,j]=seds[i].(j)
		grid_param_values[i,j]=params[i].(j)
	ENDFOR
ENDFOR
;stop

;=== loop over grid lines to compute chi2s
;we do actually store the total of grid sed for each grid sed into the grid fits file
;we could also store the total of observed seds in the observed sed save fits file (!)
facts=dblarr(Nseds)
facts[*]=1.d0
;stop
IF keyword_set(normalize) THEN BEGIN
	FOR i=0L,Nseds-1 DO BEGIN
		ref_sed=reform(grid_seds[i,*])
		;total_ref_sed1=total_seds[i]
		total_ref_sed2=total(ref_sed)
		;diff=abs(total_ref_sed1-total_ref_sed2)
		;print,diff
		;total_ref_sed=total_ref_sed1    ;This does not work. I'm not sure why, as it is supposed to be the same pre-computed ...
		total_ref_sed=total_ref_sed2    ;This does !!		
		;IF diff GT 0.1 THEN stop
		facts[i]=total(sed.stokesI)/total_ref_sed
    this_ref_sed=ref_sed*facts[i]
    this_sed=sed.stokesI
    this_var=sed.sigmaII
		;this_sed=sed.stokesI/facts[i]
		;this_var=sed.sigmaII/facts[i]^2
	  ;chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
	  chi2s[i]=total((this_sed-this_ref_sed)^2/this_var)
	ENDFOR
ENDIF ELSE BEGIN
	FOR i=0L,Nseds-1 DO BEGIN
	   chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
	ENDFOR
ENDELSE

all_chi2s=chi2s
chi2=min(chi2s)
ind=where(chi2s EQ chi2,count)

;stop
;print,facts[14406],facts[ind[0]]
;print,chi2s[14406],chi2s[ind[0]]

params=reform(grid_param_values[ind[0],*])  ;these are the best fit parameters
fact=facts[ind[0]]   ; This is the corresponding scaling factor
Nfreedom=Nfilters-Nparams-1   ;This is the degree of freedom
rchi2=chi2/Nfreedom        ;This is the reduced chi2
best_grid_sed=reform(grid_seds[ind[0],*])   ;This is the best SED fround in the grid.
best_sed=best_grid_sed*fact  ;This is the best scaled grid SED matching the input SED.

;=== compute Likelyhood weighted params values
weighted_params=dblarr(Nparams)
weights=1.d0/chi2s
tot_weights=total(weights)
FOR i=0L,Nparams-1 DO BEGIN
	weighted_params[i]=total(reform(1.d0*grid_param_values[*,i])*weights)/tot_weights
ENDFOR

;stop
;==== Check bumping to table edges.
params_hit=intarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
	diff_min=abs(params[i]-table_pmins[i])/table_pmins[i]*100.
	diff_max=abs(params[i]-table_pmaxs[i])/table_pmaxs[i]*100.
	IF diff_min LE 1. THEN params_hit[i]=-1
	IF diff_max LE 1. THEN params_hit[i]=1
  ;IF params[i] LE table_pmins[i] THEN params_hit[i]=-1
  ;IF params[i] GE table_pmaxs[i] THEN params_hit[i]=1
ENDFOR

;==== Compute parameter uncertainties
;stop
conf_level=90./100.
;print,Nfreedom,conf_level
IF Nfreedom NE 0. THEN BEGIN
	dchi2=delta_chi2(Nfreedom,conf_level)
ENDIF ELSE BEGIN
	message,'Nfreedom = 0. Beware of meaningless fit',/continue
	stop
	message,'seting Nfreedom = 1. Beware of meaningless fit',/continue
  dchi2=delta_chi2(1.,conf_level)
ENDELSE
ind=where(chi2s LE chi2+dchi2,count)
params_min=fltarr(Nparams)
params_max=fltarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
	params_min[i]=min(grid_param_values[ind[0],i])
	params_max[i]=max(grid_param_values[ind[0],i])
ENDFOR
params_uncertainties=params_max-params_min

IF keyword_set(show_sed) THEN BEGIN
	xtit='Wavelength [mic]'
	ytit='SED (model=red,data=blue)'
	;should also plot uncertainties
	minv=la_min([sed.stokesI,best_sed])
	maxv=la_max([sed.stokesI,best_sed])
	cgplot,sed.wave,sed.stokesI,psym='Filled Circle',color='blue',/ylog,/ysty,xtit=xtit,ytit=ytit,/xlog,yrange=[minv,maxv]
	cgoplot,sed.wave,best_sed,psym='Filled Square',color='red'
	if not keyword_set(nostop) THEN stop
ENDIF

the_end:

RETURN,params
END