improved to save all needed values in fits tables

Jean-Philippe Bernard
1 parent fe88f96d
Showing 4 changed files with 54 additions and 19 deletions Show diff stats
src/idl/dustem_fit_intensity_example.pro
src/idl/dustem_init.pro
src/idl/dustem_read_fits_table.pro
src/idl/dustem_write_fits_table.pro
@@ -370,7 +370,9 @@ end
 ;I think modifications will have to be applied to dustemwrap_plot...
 ;idea: create a system variable that signals if the fits file has been successfully saved and use it in the dustemwrap_plot(obj/noobj) routine(s) 
 ;to add the '(From saved FITS file)' string to the title
-;
+
+;stop
+
 IF keyword_set(fits_save_and_restore) THEN BEGIN
    message,'Writing out results structure: '+fits_save_and_restore,/info
    dustem_write_fits_table,filename=fits_save_and_restore,help=help
@@ -392,7 +394,7 @@ ENDIF
  
 message,'Finished dustem_fit_intensity_example',/info
  
-stop
+;stop
  
 the_end:
  
@@ -349,7 +349,7 @@ defsysv,&#39;!dustem_show_plot&#39;,1           ;0=show no plot
 ;== This is for the default mode
 dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'}
 !dustem_inputs=ptr_new(dustem_inputs)
-
+use_model='DEFAULT'
  
 defsysv, '!dustem_HCD', ptr_new(1.00E20) ;I do not see the reason behind using a pointer in this case.
  
@@ -357,10 +357,6 @@ defsysv, &#39;!dustem_HCD&#39;, ptr_new(1.00E20) ;I do not see the reason behind using a
 ;===This is the Desert option
  
 IF keyword_set(model) THEN BEGIN
-
-defsysv, '!dustem_model', model
-
-
   CASE model OF
     'MC10':BEGIN
       (*!dustem_inputs).grain='GRAIN_MC10.DAT'
@@ -414,8 +410,11 @@ defsysv, &#39;!dustem_model&#39;, model
       stop
     END
   ENDCASE
+  use_model=model
 ENDIF
  
+defsysv, '!dustem_model', use_model
+
 ;Initialize !dustem_params structure
 IF not keyword_set(dir) THEN BEGIN
 ;  CASE getenv('DUSTEM_WHICH') OF
@@ -60,7 +60,20 @@ str_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
 ;==== add the predicted total spectrum in extension 3
 ;str_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)
  
+used_model=strtrim(sxpar(header_predicted_SED,'MODEL'),2)
+IF used_model NE 'DEFAULT' THEN BEGIN
+	dustem_init,model=used_model
+ENDIF ELSE BEGIN
+  	dustem_init,model=0
+ENDELSE
+!dustem_model=used_model
+
+;stop
+;==== This is to create (*!dustem_data).sed in case it does not exist already
 Nsed=n_elements(str_input_SED)
+sed=dustem_initialize_internal_sed(Nsed)
+(*!dustem_data).sed=ptr_new(sed)
+
 (*(*!dustem_data).sed).instru_names=dustem_filter2instru(str_input_SED.filter)
 (*(*!dustem_data).sed).filt_names=str_input_SED.filter
 (*(*!dustem_data).sed).wav=str_input_SED.wavelength
