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

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

;=== Should marginalize table for fixed parameters, if any
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

IF keyword_set(fixed_parameters_description) THEN BEGIN
	dustem_marginalize_grid,fixed_parameters_description,fixed_parameters_values,seds,Nseds,Nparams,table_params,table_pmins,table_pmaxs
ENDIF

;=== 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)
	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 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
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.

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

;==== Check bumping to table edges.
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

;==== 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(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