PRO phangs_compare_seds_isrf,isrf_class,N_seds=N_seds

;phangs_compare_seds_isrf,15,N_seds=1000L

dir=!phangs_data_dir+'/ISRF/GRIDS/'
isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
model='DBP90'
table_name=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+isrf_class_str+'.fits'
st=mrdfits(table_name,1,h)
Nst=n_elements(st)

use_filters_names=['F360M','F0770W','F1000W','F1130W','F2100W']
filters=dustem_filter_names2filters(use_filters_names)
Nfilters=n_elements(filters)

ii=100L

Nseds=100L
IF keyword_set(N_seds) THEN Nseds=N_seds

Nparams=3
true_params=fltarr(Nseds,Nparams)
fitted_params=fltarr(Nseds,Nparams)
seed=1

FOR j=0L,Nseds-1 DO BEGIN
	;seed=1
	inn=abs(randomu(seed))
    ii=inn*Nst
    ;print,inn,ii
	sed=dustem_initialize_sed(Nfilters)
	;sed=fltarr(Nfilters)
	tags=tag_names(st)
	FOR i=0L,Nfilters-1 DO BEGIN
		ind=where(tags EQ 'I'+filters[i],count)
		sed[i].filter=filters[i]
		sed[i].stokesI=st[ii].(ind[0])
	ENDFOR
	sed.sigmaii=sed.stokesI*1.e-5
	FOR i=0L,Nparams-1 DO BEGIN
	  true_params[j,i]=st[ii].(i)
	ENDFOR

	;stop

	fit_isrf_class_str='_isrfclass0'
	table_name_2=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+fit_isrf_class_str+'.fits'

	params=dustem_brute_force_fit(sed,table_name_2,filters,fact=fact,chi2=chi2,/normalize,rchi2=rchi2,show_sed=show_sed $
			   					 ,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max $
			   					 ,fixed_parameters_description=fixed_parameters_description,fixed_parameters_values=fixed_parameters_values)
    fitted_params[j,*]=params

ENDFOR

win=0L
xtit='True parameter ('+isrf_class_str+')'
ytit='recovered parameter (Mathis ISRF)'

window,win & win=win+1
cgplot,true_params[*,0],fitted_params[*,0],tit='G0',psym=4,/xlog,/ylog,xtit=xtit,ytit=ytit
cgoplot,10^!x.crange,10^!y.crange,linestyle=2,color='red',thickness=3

window,win & win=win+1
cgplot,true_params[*,1],fitted_params[*,1],tit='Ypah',psym=4,/xlog,/ylog,xtit=xtit,ytit=ytit
cgoplot,10^!x.crange,10^!y.crange,linestyle=2,color='red',thickness=3

window,win & win=win+1
cgplot,true_params[*,2],fitted_params[*,2],tit='Yvsg',psym=4,/xlog,/ylog,xtit=xtit,ytit=ytit
cgoplot,10^!x.crange,10^!y.crange,linestyle=2,color='red',thickness=3

stop

END