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