NAME: 
  render
 PURPOSE: 
  Render a rectangular projection map to a sphere.
 DESCRIPTION:
  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
  increase.

 CATEGORY:
  Image display
 CALLING SEQUENCE:
  render,map,radius,scale,pole,selat,selon,sslat,sslon,nx,ny,image
 INPUTS:
  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.
 OPTIONAL INPUT PARAMETERS:
 KEYWORD INPUT PARAMETERS:
  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] = P(g)
                HAP2[2] = B0, Emperical backscatter factor (p. 228)
                HAP2[3] = Theta(bar), surface roughness parameter.
     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).

 OUTPUTS:
  image  - Rendered image of object.  Return is a traditional I/F value.
 KEYWORD OUTPUT PARAMETERS:
  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
           radii.
  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
           radii.
 COMMON BLOCKS:
 SIDE EFFECTS:
 RESTRICTIONS:
 PROCEDURE:
 MODIFICATION HISTORY:
   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