Changeset 125 for trunk/SRC/Interpolation/map_npoints.pro

Timestamp:

07/06/06 16:10:25 (18 years ago)

Author:

pinsard

Message:

improvements of Interpolation/*.pro header

File:

• trunk/SRC/Interpolation/map_npoints.pro (modified) (3 diffs)

trunk/SRC/Interpolation/map_npoints.pro

@@ -1 @@
-;+
-;
-; @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 map_2points with the meter parameter but for n points
-;       without do loop.
+;+ 
+; 
+; @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 map_2points 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 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 
+; @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 betwen 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 
+; @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 betwen 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 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).
+; @keyword 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.
 ;
 ; @returns
-;       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.
+; 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> 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]
@@ -53 +58 @@
 ; IDL> lon1 = [10, 60, 280]
 ; IDL> lat1 = [0, 10, 45]
-; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi)
+; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, AZIMUTH = azi)
 ; IDL> help, dist, azi
-; DIST            DOUBLE    = Array[3, 3]
-; AZI             DOUBLE    = Array[3, 3]
+; 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)
+; 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]
+; DIST DOUBLE = Array[3]
+; AZI DOUBLE = Array[3]
 ; IDL> print, dist, azi
-;       2226414.0       4957944.5       10018863.
-;       90.000000       64.494450   4.9615627e-15
+; 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
+; 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]
@@ -77 +82 @@
 ; IDL> lon1 = [10, 60, 280]
 ; IDL> lat1 = [0, 10, 45]
-; IDL> mid = map_npoints(lon0, lat0, lon1, lat1, /middle, /two_by_two)
+; 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
+; 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
+; 0.0000000 0.0000000
 ; IDL> print, (map_2points(lon0[1], lat0[1], lon1[1], lat1[1], npath = 3))[*, 1]
-;       40.000000  -1.5902773e-15
+; 40.000000 -1.5902773e-15
 ; IDL> print, (map_2points(lon0[2], lat0[2], lon1[2], lat1[2], npath = 3))[*, 1]
-;       190.00000       90.000000
+; 190.00000 90.000000
 ;
 ; @history
-;       Based on the IDL function map_2points.pro,v 1.6 2001/01/15
+; Based on the IDL function map_2points.pro,v 1.6 2001/01/15
 ; Sebastien Masson (smasson\@lodyc.jussieu.fr)
-;                  October 2003
+; 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
+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
+ COMPILE_OPT idl2, strictarrsubs
 
-  IF (N_PARAMS() LT 4) THEN $
-    MESSAGE, 'Incorrect number of arguments.'
+ 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'
+ 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.'
+ 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
+ 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 
+ 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[*])
+ 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 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 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 
+ 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
+ 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
+ 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 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
+ IF keyword_set(middle) then BEGIN
 
-    s0 = 0.5d0 * acos(cosc)
+ 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
+ coss = cos(s0)
+ sins = sin(s0)
 ;
-    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 
+ lats = asin(sinlt0 * coss + coslt0 * sins * cosaz) / k
+ lons = atan(sins * sinaz, coslt0 * coss - sinlt0 * sins * cosaz) / k
 ;
-  ENDIF
+ 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
 ;
-  if arg_present(azimuth) then begin
-    azimuth = atan(sinaz, cosaz)
-    IF k NE 1.0d0 THEN azimuth = temporary(azimuth) / k
-   ENDIF
+ 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
 ;