- Almost all my fitting programs fit a generic function using the
`CURFIT`routine from Bevington's book "*Data reduction and Error Analysis for the Physical Sciences*" (McGraw-Hill 1969)

The few exceptions are simple regression programs using a plain least square routine, some polynomial fits using Bevington's`POLFIT`and some surface fitting using a routine by S.McKechnie of WKS Leiden. `CURFIT`is part of the class of "generalized least square" or Levenberg-Marquardt methods, as described e.g. in the Numerical Recipes book, and is widely used by astronomers (try a Google search !).- The typical
`CURFIT`based program must provide a Fortran function`FUNC(X,I,A,NT)`where`Y=FUNC(X,I,A,NT)`evaluates y=f(x_{i}) where f(x) is the function to be fitted and the`NT`elements of array`A`are the parameters to be fitted. - What matters are the parameters to be fitted. If you want to parametrize
the function on other parameters but not to fit them, one handy way is to
pass them through a
`COMMON`block, or write a wrapper. A wrapper is also needed to hide a routine which has a different calling interface. - On the other hand it does not matter what are the
`X`, although often and in many places of the code`X`is the independent variable of the data to be fitted.One could call

`FUNC(W,1,A,NT)`to compute y=f(w) at a single scalar point x=w.One can also use

`X`as an arbitrary array of indices. This is customary in fitting X-ray spectra, where one is not actually fitting versus energy, but versus channel number. Here usually`X(I)=I`or`X(i)=I+K`.One is not even bound to fit unidimensional data. I fitted two-dimensional functions z=f(x,y) to image data using as

`X`a one-dimensional index built from the pixel addresses`I,J`### basic CURFIT program

