;+
;
; @file_comments
; This function returns the longitude and latitude [lon, lat] of
;a point a given arc distance (-pi <= Arc_Dist <= pi), and azimuth (Az),
;from a specified location Lon0, lat0.
;       Same as LL_ARC_DISTANCE but for n points without do loop.
;
; @categories Mapping, geography
;
; @examples 
;Result = LL_NARCS_DISTANCES(Lon, lat0, Arc_Dist, Az)
;
;    @param Lon0 {in}{required}  An array containing the longitude of the starting point.
;             Values are assumed to be in radians unless the keyword
;             DEGREES is set.
;    @param Lat0 {in}{required}  An array containing the latitude of the starting point.
;             Values are assumed to be in radians unless the keyword
;             DEGREES is set.
;    @param Arc_Dist {in}{required}  The arc distance from Lon_lat0. The value must be between
; -!PI and +!PI. To express distances in arc units, divide
;  by the radius of the globe expressed in the original units.
;  For example, if the radius of the earth is 6371 km, divide
;  the distance in km by 6371 to obtain the arc distance.    
;    @param Az {in}{required}   The azimuth from Lon_lat0. The value is assumed to be in
;  radians unless the keyword DEGREES is set.
;
; @keyword    DEGREES  Set this keyword to express all measurements and
;  results in degrees.
;
; @returns
; a (2, n) array containing the 
;       longitude / latitude of the resultings points. Values are in radians
;       unless the keyword DEGREES is set.
;
; @file_comments
;Formula from Map Projections - a working manual.  USGS paper
;1395.  Equations (5-5) and (5-6).
;
; @examples 
;Lon_lat0 = [1.0, 2.0]; Initial point specified in radians
;Arc_Dist = 2.0; Arc distance in radians
;Az = 1.0; Azimuth in radians
;Result = LL_ARC_DISTANCE(Lon_lat0, Arc_Dist, Az)
;PRINT, Result
;       2.91415    -0.622234
;
;IDL> lon0 = [-10, 20, 100]
;IDL> lat0 = [0, -10, 45]
;IDL> lon1 = [10, 60, 280]
;IDL> lat1 = [0, 10, 45]
;IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two)
;IDL> earthradius = 6378206.4d0
;IDL> res = ll_narcs_distances(lon0, lat0, dist/earthradius, azi, /degrees)
;IDL> print, reform(res[0, *])
;       10.000000       60.000000       280.00000
;IDL> print, reform(res[1, *])
;          1.1999280e-15       10.000000       45.000000
;
; @history
;       Based on the IDL function ll_arc_distance.pro,v 1.11 2003/02/03
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                  August 2005
;-

;+
; @file_comments Return the [lon, lat] of the point a given arc distance 
;(-pi <= arc_dist <= pi),
; and azimuth (az), from lon_lat0.
;-
;
FUNCTION LL_NARCS_DISTANCES, lon0, lat0, arc_dist, az, DEGREES = degs
;
  IF n_elements(lon0) NE n_elements(lat0) $
    OR n_elements(lon0) NE n_elements(arc_dist) $
    OR n_elements(lon0) NE n_elements(az) THEN return, -1

  cdist = cos(arc_dist[*])      ;Arc_Dist is always in radians.
  sdist = sin(arc_dist[*])

  if keyword_set(degs) then s = !dpi/180.0 else s = 1.0d0

  ll = lat0[*] * s              ;To radians
  sinll1 = sin(ll)
  cosll1 = cos(ll)
  azs = az[*] * s
  phi = asin(sinll1 * cdist + cosll1 * sdist * cos(azs))
  ll = lon0[*] * s              ;To radians
  lam = ll + atan(sdist * sin(azs), $
                  cosll1 * cdist - sinll1 * sdist * cos(azs))

  zero = where(arc_dist eq 0, count)
  IF count NE 0 THEN BEGIN 
    lam[zero] = lon0[zero]
    phi[zero] = lat0[zero]
  ENDIF 

  if keyword_set(degs) then return, transpose([[lam], [phi]]) / s $
  ELSE return, transpose([[lam], [phi]])

end