1 /*
2 See the top of datum.h for information on this code.
3 N7TAP
4
5 Portions Copyright (C) 2000-2019 The Xastir Group
6
7 */
8
9
10 #ifdef HAVE_CONFIG_H
11 #include "config.h"
12 #endif // HAVE_CONFIG_H
13
14 #include "snprintf.h"
15
16 #include <math.h>
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <ctype.h>
20
21 #include "xastir.h"
22 #include "datum.h"
23 #include "main.h"
24 #include "util.h"
25
26 // Must be last include file
27 #include "leak_detection.h"
28
29
30
31 // ellipsoid: index into the gEllipsoid[] array, in which
32 // a is the ellipsoid semimajor axis
33 // invf is the inverse of the ellipsoid flattening f
34 // dx, dy, dz: ellipsoid center with respect to WGS84 ellipsoid center
35 // x axis is the prime meridian
36 // y axis is 90 degrees east longitude
37 // z axis is the axis of rotation of the ellipsoid
38
39 // The following values for dx, dy and dz were extracted from the output of
40 // the GARMIN PCX5 program. The output also includes values for da and df, the
41 // difference between the reference ellipsoid and the WGS84 ellipsoid semi-
42 // major axis and flattening, respectively. These are replaced by the
43 // data contained in the structure array gEllipsoid[], which was obtained from
44 // the Defence Mapping Agency document number TR8350.2, "Department of Defense
45 // World Geodetic System 1984."
46
47 /*
48 The above are the original comments by John Waers.
49
50 Curt Mills, WE7U wrote a perl version of this code and added more datums
51 from Peter H. Dana's website.
52
53 Reference Ellipsoids:
54 http://www.Colorado.EDU/geography/gcraft/notes/datum/elist.html
55
56 Reference Ellipsoids and Datums:
57 http://www.Colorado.EDU/geography/gcraft/notes/datum/edlist.html
58
59 I've loaded the numbers from that second, newer web page.
60 N7TAP
61 */
62
63 /* Keep the enum in datum.h up to date with the order of this array */
64 const Ellipsoid gEllipsoid[] =
65 {
66 // name a 1/f
67 { "Airy 1830", 6377563.396, 299.3249646 },
68 { "Modified Airy", 6377340.189, 299.3249646 },
69 { "Australian National", 6378160.0, 298.25 },
70 { "Bessel 1841", 6377397.155, 299.1528128 },
71 { "Bessel 1841 (Namibia)", 6377483.865, 299.1528128 },
72 { "Clarke 1866", 6378206.4, 294.9786982 },
73 { "Clarke 1880", 6378249.145, 293.465 },
74 { "Everest (India 1830)", 6377276.345, 300.8017 },
75 { "Everest (India 1956)", 6377301.243, 300.8017 },
76 { "Everest (Sabah Sarawak)", 6377298.556, 300.8017 },
77 { "Everest (Malaysia 1969)", 6377295.664, 300.8017 },
78 { "Everest (Malay. & Sing)", 6377304.063, 300.8017 },
79 { "Everest (Pakistan)", 6377309.613, 300.8017 },
80 { "Fischer 1960 (Mercury)", 6378166.0, 298.3 },
81 { "Modified Fischer 1960", 6378155.0, 298.3 },
82 { "Fischer 1968", 6378150.0, 298.3 },
83 { "Helmert 1906", 6378200.0, 298.3 },
84 { "Hough 1960", 6378270.0, 297.0 },
85 { "Indonesian 1974", 6378160.0, 298.247 },
86 { "International 1924", 6378388.0, 297.0 },
87 { "Krassovsky 1940", 6378245.0, 298.3 },
88 { "GRS 67", 6378160.0, 298.247167427 },
89 { "GRS 80", 6378137.0, 298.257222101 },
90 { "South American 1969", 6378160.0, 298.25 },
91 { "WGS 60", 6378165.0, 298.3 },
92 { "WGS 66", 6378145.0, 298.25 },
93 { "WGS 72", 6378135.0, 298.26 },
94 { "WGS 84", 6378137.0, 298.257223563 }
95 };
96
97
98
99
100
101 /* Keep correct indices to commonly used datums in the enum in datum.h */
102 /* Feel free to add mnemonic indices for datums that you use */
103 const Datum gDatum[] =
104 {
105 // name ellipsoid dx dy dz
106 { "Adindan (Burkina Faso)", E_CLARKE_80, -118, -14, 218 }, // 0
107 { "Adindan (Cameroon)", E_CLARKE_80, -134, -2, 210 }, // 1
108 { "Adindan (Ethiopia)", E_CLARKE_80, -165, -11, 206 }, // 2
109 { "Adindan (Mali)", E_CLARKE_80, -123, -20, 220 }, // 3
110 { "Adindan (MEAN FOR Ethiopia; Sudan)", E_CLARKE_80, -166, -15, 204 }, // 4
111 { "Adindan (Senegal)", E_CLARKE_80, -128, -18, 224 }, // 5
112 { "Adindan (Sudan)", E_CLARKE_80, -161, -14, 205 }, // 6
113 { "Afgooye (Somalia)", E_KRASS_40, -43, -163, 45 }, // 7
114 { "Ain el Abd 1970 (Bahrain)", E_INT_24, -150, -250, -1 }, // 8
115 { "Ain el Abd 1970 (Saudi Arabia)", E_INT_24, -143, -236, 7 }, // 9
116 { "American Samoa 1962 (American Samoa Islands)", E_CLARKE_66, -115, 118, 426 }, // 10
117 { "Anna 1 Astro 1965 (Cocos Islands)", E_AUS_NAT, -491, -22, 435 }, // 11
118 { "Antigua Island Astro 1943 (Antigua (Leeward Islands))", E_CLARKE_80, -270, 13, 62 }, // 12
119 { "Arc 1950 (Botswana)", E_CLARKE_80, -138, -105, -289 }, // 13
120 { "Arc 1950 (Burundi)", E_CLARKE_80, -153, -5, -292 }, // 14
121 { "Arc 1950 (Lesotho)", E_CLARKE_80, -125, -108, -295 }, // 15
122 { "Arc 1950 (Malawi)", E_CLARKE_80, -161, -73, -317 }, // 16
123 { "Arc 1950 (MEAN FOR Botswana; Lesotho; Malawi; Swaziland; Zaire; Zambia; Zimbabwe)", E_CLARKE_80, -143, -90, -294 }, // 17
124 { "Arc 1950 (Swaziland)", E_CLARKE_80, -134, -105, -295 }, // 18
125 { "Arc 1950 (Zaire)", E_CLARKE_80, -169, -19, -278 }, // 19
126 { "Arc 1950 (Zambia)", E_CLARKE_80, -147, -74, -283 }, // 20
127 { "Arc 1950 (Zimbabwe)", E_CLARKE_80, -142, -96, -293 }, // 21
128 { "Arc 1960 (MEAN FOR Kenya; Tanzania)", E_CLARKE_80, -160, -6, -302 }, // 22
129 { "Arc 1960 (Kenya)", E_CLARKE_80, -157, -2, -299 }, // 23
130 { "Arc 1960 (Taanzania)", E_CLARKE_80, -175, -23, -303 }, // 24
131 { "Ascension Island 1958 (Ascension Island)", E_INT_24, -205, 107, 53 }, // 25
132 { "Astro Beacon E 1945 (Iwo Jima)", E_INT_24, 145, 75, -272 }, // 26
133 { "Astro DOS 71/4 (St Helena Island)", E_INT_24, -320, 550, -494 }, // 27
134 { "Astro Tern Island (FRIG) 1961 (Tern Island)", E_INT_24, 114, -116, -333 }, // 28
135 { "Astronomical Station 1952 (Marcus Island)", E_INT_24, 124, -234, -25 }, // 29
136 { "Australian Geodetic 1966 (Australia; Tasmania)", E_AUS_NAT, -133, -48, 148 }, // 30
137 { "Australian Geodetic 1984 (Australia; Tasmania)", E_AUS_NAT, -134, -48, 149 }, // 31
138 { "Ayabelle Lighthouse (Djibouti)", E_CLARKE_80, -79, -129, 145 }, // 32
139 { "Bellevue (IGN) (Efate & Erromango Islands)", E_INT_24, -127, -769, 472 }, // 33
140 { "Bermuda 1957 (Bermuda)", E_CLARKE_66, -73, 213, 296 }, // 34
141 { "Bissau (Guinea-Bissau)", E_INT_24, -173, 253, 27 }, // 35
142 { "Bogota Observatory (Colombia)", E_INT_24, 307, 304, -318 }, // 36
143 { "Bukit Rimpah (Indonesia (Bangka & Belitung Ids))", E_BESS_41, -384, 664, -48 }, // 37
144 { "Camp Area Astro (Antarctica (McMurdo Camp Area))", E_INT_24, -104, -129, 239 }, // 38
145 { "Campo Inchauspe (Argentina)", E_INT_24, -148, 136, 90 }, // 39
146 { "Canton Astro 1966 (Phoenix Islands)", E_INT_24, 298, -304, -375 }, // 40
147 { "Cape (South Africa)", E_CLARKE_80, -136, -108, -292 }, // 41
148 { "Cape Canaveral (Bahamas; Florida)", E_CLARKE_66, -2, 151, 181 }, // 42
149 { "Carthage (Tunisia)", E_CLARKE_80, -263, 6, 431 }, // 43
150 { "Chatham Island Astro 1971 (New Zealand (Chatham Island))", E_INT_24, 175, -38, 113 }, // 44
151 { "Chua Astro (Paraguay)", E_INT_24, -134, 229, -29 }, // 45
152 { "Corrego Alegre (Brazil)", E_INT_24, -206, 172, -6 }, // 46
153 { "Dabola (Guinea)", E_CLARKE_80, -83, 37, 124 }, // 47
154 { "Deception Island (Deception Island; Antarctia)", E_CLARKE_80, 260, 12, -147 }, // 48
155 { "Djakarta (Batavia) (Indonesia (Sumatra))", E_BESS_41, -377, 681, -50 }, // 49
156 { "DOS 1968 (New Georgia Islands (Gizo Island))", E_INT_24, 230, -199, -752 }, // 50
157 { "Easter Island 1967 (Easter Island)", E_INT_24, 211, 147, 111 }, // 51
158 { "Estonia; Coordinate System 1937 (Estonia)", E_BESS_41, 374, 150, 588 }, // 52
159 { "European 1950 (Cyprus)", E_INT_24, -104, -101, -140 }, // 53
160 { "European 1950 (Egypt)", E_INT_24, -130, -117, -151 }, // 54
161 { "European 1950 (England; Channel Islands; Scotland; Shetland Islands)", E_INT_24, -86, -96, -120 }, // 55
162 { "European 1950 (England; Ireland; Scotland; Shetland Islands)", E_INT_24, -86, -96, -120 }, // 56
163 { "European 1950 (Finland; Norway)", E_INT_24, -87, -95, -120 }, // 57
164 { "European 1950 (Greece)", E_INT_24, -84, -95, -130 }, // 58
165 { "European 1950 (Iran)", E_INT_24, -117, -132, -164 }, // 59
166 { "European 1950 (Italy (Sardinia))", E_INT_24, -97, -103, -120 }, // 60
167 { "European 1950 (Italy (Sicily))", E_INT_24, -97, -88, -135 }, // 61
168 { "European 1950 (Malta)", E_INT_24, -107, -88, -149 }, // 62
169 { "European 1950 (MEAN FOR Austria; Belgium; Denmark; Finland; France; W Germany; Gibraltar; Greece; Italy; Luxembourg; Netherlands; Norway; Portugal; Spain; Sweden; Switzerland)", E_INT_24, -87, -98, -121 }, // 63
170 { "European 1950 (MEAN FOR Austria; Denmark; France; W Germany; Netherlands; Switzerland)", E_INT_24, -87, -96, -120 }, // 64
171 { "European 1950 (MEAN FOR Iraq; Israel; Jordan; Lebanon; Kuwait; Saudi Arabia; Syria)", E_INT_24, -103, -106, -141 }, // 65
172 { "European 1950 (Portugal; Spain)", E_INT_24, -84, -107, -120 }, // 66
173 { "European 1950 (Tunisia)", E_INT_24, -112, -77, -145 }, // 67
174 { "European 1979 (MEAN FOR Austria; Finland; Netherlands; Norway; Spain; Sweden; Switzerland)", E_INT_24, -86, -98, -119 }, // 68
175 { "Fort Thomas 1955 (Nevis; St. Kitts (Leeward Islands))", E_CLARKE_80, -7, 215, 225 }, // 69
176 { "Gan 1970 (Republic of Maldives)", E_INT_24, -133, -321, 50 }, // 70
177 { "Geodetic Datum 1949 (New Zealand)", E_INT_24, 84, -22, 209 }, // 71
178 { "Graciosa Base SW 1948 (Azores (Faial; Graciosa; Pico; Sao Jorge; Terceira))", E_INT_24, -104, 167, -38 }, // 72
179 { "Guam 1963 (Guam)", E_CLARKE_66, -100, -248, 259 }, // 73
180 { "Gunung Segara (Indonesia (Kalimantan))", E_BESS_41, -403, 684, 41 }, // 74
181 { "GUX 1 Astro (Guadalcanal Island)", E_INT_24, 252, -209, -751 }, // 75
182 { "Herat North (Afghanistan)", E_INT_24, -333, -222, 114 }, // 76
183 { "Hermannskogel Datum (Croatia -Serbia, Bosnia-Herzegovina)", E_BESS_41_NAM, 653, -212, 449 }, // 77
184 { "Hjorsey 1955 (Iceland)", E_INT_24, -73, 46, -86 }, // 78
185 { "Hong Kong 1963 (Hong Kong)", E_INT_24, -156, -271, -189 }, // 79
186 { "Hu-Tzu-Shan (Taiwan)", E_INT_24, -637, -549, -203 }, // 80
187 { "Indian (Bangladesh)", E_EVR_IND_30, 282, 726, 254 }, // 81
188 { "Indian (India; Nepal)", E_EVR_IND_56, 295, 736, 257 }, // 82
189 { "Indian (Pakistan)", E_EVR_PAK, 283, 682, 231 }, // 83
190 { "Indian 1954 (Thailand)", E_EVR_IND_30, 217, 823, 299 }, // 84
191 { "Indian 1960 (Vietnam (Con Son Island))", E_EVR_IND_30, 182, 915, 344 }, // 85
192 { "Indian 1960 (Vietnam (Near 16N))", E_EVR_IND_30, 198, 881, 317 }, // 86
193 { "Indian 1975 (Thailand)", E_EVR_IND_30, 210, 814, 289 }, // 87
194 { "Indonesian 1974 (Indonesia)", E_IND_74, -24, -15, 5 }, // 88
195 { "Ireland 1965 (Ireland)", E_MOD_AIRY, 506, -122, 611 }, // 89
196 { "ISTS 061 Astro 1968 (South Georgia Islands)", E_INT_24, -794, 119, -298 }, // 90
197 { "ISTS 073 Astro 1969 (Diego Garcia)", E_INT_24, 208, -435, -229 }, // 91
198 { "Johnston Island 1961 (Johnston Island)", E_INT_24, 189, -79, -202 }, // 92
199 { "Kandawala (Sri Lanka)", E_EVR_IND_30, -97, 787, 86 }, // 93
200 { "Kerguelen Island 1949 (Kerguelen Island)", E_INT_24, 145, -187, 103 }, // 94
201 { "Kertau 1948 (West Malaysia & Singapore)", E_EVR_MAL_SING, -11, 851, 5 }, // 95
202 { "Kusaie Astro 1951 (Caroline Islands)", E_INT_24, 647, 1777, -1124 }, // 96
203 { "Korean Geodetic System (South Korea)", E_GRS_80, 0, 0, 0 }, // 97
204 { "L. C. 5 Astro 1961 (Cayman Brac Island)", E_CLARKE_66, 42, 124, 147 }, // 98
205 { "Leigon (Ghana)", E_CLARKE_80, -130, 29, 364 }, // 99
206 { "Liberia 1964 (Liberia)", E_CLARKE_80, -90, 40, 88 }, // 100
207 { "Luzon (Philippines (Excluding Mindanao))", E_CLARKE_66, -133, -77, -51 }, // 101
208 { "Luzon (Philippines (Mindanao))", E_CLARKE_66, -133, -79, -72 }, // 102
209 { "M'Poraloko (Gabon)", E_CLARKE_80, -74, -130, 42 }, // 103
210 { "Mahe 1971 (Mahe Island)", E_CLARKE_80, 41, -220, -134 }, // 104
211 { "Massawa (Ethiopia (Eritrea))", E_BESS_41, 639, 405, 60 }, // 105
212 { "Merchich (Morocco)", E_CLARKE_80, 31, 146, 47 }, // 106
213 { "Midway Astro 1961 (Midway Islands)", E_INT_24, 912, -58, 1227 }, // 107
214 { "Minna (Cameroon)", E_CLARKE_80, -81, -84, 115 }, // 108
215 { "Minna (Nigeria)", E_CLARKE_80, -92, -93, 122 }, // 109
216 { "Montserrat Island Astro 1958 (Montserrat (Leeward Islands))", E_CLARKE_80, 174, 359, 365 }, // 110
217 { "Nahrwan (Oman (Masirah Island))", E_CLARKE_80, -247, -148, 369 }, // 111
218 { "Nahrwan (Saudi Arabia)", E_CLARKE_80, -243, -192, 477 }, // 112
219 { "Nahrwan (United Arab Emirates)", E_CLARKE_80, -249, -156, 381 }, // 113
220 { "Naparima BWI (Trinidad & Tobago)", E_INT_24, -10, 375, 165 }, // 114
221 { "North American 1927 (Alaska (Excluding Aleutian Ids))", E_CLARKE_66, -5, 135, 172 }, // 115
222 { "North American 1927 (Alaska (Aleutian Ids East of 180W))", E_CLARKE_66, -2, 152, 149 }, // 116
223 { "North American 1927 (Alaska (Aleutian Ids West of 180W))", E_CLARKE_66, 2, 204, 105 }, // 117
224 { "North American 1927 (Bahamas (Except San Salvador Id))", E_CLARKE_66, -4, 154, 178 }, // 118
225 { "North American 1927 (Bahamas (San Salvador Island))", E_CLARKE_66, 1, 140, 165 }, // 119
226 { "North American 1927 (Canada (Alberta; British Columbia))", E_CLARKE_66, -7, 162, 188 }, // 120
227 { "North American 1927 (Canada (Manitoba; Ontario))", E_CLARKE_66, -9, 157, 184 }, // 121
228 { "North American 1927 (Canada (New Brunswick; Newfoundland; Nova Scotia; Quebec))", E_CLARKE_66, -22, 160, 190 }, // 122
229 { "North American 1927 (Canada (Northwest Territories; Saskatchewan))", E_CLARKE_66, 4, 159, 188 }, // 123
230 { "North American 1927 (Canada (Yukon))", E_CLARKE_66, -7, 139, 181 }, // 124
231 { "North American 1927 (Canal Zone)", E_CLARKE_66, 0, 125, 201 }, // 125
232 { "North American 1927 (Cuba)", E_CLARKE_66, -9, 152, 178 }, // 126
233 { "North American 1927 (Greenland (Hayes Peninsula))", E_CLARKE_66, 11, 114, 195 }, // 127
234 { "North American 1927 (MEAN FOR Antigua; Barbados; Barbuda; Caicos Islands; Cuba; Dominican Republic; Grand Cayman; Jamaica; Turks Islands)", E_CLARKE_66, -3, 142, 183 }, // 128
235 { "North American 1927 (MEAN FOR Belize; Costa Rica; El Salvador; Guatemala; Honduras; Nicaragua)", E_CLARKE_66, 0, 125, 194 }, // 129
236 { "North American 1927 (MEAN FOR Canada)", E_CLARKE_66, -10, 158, 187 }, // 130
237 { "North American 1927 (MEAN FOR CONUS)", E_CLARKE_66, -8, 160, 176 }, // 131
238 { "North American 1927 (MEAN FOR CONUS (East of Mississippi; River Including Louisiana; Missouri; Minnesota))", E_CLARKE_66, -9, 161, 179 }, // 132
239 { "North American 1927 (MEAN FOR CONUS (West of Mississippi; River Excluding Louisiana; Minnesota; Missouri))", E_CLARKE_66, -8, 159, 175 }, // 133
240 { "North American 1927 (Mexico)", E_CLARKE_66, -12, 130, 190 }, // 134
241 { "North American 1983 (Alaska (Excluding Aleutian Ids))", E_GRS_80, 0, 0, 0 }, // 135
242 { "North American 1983 (Aleutian Ids)", E_GRS_80, -2, 0, 4 }, // 136
243 { "North American 1983 (Canada)", E_GRS_80, 0, 0, 0 }, // 137
244 { "North American 1983 (CONUS)", E_GRS_80, 0, 0, 0 }, // 138
245 { "North American 1983 (Hawaii)", E_GRS_80, 1, 1, -1 }, // 139
246 { "North American 1983 (Mexico; Central America)", E_GRS_80, 0, 0, 0 }, // 140
247 { "North Sahara 1959 (Algeria)", E_CLARKE_80, -186, -93, 310 }, // 141
248 { "Observatorio Meteorologico 1939 (Azores (Corvo & Flores Islands))", E_INT_24, -425, -169, 81 }, // 142
249 { "Old Egyptian 1907 (Egypt)", E_HELM_06, -130, 110, -13 }, // 143
250 { "Old Hawaiian (Hawaii)", E_CLARKE_66, 89, -279, -183 }, // 144
251 { "Old Hawaiian (Kauai)", E_CLARKE_66, 45, -290, -172 }, // 145
252 { "Old Hawaiian (Maui)", E_CLARKE_66, 65, -290, -190 }, // 146
253 { "Old Hawaiian (MEAN FOR Hawaii; Kauai; Maui; Oahu)", E_CLARKE_66, 61, -285, -181 }, // 147
254 { "Old Hawaiian (Oahu)", E_CLARKE_66, 58, -283, -182 }, // 148
255 { "Oman (Oman)", E_CLARKE_80, -346, -1, 224 }, // 149
256 { "Ordnance Survey Great Britain 1936 (England)", E_AIRY_30, 371, -112, 434 }, // 150
257 { "Ordnance Survey Great Britain 1936 (England; Isle of Man; Wales)", E_AIRY_30, 371, -111, 434 }, // 151
258 { "Ordnance Survey Great Britain 1936 (MEAN FOR England; Isle of Man; Scotland; Shetland Islands; Wales)", E_AIRY_30, 375, -111, 431 }, // 152
259 { "Ordnance Survey Great Britain 1936 (Scotland; Shetland Islands)", E_AIRY_30, 384, -111, 425 }, // 153
260 { "Ordnance Survey Great Britain 1936 (Wales)", E_AIRY_30, 370, -108, 434 }, // 154
261 { "Pico de las Nieves (Canary Islands)", E_INT_24, -307, -92, 127 }, // 155
262 { "Pitcairn Astro 1967 (Pitcairn Island)", E_INT_24, 185, 165, 42 }, // 156
263 { "Point 58 (MEAN FOR Burkina Faso & Niger)", E_CLARKE_80, -106, -129, 165 }, // 157
264 { "Pointe Noire 1948 (Congo)", E_CLARKE_80, -148, 51, -291 }, // 158
265 { "Porto Santo 1936 (Porto Santo; Madeira Islands)", E_INT_24, -499, -249, 314 }, // 159
266 { "Provisional South American 1956 (Bolivia)", E_INT_24, -270, 188, -388 }, // 160
267 { "Provisional South American 1956 (Chile (Northern; Near 19S))", E_INT_24, -270, 183, -390 }, // 161
268 { "Provisional South American 1956 (Chile (Southern; Near 43S))", E_INT_24, -305, 243, -442 }, // 162
269 { "Provisional South American 1956 (Colombia)", E_INT_24, -282, 169, -371 }, // 163
270 { "Provisional South American 1956 (Ecuador)", E_INT_24, -278, 171, -367 }, // 164
271 { "Provisional South American 1956 (Guyana)", E_INT_24, -298, 159, -369 }, // 165
272 { "Provisional South American 1956 (MEAN FOR Bolivia; Chile; Colombia; Ecuador; Guyana; Peru; Venezuela)", E_INT_24, -288, 175, -376 }, // 166
273 { "Provisional South American 1956 (Peru)", E_INT_24, -279, 175, -379 }, // 167
274 { "Provisional South American 1956 (Venezuela)", E_INT_24, -295, 173, -371 }, // 168
275 { "Provisional South Chilean 1963 (Chile (Near 53S) (Hito XVIII))", E_INT_24, 16, 196, 93 }, // 169
276 { "Puerto Rico (Puerto Rico; Virgin Islands)", E_CLARKE_66, 11, 72, -101 }, // 170
277 { "Pulkovo 1942 (Russia)", E_KRASS_40, 28, -130, -95 }, // 171
278 { "Qatar National (Qatar)", E_INT_24, -128, -283, 22 }, // 172
279 { "Qornoq (Greenland (South))", E_INT_24, 164, 138, -189 }, // 173
280 { "Reunion (Mascarene Islands)", E_INT_24, 94, -948, -1262 }, // 174
281 { "Rome 1940 (Italy (Sardinia))", E_INT_24, -225, -65, 9 }, // 175
282 { "S-42 (Pulkovo 1942) (Hungary)", E_KRASS_40, 28, -121, -77 }, // 176
283 { "S-42 (Pulkovo 1942) (Poland)", E_KRASS_40, 23, -124, -82 }, // 177
284 { "S-42 (Pulkovo 1942) (Czechoslavakia)", E_KRASS_40, 26, -121, -78 }, // 178
285 { "S-42 (Pulkovo 1942) (Latvia)", E_KRASS_40, 24, -124, -82 }, // 179
286 { "S-42 (Pulkovo 1942) (Kazakhstan)", E_KRASS_40, 15, -130, -84 }, // 180
287 { "S-42 (Pulkovo 1942) (Albania)", E_KRASS_40, 24, -130, -92 }, // 181
288 { "S-42 (Pulkovo 1942) (Romania)", E_KRASS_40, 28, -121, -77 }, // 182
289 { "S-JTSK (Czechoslavakia (Prior 1 JAN 1993))", E_BESS_41, 589, 76, 480 }, // 183
290 { "Santo (DOS) 1965 (Espirito Santo Island)", E_INT_24, 170, 42, 84 }, // 184
291 { "Sao Braz (Azores (Sao Miguel; Santa Maria Ids))", E_INT_24, -203, 141, 53 }, // 185
292 { "Sapper Hill 1943 (East Falkland Island)", E_INT_24, -355, 21, 72 }, // 186
293 { "Schwarzeck (Namibia)", E_BESS_41_NAM, 616, 97, -251 }, // 187
294 { "Selvagem Grande 1938 (Salvage Islands)", E_INT_24, -289, -124, 60 }, // 188
295 { "Sierra Leone 1960 (Sierra Leone)", E_CLARKE_80, -88, 4, 101 }, // 189
296 { "South American 1969 (Argentina)", E_S_AMER_69, -62, -1, -37 }, // 190
297 { "South American 1969 (Bolivia)", E_S_AMER_69, -61, 2, -48 }, // 191
298 { "South American 1969 (Brazil)", E_S_AMER_69, -60, -2, -41 }, // 192
299 { "South American 1969 (Chile)", E_S_AMER_69, -75, -1, -44 }, // 193
300 { "South American 1969 (Colombia)", E_S_AMER_69, -44, 6, -36 }, // 194
301 { "South American 1969 (Ecuador)", E_S_AMER_69, -48, 3, -44 }, // 195
302 { "South American 1969 (Ecuador (Baltra; Galapagos))", E_S_AMER_69, -47, 26, -42 }, // 196
303 { "South American 1969 (Guyana)", E_S_AMER_69, -53, 3, -47 }, // 197
304 { "South American 1969 (MEAN FOR Argentina; Bolivia; Brazil; Chile; Colombia; Ecuador; Guyana; Paraguay; Peru; Trinidad & Tobago; Venezuela)", E_S_AMER_69, -57, 1, -41 }, // 198
305 { "South American 1969 (Paraguay)", E_S_AMER_69, -61, 2, -33 }, // 199
306 { "South American 1969 (Peru)", E_S_AMER_69, -58, 0, -44 }, // 200
307 { "South American 1969 (Trinidad & Tobago)", E_S_AMER_69, -45, 12, -33 }, // 201
308 { "South American 1969 (Venezuela)", E_S_AMER_69, -45, 8, -33 }, // 202
309 { "South Asia (Singapore)", E_MOD_FISCH_60, 7, -10, -26 }, // 203
310 { "Tananarive Observatory 1925 (Madagascar)", E_INT_24, -189, -242, -91 }, // 204
311 { "Timbalai 1948 (Brunei; E. Malaysia (Sabah Sarawak))", E_EVR_SAB_SAR, -679, 669, -48 }, // 205
312 { "Tokyo (Japan)", E_BESS_41, -148, 507, 685 }, // 206
313 { "Tokyo (MEAN FOR Japan; South Korea; Okinawa)", E_BESS_41, -148, 507, 685 }, // 207
314 { "Tokyo (Okinawa)", E_BESS_41, -158, 507, 676 }, // 208
315 { "Tokyo (South Korea)", E_BESS_41, -147, 506, 687 }, // 209
316 { "Tristan Astro 1968 (Tristan da Cunha)", E_INT_24, -632, 438, -609 }, // 210
317 { "Viti Levu 1916 (Fiji (Viti Levu Island))", E_CLARKE_80, 51, 391, -36 }, // 211
318 { "Voirol 1960 (Algeria)", E_CLARKE_80, -123, -206, 219 }, // 212
319 { "Wake Island Astro 1952 (Wake Atoll)", E_INT_24, 276, -57, 149 }, // 213
320 { "Wake-Eniwetok 1960 (Marshall Islands)", E_HOUGH_60, 102, 52, -38 }, // 214
321 { "WGS 1972 (Global Definition)", E_WGS_72, 0, 0, 0 }, // 215
322 { "WGS 1984 (Global Definition)", E_WGS_84, 0, 0, 0 }, // 216
323 { "Yacare (Uruguay)", E_INT_24, -155, 171, 37 }, // 217
324 { "Zanderij (Suriname)", E_INT_24, -265, 120, -358 } // 218
325 };
326
327
328
329 static const double PI = 3.14159265358979323846;
330
331
332
333
334
335 /* As you can see this little function is just a 2 step datum shift, going through WGS84. */
datum_shift(double * latitude,double * longitude,short fromDatumID,short toDatumID)336 void datum_shift(double *latitude, double *longitude, short fromDatumID, short toDatumID)
337 {
338 wgs84_datum_shift(TO_WGS_84, latitude, longitude, fromDatumID);
339 wgs84_datum_shift(FROM_WGS_84, latitude, longitude, toDatumID);
340 }
341
342
343
344
345
346 /*
347 Function to convert latitude and longitude in decimal degrees from WGS84 to
348 another datum or from another datum to WGS84. The arguments to this function
349 include a direction flag 'fromWGS84', pointers to double precision latitude
350 and longitude, and an index to the gDatum[] array.
351 */
wgs84_datum_shift(short fromWGS84,double * latitude,double * longitude,short datumID)352 void wgs84_datum_shift(short fromWGS84, double *latitude, double *longitude, short datumID)
353 {
354 double dx = gDatum[datumID].dx;
355 double dy = gDatum[datumID].dy;
356 double dz = gDatum[datumID].dz;
357
358 double phi = *latitude * PI / 180.0;
359 double lambda = *longitude * PI / 180.0;
360 double a0, b0, es0, f0; /* Reference ellipsoid of input data */
361 // a1 and b1 are never actually used, so don't declare them and set
362 // them (gcc warns about set-but-unused vars)
363 //double a1, b1, es1, f1; /* Reference ellipsoid of output data */
364 double es1, f1; /* Reference ellipsoid of output data */
365 double psi; /* geocentric latitude */
366 double x, y, z; /* 3D coordinates with respect to original datum */
367 double psi1; /* transformed geocentric latitude */
368
369 if (datumID == D_WGS_84) // do nothing if current datum is WGS84
370 {
371 return;
372 }
373
374 if (fromWGS84) /* convert from WGS84 to new datum */
375 {
376 a0 = gEllipsoid[E_WGS_84].a; /* WGS84 semimajor axis */
377 f0 = 1.0 / gEllipsoid[E_WGS_84].invf; /* WGS84 flattening */
378 // a1 is never used except to set b1, which itself is never used
379 // a1 = gEllipsoid[gDatum[datumID].ellipsoid].a;
380 f1 = 1.0 / gEllipsoid[gDatum[datumID].ellipsoid].invf;
381 }
382 else /* convert from datum to WGS84 */
383 {
384 a0 = gEllipsoid[gDatum[datumID].ellipsoid].a; /* semimajor axis */
385 f0 = 1.0 / gEllipsoid[gDatum[datumID].ellipsoid].invf; /* flattening */
386 // a1 is never used except to set b1, which is never used.
387 // a1 = gEllipsoid[E_WGS_84].a; /* WGS84 semimajor axis */
388 f1 = 1.0 / gEllipsoid[E_WGS_84].invf; /* WGS84 flattening */
389 dx = -dx;
390 dy = -dy;
391 dz = -dz;
392 }
393
394 b0 = a0 * (1 - f0); /* semiminor axis for input datum */
395 es0 = 2 * f0 - f0*f0; /* eccentricity^2 */
396
397 // b1 is never used
398 // b1 = a1 * (1 - f1); /* semiminor axis for output datum */
399 es1 = 2 * f1 - f1*f1; /* eccentricity^2 */
400
401 /* Convert geodedic latitude to geocentric latitude, psi */
402 if (*latitude == 0.0 || *latitude == 90.0 || *latitude == -90.0)
403 {
404 psi = phi;
405 }
406 else
407 {
408 psi = atan((1 - es0) * tan(phi));
409 }
410
411 /* Calculate x and y axis coordinates with respect to the original ellipsoid */
412 if (*longitude == 90.0 || *longitude == -90.0)
413 {
414 x = 0.0;
415 y = fabs(a0 * b0 / sqrt(b0*b0 + a0*a0*pow(tan(psi), 2.0)));
416 }
417 else
418 {
419 x = fabs((a0 * b0) /
420 sqrt((1 + pow(tan(lambda), 2.0)) * (b0*b0 + a0*a0 * pow(tan(psi), 2.0))));
421 y = fabs(x * tan(lambda));
422 }
423
424 if (*longitude < -90.0 || *longitude > 90.0)
425 {
426 x = -x;
427 }
428 if (*longitude < 0.0)
429 {
430 y = -y;
431 }
432
433 /* Calculate z axis coordinate with respect to the original ellipsoid */
434 if (*latitude == 90.0)
435 {
436 z = b0;
437 }
438 else if (*latitude == -90.0)
439 {
440 z = -b0;
441 }
442 else
443 {
444 z = tan(psi) * sqrt((a0*a0 * b0*b0) / (b0*b0 + a0*a0 * pow(tan(psi), 2.0)));
445 }
446
447 /* Calculate the geocentric latitude with respect to the new ellipsoid */
448 psi1 = atan((z - dz) / sqrt((x - dx)*(x - dx) + (y - dy)*(y - dy)));
449
450 /* Convert to geocentric latitude and save return value */
451 *latitude = atan(tan(psi1) / (1 - es1)) * 180.0 / PI;
452
453 /* Calculate the longitude with respect to the new ellipsoid */
454 *longitude = atan((y - dy) / (x - dx)) * 180.0 / PI;
455
456 /* Correct the resultant for negative x values */
457 if (x-dx < 0.0)
458 {
459 if (y-dy > 0.0)
460 {
461 *longitude = 180.0 + *longitude;
462 }
463 else
464 {
465 *longitude = -180.0 + *longitude;
466 }
467 }
468 }
469
470
471
472
473
474 #define deg2rad (PI / 180)
475 #define rad2deg (180.0 / PI)
476
477
478
479
480
481 /*
482 Source
483 Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
484 1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency
485 */
486 //
487 // Convert lat/long to UTM/UPS coordinates
ll_to_utm_ups(short ellipsoidID,const double lat,const double lon,double * utmNorthing,double * utmEasting,char * utmZone,int utmZoneLength)488 void ll_to_utm_ups(short ellipsoidID, const double lat, const double lon,
489 double *utmNorthing, double *utmEasting, char* utmZone, int utmZoneLength)
490 {
491 //converts lat/long to UTM coords. Equations from USGS Bulletin 1532
492 //East Longitudes are positive, West longitudes are negative.
493 //North latitudes are positive, South latitudes are negative
494 //Lat and Long are in decimal degrees
495 //Written by Chuck Gantz- chuck.gantz@globalstar.com
496
497 double a = gEllipsoid[ellipsoidID].a;
498 double f = 1.0 / gEllipsoid[ellipsoidID].invf;
499 double eccSquared = (2 * f) - (f * f);
500 double k0 = 0.9996;
501
502 double LongOrigin;
503 double eccPrimeSquared;
504 double N, T, C, A, M;
505
506 //Make sure the longitude is between -180.00 .. 179.9
507 double LongTemp = (lon+180)-(int)((lon+180)/360)*360-180; // -180.00 .. 179.9;
508
509 double LatRad = lat*deg2rad;
510 double LongRad = LongTemp*deg2rad;
511 double LongOriginRad;
512 int ZoneNumber;
513
514 ZoneNumber = (int)((LongTemp + 180)/6) + 1;
515
516 if (coordinate_system == USE_UTM_SPECIAL
517 || coordinate_system == USE_MGRS)
518 {
519
520 // Special zone for southern Norway. Used for military
521 // version of UTM (MGRS) only.
522 if ( lat >= 56.0 && lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
523 {
524 ZoneNumber = 32;
525 }
526
527 // Handle the special zones for Svalbard. Used for military
528 // version of UTM (MGRS) only.
529 if (lat >= 72.0 && lat < 84.0)
530 {
531 if (LongTemp >= 0.0 && LongTemp < 9.0)
532 {
533 ZoneNumber = 31;
534 }
535 else if (LongTemp >= 9.0 && LongTemp < 21.0)
536 {
537 ZoneNumber = 33;
538 }
539 else if (LongTemp >= 21.0 && LongTemp < 33.0)
540 {
541 ZoneNumber = 35;
542 }
543 else if (LongTemp >= 33.0 && LongTemp < 42.0)
544 {
545 ZoneNumber = 37;
546 }
547 }
548 }
549
550 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
551 LongOriginRad = LongOrigin * deg2rad;
552
553 if (lat > 84.0 || lat < -80.0)
554 {
555 // We're in the UPS areas (near the poles). ZoneNumber
556 // should not be printed in this case.
557 xastir_snprintf(utmZone,
558 utmZoneLength,
559 "%c",
560 utm_letter_designator(lat, lon));
561 }
562 else // We're in the UTM areas (not near the poles).
563 {
564 //compute the UTM Zone from the latitude and longitude
565 xastir_snprintf(utmZone,
566 utmZoneLength,
567 "%d%c",
568 ZoneNumber,
569 utm_letter_designator(lat, lon));
570 }
571
572 eccPrimeSquared = (eccSquared)/(1-eccSquared);
573
574 if (lat > 84.0 || lat < -80.0)
575 {
576 //
577 // We're dealing with UPS coordinates (near the poles)
578 //
579 // The following piece of code which implements UPS
580 // conversion is derived from code that John Waers
581 // <jfwaers@csn.net> placed in the public domain. It's from
582 // his program "MacGPS45".
583
584 double t, e, rho;
585 const double k0 = 0.994;
586 double lambda = lon * (PI/180.0);
587 double phi = fabs(lat * (PI/180.0) );
588
589 e = sqrt(eccSquared);
590 t = tan(PI/4.0 - phi/2.0) / pow( (1.0 - e * sin(phi)) / (1.0 + e * sin(phi)), (e/2.0) );
591 rho = 2.0 * a * k0 * t / sqrt(pow(1.0+e, 1.0+e) * pow(1.0-e, 1.0-e));
592 *utmEasting = rho * sin(lambda);
593 *utmNorthing = rho * cos(lambda);
594
595 if (lat > 0.0) // Northern hemisphere
596 {
597 *utmNorthing = -(*utmNorthing);
598 }
599
600 *utmEasting += 2.0e6; // Add in false easting and northing
601 *utmNorthing += 2.0e6;
602 }
603 else
604 {
605 //
606 // We're dealing with UTM coordinates
607 //
608 N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
609 T = tan(LatRad)*tan(LatRad);
610 C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
611 A = cos(LatRad)*(LongRad-LongOriginRad);
612
613 M = a*((1 -
614 eccSquared/4 -
615 3*eccSquared*eccSquared/64 -
616 5*eccSquared*eccSquared*eccSquared/256) * LatRad -
617 (3*eccSquared/8 +
618 3*eccSquared*eccSquared/32 +
619 45*eccSquared*eccSquared*eccSquared/1024) * sin(2*LatRad) +
620 (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024) * sin(4*LatRad) -
621 (35*eccSquared*eccSquared*eccSquared/3072) * sin(6*LatRad));
622
623 *utmEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
624 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
625 + 500000.0);
626
627 *utmNorthing = (double)(k0*(M+N*tan(LatRad)*
628 (A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
629 + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
630
631 if (lat < 0)
632 {
633 *utmNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
634 }
635 }
636 }
637
638
639
640
641
642 // Handles UPS/UTM coordinates equally well!
643 //
utm_letter_designator(double lat,double lon)644 char utm_letter_designator(double lat, double lon)
645 {
646 // This routine determines the correct UTM/UPS letter designator
647 // for the given latitude. Originally written by Chuck Gantz-
648 // chuck.gantz@globalstar.com
649 // Modified to handle UPS zones. --we7u
650 char LetterDesignator;
651
652 if ((84 >= lat) && (lat >= 72))
653 {
654 LetterDesignator = 'X';
655 }
656 else if ((72 > lat) && (lat >= 64))
657 {
658 LetterDesignator = 'W';
659 }
660 else if ((64 > lat) && (lat >= 56))
661 {
662 LetterDesignator = 'V';
663 }
664 else if ((56 > lat) && (lat >= 48))
665 {
666 LetterDesignator = 'U';
667 }
668 else if ((48 > lat) && (lat >= 40))
669 {
670 LetterDesignator = 'T';
671 }
672 else if ((40 > lat) && (lat >= 32))
673 {
674 LetterDesignator = 'S';
675 }
676 else if ((32 > lat) && (lat >= 24))
677 {
678 LetterDesignator = 'R';
679 }
680 else if ((24 > lat) && (lat >= 16))
681 {
682 LetterDesignator = 'Q';
683 }
684 else if ((16 > lat) && (lat >= 8))
685 {
686 LetterDesignator = 'P';
687 }
688 else if (( 8 > lat) && (lat >= 0))
689 {
690 LetterDesignator = 'N';
691 }
692 else if (( 0 > lat) && (lat >= -8))
693 {
694 LetterDesignator = 'M';
695 }
696 else if ((-8 > lat) && (lat >= -16))
697 {
698 LetterDesignator = 'L';
699 }
700 else if ((-16 > lat) && (lat >= -24))
701 {
702 LetterDesignator = 'K';
703 }
704 else if ((-24 > lat) && (lat >= -32))
705 {
706 LetterDesignator = 'J';
707 }
708 else if ((-32 > lat) && (lat >= -40))
709 {
710 LetterDesignator = 'H';
711 }
712 else if ((-40 > lat) && (lat >= -48))
713 {
714 LetterDesignator = 'G';
715 }
716 else if ((-48 > lat) && (lat >= -56))
717 {
718 LetterDesignator = 'F';
719 }
720 else if ((-56 > lat) && (lat >= -64))
721 {
722 LetterDesignator = 'E';
723 }
724 else if ((-64 > lat) && (lat >= -72))
725 {
726 LetterDesignator = 'D';
727 }
728 else if ((-72 > lat) && (lat >= -80))
729 {
730 LetterDesignator = 'C';
731 }
732 else
733 {
734 //
735 // We're dealing with UPS (N/S Pole) coordinates, not UTM
736 //
737 if (lat > 84) // North Pole, Y/Z zones
738 {
739 if ((0 <= lon) && (lon <= 180))
740 {
741 LetterDesignator = 'Z'; // E or + longitude
742 }
743 else
744 {
745 LetterDesignator = 'Y'; // W or - longitude
746 }
747 }
748 else // Lat < 80S, South Pole, A/B zones
749 {
750 if ((0 <= lon) && (lon <= 180))
751 {
752 LetterDesignator = 'B'; // E or + longitude
753 }
754 else
755 {
756 LetterDesignator = 'A'; // W or - longitude
757 }
758 }
759 }
760 return LetterDesignator;
761 }
762
763
764
765
766
767 // The following piece of code which implements UPS conversion is
768 // derived from code that John Waers <jfwaers@csn.net> placed in
769 // the public domain. It's from his program "MacGPS45".
770 //
calcPhi(double * phi,double e,double t)771 static void calcPhi(double *phi, double e, double t)
772 {
773 double old = PI/2.0 - 2.0 * atan(t);
774 short maxIterations = 20;
775
776 while ( (fabs((*phi - old) / *phi) > 1.0e-8) && maxIterations-- )
777 {
778 old = *phi;
779 *phi = PI/ 2.0 - 2.0 * atan( t * pow((1.0 - e * sin(*phi)) / ((1.0 + e * sin(*phi))), (e / 2.0)) );
780 }
781 }
782
783
784
785
786
787 // Converts from UTM/UPS coordinates to Lat/Long coordinates.
788 //
utm_ups_to_ll(short ellipsoidID,const double utmNorthing,const double utmEasting,const char * utmZone,double * lat,double * lon)789 void utm_ups_to_ll(short ellipsoidID, const double utmNorthing, const double utmEasting,
790 const char* utmZone, double *lat, double *lon)
791 {
792 // Converts UTM coords to lat/long. Equations from USGS
793 // Bulletin 1532. East Longitudes are positive, West longitudes
794 // are negative. North latitudes are positive, South latitudes
795 // are negative Lat and Long are in decimal degrees.
796 // Written by Chuck Gantz- chuck.gantz@globalstar.com
797 // Modified by WE7U to add UPS support.
798
799 double k0 = 0.9996;
800 double a = gEllipsoid[ellipsoidID].a;
801 double f = 1.0 / gEllipsoid[ellipsoidID].invf;
802 double eccSquared = (2 * f) - (f * f);
803 double eccPrimeSquared;
804 double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
805 double N1, T1, C1, R1, D, M;
806 double LongOrigin;
807 // phi1 is never used, but is set. Don't make gcc warn us
808 // double mu, phi1, phi1Rad;
809 double mu, phi1Rad;
810 double x, y;
811 int ZoneNumber;
812 char* ZoneLetter;
813 // Unused variable
814 // int NorthernHemisphere; // 1=northern hemisphere, 0=southern
815
816
817 //fprintf(stderr,"%s %f %f\n",
818 // utmZone,
819 // utmEasting,
820 // utmNorthing);
821
822 x = utmEasting;
823 y = utmNorthing;
824
825 ZoneNumber = strtoul(utmZone, &ZoneLetter, 10);
826
827 // Remove any possible leading spaces
828 remove_leading_spaces(ZoneLetter);
829
830 // Make sure the zone letter is upper-case
831 *ZoneLetter = toupper(*ZoneLetter);
832
833 //fprintf(stderr,"ZoneLetter: %s\n", ZoneLetter);
834
835 if ( *ZoneLetter == 'Y' // North Pole
836 || *ZoneLetter == 'Z' // North Pole
837 || *ZoneLetter == 'A' // South Pole
838 || *ZoneLetter == 'B') // South Pole
839 {
840
841 // The following piece of code which implements UPS
842 // conversion is derived from code that John Waers
843 // <jfwaers@csn.net> placed in the public domain. It's from
844 // his program "MacGPS45".
845
846 //
847 // We're dealing with a UPS coordinate (near the poles)
848 // instead of a UTM coordinate. We need to do entirely
849 // different calculations for UPS.
850 //
851 double e, t, rho;
852 const double k0 = 0.994;
853
854
855 //fprintf(stderr,"UPS Coordinates\n");
856
857
858 e = sqrt(eccSquared);
859
860 x -= 2.0e6; // Remove false easting and northing
861 y -= 2.0e6;
862
863 rho = sqrt(x*x + y*y);
864 t = rho * sqrt(pow(1.0+e, 1.0+e) * pow(1.0-e, 1.0-e)) / (2.0 * a * k0);
865
866 calcPhi(lat, e, t);
867
868 *lat /= (PI/180.0);
869
870 // This appears to be necessary in order to get proper
871 // positions in the south polar region
872 if (*ZoneLetter == 'A' || *ZoneLetter == 'B')
873 {
874 *lat = -*lat;
875 }
876
877 if (y != 0.0)
878 {
879 t = atan(fabs(x/y));
880 }
881 else
882 {
883 t = PI / 2.0;
884 if (x < 0.0)
885 {
886 t = -t;
887 }
888 }
889
890 if (*ZoneLetter == 'Z' || *ZoneLetter == 'Y')
891 {
892 y = -y; // Northern hemisphere
893 }
894
895 if (y < 0.0)
896 {
897 t = PI - t;
898 }
899
900 if (x < 0.0)
901 {
902 t = -t;
903 }
904
905 *lon = t / (PI/180.0);
906
907 /*
908 fprintf(stderr,"datum.c:utm_ups_to_ll(): Found UPS Coordinate: %s %f %f\n",
909 utmZone,
910 utmEasting,
911 utmNorthing);
912 */
913 return; // Done computing UPS coordinates
914 }
915
916
917 // If we make it here, we're working on UTM coordinates (not
918 // UPS coordinates).
919
920
921 x = utmEasting - 500000.0; //remove 500,000 meter offset for longitude
922 y = utmNorthing;
923
924
925 if ((*ZoneLetter - 'N') >= 0)
926 {
927 // We never use this variable
928 // NorthernHemisphere = 1;//point is in northern hemisphere
929 }
930 else
931 {
932 // we never use NorthernHemisphere
933 // NorthernHemisphere = 0;//point is in southern hemisphere
934 y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
935 }
936
937 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
938
939 eccPrimeSquared = (eccSquared)/(1-eccSquared);
940
941 M = y / k0;
942 mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
943
944 phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
945 + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
946 + (151*e1*e1*e1/96)*sin(6*mu);
947 // This variable is never used, it is just phi1Rad converted to degrees
948 // phi1 = phi1Rad*rad2deg;
949
950 N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
951 T1 = tan(phi1Rad)*tan(phi1Rad);
952 C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
953 R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
954 D = x/(N1*k0);
955
956 *lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
957 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
958 *lat *= rad2deg;
959
960 *lon = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
961 *D*D*D*D*D/120)/cos(phi1Rad);
962 *lon = LongOrigin + (*lon) * rad2deg;
963 }
964
965
966