- The basic CURFIT program assumes that
`FUNC`has such name and calls it in the finite difference derivative routine, and in the chi-square computation routine (while`FUNC`is passed as an`EXTERNAL`argument to`CURFIT`itself.If desired one can replace the derivative routine with a specific routine computing the derivative analytically (for speed), or with some other optimized variant (see below).

If desired one could propagate the named

`EXTERNAL`argument to the other routines as well (I've never done this). - The
`CURFIT`routine does not do any fit ! It simply performs a single fit iteration which leads from an*initial guess*of the parameters to a better one which has a smaller chi-square.

(Actually there is a loop, potentially an infinite loop, inside`CURFIT`in case the adjustment does not lead to a smaller chi-square, changing a tuning variable`FLAMBDA`) - Therefore the typical, plain fitting program must provide something like
the following:
- code reading
`NPTS`data points`X,Y`and eventual errors`SY` - code asking (or computing) the initial guess for the
`NT`parameters to be fitted`A` - initializing the finite difference increments
`DA`used for the computation of the partial derivatives of FUNC - initializing
`FLAMBDA`and`MODE` - DO loop on maximum number of allowed iterations
`CALL CURFIT(X,Y,SY,NPTS,NT,MODE,A,DA,SA,FLAMBDA,YFIT,CHI)`- eventually displaying partial results
- leaving the loop if a convergence criterion is achieved
- END DO

About

`MODE`and`FLAMBDA`see Bevington's code.

The convergence criterion is care of the user. Generally one checks that the difference between the chi-square`CHI`of the current iteration and the one (saved) of the previous one is below a predefined threshold. But one could also check on the parameter values themselves to converge.

On exit`A`contains the fitted parameters,`SA`a rough estimate of their errors, and`YFIT`the values of fitted function y=f(x). - code reading
- I tend not to trust errors computed by
`CURFIT`this way, but prefer to use confidence limits computed according to standard prescriptions ( Lampton, Margon & Bowyer 1976 or Avni 1976)selecting 1, 2 or more

*interesting*parameters*stepping*them over an equispaced grid

for each step do a fit on all__remaining__parameters, and save the chisquare

look for the points where the chi square differs from the minimum chi square (the best fit) by less than a given delta chi square (which depends on the confidence level and number of interesting parameters) - Several of my ad-hoc programs, all the times I need to
fit a predefined model and am just interested in the parameters, are written
as the basic programs described above, i.e. use an ad-hoc
`FUNC`routine, and call`CURFIT`within an iteration loop. - A typical variant may include the selection of a range of data to be fitted,
narrower than the original data array.

Another variant may use an ad-hoc optimized version of the`FDER`or`FDERIV`derivative routine.### Spectral fitting

- By
*spectrum*I mean any distribution of spectral flux vs an energy variable (actually this can be a wavelength, a frequency or an energy). In practice my library of spectral models supports only three "standard" representations :*Flambda vs lambda*: erg/cm^{2}/s/Å vs Å*dN/dE vs E*: photon/cm^{2}/s/keV vs keV*Fnu vs nu*: erg/cm^{2}/s/Hz vs Hz

The first is used by the UV fitting programs, the second by the X-ray fitting programs and the third is currently not used.

- The library (see report
**R7**, entry 85 for details) defines a number of spectral functions allowing to select, setting appropriate`COMMON`flags via dedicated calls :- the standard unit representation (one of the three above)
- the unit for temperature parameters (when applicable) : keV or Kelvin
- whether the user supplied normalization is comprehensive of
any physical constants (c, h, pi) or not. Current programs assume
that physical constants are
*included*and display in their log file the source-specific part of the normalization when applicable (e.g. the emission measure or the part proportional to r²/d² for a Bremsstrahlung or blackbody)

The conversions (when applicable) are either taken care by branches in the code of each individual spectral form routine, or, when the model is natively computed in dN/dE, by a conversion routine called

`MOAMED`(according to the saying "*if the mountain will not come to Mohammed, Mohammed will go to the mountain*")As a technicality, all spectral routines are defined with 2 (or 3) entry points as follows :

- a numbered function entry point
`Y=SPE`*nn*(X,A) - an optional alias entry point
`Y=`*alias*(X,A) - an info entry point
`CALL INFO`*nn*(M,NAME,PARA)

where the former ones compute the flux

`Y`at energy/lambda/frequency X, while the latter information call*returns*the number of parameters`M`. Array`A(M)`contains the user supplied values of the parameters, character array`PARA(M)`contains their names (and`NAME`the spectral form name). - One can see that the above function has
*not*the calling sequence expected for`FUNC`. In fact this is not the function to be fitted, and there are several reasons for this :- the above calling sequence is handier for non-fitting usages of the spectral model
- even with single component models, the actual spectrum is usually the
product of an
*unsabsorbed*model by (interstellar) absorption, be it UV*reddening*or X-ray*photoelectric absorption* - in most cases one desires a multiple component model (the present
version of the software supports only
*additive components*, i.e. the sum of several models, all absorbed through a common absorption (a former version supported separate by-model absorptions) ) - finally not all parameters of a model are necessarily to be fitted
- and last but not least one is not necessarily fitting a model (for
X-rays one is fitting the
*convolution*of a model with a response matrix)

- The library provides a set of
*single model*wrapper routines (not used for fitting) namely`SPECUV`and`SPECSX`(which handle the UV reddening or the X-ray absorption, and call a common dispatcher`SPECS`which uses a computed GOTO to call the desired spectral model; and a set of*multiple model*wrapper routines (namely`SPEUV2`and`SPEX2`which handle absorption and call a common dispatcher`SPESUM`). - The dispatcher
`SPESUM`takes care of calling each model with its array of parameters, but first builds the array of parameters mixing in the correct way (using flags) the parameters fixed by the user with those to be fitted.

The calling wrappers act in a similar way for the absorption parameter (A_{V}or N_{H}) which may or may not be fitted. - The UV form of
`FUNC`offers the choice of fitting plain flux vs wavelength (i.e. uses "raw"`SPEUV2`at x(i)) or of fitting the integral of the model over the data bin (by numeric integration of`SPEUV2`on 20 steps over the bin). In the latter the integral is computed once for all bins, and the stored value returned at later calls.Similarly the derivative

`FDER`routines also computes derivatives once at each parameter change, and returns the stored value at other calls. - the X-ray form of
`FUNC`also computes the model once at parameter change only, using the stored values at other calls. However it performs an important additional task, which is the*convolution*of the spectrum with the response matrix for the`I`-th channel.The

`FDER`routine in the`melib`library is also used, which computes derivatives only at fitted parameter change. - One technical adjustment done quite recently to the X-ray fitting program
to speed them up a bit, is that, instead of reading the entire response
matrix (which in XAS is structured as an image), only the values pre-analysed
to be above a given threshold are saved internally and used for convolution.
- A further, more recent addition is the
*native support*to*cross-normalizations*when fitting data from more than one instrument. In XAS (unofficial add-on), inspired in this by my Exosat software, spectra (and matrices) from different instruments are physically joined in new files.

In the case of multi-instrument files, one can fit the cross-normalization of all instruments (except one which is kept as reference). This is done saving for each instrument :- a normalization value
- a normalization flag (0=not fitted 1=fitted 2=stepped)
- a normalization index (N means that it is fitted parameter A(N))
- a couple of arrays with the start and end channel of each instrument

The

`FUNC`code takes care of this and recovers the normalizations of individual spectral bins from the relevant fitted value (otherwise they are initialized to the user-supplied fixed cross-normalization) and divides the computed count spectrum by them. - This is similar to but separate from what is done to decide whether a parameter
is fitted, stepped or fixed. There are the following arrays :
- a matrix (indexed by component and parameter) of values, which record the user supplied value of the j-th parameter of the i-th component
- a matrix of flags which tell whether the j-th parameter of the i-th component is fixed (0), fitted (1) or stepped (2)
- a matrix of indices (N means that the j-th parameter of the i-th component is fitted parameter A(N))

There is a further index array for the stepped case. It has two couple of elements (for the first and second stepped parameter). The p-th location contains I,J to indicate that the p-th stepped parameter is the j-th parameter of the i-th component.

The absorption parameter (A

_{V}or N_{H}) is handled similarly, being flagged as -1,-1 in the index array, while it does not use the matrix of values, flags and indices but a separate scalar value, flag and index.