; docformat = 'rst'
;+
;
; This is an example of testing the stability of fitting an
; astronomical spectrum by iteratively removing the most contributing
; PAH species from a subsequent fit, built around the functionality
; provided by the AmesPAHdbIDLSuite and should help confirm that the
; it has been properly installed. The source code is annotated to
; guide users and developers in the inner workings of the suite.
;
; Updated versions of the NASA Ames PAH IR Spectroscopic Database and
; more information can be found at: `www.astrochemistry.org/pahdb <https://www.astrochemistry.org/pahdb>`.
;
; :Examples:
; Call the procedure directly::
;
; IDL> stability_of_fit_a_spectrum
;
; :Author:
; Dr. Christiaan Boersma
;
; :Copyright:
; BSD licensed
;
; :History:
; Changes::
;
; 08-19-2019
; Documentation added. Christiaan Boersma.
;-
;+
; Procedure testing the stability of a spectroscopic fit.
;
; :Categories:
; Example
;-
PRO STABILITY_OF_FIT_A_SPECTRUM
COMPILE_OPT IDL2
; the Spitzer IRS/SL 10 - 15 micron spectrum of NGC7023
file = 'ngc7023.dat'
; read observations into AmesPAHdbIDLSuite_Observation
observation = OBJ_NEW('AmesPAHdbIDLSuite_Observation', $
file, $
Units=AmesPAHdbIDLSuite_CREATE_OBSERVATION_UNITS_S(AUNIT=3, OUNIT=1))
; turn wavelength into frequency
observation->AbscissaUnitsTo,1
; rebin to uniform frequency grid
observation->Rebin,5D,/Uniform
; read in the default database defined by the environement variable
; !AMESPAHDEFAULTDB or the system variable AMESPAHDEFAULTDB. use the
; keyword FILENAME if these have not been set
pahdb = OBJ_NEW('AmesPAHdbIDLSuite')
; retrieve all the transitions from the database
transitions = pahdb->getTransitionsByUID( -1 )
; have every PAH absorb 6 eV (CGS units) and include the temperature
; cascade
transitions->Cascade,6D * 1.6021765D-12
; shift data 15 wavenumber to the red
transitions->Shift,-15D
; convolve the transitions into a spectrum
spectrum = transitions->Convolve( $
Grid=observation->getGrid(), $
FWHM=15D, $
/Gaussian)
OBJ_DESTROY,transitions
; run stability test
niterations = 11
data = REPLICATE({norm: 0D, $
anion: 0D, $
neutral: 0D, $
cation:0D, $
small:0D, $
large:0D, $
pure:0D, $
nitrogen: 0D,$
solo: 0L, $
duo: 0L, $
trio: 0L, $
quartet: 0L, $
quintet: 0L}, niterations)
FOR i = 0, niterations - 1 DO BEGIN
; fit the spectrum
fit = spectrum->Fit(observation)
; store breakdown
d = data[i]
STRUCT_ASSIGN,fit->getBreakdown(), d
data[i] = d
data[i].norm = fit->getNorm()
; sort the spectra based on their total contribution to the flux
fit->Sort,/Flux
; remove most contributing species from spectrum
spectrum->Intersect,(fit->getUids())[1:*]
; clean up fit
OBJ_DESTROY,fit
ENDFOR
PRINT
; clean up
OBJ_DESTROY,[spectrum, pahdb]
; plot ionized fraction
PLOT,INDGEN(niterations), data.cation / (data.cation + data.neutral),XTITLE='iteration [#] / species removed [#]',YTITLE='ionized fraction [ratio]',YRANGE=[0,1],/YSTYLE
key = ''
IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."
; plot large fraction
PLOT,INDGEN(niterations), data.large / (data.large + data.small),XTITLE='iteration [#] / species removed [#]',YTITLE='large fraction [ratio]',YRANGE=[0,1],/YSTYLE
IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."
; plot nitrogen fraction
PLOT,INDGEN(niterations), data.nitrogen / (data.nitrogen + data.pure),XTITLE='iteration [#] / species removed [#]',YTITLE='nitrogen fraction [ratio]',YRANGE=[0,1],/YSTYLE
IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."
; plot solo/(duo+trio)
PLOT,INDGEN(niterations), data.solo / (data.duo + data.trio),XTITLE='iteration [#] / species removed [#]',YTITLE='solo/duo+trio [ratio]',YRANGE=[0,1],/YSTYLE
END