source: trunk/SRC/Interpolation/ll_narcs_distances.pro @ 242

;+
;
; @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 <proidl>LL_ARC_DISTANCE</proidl> but for n points without do loop.
;
; Formula from Map Projections - a working manual.  USGS paper
; 1395. Equations (5-5) and (5-6).
;
; @categories
; Mapping, geography
;
; @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 resulting points.
; Values are in radians unless the keyword DEGREES is set.
;
; @examples
; IDL> Lon_lat0 = [1.0, 2.0]; Initial point specified in radians
; IDL> Arc_Dist = 2.0; Arc distance in radians
; IDL> Az = 1.0; Azimuth in radians
; IDL> Result = LL_ARC_DISTANCE(Lon_lat0, Arc_Dist, Az)
; IDL> 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
;
; @version
; $Id$
;
;-
;
FUNCTION ll_narcs_distances, lon0, lat0, arc_dist, az, DEGREES = degs
;
  compile_opt idl2, strictarrsubs
;
;
  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