Render a rectangular projection map to a sphere.
  The most confusing aspect of this program is by far the image orientation.
  I can't make the conventions natural for everyone (not even myself), but
  here's what this program does:

  In the output image there is a natural coordinate system, (x,y).  Each
  row, ie., image[*,i] corresponds to a fixed value of y.  image[*,0]
  is the row with the LOWEST value of y.  This may seem strange, but it
  permits a "normal" orientation of the image if displayed with !order=0
  (the usual IDL default).  The x direction (image[i,*]) works in the
  obvious way, image[0,*] is the column with the LOWEST value of x.  So,
  viewed with !ORDER=0, this is an image that has up at the top.  The
  top is also considered NORTH when dealing with images that one sees in
  the sky (or relative to your point of view).

  The position angle of the pole sets the rotation of the sphere in the
  plane of the sky (image).  North is to the top (+y) and East is to the
  left (-x).  The position angle (POLE) is the angle of the apparent
  North pole of the sphere relative to up (North) in the image, measured
  From up (North) counter-clockwise (also known as eastward from North).

  Visualizing the indexing of the input map is equally confusing.  In
  the coordinate system containing the map, North is +y, East is -x, as
  before.  These coordinates are plane-of-sky on the object.  The map
  has its own coordinates, latitude and longitude.  The first index into
  the map is treated like x but is actually East Longitude.  The left edge
  of image[0,*] is precisely at 0 degrees longitude.  As the first index
  increases (increasing x?), the longitude increases.  The right edge of
  the last column in the map is 360 degrees.

  Ok, here is the weird part.  The first row in the map is the NORTHERN-most
  set of pixels.  The last row is the SOUTHERN-most set.  So, map[0,0]
  touches the prime meridian and the north pole.
      Thus map[i,j]:
           longitude = (i+0.5)/nl * 360.0
           latitude  = 90.0 - (j+0.5)/nt * 180.0
        where nl is the width of map and nt is the height of the map.
      These formulas give the lat,lon of the CENTER of the pixel in degrees.

  Now, understanding the sub-"earth" and sub-solar longitudes becomes
  easy.  The sub-"earth" point is the lat,lon of the map that is in the
  center of the projected disk.  The sub-solar point is the lat,lon of
  the map that is nearest to the sun (normal solar illumination).  Note,
  that as the object rotates, the sub-earth longitude will decrease, not

  Image display
  map - Array containing a full map of surface
  radius - Radius of object (eg., in km)
  scale  - Scale of image (eg., km/pixel)
  pole   - Position angle of pole, east from north (degrees)
  selat  - Sub-earth latitude (degrees)
  selon  - Sub-earth longitude (degrees)
  sslat  - Sub-solar latitude (degrees).  Only used for Hapke functions.
  sslon  - Sub-solar longitude (degrees).  Only used for Hapke functions.
  nx     - X size of output image in pixels.
  ny     - Y xize of output image in pixels.
  SILENT - Flag, if set suppresses all printed output
  NODISPLAY - Flag, if set suppresses graphical output
  GEOM - Undefined (default), no special action.
         If defined but not a structure, then upon return geom will contain
            all the geometric information needed to render the image with a
            new map.
         If defined and is a structure, then the contents are taken to be
            that needed to do the final image calculation without redoing
            all the geometric computation.
  XOFFSET - x offset, of object from center of image (in pixels) default=0
  YOFFSET - y offset, of object from center of image (in pixels) default=0

  Limb darkening model to use (set one of these, at most):
  LIN   - Linear limb-darkening coeff, map is normal albedo
  MIN   - Minnaert limb-darkening coeff, map is normal albedo
  HAPKE - Hapke scattering model, map is single scattern albedo
             HAPKE[0] = h (old style, circa 1981)
             HAPKE[1] = P(0)
  HAP2  - 2nd generation Hapke scattering model, 1986 and the book,
            "Theory of Reflectance and Emittance Spectroscopy"
             HAP2[0] = h (new style, circa 1986, p. 226 in book)
             HAP2[1:PPARMS] = P(g), or HG function parameters
             HAP2[PPARMS+1] = B0, Emperical backscatter factor (p. 228)
             HAP2[PPARMS+2] = Theta(bar), surface roughness
                                 parameter (radians).
           PPARMS is deduced from the length of this vector.  The minimum
             length is 4, for which PPARMS=1.  Only PPARMS=1 is allowed
             if you are not using a single-particle scattering function.
  H93   - When used with HAP2, will select which approximation to the
             H-function is used.  Default is the simplest (and fastest)
             formula from the 1981 Hapke paper.  H93, when set, will
             force the use of eqn. 8.57 from the Hapke book (from 1993).
  Pfn   - Specify function to use for P(g) instead of the default constant.
            The function must be a procedure taking arguments g,a,F,/radians
              g  phase angle in radians (with keyword /radians set)
              a  an array of Pparms parameters
              F  phase function evaluated at phase angles g
            Use "" as a model.
          The default is to use the scalar value in HAP2 for the value of
          the phase function.  This is ignored if HAP2 not being used.

  image  - Rendered image of object.  Return is a traditional I/F value.
  XARR - Optional return of the array of plane-of-sky x positions for each
           pixel in the image.  These values are in units of the object
  YARR - Optional return of the array of plane-of-sky y positions for each
           pixel in the image.  These values are in units of the object
   95/07/24, Written by Marc W. Buie, Lowell Observatory
   96/04/11, MWB, fixed numerous bugs, now works for sun and viewpoint
                    being different
   97/08/21, MWB, added HAP2 keyword
   97/08/27, MWB, added GEOM keyword
   97/09/18, MWB, repaired incorrect mapping of map to surface.
   2003/04/03, MWB, added XOFFSET and YOFFSET keywords
   2009/02/01, MWB, improved parameter validation and some documentation
   2012/01/29, MWB, workaround added for IDL v8.1, was crashing IDL if
                     xoffset/yoffset supplied as double precision
   2015/05/13, MWB, added support for single-particle scattering functions