README_oclc_epq03

These routines implement the model of Elliot, Person and Qu 2003.
Analysis of stellar occultation data for planetary atmospheres. II - 
Inversion, with application to Pluto and Triton.  Ap J 126: 1041-1079.

The routines are called oclc_epq03_*.  oclc stands for OCcultation
LightCurve.  EPQ03 refers to the paper cited above.
Other occultation lightcurves can be implemented, and will be prefixed
oclc_<other>_.pro.

---------------------------------------------------------------------------
OVERVIEW
---------------------------------------------------------------------------

This formulation is a numerical implementation of the Abel transform,
integrating the bending angle along the line of sight through the
atmosphere. 

It is a challenge to keep track of the shells on which the variables are
defines, since some are tabulated at shell centers and others at shell
boundaries.

We begin by reproducing Figure 2 and caption from EPQ03:

======================================================================================================================
          .
 Boundary {############################## r ^#######                     y ^
 Integral {#################################|#######                       |
   Region {#################################|#######                       |
           _ r_b = r_{i_b-1/2}  _ _ _ _ _ _ ________________...............|
          {                    {            |               \_______________ y_b = y_{i_b-1/2} _ _
          {                    {            |                              |                      }
          {       Delta r_i_b  {      r_i_b +                              |                      }
          {                    {            |                              |                      }
Inversion {                    {            |                              + y_i_b                } Delta y_i_b
Summation {  r_b = r_{i_b+1/2}  _ _ _ _ _ _ ___________....................|                      }
   Region {                    {            |          \____theta_{i_b+1/2}|                      }
          {                    {            |                    \__________ y_{i_b+1/2} _ _ _ _ _
          {    Delta r_{i_b+1} {  r_i_{b+1} +                              |                      }
          {                    {            |                              |                      }
          {                    {            |                              |                      }
          {  r_b = r_{i_b+3/2}  _ _ _ _ _ _ _______........................+ y_{i_b+1}            } Delta y_{i_b+3/2}
          {                                 |      \_____   theta_{i_b+3/2}|                      }
          {                                 |            \_____            |                      }
          {                                 |                  \_____      |                      }
           _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _|_ _ _ _ _ _ _ _ _ _ _ _ \______ y_i_{b+3/2} _ _ _ _ _
                                            |                              |
                                           //                             //
                                        Body Plane                Observer Plane
                                           //                             //
                                           |                               |
                                     r = 0 +                         y = 0 +
                                     
                                           <--------------- D ------------->
                                           
EPQ03 Figure 2. Shells used in inversion summations.  This magnified
view of the body section of Fig 1 shows our scheme for the inversion
shells.  Note that the integral shell indeices, i, occur at shell
midpoints, while the index of i-1/2 denotes the upper boundary of the
shell i and i+1/2 denotes its lowerboundary.  Corresponding values on
the y-axis are determined by the ray paths for the refracted
starlight, and the shaded area above is the inersion boundary integral
region. The first shell of the inversion region ("boundary shell")
begins at radius r=r_b and has an index i = i_b

======================================================================================================================

To maintain the strictest 1-to-1 correspondence with the cited paper, we
will use the shell index scheme so defined.  

The scheme tabulates some values on shell centers and others on shell boundarys.

shell center (integer i): y, Delta y, phi, sigma(phi), Delta r, Delta theta
shell edge (half-integer i): r, theta, nu, n, p, T, H

---------------------------------------------------------------------------
VARIABLES
---------------------------------------------------------------------------

dist       
EPQ03 notation: D 
Distance between target and observer
scalar
cm

dtheta_dr     
EPQ03 notation: d theta/dr 
Derivative of the bending angle with respect to planet radius,
positive for an isothermal atmosphere.
scalar
radian/cm

deltatheta
EPQ03 notation: Delta theta or Delta theta_j
change in theta over one shell
negative for an isothermal atmosphere.

delta_y
EPQ03 notation: Delta y_i
width of the shadow radius shell, defined as y_{i+1/2} - y_{i-1/2}
so delta_y is negative for an isothermal atmosphere.

i
EPQ03 notation: i
index of shells ranging from 1 (topmost) to imax

ib
EPQ03 notation: i_b
index of first (highest) shell for the inversion.
1 <= i < ib      shells of boundary integral region
ib <= i <= imax  shells of inversion
scalar
unitless

iminus
EPQ03 notation: i - 1/2
index of top of shell i
unitless

imax
largest shell index (closest to surface)
scalar
unitless

j
EPQ03 notation: j
summation shell index for calculating r_i, theta_i, nu_i, p_i, 
ranging from ib to i.
scalar
unitless  

nu
Refractivity. 
Unitless

phi
EPQ03 notation: phi
normalized stellar flux.
unitless.

phiarr
EPQ03 notation: phi_i
normalized stellar flux array, indexed from 1 to imax
unitless.

r
EPQ03 notation: r
planet radius (distance from planet center).
cm

rarr
EPQ03 notation: r_{i + 1/2}
planet radius array, indexed form i=1/2 to i = imax + 1/2
cm

rb

y
shadow radius (distance from shadow center), negative for
far-limb contribution.
[same as rho in oclc_ey92]
scalar.
cm.

ymid
EPQ03 notation: y_i (eg eq. 34)
shadow radius array tabulted at shell mid, 
indexed from 0 to imax-1,
cm

y_i
EPQ03 notation: y_i
shadow radius at shell center

y_iminus
EPQ03 notation: y_{i - 1/2}
shadow radius at shell top

yval
EPQ03 notation: y_i
shadow radius array, indexed from 1 to imax
cm

yb
shadow radius at beginning of inversion
equivalent to y_{i_b - 1/2}

theta
bending angle, negative for an isothermal atmosphere.
radians.

---------------------------------------------------------------------------
LOW-LEVEL FUNCTIONS
---------------------------------------------------------------------------

vali = oclc_epq03_valget(arr,i)

oclc_epq03_valset, arr, i, vali

rarr = oclc_eqp03_r(yarr, delta_yarr, phiarr, ib, imax)