;+ ; 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