source: trunk/Interpolation/map_npoints.pro @ 59

;+
; NAME:
;       Map_nPoints
;
; PURPOSE:
;       Return the distance in meter between all np0 points P0 and all
;       np1 points P1 on a sphere. If keyword /TWO_BY_TWO is given then
;       returns the distances between number n of P0 points and number
;       n of P1 points (in that case, np0 and np1 must be equal).
;       Same as map_2points with the meter parameter but for n points
;       without do loop.
;
; CATEGORY:
;       Maps.
;
; CALLING SEQUENCE:
;       Result = Map_nPoints(lon0, lat0, lon1, lat1)
;
; INPUTS:
;       Lon0, Lat0 = np0 elements vector. longitudes and latitudes of
;       np0 points P0 
;       Lon1, Lat1 = np1 elements vector. longitude and latitude of
;       np1 points P1 
;
; KEYWORD PARAMETERS:
;
;   AZIMUTH: A named variable that will receive the azimuth of the great
;       circle  connecting the two points, P0 to P1
;   /MIDDLE: to get the longitude/latitude of the middle point betwen P0 and P1.
;   RADIANS = if set, inputs and angular outputs are in radians, otherwise
;       degrees.
;   RADIUS: If given, return the distance between the two points
;       calculated using the given radius.
;       Default value is the earth radius : 6378206.4d0
;   TWO_BY_TWO:If given,then Map_nPoints returns the distances between
;       number n of P0 points and number n of P1 points (in that case,
;       np0 and np1 must be equal).
;
; OUTPUTS:
;       An (np0,np1) array giving the distance in meter between np0
;       points P0 and np1 points P1. Element (i,j) of the ouput is the
;       distance between element P0[i] and P1[j].
;       If keyword /TWO_BY_TWO is given then Map_nPoints returns
;       an np-element vector giving the distance in meter between P0[i]
;       and P1[i] (in that case, we have np0 = np1 = np)
;       if /MIDDLE see this keyword.
;
; EXAMPLES:
;       IDL> print, $
;       map_npoints([-105.15,1],[40.02,1],[-0.07,100,50],[51.30,20,0])
;              7551369.3       5600334.8
;              12864354.       10921254.
;              14919237.       5455558.8
;       
;       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)
;       IDL> help, dist, azi
;       DIST            DOUBLE    = Array[3, 3]
;       AZI             DOUBLE    = Array[3, 3]
;       IDL> print, dist[4*lindgen(3)], azi[4*lindgen(3)]
;              2226414.0       4957944.5       10018863.
;              90.000000       64.494450   4.9615627e-15
;       IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two)
;       IDL> help, dist, azi
;       DIST            DOUBLE    = Array[3]
;       AZI             DOUBLE    = Array[3]
;       IDL> print, dist, azi
;              2226414.0       4957944.5       10018863.
;              90.000000       64.494450   4.9615627e-15
;       IDL> print, map_2points(lon0[0], lat0[0], lon1[0], lat1[0])
;              20.000000       90.000000
;       IDL> print, map_npoints(lon0[0], lat0[0], lon1[0], lat1[0], azi=azi)/6378206.4d0 / !dtor, azi
;              20.000000
;              90.000000
;
;       IDL> lon0 = [-10, 20, 100]
;       IDL> lat0 = [0, -10, 45]
;       IDL> lon1 = [10, 60, 280]
;       IDL> lat1 = [0, 10, 45]
;       IDL> mid = map_npoints(lon0, lat0, lon1, lat1, /middle, /two_by_two)
;       IDL> print, reform(mid[0,*]), reform(mid[1,*])
;              0.0000000       40.000000       190.00000
;              0.0000000  -1.5902773e-15       90.000000
;       IDL> print, (map_2points(lon0[0], lat0[0], lon1[0], lat1[0], npath = 3))[*, 1]
;              0.0000000       0.0000000
;       IDL> print, (map_2points(lon0[1], lat0[1], lon1[1], lat1[1], npath = 3))[*, 1]
;              40.000000  -1.5902773e-15
;       IDL> print, (map_2points(lon0[2], lat0[2], lon1[2], lat1[2], npath = 3))[*, 1]
;              190.00000       90.000000
;       
; MODIFICATION HISTORY:
;       Based on the IDL function map_2points.pro,v 1.6 2001/01/15
;       Sebastien Masson (smasson@lodyc.jussieu.fr)
;                  October 2003
;-
Function Map_npoints, lon0, lat0, lon1, lat1, azimuth = azimuth $
  ,  RADIANS = radians, RADIUS = radius, MIDDLE = middle, TWO_BY_TWO = two_by_two

  COMPILE_OPT idl2
  ON_ERROR, 2                   ; return to caller

  IF (N_PARAMS() LT 4) THEN $
    MESSAGE, 'Incorrect number of arguments.'

  np0 = n_elements(lon0) 
  IF n_elements(lat0) NE np0 THEN $
    MESSAGE, 'lon0 and lat0 must have the same number of elements'
  np1 = n_elements(lon1) 
  IF n_elements(lat1) NE np1 THEN $
    MESSAGE, 'lon1 and lat1 must have the same number of elements'
  if keyword_set(two_by_two) AND np0 NE np1 then $
    MESSAGE, 'When using two_by_two keyword, P0 and P1 must have the same number of elements'

  mx = MAX(ABS([lat0, lat1]))
  pi2 = !dpi/2
  IF (mx GT (KEYWORD_SET(radians) ? pi2 : 90)) THEN $
    MESSAGE, 'Value of Latitude is out of allowed range.'

  k = KEYWORD_SET(radians) ? 1.0d0 : !dpi/180.0
