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