PRO dustem_add_linear_params2grid,input_table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=out_filename,plog=plog,show_seds=show_seds,help=help

;+
; NAME:
;    dustem_add_linear_params2grid
; PURPOSE:
;    Adds entries for linear parameters to a Dustemwrap grid fits file
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_add_linear_params2grid,input_table_name,pd,iv_min,iv_max,iv_Nvalues[,filename=][,plog=][,/show_seds][,/help]
; INPUTS:
;    input_table_name       : dustemwrap model name to be used
;    pd                     : dustemwrap parameter description array. Can only have one element for now
;    iv_min                 : minimum values for parameter in pd
;    iv_max                 : maximum values for parameter in pd
;    iv_Nvalues             : number of parameter values for parameter in pd
; OPTIONAL INPUT PARAMETERS:
;    plog                   : array indicating if the parameter is to be sampled in linear (0) or log scale (1)
;    out_filename           : output fits file name for the table (default ='/tmp/dustem_seds.fits')
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    A file is written
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
;    
; EXAMPLES
;	dustem_init,model='DL07'
;	input_table_name='/Volumes/PILOT_FLIGHT1/PHANGS//ISRF/GRIDS/TEST_DL07_MuseISRF_JWST_G0_4Phangs_noionis_isrfclass15.fits'
;	pd=['(*!dustem_params).grains[0].MDUST_O_MH']
;	default_value=(dustem_get_param_values_default(pd))[0]
;	iv_min=[default_value] & iv_max=[default_value*3.] & iv_Nvalues=[3] & plog=[0]
;	output_table_name='/Volumes/PILOT_FLIGHT1/PHANGS//ISRF/GRIDS/TEST_DL07_MuseISRF_JWST_G0_YPAH_4Phangs_noionis_isrfclass15_essai.fits'
;	dustem_add_linear_params2grid,input_table_name,pd,iv_min,iv_max,iv_Nvalues,filename=output_table_name,plog=plog,/show_seds
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (May 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_add_linear_params2grid'
  goto,the_end
ENDIF

use_plog=0
if keyword_set(plog) THEN use_plog=1
parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=use_plog)

;put the requested grid into !dustem_grid
st0=mrdfits(input_table_name,0,h0)   ;This gets the main header
Ngrains=sxpar(h0,'NGRAIN')
model=strtrim(sxpar(h0,'MODEL'),2)
old_Ncomb=sxpar(h0,'NCOMB')
old_Nparams=sxpar(h0,'NPARAM')
old_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')

old_parameter_description=strarr(old_Nparams)
old_parameter_pmin=fltarr(old_Nparams)
old_parameter_pmax=fltarr(old_Nparams)
old_parameter_pnval=lonarr(old_Nparams)
old_parameter_plog=intarr(old_Nparams)
old_fparameter_description=strarr(old_NFparams)
old_fparameter_value=fltarr(old_NFparams)
FOR i=0L,old_Nparams-1 DO BEGIN
	old_parameter_description[i]=sxpar(h0,'PNAME'+strtrim(i+1,2))
	old_parameter_pmin[i]=sxpar(h0,'PMIN'+strtrim(i+1,2))
	old_parameter_pmax[i]=sxpar(h0,'PMAX'+strtrim(i+1,2))
	old_parameter_pnval[i]=sxpar(h0,'PNVAL'+strtrim(i+1,2))
	old_parameter_plog[i]=sxpar(h0,'PLOG'+strtrim(i+1,2))
ENDFOR
FOR i=0L,old_NFparams-1 DO BEGIN
	old_fparameter_description[i]=sxpar(h0,'FPNAME'+strtrim(i+1,2))
	old_fparameter_value[i]=sxpar(h0,'FPVALUE'+strtrim(i+1,2))
ENDFOR

new_parameter_description=[old_parameter_description,pd]
new_parameter_pmin=[old_parameter_pmin,iv_min]
new_parameter_pmax=[old_parameter_pmax,iv_max]
new_parameter_pnval=[old_parameter_pnval,iv_Nvalues]
new_parameter_plog=[old_parameter_plog,use_plog]

new_parameter_values=dustem_param_range2param_values(new_parameter_pmin,new_parameter_pmax,new_parameter_pnval,Nc=Nc,log=new_parameter_plog)

;=== read parameter values
old_param_values=mrdfits(input_table_name,1,h1)
old_Nparams=N_tags(old_param_values)
;=== read seds
old_seds=mrdfits(input_table_name,2,h2)
filters=tag_names(old_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)
old_I_tot=mrdfits(input_table_name,4,h4)
old_ext_tot=mrdfits(input_table_name,5,h5)
;Nwavs_ir=n_elements(wavs.wav)
;Nwavs_uv=n_tags(old_ext_tot[0])
Nwavs=n_elements(wavs.wav)

