---

Spectroscopic database fitting

The first step in fitting astronomical observations is loading astronomical observations.

Dealing with astronomical observations

Astronomical observations can be handled by the 'AmesPAHdbIDLSuite_Observation'-object, which is able to read text, ISO-SWS, and Spitzer-IRS-files. A convenience routine is available to manage units; AmesPAHdbIDLSuite_CREATE_OBSERVATION_UNITS_S.

observation = OBJ_NEW('AmesPAHdbIDLSuite_Observation', $
                     'myObservationFile', $
              Units=AmesPAHdbIDLSuite_CREATE_OBSERVATION_UNITS_S())

In addition, an 'AmesPAHdbIDLSuite_Observation'-object can be created using keyword-initializers.

observation = OBJ_NEW('AmesPAHdbIDLSuite_Observation', X=frequency, $
                                                       Y=intensity, $
                                                       ErrY=ystdev, $
              Units=AmesPAHdbIDLSuite_CREATE_OBSERVATION_UNITS_S())

The 'AmesPAHdbIDLSuite_Observation'-object exposes the observation and provides the 'Plot', and 'Write'-methods for output. The 'Plot'-method will display the observation and accepts the 'Oplot', and 'Color'-keywords to control overplotting and color, respectively. Through IDL's keyword inheritance mechanism additional keywords accepted by IDL's 'PLOT'-procedure can be passed.

observation->Plot,XRANGE=[2.5,15],/XSTYLE

The 'Write'-method will write the observation to a single text (.txt) file. Optionally, a prefix can be given that will be prepended to the filename.

observation->Write,'myPrefix'

NB Text-files can have up to five columns, organized as follows. Column 1: abscissa, Column 2: ordinate, Column 3: continuum, Column 4: uncertainty in ordinate, and Column 5: uncertainty in abscissa.

The 'AmesPAHdbIDLSuite_Observation'-object's 'Rebin', 'AbscissaUnitsTo', 'SetGridRange', and 'GetGrid'-methods can rebin the observation onto a specified grid or, with the 'Uniform'-Keyword set, onto a uniform created grid with specified sampling; convert the units associated with the abscissa; change the grid range; and return the grid, respectively.

observation->Rebin,myGrid

observation->Rebin,5D,/Uniform

observation->AbsciccaUnitsTo

observation->SetGridRange,500,2000

grid = observation->GetGrid()

Database fitting

Spectroscopic database fitting is handled by the 'AmesPAHdbIDLSuite_Spectrum'-object and it either accepts an 'AmesPAHdbIDLSuite_Observation'-object, or a simple array of ordinates with an optional array of ordinate uncertainties. Whether ordinate uncertainties are provided or not, the 'AmesPAHdbIDLSuite_Spectrum'-object's 'Fit'-method will perform a non-negative least-chi-square or non-negative least-square fit, respectively, and return an 'AmesPAHdbIDLSuite_Fitted_Spectrum'- object.

fit = spectrum->Fit(observation)

Optionally, a callback-procedure can be specified with the 'CALLBACK_NNLS'-keyword that is called after each iteration and can be used to, for example, monitor progress.

PRO CB,A,b

  COMMON COM, observation

  observation->Plot

  dim = SIZE(A, /DIMENSIONS)

  yfit = TOTAL(REBIN(b, dim[0], dim[1]) * b, 1)

  OPLOT,observation->GetGrid(),yfit,COLOR=2

END

fit = spectrum->Fit(intensity, uncertainty, CALLBACK_NNLS='CB')

In addition, optionally the 'TOLERANCE_NNLS'- and 'MAXITER_NNLS'-keywords can be used to set the tolerance for convergence and the maximum number of iterations, respectively.

fit = spectrum->Fit(observation, TOLERANCE_NNLS=1D-20, MAXITER_NNLS=256L)

The 'AmesPAHdbIDLSuite_Fitted_Spectrum'-object exposes the fit and provides the 'Plot', and 'Write'-methods for output. The 'Plot'-method accepts the 'Residual', 'Size', 'Charge', 'Composition' and 'DistributionSize-'keywords, which selectively display the residual of the fit, or either the size, charge, compositional breakdown, the size distribution. Without these keywords the fit itself is displayed.

fit->Plot,/Charge

Optionally, the 'Wavelength', 'Stick', 'Oplot', 'Legend', and 'Color'-keywords can be given to the 'Plot'-method to control the abscissa, stick representation, overplotting, legend and color, respectively. Through IDL's keyword inheritance mechanism additional keywords accepted by IDL's 'PLOT'-procedure can be passed. Furthermore, the 'Small'- and 'Medium-keywords can be used to set the small and medium cutoff size, defaulting to 50 and 70, respectively, when plotting the size breakdown.

fit->Plot,/Size,Small=20L,Medium=50L,/Wavelength,XTITLE=[2.5,15],/XSTYLE

The 'AmesPAHdbIDLSuite_Fitted_Spectrum'-object's 'Write'-method will write the fit to a single text (.txt) file. Optionally, a prefix can be given that will be prepended to the filename.

fit->Write,'myPrefix'

The 'AmesPAHdbIDLSuite_Fitted_Spectrum'-object's 'GetClasses', and 'GetBreakdown'-methods return the fit broken down by charge, size, and composition, where the first provides the spectrum for each component and the latter its relative contribution.

classes = fit->GetClasses()

breakdown = fit->GetBreakdown()

Optionally the 'Small'- and 'Medium'- keywords can be set, defaulting to 50 and 90, respectively, which control the small and medium cutoff size in number of carbon atoms.

classes = fit->GetClasses(Small=20L, Medium=50L)

The 'GetBreakdown'-method also accepts the 'Flux'-keyword, which controls whether the relative breakdown should be reported based on fitted weight or integrated flux.

breakdown = fit->GetBreakdown(/Flux)

The quality of a fit can be assessed using the 'AmesPAHdbIDLSuite_Fitted_Spectrum'-object's 'GetError'-method, which returns a struct with the integral of the absolute residual over the integral of the observations for the entire spectrum and piecewise for the different PAH bands.

error =fit->GetError()

Monte Carlo

Neither non-negative least-chi-square and non-negative least-square fits provide uncertainties on their own. However, in combination with a Monte Carlo approach that varies the observed spectrum within its uncertainties they can.

Monte Carlo fitting is handled by the 'AmesPAHdbIDLSuite_Spectrum'-object and it either accepts an 'AmesPAHdbIDLSuite_Observation'-object, or a simple array of ordinates with a required array of ordinate uncertainties, and the number of samples. The 'AmesPAHdbIDLSuite_Spectrum'-object's 'MCFit'-method will perform a non-negative least-chi-square fit and return an 'AmesPAHdbIDLSuite_MCFitted_Spectrum'-object.

mcfit = spectrum->MCFit(observation, 5000L)

By default a normal distribution is used to draw from. However, the 'Uniform'-keyword can be set to draw from a uniform distribution instead.

mcfit = spectrum->MCFit(observation, 5000L, /Uniform)

The 'AmesPAHdbIDLSuite_MCFitted_Spectrum'-object offers mostly the same set of methods as the regular 'AmesPAHdbIDLSuite_Fitted_Spectrum'-object does, but transparently handles the necessary Monte Carlo statistics. For example, the 'AmesPAHdbIDLSuite_MCFitted_Spectrum'-object's 'GetClasses' returns spectra of the fit broken down by charge, size, and composition. The statistics are as returned by IDL's 'MOMENT'-function, i.e., a four elements containing the mean, variance, skewness, and kurtosis.

mcclasses = mcfit->GetClasses()