; docformat = 'rst'

;+
;
; This is an example of fitting an astronomical spectrum using a Monte Carlo
; approach and should confirm that the AmesPAHdbIDLSuite has been correctly
; 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> mc_fit_a_spectrum
;
; :Author:
;   Dr. Christiaan Boersma
;
; :Copyright:
;   BSD licensed
;
; :History:
;   Changes::
;
;     04-15-2023
;     Remove redundant code and don't destroy observation. Christiaan Boersma.
;     05-03-2022
;     Add plotting size distribution. Christiaan Boersma.
;     04-28-2022
;     Adapt for piecewise errors. Christiaan Boersma.
;     04-27-2022
;     Use MCFit instead of Fit. Christiaan Boersma.
;     07-06-2021
;     Cleaned up progress bar. Christiaan Boersma.
;     04-30-2021
;     First version. Christiaan Boersma.
;-

;+
; Procedure demonstrating fitting an astronomical spectrum using the
; AmesPAHdbIDLSuite and establishing uncertainties using a Monte-Carlo
; approach.
;
; :Categories:
;   Example
;-
PRO MC_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

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

  ; search for the UIDs for a subset of PAHs
  uids = pahdb->Search("magnesium=0 oxygen=0 iron=0 silicium=0 chx=0 ch2=0 c>20 h>0")

  ; put back fullerenes, which have h=0
  fullerenes = [717, 720, 723, 735, 736, 737]

  uids = [uids, fullerenes]

  ; retrieve the transitions for the subset of PAHs
  transitions = pahdb->getTransitionsByUID(uids)

  ; apply full temperature cascade emission model with an 8 eV photon
  transitions->Cascade,8D*1.602D-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)

  ; clean up transitions
  OBJ_DESTROY,[transitions]

  ; fit the spectrum using Monte Carlo approach
  mcfit = spectrum->MCFit(observation, 1024)

  ; clean up
  OBJ_DESTROY,[spectrum]

  ; print results
  mcbd = mcfit->GetBreakdown()

  tags = TAG_NAMES(mcbd)

  ntags = N_TAGS(mcbd)

  PRINT,FORMAT='(A8,X,A6,X,A6)','','MEAN','STDEV'
  FOR i = 0, ntags - 1 DO $
    PRINT,FORMAT='(A8,X,F6.2,X,F6.3)',tags[i],mcbd.(i)[0],SQRT(mcbd.(i)[1])
  mcerr = mcfit->GetError()

  tags = TAG_NAMES(mcerr)

  ntags = N_TAGS(mcerr)

  FOR i = 0, ntags - 1 DO $
    PRINT,FORMAT='(A8,X,F6.2,X,F6.3)',tags[i],mcerr.(i)[0],mcerr.(i)[1]

  mcfit->Plot,/Wavelength

  key = ''
  IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."

  mcfit->Plot,/Wavelength,/Size

  IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."

  mcfit->Plot,/Wavelength,/Charge

  IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."

  mcfit->Plot,/Wavelength,/Composition

  IF !D.NAME EQ 'X' THEN READ,key,PROMPT="Press <enter> to continue..."

  mcfit->Plot,/DistributionSize,NBins=10L,Min=20,Max=200

  OBJ_DESTROY,[mcfit, pahdb, observation]

END