NAME:
  airmass
 PURPOSE: (one line)
  Compute airmass for one or more times.
 DESCRIPTION:
  This is should be a pretty good function for computing the air mass factor.
  The default is to use the cosine based formula derived by David Tholen
  but the older secant based formula from Hardie is still available.  The
  zenith angle is corrected for refraction before using either formula (see
  REFRAC).  The defaults on the atmospheric conditions are STP.  This function
  should be quite good up to 5 airmasses.  This formula will work up to
  a zenith angle of 80 degrees after which the computation id not done.
 CATEGORY:
  Astronomy
 CALLING SEQUENCE:
  am = airmass(jd,ra,dec,lat,lon,wave,pressure,temp,relhum)
 INPUTS:
  jd  - Julian date (must be double precision to get nearest second).
  ra  - Right ascension (of date) in radians.
  dec - Declination (of date) in radians.
  lat - Latitude of observatory in radians.
  lon - West longitude of observatory in radians.
 OPTIONAL INPUT PARAMETERS:
  wave     - wavelength of light, in microns (default=0.56)
  pressure - atmospheric pressure in mm of Hg (default=760.0)
  temp     - atmospheric temperature in degrees C (default=0.0)
  relhum   - Relative humidity (in percent) (default=0.0)
 KEYWORD INPUT PARAMETERS:
  UT  - Time, in hours to add to JD to get the correct Universal Time.
  HARDIE - Flag, if set causes Hardie formula to be used.
 KEYWORD OUTPUT PARAMETERS:
  ALT - Optional return of the altitude for each airmass.
  LHA - Optional return of the local hour angle.
  LST - Optional return of the local sidereal time.
  AZI - Optional return of the azimuth (west from south).
 OUTPUTS:
  Return value is the airmass in single precision.
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
 RESTRICTIONS:
  Any input may be a vector.  If more than one is a vector then the
  lengths must match.  The return will have the same dimensions as
  the input.
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 1992 March, by Marc W. Buie, Lowell Observatory
  94/05/05 - MWB, modified to split out the LST calculation, added LST and
                LHA optional keyword outputs and UT keyword input.'
  97/03/03 - MWB, added the Tholen airmass equation as the default.  Hardie
                is still included as an option.  Also added refraction
                to the calculation which added numerous inputs.
  97/10/23 - MWB, added the AZI keyword.
  2002/01/08, MWB, now allows 2-d inputs on jd,ra,dec