old_I_grains=replicate(old_I_tot[0],[Ngrains,old_Ncomb]) ;in MJy/sr for NH=1e20 H/cm2
old_ext_grains=replicate(old_ext_tot[0],[Ngrains,old_Ncomb])

iu=6    ;first extension with per-grain emission spectrum
FOR i=0L,Ngrains-1 DO BEGIN
	aa=mrdfits(input_table_name,iu,hg)    ;in MJy/sr for NH=1e20 H/cm2
	old_I_grains[i,*]=aa
	iu=iu+1
ENDFOR
FOR i=0L,Ngrains-1 DO BEGIN
	aa=mrdfits(input_table_name,iu,hg)
	old_ext_grains[i,*]=aa
	iu=iu+1
ENDFOR

New_Ncomb=old_Ncomb*iv_Nvalues[0]
New_Nparams=old_Nparams+1

new_seds=replicate(old_seds[0],New_Ncomb)
new_I_tot=replicate(old_I_tot[0],New_Ncomb)
new_ext_tot=replicate(old_ext_tot[0],New_Ncomb)

new_I_grains=replicate(old_I_tot[0],[Ngrains,new_Ncomb])
new_ext_grains=replicate(old_ext_tot[0],[Ngrains,new_Ncomb])

;stop

;=== define new quantities
;=== do the new parameter value structure
;format_template=old_param_values.P1
format_template='0.E0'
IF use_double THEN BEGIN
	format_template='0.D0'
ENDIF
str='one_param_st={'
FOR i=0L,new_Nparams-1 DO BEGIN
  str=str+'P'+strtrim(i+1,2)+':'+format_template
  IF i NE new_Nparams-1 THEN str=str+','
ENDFOR
str=str+'}'
toto=execute(str)
new_param_values=replicate(one_param_st,New_Ncomb)
;FOR i=0L,old_Nparams-1 DO BEGIN
;  new_param_values[*].(i)=old_param_values[].(i)
;ENDFOR
;=== fill in the new parameter values structure
ii=0L
FOR i=0L,old_Ncomb-1 DO BEGIN
	FOR k=0L,iv_Nvalues[0]-1 DO BEGIN
		FOR j=0L,old_Nparams-1 DO BEGIN
			new_param_values[ii].(j)=old_param_values[i].(j)
		ENDFOR
		new_param_values[ii].(new_Nparams-1)=*parameter_values[k]
		ii=ii+1
	ENDFOR
ENDFOR

;stop
new_seds=ptrarr(New_Ncomb)
em_spectra=ptrarr(New_Ncomb)
ext_spectra=ptrarr(New_Ncomb)

xrange=[1.,1.e3]
colors=long(range_gen(New_Ncomb,[0,255]))
loadct,13
input_spec=fltarr(Nwavs)   ;input spectrum in MJy/sr

