1 | ;+ |
---|
2 | ; |
---|
3 | ; @file_comments |
---|
4 | ;Return the distance in meter between all np0 points P0 and all |
---|
5 | ; np1 points P1 on a sphere. If keyword /TWO_BY_TWO is given then |
---|
6 | ; returns the distances between number n of P0 points and number |
---|
7 | ; n of P1 points (in that case, np0 and np1 must be equal). |
---|
8 | ; Same as map_2points with the meter parameter but for n points |
---|
9 | ; without do loop. |
---|
10 | ; |
---|
11 | ; @categories Maps |
---|
12 | ; |
---|
13 | ; @examples |
---|
14 | ;Result = Map_nPoints(lon0, lat0, lon1, lat1) |
---|
15 | ; |
---|
16 | ;@param Lon0 Lat0 {in}{required} np0 elements vector. longitudes and latitudes of np0 points P0 |
---|
17 | ;@param Lon1 Lat1 {in}{required} np1 elements vector. longitude and latitude of np1 points P1 |
---|
18 | ; |
---|
19 | ; @keyword AZIMUTH A named variable that will receive the azimuth of the great |
---|
20 | ; circle connecting the two points, P0 to P1 |
---|
21 | ; @keyword /MIDDLE to get the longitude/latitude of the middle point betwen P0 and P1. |
---|
22 | ; @keyword RADIANS = if set, inputs and angular outputs are in radians, otherwise |
---|
23 | ;degrees. |
---|
24 | ; @keyword RADIUS If given, return the distance between the two points |
---|
25 | ;calculated using the given radius. |
---|
26 | ; Default value is the earth radius : 6378206.4d0 |
---|
27 | ; @keyword TWO_BY_TWO:If given,then Map_nPoints returns the distances between |
---|
28 | ; number n of P0 points and number n of P1 points (in that case, |
---|
29 | ; np0 and np1 must be equal). |
---|
30 | ; |
---|
31 | ; @returns |
---|
32 | ; An (np0,np1) array giving the distance in meter between np0 |
---|
33 | ; points P0 and np1 points P1. Element (i,j) of the ouput is the |
---|
34 | ; distance between element P0[i] and P1[j]. |
---|
35 | ; If keyword /TWO_BY_TWO is given then Map_nPoints returns |
---|
36 | ; an np-element vector giving the distance in meter between P0[i] |
---|
37 | ; and P1[i] (in that case, we have np0 = np1 = np) |
---|
38 | ; if /MIDDLE see this keyword. |
---|
39 | ; |
---|
40 | ; @examples |
---|
41 | ;IDL> print, $ |
---|
42 | ;map_npoints([-105.15,1],[40.02,1],[-0.07,100,50],[51.30,20,0]) |
---|
43 | ; 7551369.3 5600334.8 |
---|
44 | ; 12864354. 10921254. |
---|
45 | ; 14919237. 5455558.8 |
---|
46 | ; |
---|
47 | ;IDL> lon0 = [-10, 20, 100] |
---|
48 | ;IDL> lat0 = [0, -10, 45] |
---|
49 | ;IDL> lon1 = [10, 60, 280] |
---|
50 | ;IDL> lat1 = [0, 10, 45] |
---|
51 | ;IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi) |
---|
52 | ;IDL> help, dist, azi |
---|
53 | ;DIST DOUBLE = Array[3, 3] |
---|
54 | ;AZI DOUBLE = Array[3, 3] |
---|
55 | ;IDL> print, dist[4*lindgen(3)], azi[4*lindgen(3)] |
---|
56 | ; 2226414.0 4957944.5 10018863. |
---|
57 | ; 90.000000 64.494450 4.9615627e-15 |
---|
58 | ;IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two) |
---|
59 | ;IDL> help, dist, azi |
---|
60 | ;DIST DOUBLE = Array[3] |
---|
61 | ;AZI DOUBLE = Array[3] |
---|
62 | ;IDL> print, dist, azi |
---|
63 | ; 2226414.0 4957944.5 10018863. |
---|
64 | ; 90.000000 64.494450 4.9615627e-15 |
---|
65 | ;IDL> print, map_2points(lon0[0], lat0[0], lon1[0], lat1[0]) |
---|
66 | ; 20.000000 90.000000 |
---|
67 | ;IDL> print, map_npoints(lon0[0], lat0[0], lon1[0], lat1[0], azi=azi)/6378206.4d0 / !dtor, azi |
---|
68 | ; 20.000000 |
---|
69 | ; 90.000000 |
---|
70 | ; |
---|
71 | ;IDL> lon0 = [-10, 20, 100] |
---|
72 | ;IDL> lat0 = [0, -10, 45] |
---|
73 | ;IDL> lon1 = [10, 60, 280] |
---|
74 | ;IDL> lat1 = [0, 10, 45] |
---|
75 | ;IDL> mid = map_npoints(lon0, lat0, lon1, lat1, /middle, /two_by_two) |
---|
76 | ;IDL> print, reform(mid[0,*]), reform(mid[1,*]) |
---|
77 | ; 0.0000000 40.000000 190.00000 |
---|
78 | ; 0.0000000 -1.5902773e-15 90.000000 |
---|
79 | ;IDL> print, (map_2points(lon0[0], lat0[0], lon1[0], lat1[0], npath = 3))[*, 1] |
---|
80 | ; 0.0000000 0.0000000 |
---|
81 | ;IDL> print, (map_2points(lon0[1], lat0[1], lon1[1], lat1[1], npath = 3))[*, 1] |
---|
82 | ; 40.000000 -1.5902773e-15 |
---|
83 | ;IDL> print, (map_2points(lon0[2], lat0[2], lon1[2], lat1[2], npath = 3))[*, 1] |
---|
84 | ; 190.00000 90.000000 |
---|
85 | ; |
---|
86 | ; @history |
---|
87 | ; Based on the IDL function map_2points.pro,v 1.6 2001/01/15 |
---|
88 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
89 | ; October 2003 |
---|
90 | ;- |
---|
91 | Function Map_npoints, lon0, lat0, lon1, lat1, azimuth = azimuth $ |
---|
92 | , RADIANS = radians, RADIUS = radius, MIDDLE = middle, TWO_BY_TWO = two_by_two |
---|
93 | |
---|
94 | COMPILE_OPT idl2 |
---|
95 | ON_ERROR, 2 ; return to caller |
---|
96 | |
---|
97 | IF (N_PARAMS() LT 4) THEN $ |
---|
98 | MESSAGE, 'Incorrect number of arguments.' |
---|
99 | |
---|
100 | np0 = n_elements(lon0) |
---|
101 | IF n_elements(lat0) NE np0 THEN $ |
---|
102 | MESSAGE, 'lon0 and lat0 must have the same number of elements' |
---|
103 | np1 = n_elements(lon1) |
---|
104 | IF n_elements(lat1) NE np1 THEN $ |
---|
105 | MESSAGE, 'lon1 and lat1 must have the same number of elements' |
---|
106 | if keyword_set(two_by_two) AND np0 NE np1 then $ |
---|
107 | MESSAGE, 'When using two_by_two keyword, P0 and P1 must have the same number of elements' |
---|
108 | |
---|
109 | mx = MAX(ABS([lat0, lat1])) |
---|
110 | pi2 = !dpi/2 |
---|
111 | IF (mx GT (KEYWORD_SET(radians) ? pi2 : 90)) THEN $ |
---|
112 | MESSAGE, 'Value of Latitude is out of allowed range.' |
---|
113 | |
---|
114 | k = KEYWORD_SET(radians) ? 1.0d0 : !dpi/180.0 |
---|
115 | ;Earth equatorial radius, meters, Clarke 1866 ellipsoid |
---|
116 | r_sphere = n_elements(RADIUS) NE 0 ? RADIUS : 6378206.4d0 |
---|
117 | ; |
---|
118 | coslt1 = cos(k*lat1[*]) |
---|
119 | sinlt1 = sin(k*lat1[*]) |
---|
120 | coslt0 = cos(k*lat0[*]) |
---|
121 | sinlt0 = sin(k*lat0[*]) |
---|
122 | ; |
---|
123 | IF np0 EQ np1 AND np1 EQ 1 THEN two_by_two = 1 |
---|
124 | ; |
---|
125 | if NOT keyword_set(two_by_two) THEN BEGIN |
---|
126 | coslt1 = replicate(1.0d0, np0)#temporary(coslt1) |
---|
127 | sinlt1 = replicate(1.0d0, np0)#temporary(sinlt1) |
---|
128 | coslt0 = temporary(coslt0)#replicate(1.0d0, np1) |
---|
129 | sinlt0 = temporary(sinlt0)#replicate(1.0d0, np1) |
---|
130 | ENDIF |
---|
131 | ; |
---|
132 | if keyword_set(two_by_two) THEN BEGIN |
---|
133 | cosl0l1 = cos(k*(lon1[*]-lon0[*])) |
---|
134 | sinl0l1 = sin(k*(lon1[*]-lon0[*])) |
---|
135 | ENDIF ELSE BEGIN |
---|
136 | cosl0l1 = cos(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1))) |
---|
137 | sinl0l1 = sin(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1))) |
---|
138 | ENDELSE |
---|
139 | |
---|
140 | cosc = sinlt0 * sinlt1 + coslt0 * coslt1 * cosl0l1 ;Cos of angle between pnts |
---|
141 | ; Avoid roundoff problems by clamping cosine range to [-1,1]. |
---|
142 | cosc = -1.0d0 > cosc < 1.0d0 |
---|
143 | ; |
---|
144 | if arg_present(azimuth) OR keyword_set(middle) then begin |
---|
145 | sinc = sqrt(1.0d0 - cosc*cosc) |
---|
146 | bad = where(abs(sinc) le 1.0e-7) |
---|
147 | IF bad[0] NE -1 THEN sinc[bad] = 1 |
---|
148 | cosaz = (coslt0 * sinlt1 - sinlt0*coslt1*cosl0l1) / sinc |
---|
149 | sinaz = sinl0l1*coslt1/sinc |
---|
150 | IF bad[0] NE -1 THEN BEGIN |
---|
151 | sinc[bad] = 0.0d0 |
---|
152 | sinaz[bad] = 0.0d0 |
---|
153 | cosaz[bad] = 1.0d0 |
---|
154 | ENDIF |
---|
155 | ENDIF |
---|
156 | ; |
---|
157 | IF keyword_set(middle) then BEGIN |
---|
158 | |
---|
159 | s0 = 0.5d0 * acos(cosc) |
---|
160 | ; |
---|
161 | coss = cos(s0) |
---|
162 | sins = sin(s0) |
---|
163 | ; |
---|
164 | lats = asin(sinlt0 * coss + coslt0 * sins * cosaz) / k |
---|
165 | lons = atan(sins * sinaz, coslt0 * coss - sinlt0 * sins * cosaz) / k |
---|
166 | ; |
---|
167 | if keyword_set(two_by_two) THEN BEGIN |
---|
168 | return, transpose([[lon0[*] + lons], [lats]]) |
---|
169 | ENDIF ELSE BEGIN |
---|
170 | return, [ [[lon0[*]#replicate(1.0d0, np1) + lons]], [[lats]] ] |
---|
171 | ENDELSE |
---|
172 | ; |
---|
173 | ENDIF |
---|
174 | ; |
---|
175 | if arg_present(azimuth) then begin |
---|
176 | azimuth = atan(sinaz, cosaz) |
---|
177 | IF k NE 1.0d0 THEN azimuth = temporary(azimuth) / k |
---|
178 | ENDIF |
---|
179 | return, acos(cosc) * r_sphere |
---|
180 | ; |
---|
181 | end |
---|