;+
;
; @file_comments
; 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 <proidl>MAP_2POINTS</proidl> with the meter parameter but for n
; points without do loop.
;
; @categories
; Maps
;
; @param Lon0 {in}{required}
; @param Lat0 {in}{required}
; np0 elements vector. longitudes and latitudes of np0 points P0
;
; @param Lon1 {in}{required}
; @param Lat1 {in}{required}
; np1 elements vector. longitude and latitude of np1 points P1
;
; @keyword AZIMUTH
; A named variable that will receive the azimuth of the great
; circle connecting the two points, P0 to P1
;
; @keyword MIDDLE
; to get the longitude/latitude of the middle point between P0 and P1.
;
; @keyword RADIANS
; if set, inputs and angular outputs are in radians, otherwise degrees.
;
; @keyword RADIUS {default=6378206.4d0}
; If given, return the distance between the two points calculated using the
; given radius.
; Default value is the Earth radius.
;
; @keyword TWO_BY_TWO
; If given, then <pro>map_npoints</pro> returns the distances between 
; number n of P0 points and number n of P1 pointsi.
; In that case, np0 and np1 must be equal.
;
; @returns
; An (np0,np1) array giving the distance in meter between np0
; points P0 and np1 points P1. Element (i,j) of the output is the
; distance between element P0[i] and P1[j].
; If keyword /TWO_BY_TWO is given then <pro>map_npoints</pro> returns
; an np-elements 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, $
;   IDL> 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], AZIMUTH=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
;
; @history
; Based on the IDL function map_2points.pro,v 1.6 2001/01/15
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; October 2003
;
; @version
; $Id$
;
;-
FUNCTION map_npoints, lon0, lat0, lon1, lat1, AZIMUTH=azimuth $
                    , RADIANS=radians, RADIUS=radius, MIDDLE=middle $
                    , TWO_BY_TWO=two_by_two
;
 compile_opt idl2, strictarrsubs
;
 IF (N_PARAMS() LT 4) THEN $
 ras = report('Incorrect number of arguments.')

 np0 = n_elements(lon0)
 IF n_elements(lat0) NE np0 THEN $
 ras = report('lon0 and lat0 must have the same number of elements')
 np1 = n_elements(lon1)
 IF n_elements(lat1) NE np1 THEN $
 ras = report('lon1 and lat1 must have the same number of elements')
 if keyword_set(two_by_two) AND np0 NE np1 then $
 ras = report('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 $
 ras = report('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