;Earth equatorial radius, meters, Clarke 1866 ellipsoid
  r_sphere =  n_elements(RADIUS) NE 0 ? RADIUS : 6378206.4d0 
;
  coslt1 = cos(k*lat1[*])
  sinlt1 = sin(k*lat1[*])
  coslt0 = cos(k*lat0[*])
  sinlt0 = sin(k*lat0[*])
;
  IF np0 EQ np1 AND np1 EQ 1 THEN two_by_two = 1
;
  if NOT keyword_set(two_by_two) THEN BEGIN 
    coslt1 = replicate(1.0d0, np0)#temporary(coslt1)
    sinlt1 = replicate(1.0d0, np0)#temporary(sinlt1)
    coslt0 = temporary(coslt0)#replicate(1.0d0, np1)
    sinlt0 = temporary(sinlt0)#replicate(1.0d0, np1)
  ENDIF 
;
  if keyword_set(two_by_two) THEN BEGIN 
    cosl0l1 = cos(k*(lon1[*]-lon0[*]))
    sinl0l1 = sin(k*(lon1[*]-lon0[*]))
  ENDIF ELSE BEGIN 
    cosl0l1 = cos(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1)))
    sinl0l1 = sin(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1)))
  ENDELSE 

  cosc = sinlt0 * sinlt1 + coslt0 * coslt1 * cosl0l1 ;Cos of angle between pnts
; Avoid roundoff problems by clamping cosine range to [-1,1].
  cosc = -1.0d0 > cosc < 1.0d0
;
  if arg_present(azimuth) OR keyword_set(middle) then begin
    sinc = sqrt(1.0d0 - cosc*cosc)
    bad = where(abs(sinc) le 1.0e-7)
    IF bad[0] NE -1 THEN sinc[bad] = 1
    cosaz = (coslt0 * sinlt1 - sinlt0*coslt1*cosl0l1) / sinc 
    sinaz = sinl0l1*coslt1/sinc
    IF bad[0] NE -1 THEN BEGIN 
      sinc[bad] = 0.0d0
      sinaz[bad] = 0.0d0
      cosaz[bad] = 1.0d0
    ENDIF
  ENDIF
;
  IF keyword_set(middle) then BEGIN

    s0 = 0.5d0 * acos(cosc)
 ;
    coss = cos(s0)
    sins = sin(s0)
;    
    lats = asin(sinlt0 * coss + coslt0 * sins * cosaz) / k
    lons = atan(sins * sinaz, coslt0 * coss - sinlt0 * sins * cosaz) / k
;
    if keyword_set(two_by_two) THEN BEGIN 
      return, transpose([[lon0[*] + lons], [lats]])
    ENDIF ELSE BEGIN 
      return, [ [[lon0[*]#replicate(1.0d0, np1) + lons]], [[lats]] ]
    ENDELSE 
;
  ENDIF
;
  if arg_present(azimuth) then begin
    azimuth = atan(sinaz, cosaz)
    IF k NE 1.0d0 THEN azimuth = temporary(azimuth) / k
   ENDIF
 return, acos(cosc) * r_sphere
;
end