FUNCTION dustem_brute_force_fit,sed $
							   ,table_name $
							   ,filters $
							   ,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

;+
; 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:

;    None

; 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

; 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
ENDIF

seds=!dustem_grid.seds
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
;dustem_grid={Nsed:Nsed,Nfilters:Nfilters,Nparams:Ntparams,filters:filters,params:table_params,seds:st,pmin_values:table_pmins,pmax_values:table_pmaxs}

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)
	ENDFOR

	total_seds[i]=seds[i].total

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)
	ENDFOR
ENDFOR
;stop

;=== loop over grid lines to compute chi2s

;we should 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
IF keyword_set(normalize) THEN BEGIN
	FOR i=0L,Nseds-1 DO BEGIN

		;facts[i]=total(sed.stokesI)/total(grid_seds[i,*])
		facts[i]=total(sed.stokesI)/total_seds[i]

		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-grid_seds[i,*])^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
;stop

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

params=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.

;==== should check bumping to table edges there

params_hit=intarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
  IF params[i] EQ table_pmins[i] THEN params_hit[i]=-1
  IF params[i] EQ table_pmaxs[i] THEN params_hit[i]=1
ENDFOR

;==== should compute parameter uncertainties there

;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(params[ind,i])
	params_max[i]=max(params[ind,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