@@ -77,29 +90,34 @@ best_fit_chi2=sxpar(header_predicted_SED,&#39;CHI2&#39;)
 best_fit_rchi2=sxpar(header_predicted_SED,'RCHI2')
  
 ;===== Fill in header with parameter values
-parameters_desc=*(*!dustem_fit).PARAM_DESCS
-parameters_init_values=*(*!dustem_fit).PARAM_INIT_VALUES
-best_parameters=*(*!dustem_fit).CURRENT_PARAM_VALUES
-best_parameters_uncertainties=*(*!dustem_fit).CURRENT_PARAM_ERRORS
-best_fit_chi2=(*!dustem_fit).CHI2
-best_fit_rchi2=(*!dustem_fit).RCHI2
+;parameters_desc=*(*!dustem_fit).PARAM_DESCS
+;parameters_init_values=*(*!dustem_fit).PARAM_INIT_VALUES
+;best_parameters=*(*!dustem_fit).CURRENT_PARAM_VALUES
+;best_parameters_uncertainties=*(*!dustem_fit).CURRENT_PARAM_ERRORS
+;best_fit_chi2=(*!dustem_fit).CHI2
+;best_fit_rchi2=(*!dustem_fit).RCHI2
  
 Nparams=sxpar(header_predicted_SED,'NPARAMS')
 parameters_desc=strarr(Nparams)
 best_parameters=dblarr(Nparams)
 best_parameters_init_values=dblarr(Nparams)
 best_parameters_uncertainties=dblarr(Nparams)
+parameters_func_values=dblarr(Nparams)
 FOR i=0L,Nparams-1 DO BEGIN
 	parameters_desc[i]=sxpar(header_predicted_SED,'PARN'+strtrim(i+1,2))
 	best_parameters[i]=sxpar(header_predicted_SED,'PARV'+strtrim(i+1,2))
 	best_parameters_uncertainties[i]=sxpar(header_predicted_SED,'PARU'+strtrim(i+1,2))
 	best_parameters_init_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
+	parameters_func_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
 ENDFOR
  
-*(*!dustem_fit).PARAM_DESCS=parameters_desc
-*(*!dustem_fit).CURRENT_PARAM_VALUES=best_parameters
-*(*!dustem_fit).CURRENT_PARAM_ERRORS=best_parameters_uncertainties
-*(*!dustem_fit).PARAM_INIT_VALUES=best_parameters_init_values
+(*!dustem_fit).PARAM_DESCS=ptr_new(parameters_desc)
+(*!dustem_fit).PARAM_INIT_VALUES=ptr_new(best_parameters_init_values)
+(*!dustem_fit).CURRENT_PARAM_VALUES=ptr_new(best_parameters)
+(*!dustem_fit).CURRENT_PARAM_ERRORS=ptr_new(best_parameters_uncertainties)
+(*!dustem_fit).CHI2=best_fit_chi2
+(*!dustem_fit).RCHI2=best_fit_rchi2
+(*!dustem_fit).PARAM_FUNC=ptr_new(parameters_func_values)
  
 ;=== Below is only to get the same output form as dustem_compute_sed
 ;Note: If the same model is used, the outputs should be the same as what is in the fits file ....
@@ -47,6 +47,8 @@ IF not keyword_set(function_name) THEN BEGIN
 ENDIF ELSE BEGIN
 ENDELSE
  
+used_model=!dustem_model
+
 Ngrains=n_tags(dustem_spectra_st.sed)-2
  
 ;===== define structures containing results
@@ -89,6 +91,7 @@ best_parameters=*(*!dustem_fit).CURRENT_PARAM_VALUES
 best_parameters_uncertainties=*(*!dustem_fit).CURRENT_PARAM_ERRORS
 best_fit_chi2=(*!dustem_fit).CHI2
 best_fit_rchi2=(*!dustem_fit).RCHI2
+parameters_func_values=*(*!dustem_fit).PARAM_FUNC
  
 Nparams=n_elements(parameters_desc)
  
@@ -169,6 +172,12 @@ FOR i=0L,Nparams-1 DO BEGIN
 	sxaddpar,header_predicted_SED,'PARN'+strtrim(i+1,2),parameters_desc[i],' DustemWrap parameter name'
 ENDFOR
 sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARN'+strtrim(i,2)
+
+;====== write dustem wrap parameter best fit values
+sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap dust model used ======',after='COMMENT'
+sxaddpar,header_predicted_SED,'MODEL',used_model,'Dust model used'
+sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='MODEL'
+
 ;====== write dustem wrap parameter best fit values
 sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter best fit values ======',after='COMMENT'
 FOR i=0L,Nparams-1 DO BEGIN
@@ -181,14 +190,21 @@ sxaddpar,header_predicted_SED,&#39;COMMENT&#39;,&#39;======= Dustemwrap best fit parameter u
 FOR i=0L,Nparams-1 DO BEGIN
 	sxaddpar,header_predicted_SED,'PARU'+strtrim(i+1,2),best_parameters_uncertainties[i],'Param. '+parameters_desc[i]
 ENDFOR
-sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARIU'+strtrim(i,2)
+sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARU'+strtrim(i,2)
  
 ;====== write dustem wrap parameter initial fit values
 sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter initial fit values ======',after='COMMENT'
 FOR i=0L,Nparams-1 DO BEGIN
 	sxaddpar,header_predicted_SED,'PARI'+strtrim(i+1,2),parameters_init_values[i],'Param. '+parameters_desc[i]
 ENDFOR
-sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARIV'+strtrim(i,2)
+sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARI'+strtrim(i,2)
+
+;====== write dustem wrap parameter func values (not sure exactly what this is)
+sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter func values ======',after='COMMENT'
+FOR i=0L,Nparams-1 DO BEGIN
+	sxaddpar,header_predicted_SED,'PARF'+strtrim(i+1,2),parameters_func_values[i],'Param. '+parameters_desc[i]
+ENDFOR
+sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARF'+strtrim(i,2)
  
 ;===== complement extension header for predicted spectrum total
 sxaddpar,header_predicted_spectrum_tot,'TITLE','OUTPUT BEST FIT TOTAL SPECTRUM BY DUSTEMWRAP'