1 | ;+ |
---|
2 | ; |
---|
3 | ; @file_comments |
---|
4 | ; This function returns the longitude and latitude [lon, lat] of |
---|
5 | ; a point a given arc distance (-pi <= Arc_Dist <= pi), and azimuth (Az), |
---|
6 | ; from a specified location Lon0, Lat0. |
---|
7 | ; Same as <proidl>LL_ARC_DISTANCE</proidl> but for n points without do loop. |
---|
8 | ; |
---|
9 | ; Formula from Map Projections - a working manual. USGS paper |
---|
10 | ; 1395. Equations (5-5) and (5-6). |
---|
11 | ; |
---|
12 | ; @categories |
---|
13 | ; Mapping, geography |
---|
14 | ; |
---|
15 | ; @param Lon0 {in}{required} |
---|
16 | ; An array containing the longitude of the starting point. |
---|
17 | ; Values are assumed to be in radians unless the keyword DEGREES is set. |
---|
18 | ; |
---|
19 | ; @param Lat0 {in}{required} |
---|
20 | ; An array containing the latitude of the starting point. |
---|
21 | ; Values are assumed to be in radians unless the keyword DEGREES is set. |
---|
22 | ; |
---|
23 | ; @param Arc_Dist {in}{required} |
---|
24 | ; The arc distance from Lon_lat0. The value must be between |
---|
25 | ; -!PI and +!PI. To express distances in arc units, divide |
---|
26 | ; by the radius of the globe expressed in the original units. |
---|
27 | ; For example, if the radius of the earth is 6371 km, divide |
---|
28 | ; the distance in km by 6371 to obtain the arc distance. |
---|
29 | ; |
---|
30 | ; @param Az {in}{required} |
---|
31 | ; The azimuth from Lon_lat0. The value is assumed to be in |
---|
32 | ; radians unless the keyword DEGREES is set. |
---|
33 | ; |
---|
34 | ; @keyword DEGREES |
---|
35 | ; Set this keyword to express all measurements and results in degrees. |
---|
36 | ; |
---|
37 | ; @returns |
---|
38 | ; a (2,n) array containing the longitude/latitude of the resulting points. |
---|
39 | ; Values are in radians unless the keyword DEGREES is set. |
---|
40 | ; |
---|
41 | ; @examples |
---|
42 | ; |
---|
43 | ; IDL> Lon_lat0 = [1.0, 2.0]; Initial point specified in radians |
---|
44 | ; IDL> Arc_Dist = 2.0; Arc distance in radians |
---|
45 | ; IDL> Az = 1.0; Azimuth in radians |
---|
46 | ; IDL> Result = LL_ARC_DISTANCE(Lon_lat0, Arc_Dist, Az) |
---|
47 | ; IDL> PRINT, Result |
---|
48 | ; 2.91415 -0.622234 |
---|
49 | ; |
---|
50 | ; IDL> lon0 = [-10, 20, 100] |
---|
51 | ; IDL> lat0 = [0, -10, 45] |
---|
52 | ; IDL> lon1 = [10, 60, 280] |
---|
53 | ; IDL> lat1 = [0, 10, 45] |
---|
54 | ; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two) |
---|
55 | ; IDL> earthradius = 6378206.4d0 |
---|
56 | ; IDL> res = ll_narcs_distances(lon0, lat0, dist/earthradius, azi, /degrees) |
---|
57 | ; IDL> print, reform(res[0, *]) |
---|
58 | ; 10.000000 60.000000 280.00000 |
---|
59 | ; IDL> print, reform(res[1, *]) |
---|
60 | ; 1.1999280e-15 10.000000 45.000000 |
---|
61 | ; |
---|
62 | ; @history |
---|
63 | ; Based on the IDL function ll_arc_distance.pro,v 1.11 2003/02/03 |
---|
64 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
65 | ; August 2005 |
---|
66 | ; |
---|
67 | ; @version |
---|
68 | ; $Id$ |
---|
69 | ; |
---|
70 | ;- |
---|
71 | FUNCTION ll_narcs_distances, lon0, lat0, arc_dist, az, DEGREES=degs |
---|
72 | ; |
---|
73 | compile_opt idl2, strictarrsubs |
---|
74 | ; |
---|
75 | ; |
---|
76 | IF n_elements(lon0) NE n_elements(lat0) $ |
---|
77 | OR n_elements(lon0) NE n_elements(arc_dist) $ |
---|
78 | OR n_elements(lon0) NE n_elements(az) THEN return, -1 |
---|
79 | |
---|
80 | cdist = cos(arc_dist[*]) ;Arc_Dist is always in radians. |
---|
81 | sdist = sin(arc_dist[*]) |
---|
82 | |
---|
83 | if keyword_set(degs) then s = !dpi/180.0 else s = 1.0d0 |
---|
84 | |
---|
85 | ll = lat0[*] * s ;To radians |
---|
86 | sinll1 = sin(ll) |
---|
87 | cosll1 = cos(ll) |
---|
88 | azs = az[*] * s |
---|
89 | phi = asin(sinll1 * cdist + cosll1 * sdist * cos(azs)) |
---|
90 | ll = lon0[*] * s ;To radians |
---|
91 | lam = ll + atan(sdist * sin(azs), $ |
---|
92 | cosll1 * cdist - sinll1 * sdist * cos(azs)) |
---|
93 | |
---|
94 | zero = where(arc_dist eq 0, count) |
---|
95 | IF count NE 0 THEN BEGIN |
---|
96 | lam[zero] = lon0[zero] |
---|
97 | phi[zero] = lat0[zero] |
---|
98 | ENDIF |
---|
99 | |
---|
100 | if keyword_set(degs) then return, transpose([[lam], [phi]]) / s $ |
---|
101 | ELSE return, transpose([[lam], [phi]]) |
---|
102 | |
---|
103 | end |
---|