108 {
109
110
111
112
113
114
115 const double PI = 3.14159265;
116
117 const double deg2rad = PI / 180;
118
119
122 double k0 = 0.9996;
123
124 double LongOrigin;
125 double eccPrimeSquared;
126 double N, T, C, A, M;
127
128
129 double LongTemp = (Long + 180) - int((Long + 180) / 360)*360 - 180;
130
131 double LatRad = Lat*deg2rad;
132 double LongRad = LongTemp*deg2rad;
133 double LongOriginRad;
134 int ZoneNumber;
135
136 ZoneNumber = int((LongTemp + 180) / 6) + 1;
137
138 if (Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0)
139 ZoneNumber = 32;
140
141
142 if (Lat >= 72.0 && Lat < 84.0) {
143 if (LongTemp >= 0.0 && LongTemp < 9.0) ZoneNumber = 31;
144 else if (LongTemp >= 9.0 && LongTemp < 21.0) ZoneNumber = 33;
145 else if (LongTemp >= 21.0 && LongTemp < 33.0) ZoneNumber = 35;
146 else if (LongTemp >= 33.0 && LongTemp < 42.0) ZoneNumber = 37;
147 }
148 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;
149 LongOriginRad = LongOrigin * deg2rad;
150
151
153
154 eccPrimeSquared = (eccSquared) / (1 - eccSquared);
155
156 N =
a / sqrt(1 - eccSquared * sin(LatRad) * sin(LatRad));
157 T = tan(LatRad) * tan(LatRad);
158 C = eccPrimeSquared * cos(LatRad) * cos(LatRad);
159 A = cos(LatRad)*(LongRad - LongOriginRad);
160
161 M =
a * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad
162 - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(2 * LatRad)
163 + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(4 * LatRad)
164 - (35 * eccSquared * eccSquared * eccSquared / 3072) * sin(6 * LatRad));
165
166 UTMEasting = (double) (k0 * N * (A + (1 - T + C) * A * A * A / 6
167 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120)
168 + 500000.0);
169
170 UTMNorthing = (double) (k0 * (M + N * tan(LatRad)*(A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
171 + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720)));
172 if (Lat < 0)
173 UTMNorthing += 10000000.0;
174}
static Ellipsoid ellipsoid[]
char UTMLetterDesignator(double Lat)
double eccentricitySquared