PRO dustem_read_grid_table2arrays,input_table_name, $
								  Ngrains,model,Ncomb,NFparams,INCSP,use_polar,use_double, $
								  parameter_description,parameter_pmin,parameter_pmax,parameter_pnval,parameter_plog,fparameter_description,fparameter_value, $
								  param_values,Nparams,seds,filters,Nfilters,wavs,Nwavs, $
								  I_tot_vec,ext_tot_vec,I_grains_vec,ext_grains_vec,grain_keywords=grain_keywords

fst=file_info(input_table_name)
IF fst.exists NE 1 THEN BEGIN
	message,'input table '+input_table_name+' not found',/continue
	stop
ENDIF

st0=mrdfits(input_table_name,0,h0)   ;This gets the main header
Ngrains=sxpar(h0,'NGRAIN')
grain_keywords=strarr(Ngrains)
FOR i=0L,Ngrains-1 DO BEGIN
	grain_keywords[i]=sxpar(h0,'KEYWG'+strtrim(i+1,2),count=count)
	IF count EQ 0 THEN grain_keywords[i]='?'
ENDFOR
model=strtrim(sxpar(h0,'MODEL'),2)
Ncomb=sxpar(h0,'NCOMB')
Nparams=sxpar(h0,'NPARAM')
NFparams=sxpar(h0,'NFPARAM')
INCSP=sxpar(h0,'INCSP')    ;weither grid includes emission and extinction spectra
use_polar=sxpar(h0,'POLAR')
use_double=sxpar(h0,'DOUBLE')

parameter_description=strarr(Nparams)
parameter_pmin=fltarr(Nparams)
parameter_pmax=fltarr(Nparams)
parameter_pnval=lonarr(Nparams)
parameter_plog=intarr(Nparams)
fparameter_description=strarr(NFparams)
fparameter_value=fltarr(NFparams)
FOR i=0L,Nparams-1 DO BEGIN
	parameter_description[i]=sxpar(h0,'PNAME'+strtrim(i+1,2))
	parameter_pmin[i]=sxpar(h0,'PMIN'+strtrim(i+1,2))
	parameter_pmax[i]=sxpar(h0,'PMAX'+strtrim(i+1,2))
	parameter_pnval[i]=sxpar(h0,'PNVAL'+strtrim(i+1,2))
	parameter_plog[i]=sxpar(h0,'PLOG'+strtrim(i+1,2))
ENDFOR
FOR i=0L,NFparams-1 DO BEGIN
	fparameter_description[i]=sxpar(h0,'FPNAME'+strtrim(i+1,2))
	fparameter_value[i]=sxpar(h0,'FPVALUE'+strtrim(i+1,2))
ENDFOR

;=== read parameter values
param_values=mrdfits(input_table_name,1,h1)
Nparams=N_tags(param_values)
;=== read seds
seds=mrdfits(input_table_name,2,h2)
filters=tag_names(seds[0])
Nfilters=n_elements(filters)-1    ;removing the TOTAL column
filters=strmid(filters[0:Nfilters-1],1,10000)
wavs=mrdfits(input_table_name,3,h3)
I_tot=mrdfits(input_table_name,4,h4)   ;in MJy/sr
ext_tot=mrdfits(input_table_name,5,h5)   ;in MJy/sr

Nwavs=n_elements(wavs.wav)
I_tot_vec=fltarr(Nwavs,Ncomb)
ext_tot_vec=fltarr(Nwavs,Ncomb)
;Here we assume igrain=0, only to create the emission and extinction arrays
I_grains=replicate(I_tot[0],[Ngrains,Ncomb]) ;in MJy/sr for NH=1e20 H/cm2
ext_grains=replicate(ext_tot[0],[Ngrains,Ncomb])

I_grains_vec=fltarr(Nwavs,Ngrains,Ncomb) ;in MJy/sr for NH=1e20 H/cm2
ext_grains_vec=fltarr(Nwavs,Ngrains,Ncomb)

; search for iu, the first fits extension of the grid with per-grain emission spectrum
found_ext=0
kkk=1L
WHILE found_ext EQ 0 DO BEGIN
	str=sxpar(h0,'UNIT'+strtrim(kkk,2))
	IF str EQ 'Grain1 Emission Spectrum [MJy/sr for NH=1e20 H/cm2]' THEN BEGIN
		found_ext=1
		iu=kkk
	ENDIF ELSE BEGIN
		kkk=kkk+1
  ENDELSE
ENDWHILE
;stop
FOR i=0L,Ngrains-1 DO BEGIN
	aa=mrdfits(input_table_name,iu,hg)    ;in MJy/sr for NH=1e20 H/cm2
	I_grains[i,*]=aa
	iu=iu+1
ENDFOR
;=== This actually assumes that extinction spectra will always be just after emission spectra
FOR i=0L,Ngrains-1 DO BEGIN
	aa=mrdfits(input_table_name,iu,hg)
	ext_grains[i,*]=aa
	iu=iu+1
ENDFOR

;fill in old wavelengths vectors from structures in fits file
FOR k=0L,Nwavs-1 DO BEGIN
  	I_tot_vec[k,*]=I_tot.(k)
  	ext_tot_vec[k,*]=ext_tot.(k)
  	I_grains_vec[k,*,*]=I_grains.(k)
 	ext_grains_vec[k,*,*]=ext_grains.(k)
ENDFOR

;removing un-needed old variables
I_grains=0
ext_grains=0
ext_tot=0
I_tot=0
;seds=0

END