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

;JPB note: If fact=1 which happens, we could do nothing (to save time) and just copy that entry of the original table to the new one.
;I started coding this but gave up, as the old seds come out in a format different than needed.
;JPB note: We could also do without saving intermediate files.

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=10L   ;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])
  use_input_table_name=temporary_table_name
  use_output_table_name=temporary_table_name 	
  IF ipd EQ 0 THEN BEGIN
  	use_input_table_name=input_table_name
  ENDIF ;ELSE BEGIN
  IF ipd EQ Npd-1 THEN BEGIN
	  	use_output_table_name=use_out_filename
	ENDIF
	;===== read the requested grid fits file
	dustem_read_grid_table2arrays,use_input_table_name $
								  ,Ngrains,model,old_Ncomb,NFparams,INCSP,use_polar,use_double $
								  ,old_parameter_description,old_parameter_pmin,old_parameter_pmax,old_parameter_pnval,old_parameter_plog,old_fparameter_description,old_fparameter_value $
									,old_param_values,old_Nparams,old_seds,filters,Nfilters,wavs,Nwavs $
								  ,old_I_tot_vec,old_ext_tot_vec,old_I_grains_vec,old_ext_grains_vec
	
	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)
  filter_wavs=dustem_filter2wav(filters)
  sed_index=where(filters NE 'SPECTRUM' and filter_wavs gt 0.,count_sed)
  ;stop

	empty_sed=dustem_initialize_sed(Nfilters)
	empty_sed.filter=filters
	empty_sed.wave=filter_wavs

  ;==== remove temporary fits file
	fst=file_info(temporary_table_name)
	IF fst.exists EQ 1 THEN BEGIN
		message,'removing '+temporary_table_name,/continue
    str='rm '+temporary_table_name
    spawn,str
		;stop
	ENDIF

	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_vec=fltarr(Nwavs,New_Ncomb)
	new_I_tot_vec=fltarr(Nwavs,New_Ncomb)
	new_ext_tot_vec=fltarr(Nwavs,New_Ncomb)

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

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

	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')

	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
		st=dustem_make_empty_fortran_st(Ngrains)
		;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 OR i EQ old_Ncomb-1 THEN BEGIN
				message,'Doing initial line '+strtrim(i,2)+' of '+strtrim(old_Ncomb-1,2)+' ('+strtrim(1.*i/(old_Ncomb-1)*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 !
        ;message,'fact='+strtrim(fact,2),/continue
;stop
        ;IF abs(fact-1.) LE 1.e-5 THEN BEGIN
        ;	message,'Doing nothing ...',/continue
        ;	stop
        ;	new_seds[ii]=ptr_new(old_seds[i])
        ;	FOR kk=0L,Ngrains-1 DO BEGIN
				;		st.sed[*].(kk+1)=old_I_grains_vec[*,kk,ii]/fact_MJy
				;		;st.ext[*].ext_tot=old_ext_grains_vec[*,kk,ii]
				;	ENDFOR
				;	st.sed[*].em_tot=old_I_tot_vec[*,ii]
				;	st.ext[*].ext_tot=old_ext_tot_vec[*,ii]
        ;	em_spectra[ii]=ptr_new(st.sed)
        ;	ext_spectra[ii]=ptr_new(st.ext)
        ;	;stop
        ;	goto,nothing_to_do
        ;ENDIF
        ;transfer per-grain emission to new array
				new_I_grains_vec[*,*,ii]=old_I_grains_vec[*,*,i]
				;modifiy per-grain emission for grain type to be modified
				new_I_grains_vec[*,igrain,ii]=old_I_grains_vec[*,igrain,i]*fact
        ;modify total emission spectra accordingly
				new_I_tot_vec[*,ii]=old_I_tot_vec[*,i]-old_I_grains_vec[*,igrain,i]+new_I_grains_vec[*,igrain,ii]
        ;modifiy per-grain extinction for grain type to be modified
				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,*]
				new_ext_tot_vec[*,ii]=new_ext_tot_vec[*,ii]-old_ext_grains_vec[*,igrain,i]+new_ext_grains_vec[*,igrain,ii]
				input_spec=new_I_tot_vec[*,ii]
				;=== 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
				;	st=dustem_make_empty_fortran_st(Ngrains)
				;ENDIF
				;=== fill in st with the new emission and extinction spectrum
				;=== reform is needed here because of the structure aspect of st
				FOR kk=0L,Ngrains-1 DO BEGIN
					st.sed[*].(kk+1)=new_I_grains_vec[*,kk,ii]/fact_MJy
					;st.ext[*].ext_tot=new_ext_grains_vec[*,kk,ii]
				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
				st.ext[*].ext_tot=new_ext_tot_vec[*,ii]
				;=== 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,filters=filters)   ;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,filters=filters,sed_index=sed_index,/no_sort)   ;finally recompute the SED
				new_sed=empty_sed
				new_sed.stokesI=dustem_Ised
				new_sed.sigmaII=0.01*dustem_Ised
				new_seds[ii]=ptr_new(new_sed)

				nothing_to_do:

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

  message,'Writing '+use_output_table_name,/continue

	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