User:Wwheaton/test9

Discretized NNLS Radon Transform for Earth Occultation Mapping

This approach is based on the Lawson Non-Negative Least Squares algorithm, nnls, now available in IDL, C, Fortran 90, and Mathlab, as well as the original F77. This algorithm proved spectacularly successful in stabilizing inverse imaging problems in gamma-ray and radio astronomy.

The basic idea is is simply to regard each datum, di in data bin i, as a single linear equation, i = 1,2,...,I:

di ≈ E[Σj (Mi,j Fj ) + Bi ]                                    (1)

where Fj is the flux from the j-th pixel in the sky, Mi,j is the matrix element connecting i and j, and E[-] is the statistical expectation operator. Note that F, B, and M are necessarily non-negative for physical reasons. The non-gamma-ray background B is determined for each i by numerous non-negative terms involving orbital and other parameters, which need to be estimated along with the Fj. Then in principle, we simply solve the system of I equations in J unknowns to find the best linear least squares approximation. Because the number of equations I is very large, the problem could be solvable by least squares.

Yet it is a good example of the kind of inverse problem that is classically considered ill-posed, and either singular or terribly badly conditioned, making such direct methods useless in practice. We will find, however, that the careful and systematic imposition of the non-negativity constraint, using the Lawson nnls algorithm, powerfully regularizes the problem, making it solvable in practice as well as in principle. There are a number of detailed problems and complications, which are treated below.

Unknowns
The number of unknowns J is large, proportional to the spatial resolution we can achieve. Based on the sharpness of the Earth limb (~0.2 degrees) this could in principle be on the order of 41253*((1/0.2)2) ~ 106. While this number of unknowns is very large, effective means of drastically reducing it, to the order of the number of detected sources, will be discussed below.

Detector response
The basic equation (1) neglects the fact that the detector response, D, depends on angles and energy. The form of the dependence is D(phi,theta,E,k), where phi and theta are the angles from the detector to the source in detector co-ordinates, E is photon energy, and k is one of the 16 pulse height channels, 0-15. All the D are non-negative, and may be taken to be in units of area per keV for each channel k. D can be estimated fairly accurately by detector simulation Monte Carlo programs. Furthermore, the (E,k) dependence is triangular, in the sense that a given E only has non-zero matrix elements for pulse-height channels k that correspond to lower energies. It is also almost certain that D can be factored, at least to some useful approximation, into angular and energy factors, so

M = T D = T A0 R(phi,theta) C(E,k),                                    (2)

where T is the Earth atmospheric transmission function, A0 is the geometric area of the detector, 0<R<1 is the angular response function, and C is the energy matrix, which is mostly dominated in its off-diagonal elements by Compton scattering, with an important contribution due to the broad energy resolution of the detectors, and some contribution due to pair-production.

The triangular character of the Compton energy matrix C, together with its non-negativity, means that the highest channel K can only be due to contributions from the highest energy, CK. Once the flux at the highest energy has been determined, the next lower flux can be found without reference to unknowns of lower energy. And so on, working downwards, until the entire spectrum is determined. This procedure is feasible because both the fluxes and the matrix elements are non-negative.

Transmission functions
The Earth transmission function values Ti,j are almost all 0.0 (~ 1/3) or 1.0 (~2/3), but a small fraction (<1%, at the limb of the Earth's atmosphere) will be between 0 and 1, and these will be important for the ultimate limit in source positioning accuracy. The problem of computing the air-mass x [gm/cmsq] through the atmosphere to a given position on the sky for the Earth's atmosphere and oblateness has been well solved for BATSE; this then gives the transmission as a function of energy E when combined with the gamma-ray absorption mu(E).

Atmospheric Compton scattering
This is essentially, logically, part of the detector response, considering the Earth as part of the collimating system of the instrument. It can be computed from first principles, and the computation can be checked by calibration against the amplitudes of slower pulsars like Vela X-1. For efficiency, the computation should probably done by constructing a large table, once and for all, which would allow us to look up the atmospheric Compton scattering matrix contribution to the source response. The table would be scattered flux (for a source of unit intensity) as a function of the elevation of the source (measured from the nadir, say), the scattering azimuth, the incident photon energy, and the scattered photon energy. It should be a fairly smooth function of angles and energy, and so need not be evaluated at an enormous number of points.

Reduction of number of unknowns to a smaller set of decoupled problems
Our primary interest here is in discovering point sources, and locating them well enough that they can be identified with optical or IR objects. Then the number of unknowns can be reduced to the number of sources, Js. To accomplish this, I propose starting with a coarser set of unknowns. Our 4X4 degree set of patches would give us ~2500 unknowns, which ought to be viable computationally, and give us a coarse resolution view of the whole sky. Then patches showing evidence of positive flux I would sub-divide into smaller pixel unknowns, while keeping the sparse (mostly high-galactic-latitude) regions with a coarse grid. For the ultimate in position resolution, I would zoom in on unknown sources one at a time, in order of significance, going all the way down to 0.2 degree pixels in the immediate region around the source. Once the source position has been accurately determined, we put that source into the catalog at that position, for all energies, and restore the coarse pixel for the sky background in that region, giving 2 unknowns. Continue in this way until all significant sky excesses have been resolved.

Background models
Essentially this is just the EBOP multi-component linear model we have, with some possible improvements that are pending.