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 descriptions array
;    iv_min                 : minimum values for parameters in pd
;    iv_max                 : maximum values for parameters in pd
;    iv_Nvalues             : number of parameter values for parameters in pd
; OPTIONAL INPUT PARAMETERS:
;    plog                   : array indicating if the parameters are to be sampled in linear (0) or log scale (1)
;    out_filename           : output fits file name for the table (default ='/tmp/dustem_add_linear_params2grid_final_file.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:
;    Note: A fits file is created at each parameter addition
;    Note: grids do not retain scattering and absorption spectra separately, only totalm extinction. So these are not handled here.
; 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

Npd=n_elements(pd)

use_plog=lonarr(Npd)
IF keyword_set(plog) THEN use_plog=plog
;This is the combination values of parameters to be added
parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=use_plog)

temporary_table_name='/tmp/dustem_add_linear_params2grid_temporary_file.fits'
use_out_filename='/tmp/dustem_add_linear_params2grid_final_file.fits'
IF keyword_set(out_filename) THEN use_out_filename=out_filename

IF keyword_set(show_seds) THEN BEGIN
	xrange=[1.,1.e3]
	loadct,13
ENDIF

fact_MJy=dustem_get_emission_conversion_factor()
eachi=100L   ;will print progression at eachi

FOR ipd=0L,Npd-1 DO BEGIN
	;==== These are the parameter values of the parameter to be added
	this_parameter_values=dustem_param_range2param_values(iv_min[ipd],iv_max[ipd],iv_Nvalues[ipd],Nc=Nc,log=use_plog[ipd])
  IF ipd EQ 0 THEN BEGIN
  	use_input_table_name=input_table_name
  	use_output_table_name=temporary_table_name
  ENDIF ELSE BEGIN
  	IF ipd EQ Npd-1 THEN BEGIN
	  	use_input_table_name=temporary_table_name
	  	use_output_table_name=use_out_filename
	  ENDIF ELSE BEGIN
  		use_input_table_name=temporary_table_name
  		use_output_table_name=temporary_table_name 	
	  ENDELSE
  ENDELSE
	;===== read the requested grid fits file
	fst=file_info(use_input_table_name)
	IF fst.exists NE 1 THEN BEGIN
		message,'input table '+use_input_table_name+' not found',/continue
		stop
	ENDIF
	st0=mrdfits(use_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[ipd]]
	new_parameter_pmin=[old_parameter_pmin,iv_min[ipd]]
	new_parameter_pmax=[old_parameter_pmax,iv_max[ipd]]
	new_parameter_pnval=[old_parameter_pnval,iv_Nvalues[ipd]]
	new_parameter_plog=[old_parameter_plog,use_plog[ipd]]

  ;This is the new parameter values
  ;It is wrong and used only to save.
  ;This should not be used in that order to write the fits table (see below)
	;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(use_input_table_name,1,h1)
	old_Nparams=N_tags(old_param_values)
	;=== read seds
	old_seds=mrdfits(use_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(use_input_table_name,3,h3)
	old_I_tot=mrdfits(use_input_table_name,4,h4)   ;in MJy/sr
	old_ext_tot=mrdfits(use_input_table_name,5,h5)   ;in MJy/sr
	Nwavs=n_elements(wavs.wav)
  old_I_tot_vec=fltarr(old_Ncomb,Nwavs)
  old_ext_tot_vec=fltarr(old_Ncomb,Nwavs)
  ;Here we assume igrain=0, only to create the emission and extinction arrays
	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])

	old_I_grains_vec=fltarr(Ngrains,old_Ncomb,Nwavs) ;in MJy/sr for NH=1e20 H/cm2
	old_ext_grains_vec=fltarr(Ngrains,old_Ncomb,Nwavs)
	;old_seds_vec=fltarr(Nfilters,Nwavs)

	; 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(use_input_table_name,iu,hg)    ;in MJy/sr for NH=1e20 H/cm2
		old_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(use_input_table_name,iu,hg)
		old_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
  	old_I_tot_vec[*,k]=old_I_tot.(k)
  	old_ext_tot_vec[*,k]=old_ext_tot.(k)
  	old_I_grains_vec[*,*,k]=old_I_grains.(k)
 		old_ext_grains_vec[*,*,k]=old_ext_grains.(k)
		;old_seds_vec[*,k]=old_seds.(k)
  ENDFOR

	New_Ncomb=old_Ncomb*iv_Nvalues[ipd]     ;This is the new number of parameter values combinations
	New_Nparams=old_Nparams+1

  ;=== create new arrays for seds, total and partial emission and extinction
	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_seds_vec=fltarr(New_Ncomb,Nwavs)
	new_I_tot_vec=fltarr(New_Ncomb,Nwavs)
	new_ext_tot_vec=fltarr(New_Ncomb,Nwavs)

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

	new_I_grains_vec=fltarr(Ngrains,new_Ncomb,Nwavs)
	new_ext_grains_vec=fltarr(Ngrains,new_Ncomb,Nwavs)

	;=== define new quantities
	;=== do the new parameter value structure
	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)
	;=== fill in the new parameter values structure
	ii=0L
  message,'========= Old parameter values =========',/continue
	FOR i=0L,old_Ncomb-1 DO BEGIN
		print,old_param_values[i]
		FOR j=0L,iv_Nvalues[ipd]-1 DO BEGIN
			FOR k=0L,old_Nparams-1 DO BEGIN
				;IF ipd EQ 1 THEN stop
				new_param_values[ii].(k)=old_param_values[i].(k)   ;This fills in with the content of the fits table
			ENDFOR
			new_param_values[ii].(new_Nparams-1)=(*this_parameter_values[j])[0] ;this adds the new values for the new parameter
			ii=ii+1
		ENDFOR
	ENDFOR
  message,'========= New parameter values =========',/continue
  FOR kki=0L,New_Ncomb-1 DO print,new_param_values[kki]
	;stop

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

	input_spec=fltarr(Nwavs)   ;input spectrum in MJy/sr

	this_pd=pd[ipd]
	;===== searching for parameter descriptions like '(*!dustem_params).grains[0].MDUST_O_MH'
	pp=strpos(this_pd,'MDUST_O_MH')

  old_spec=fltarr(Nwavs)

	IF pp NE -1 THEN BEGIN
		;determine which grain type is adressed
		p1=strpos(this_pd,'[')
	  p2=strpos(this_pd,']')
	  igrain=long(strmid(this_pd,p1+1,p2-p1-1))
		IF incsp NE 1 THEN BEGIN
			message,'Cannot compute new table from grid with no spectra included',/continue
			stop
		ENDIF
		ii=0L   ;global counter
		;get the default parameter value for the model in use
		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(this_pd,param_current_values=param_current_values))[0]
		FOR i=0L,old_Ncomb-1 DO BEGIN
			IF i mod eachi EQ 0 THEN BEGIN
				message,'Doing initial line '+strtrim(i,2)+' of '+strtrim(old_Ncomb,2)+' ('+strtrim(1.*i/old_Ncomb*100,2)+' %)',/continue
			ENDIF
			FOR j=0L,iv_Nvalues[ipd]-1 DO BEGIN
				old_spec=reform(old_I_tot_vec[i,*])
				;stop
				fact=(*this_parameter_values[j])[0]/default_value    ;Caution, this assumes that the default value was not fixed to a different value when running the grid !
				new_I_tot_vec[ii,*]=old_I_tot_vec[i,*]
				new_I_grains_vec[*,ii,*]=old_I_grains_vec[*,i,*]
				new_I_grains_vec[igrain,ii,*]=reform(old_I_grains_vec[igrain,i,*])*fact
				new_I_tot_vec[ii,*]=new_I_tot_vec[ii,*]-old_I_grains_vec[igrain,i,*]+new_I_grains_vec[igrain,ii,*]
				new_ext_grains_vec[igrain,ii,*]=old_ext_grains_vec[igrain,i,*]*fact
				new_ext_tot_vec[ii,*]=new_ext_tot_vec[ii,*]-old_ext_grains_vec[igrain,i,*]+new_ext_grains_vec[igrain,ii,*]
				input_spec=reform(new_I_tot_vec[ii,*])
        ;stop

				;==== Could do that directly if 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.
				;=== Note: The st structure is useful to construct, even if not used in dustem_compute_sed, to get the right shape for the spectra
				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,[this_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 new emission and extinction spectrum
				FOR kk=0L,Ngrains-1 DO BEGIN
					aa=reform(new_I_grains_vec[kk,ii,*])
					st.sed[*].(kk+1)=aa/fact_MJy
					aa=reform(new_ext_grains_vec[kk,ii,*])
					st.ext[*].ext_tot=aa
				ENDFOR
				st.sed[*].em_tot=st.sed[*].(1) ;This is the emission of grain type 1
				FOR kk=1L,Ngrains-1 DO BEGIN
					st.sed[*].em_tot=st.sed[*].em_tot+st.sed[*].(kk+1)
				ENDFOR
				;=== This below is why it is useful to get st. Caution, st.sed must be in fortran units, not MJy/sr
				em_spectra[ii]=ptr_new(st.sed)
				ext_spectra[ii]=ptr_new(st.ext)
				;At that stage, st.sed.em_tot should equal input_spec, which is not the case
				diff_spec=input_spec-st.sed.em_tot*fact_MJy
				dustem_Ised=dustem_compute_sed([0.],st=st)   ;finally recompute the SED (should give the same result as below)
				;The line below gives the same result as the above, which is a confirmation that st and input_spec are consistent
				;dustem_Ised=dustem_compute_sed([0.],input_spec=input_spec)   ;finally recompute the SED
				new_sed=sed
				new_sed.stokesI=dustem_Ised
				new_seds[ii]=ptr_new(new_sed)
				IF keyword_set(show_seds) THEN BEGIN
				  colors=long(range_gen(New_Ncomb,[0,255]))
					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]'
				    	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,igrain,ii,ipd,i,j,fact,max(old_spec),max(st.sed.em_tot*fact_MJy),max(diff_spec)
			    	;print,new_param_values[ii]
				ENDIF
				ii=ii+1
			ENDFOR
		ENDFOR
	ENDIF ELSE BEGIN
		message,'parameter description '+this_pd,' not processed',/continue
		stop
	ENDELSE
	;===== write the output grid
	;help,model,new_parameter_description,new_parameter_pmin,new_parameter_pmax,new_parameter_pnval,new_parameter_plog,new_param_values,new_seds

  ;recompute the pointer to new_parameter_values from effective parameter values
  ;This is needed to have the grid saved actually be saved with the right parameter order ...
  message,'========= New parameter values written to table =========',/continue
  new_parameter_values=ptrarr(New_Ncomb)
  new_Nparams=old_Nparams+1
  vec=dblarr(new_Nparams)
  FOR ii=0L,New_Ncomb-1 DO BEGIN
  	FOR kk=0L,new_Nparams-1 DO BEGIN
  		vec[kk]=new_param_values[ii].(kk)
  	ENDFOR
  	new_parameter_values[ii]=ptr_new(vec)
  	;print,*new_parameter_values[ii]
  ENDFOR
 	;stop

	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_output_table_name, $
						 use_polarisation=use_polarisation, $
						 double=use_double
	 ;stop
ENDFOR

the_end:

END