CASE pd OF
	'(*!dustem_params).grains[0].MDUST_O_MH':BEGIN
		IF incsp NE 1 THEN BEGIN
			message,'Cannot compute new table from grid with no spectra included',/continue
			stop
		ENDIF
		igrain=0   ;This is the grain type corresponding to pd
		ii=0L
		;get the default parameter value
		dustem_init,model=model
		;======= compute the scaling factor
		;== This is the default model value for the parameter to be changed (ie MDUST_O_MH)
		default_value=(dustem_get_param_values_default(pd,param_current_values=param_current_values))[0]
		;fact=(*parameter_values[0])/default_value
		FOR i=0L,old_Ncomb-1 DO BEGIN
			FOR j=0L,iv_Nvalues[0]-1 DO BEGIN
				;param_fact=*parameter_values[j]
				fact=(*parameter_values[j])[0]/default_value    ;Caution, this assumes that the default value was not fixed to a different value when running the grid !
				FOR k=0L,Nwavs-1 DO BEGIN
					new_I_grains[igrain,ii].(k)=old_I_grains[igrain,i].(k)*fact
					new_I_tot[ii].(k)=new_I_tot[ii].(k)-old_I_grains[igrain,i].(k)+new_I_grains[igrain,ii].(k)
					new_ext_grains[igrain,ii].(k)=old_ext_grains[igrain,i].(k)*fact
					new_ext_tot[ii].(k)=new_ext_tot[ii].(k)-old_ext_grains[igrain,i].(k)+new_ext_grains[igrain,ii].(k)
					input_spec[k]=new_I_tot[ii].(k)
				ENDFOR
				;==== Could do that is fucking dustem_compute_sed did not need data
				;dustem_Ised=dustem_compute_sed([0.],input_spec=input_spec)
				;stop
				;=== run dustem once, just to get the st structure
				IF ii EQ 0 THEN BEGIN
					;stop
					;====recompute SED
					;use_next=1
					;Nfilt=1
					;ext=dustem_initialize_ext(use_next)
					sed=dustem_initialize_sed(Nfilters)
					sed.filter=filters
					sed.wave=dustem_filter2wav(filters)
					sed.stokesI=0.1
					sed.sigmaII=0.01
					dustem_set_data,m_fit=sed,m_show=sed
					dustem_init_params,model,[pd],[default_value],fpd=fpd,fiv=fiv,pol=use_polarisation,/force_reset
					st=0
					dustem_Ised=dustem_compute_sed([1.],st=st)   ;This is to compute st for the first time
					;stop
				ENDIF
				;fill in st with the spectrum from the grid
				;goto,jump_there
				sst=st.sed.(igrain+1)
				FOR k=0L,Nwavs-1 DO BEGIN
					aa=new_I_grains[0,ii].(k)
					sst[k]=aa
					;(st.sed.(igrain+1))[k]=aa    ;This does produce an error for some reason
				ENDFOR
				st.sed.(igrain+1)=sst   ;This is the new spectrum with modified PAH abundance
				;recompute total emission spectrum
				st.sed.em_tot=st.sed.(1)
				FOR kk=1L,Ngrains-1 DO BEGIN
					st.sed.em_tot=st.sed.em_tot+st.sed.(kk+1)
				ENDFOR
				;recompute total extinction spectrum
				;IDL> help,st.ext.abs_grain,/str
				;<Expression>    FLOAT     = Array[5, 800]
				st.ext.ext_tot=st.ext.abs_grain[1]+st.ext.sca_grain[1]
				FOR kk=1L,Ngrains-1 DO BEGIN
					st.ext.ext_tot=st.ext.ext_tot+st.ext.abs_grain[kk]+st.ext.sca_grain[kk]
				ENDFOR
				em_spectra[ii]=ptr_new(st.sed)
				ext_spectra[ii]=ptr_new(st.ext)
				;jump_there:
				;input_spec=st.sed.em_tot
				;dustem_Ised=dustem_compute_sed([0.],st=st)   ;finally recompute the SED
				dustem_Ised=dustem_compute_sed([0.],input_spec=input_spec)   ;finally recompute the SED
				;print,st.sed[100].em_tot,input_spec[100],dustem_Ised[10]
				new_sed=sed
				new_sed.stokesI=dustem_Ised
				new_seds[ii]=ptr_new(new_sed)
				IF keyword_set(show_seds) THEN BEGIN
					wavs=(*new_seds[ii]).wave
					Is=(*new_seds[ii]).stokesI
	    			IF ii EQ 0 THEN BEGIN
					    tit='Grid SEDs model '+model
					    xtit='Wavelength [mic]'
					    ytit='SED [MJy/sr for NH=1.e20 H/cm^2]'
					    ;stop
				    	cgplot,wavs,Is,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog,xrange=xrange,/xsty
			    	ENDIF ELSE BEGIN
			    		cgoplot,wavs,Is,psym=-4,color=colors[ii]
			    	ENDELSE
			    	ishow=100
			    	print,ii,i,j,fact
			    	print,new_I_grains[0,ii].W55,(st.sed.EM_GRAIN_1)[ishow],(st.sed.EM_tot)[ishow]
			    	stop
				ENDIF
				ii=ii+1
			ENDFOR
		ENDFOR
	END
ENDCASE

;stop

;===== write the output grid
use_out_filename='/tmp/dustem_grid_with_added_params.fits'
IF keyword_set(out_filename) THEN use_out_filename=out_filename

;parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=plog)
;seds=ptrarr(Nc)

dustem_write_sed_grid,model,new_parameter_description, $
					 new_parameter_pmin,new_parameter_pmax,new_parameter_pnval,new_parameter_plog, $
					 new_parameter_values,new_seds, $
					 em_spectra=em_spectra, $
					 ext_spectra=ext_spectra, $
					 fpd=old_fparameter_description, $
					 fiv=old_fparameter_value, $
					 filename=use_out_filename, $
					 use_polarisation=use_polarisation, $
					 double=use_double


;stop
the_end:

END