Rev Author Line No. Line
1212 kakl 1 /************************************************************************
2 *
3 * File: Utm.cpp
4 * RCS: $Header: /cvsroot/stelvio/stelvio/NavStar/Utm.cpp,v 1.1 2001/03/18 20:07:03 steve_l Exp $
5 * Author: Steve Loughran
6 * Created: 2001
7 * Language: C++
8 * Package:
9 * Status: Experimental
10 * @doc
11 *
12 ************************************************************************/
13  
14 /*
15 This is code to do UTM conversion.
16  
17 I took this code from Jason Bevins' GPS thing which blagged the VB algorithms
18 from the Mapping Datum Transformation Software (MADTRAN) program,
19 written in PowerBasic. To get the source code for MADTRAN, go to:
20  
21 http://164.214.2.59/publications/guides/MADTRAN/index.html
22  
23 this version retains the core algorithms as static functions
24  
25 */
26  
1269 kakl 27 //#include "stdafx.h"
28 //#include "Common.h"
29 //#include "osgb.h"
1212 kakl 30 #include "utm.h"
31 #include <math.h>
32  
33 // Some constants used by these functions.
34 static const double fe = 500000.0;
35 static const double ok = 0.9996;
36  
37 // An array containing each vertical UTM zone.
38 static char cArray[] = "CDEFGHJKLMNPQRSTUVWX";
39  
40  
41 /////////////////////////////////////////////////////////////////////////////
42 // Miscellaneous functions for these UTM conversion formulas.
43  
44 double CalculateESquared (double a, double b)
45 {
46 return ((a * a) - (b * b)) / (a * a);
47 }
48  
49  
50 double CalculateE2Squared (double a, double b)
51 {
52 return ((a * a) - (b * b)) / (b * b);
53 }
54  
55  
56 double denom (double es, double sphi)
57 {
58 double sinSphi = sin (sphi);
59 return sqrt (1.0 - es * (sinSphi * sinSphi));
60 }
61  
62  
63 double sphsr (double a, double es, double sphi)
64 {
65 double dn = denom (es, sphi);
66 return a * (1.0 - es) / (dn * dn * dn);
67 }
68  
69  
70 double sphsn (double a, double es, double sphi)
71 {
72 double sinSphi = sin (sphi);
73 return a / sqrt (1.0 - es * (sinSphi * sinSphi));
74 }
75  
76  
77 double sphtmd (double ap, double bp, double cp, double dp, double ep,
78 double sphi)
79 {
80 return (ap * sphi) - (bp * sin (2.0 * sphi)) + (cp * sin (4.0 * sphi))
81 - (dp * sin (6.0 * sphi)) + (ep * sin (8.0 * sphi));
82 }
83  
84  
85 //=======================================================================
86 // Purpose:
87 // This function converts the specified lat/lon coordinate to a UTM
88 // coordinate.
89 // Parameters:
90 // double a:
91 // Ellipsoid semi-major axis, in meters. (For WGS84 datum, use 6378137.0)
92 // double f:
93 // Ellipsoid flattening. (For WGS84 datum, use 1 / 298.257223563)
94 // int& utmXZone:
95 // Upon exit, this parameter will contain the hotizontal zone number of
96 // the UTM coordinate. The returned value for this parameter is a number
97 // within the range 1 to 60, inclusive.
98 // char& utmYZone:
99 // Upon exit, this parameter will contain the zone letter of the UTM
100 // coordinate. The returned value for this parameter will be one of:
101 // CDEFGHJKLMNPQRSTUVWX.
102 // double& easting:
103 // Upon exit, this parameter will contain the UTM easting, in meters.
104 // double& northing:
105 // Upon exit, this parameter will contain the UTM northing, in meters.
106 // double lat, double lon:
107 // The lat/lon coordinate to convert.
108 // Notes:
109 // - The code in this function is a C conversion of some of the source code
110 // from the Mapping Datum Transformation Software (MADTRAN) program,
111 // written in PowerBasic. To get the source code for MADTRAN, go to:
112 //
113 // http://164.214.2.59/publications/guides/MADTRAN/index.html
114 //
115 // and download MADTRAN.ZIP
116 // - If the UTM zone is out of range, the y-zone character is set to the
117 // asterisk character ('*').
118 //=======================================================================
119  
1269 kakl 120 void LatLonToUtm (double a, double f, int* utmXZone, char* utmYZone,
121 double* easting, double* northing, double lat, double lon)
1212 kakl 122 {
123 double recf;
124 double b;
125 double eSquared;
126 double e2Squared;
127 double tn;
128 double ap;
129 double bp;
130 double cp;
131 double dp;
132 double ep;
133 double olam;
134 double dlam;
135 double s;
136 double c;
137 double t;
138 double eta;
139 double sn;
140 double tmd;
141 double t1, t2, t3, t6, t7;
142 double nfn;
143  
144 if (lon <= 0.0) {
1269 kakl 145 *utmXZone = 30 + (int)(lon / 6.0);
1212 kakl 146 } else {
1269 kakl 147 *utmXZone = 31 + (int)(lon / 6.0);
1212 kakl 148 }
149 if (lat < 84.0 && lat >= 72.0) {
150 // Special case: zone X is 12 degrees from north to south, not 8.
1269 kakl 151 *utmYZone = cArray[19];
1212 kakl 152 } else {
1269 kakl 153 *utmYZone = cArray[(int)((lat + 80.0) / 8.0)];
1212 kakl 154 }
155 if (lat >= 84.0 || lat < -80.0) {
156 // Invalid coordinate; the vertical zone is set to the invalid
157 // character.
1269 kakl 158 *utmYZone = '*';
1212 kakl 159 }
160  
1269 kakl 161 double latRad = lat * (M_PI/180);
162 double lonRad = lon * (M_PI/180);
163 // double latRad = lat * deg2rad;
164 // double lonRad = lon * deg2rad;
1212 kakl 165 recf = 1.0 / f;
166 b = a * (recf - 1.0) / recf;
167 eSquared = CalculateESquared (a, b);
168 e2Squared = CalculateE2Squared (a, b);
169 tn = (a - b) / (a + b);
170 ap = a * (1.0 - tn + 5.0 * ((tn * tn) - (tn * tn * tn)) / 4.0 + 81.0 *
171 ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 64.0);
172 bp = 3.0 * a * (tn - (tn * tn) + 7.0 * ((tn * tn * tn)
173 - (tn * tn * tn * tn)) / 8.0 + 55.0 * (tn * tn * tn * tn * tn) / 64.0)
174 / 2.0;
175 cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn)
176 - (tn * tn * tn * tn * tn)) / 4.0) / 16.0;
177 dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0
178 * (tn * tn * tn * tn * tn) / 16.0) / 48.0;
179 ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0;
1269 kakl 180 olam = (*utmXZone * 6 - 183) * (M_PI/180);
1212 kakl 181 dlam = lonRad - olam;
182 s = sin (latRad);
183 c = cos (latRad);
184 t = s / c;
185 eta = e2Squared * (c * c);
186 sn = sphsn (a, eSquared, latRad);
187 tmd = sphtmd (ap, bp, cp, dp, ep, latRad);
188 t1 = tmd * ok;
189 t2 = sn * s * c * ok / 2.0;
190 t3 = sn * s * (c * c * c) * ok * (5.0 - (t * t) + 9.0 * eta + 4.0
191 * (eta * eta)) / 24.0;
192 if (latRad < 0.0) nfn = 10000000.0; else nfn = 0;
1269 kakl 193 *northing = nfn + t1 + (dlam * dlam) * t2 + (dlam * dlam * dlam
1212 kakl 194 * dlam) * t3 + (dlam * dlam * dlam * dlam * dlam * dlam) + 0.5;
195 t6 = sn * c * ok;
196 t7 = sn * (c * c * c) * (1.0 - (t * t) + eta) / 6.0;
1269 kakl 197 *easting = fe + dlam * t6 + (dlam * dlam * dlam) * t7 + 0.5;
198 if (*northing >= 9999999.0) *northing = 9999999.0;
1212 kakl 199 }
200  
201 //=======================================================================
202 // Purpose:
203 // This function converts the specified lat/lon coordinate to a UTM
204 // coordinate in the WGS84 datum. (See the comment block for the
205 // LatLonToUtm() member function.)
206 //=======================================================================
207  
1269 kakl 208 void LatLonToUtmWGS84 (int* utmXZone, char* utmYZone,
209 double* easting, double* northing, double lat, double lon)
1212 kakl 210 {
211 LatLonToUtm (6378137.0, 1 / 298.257223563, utmXZone, utmYZone,
212 easting, northing, lat, lon);
213 }