/Designs/skrysohledac2/SW/utm.c
24,9 → 24,9
 
*/
 
#include "stdafx.h"
#include "Common.h"
#include "osgb.h"
//#include "stdafx.h"
//#include "Common.h"
//#include "osgb.h"
#include "utm.h"
#include <math.h>
 
117,8 → 117,8
// asterisk character ('*').
//=======================================================================
 
void LatLonToUtm (double a, double f, int& utmXZone, char& utmYZone,
double& easting, double& northing, double lat, double lon)
void LatLonToUtm (double a, double f, int* utmXZone, char* utmYZone,
double* easting, double* northing, double lat, double lon)
{
double recf;
double b;
142,24 → 142,26
double nfn;
 
if (lon <= 0.0) {
utmXZone = 30 + (int)(lon / 6.0);
*utmXZone = 30 + (int)(lon / 6.0);
} else {
utmXZone = 31 + (int)(lon / 6.0);
*utmXZone = 31 + (int)(lon / 6.0);
}
if (lat < 84.0 && lat >= 72.0) {
// Special case: zone X is 12 degrees from north to south, not 8.
utmYZone = cArray[19];
*utmYZone = cArray[19];
} else {
utmYZone = cArray[(int)((lat + 80.0) / 8.0)];
*utmYZone = cArray[(int)((lat + 80.0) / 8.0)];
}
if (lat >= 84.0 || lat < -80.0) {
// Invalid coordinate; the vertical zone is set to the invalid
// character.
utmYZone = '*';
*utmYZone = '*';
}
 
double latRad = lat * deg2rad;
double lonRad = lon * deg2rad;
double latRad = lat * (M_PI/180);
double lonRad = lon * (M_PI/180);
// double latRad = lat * deg2rad;
// double lonRad = lon * deg2rad;
recf = 1.0 / f;
b = a * (recf - 1.0) / recf;
eSquared = CalculateESquared (a, b);
175,7 → 177,7
dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0
* (tn * tn * tn * tn * tn) / 16.0) / 48.0;
ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0;
olam = (utmXZone * 6 - 183) * deg2rad;
olam = (*utmXZone * 6 - 183) * (M_PI/180);
dlam = lonRad - olam;
s = sin (latRad);
c = cos (latRad);
188,12 → 190,12
t3 = sn * s * (c * c * c) * ok * (5.0 - (t * t) + 9.0 * eta + 4.0
* (eta * eta)) / 24.0;
if (latRad < 0.0) nfn = 10000000.0; else nfn = 0;
northing = nfn + t1 + (dlam * dlam) * t2 + (dlam * dlam * dlam
*northing = nfn + t1 + (dlam * dlam) * t2 + (dlam * dlam * dlam
* dlam) * t3 + (dlam * dlam * dlam * dlam * dlam * dlam) + 0.5;
t6 = sn * c * ok;
t7 = sn * (c * c * c) * (1.0 - (t * t) + eta) / 6.0;
easting = fe + dlam * t6 + (dlam * dlam * dlam) * t7 + 0.5;
if (northing >= 9999999.0) northing = 9999999.0;
*easting = fe + dlam * t6 + (dlam * dlam * dlam) * t7 + 0.5;
if (*northing >= 9999999.0) *northing = 9999999.0;
}
 
//=======================================================================
203,187 → 205,9
// LatLonToUtm() member function.)
//=======================================================================
 
void LatLonToUtmWGS84 (int& utmXZone, char& utmYZone,
double& easting, double& northing, double lat, double lon)
void LatLonToUtmWGS84 (int* utmXZone, char* utmYZone,
double* easting, double* northing, double lat, double lon)
{
LatLonToUtm (6378137.0, 1 / 298.257223563, utmXZone, utmYZone,
easting, northing, lat, lon);
}
 
 
 
 
//=======================================================================
// Purpose:
// This function converts the specified UTM coordinate to a lat/lon
// coordinate.
// Pre:
// - utmXZone must be between 1 and 60, inclusive.
// - utmYZone must be one of: CDEFGHJKLMNPQRSTUVWX
// Parameters:
// double a:
// Ellipsoid semi-major axis, in meters. (For WGS84 datum, use 6378137.0)
// double f:
// Ellipsoid flattening. (For WGS84 datum, use 1 / 298.257223563)
// int utmXZone:
// The horizontal zone number of the UTM coordinate.
// char utmYZone:
// The vertical zone letter of the UTM coordinate.
// double easting, double northing:
// The UTM coordinate to convert.
// double& lat:
// Upon exit, lat contains the latitude.
// double& lon:
// Upon exit, lon contains the longitude.
// Notes:
// The code in this function is a C conversion of some of the source code
// from the Mapping Datum Transformation Software (MADTRAN) program, written
// in PowerBasic. To get the source code for MADTRAN, go to:
//
// http://164.214.2.59/publications/guides/MADTRAN/index.html
//
// and download MADTRAN.ZIP
//=======================================================================
 
void UtmToLatLon (double a, double f, int utmXZone, char utmYZone,
double easting, double northing, double& lat, double& lon)
{
double recf;
double b;
double eSquared;
double e2Squared;
double tn;
double ap;
double bp;
double cp;
double dp;
double ep;
double nfn;
double tmd;
double sr;
double sn;
double ftphi;
double s;
double c;
double t;
double eta;
double de;
double dlam;
double olam;
 
recf = 1.0 / f;
b = a * (recf - 1) / recf;
eSquared = CalculateESquared (a, b);
e2Squared = CalculateE2Squared (a, b);
tn = (a - b) / (a + b);
ap = a * (1.0 - tn + 5.0 * ((tn * tn) - (tn * tn * tn)) / 4.0 + 81.0 *
((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 64.0);
bp = 3.0 * a * (tn - (tn * tn) + 7.0 * ((tn * tn * tn)
- (tn * tn * tn * tn)) / 8.0 + 55.0 * (tn * tn * tn * tn * tn) / 64.0)
/ 2.0;
cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn)
- (tn * tn * tn * tn * tn)) / 4.0) / 16.0;
dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0
* (tn * tn * tn * tn * tn) / 16.0) / 48.0;
ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0;
if ((utmYZone <= 'M' && utmYZone >= 'C')
|| (utmYZone <= 'm' && utmYZone >= 'c')) {
nfn = 10000000.0;
} else {
nfn = 0;
}
tmd = (northing - nfn) / ok;
sr = sphsr (a, eSquared, 0.0);
ftphi = tmd / sr;
double t10, t11, t14, t15;
for (int i = 0; i < 5; i++) {
t10 = sphtmd (ap, bp, cp, dp, ep, ftphi);
sr = sphsr (a, eSquared, ftphi);
ftphi = ftphi + (tmd - t10) / sr;
}
sr = sphsr (a, eSquared, ftphi);
sn = sphsn (a, eSquared, ftphi);
s = sin (ftphi);
c = cos (ftphi);
t = s / c;
eta = e2Squared * (c * c);
de = easting - fe;
t10 = t / (2.0 * sr * sn * (ok * ok));
t11 = t * (5.0 + 3.0 * (t * t) + eta - 4.0 * (eta * eta) - 9.0 * (t * t)
* eta) / (24.0 * sr * (sn * sn * sn) * (ok * ok * ok * ok));
lat = ftphi - (de * de) * t10 + (de * de * de * de) * t11;
t14 = 1.0 / (sn * c * ok);
t15 = (1.0 + 2.0 * (t * t) + eta) / (6 * (sn * sn * sn) * c
* (ok * ok * ok));
dlam = de * t14 - (de * de * de) * t15;
olam = (utmXZone * 6 - 183.0) * deg2rad;
lon = olam + dlam;
lon *= rad2deg;
lat *= rad2deg;
}
 
//=======================================================================
// Purpose:
// This function converts the specified UTM coordinate to a lat/lon
// coordinate in the WGS84 datum. (See the comment block for the
// UtmToLatLon() member function.
//=======================================================================
 
void UtmToLatLonWGS84 (int utmXZone, char utmYZone, double easting,
double northing, double& lat, double& lon)
{
UtmToLatLon (6378137.0, 1 / 298.257223563, utmXZone, utmYZone,
easting, northing, lat, lon);
}
 
//=======================================================================
/**
@func Build a position string
@parm target. must be 30 characters or longer.
*/
//=======================================================================
 
void CUtmPoint::GetString(TCHAR *position) const
{
_stprintf(position,
_T("%02d%c %06d %07d"),
m_xzone, m_yzone,
(int)m_easting,
(int)m_northing);
}
 
//=======================================================================
/**
@func get the position of a UTM point
@parm point out
*/
//=======================================================================
 
void CUtmPoint::ToPosition(CPosition &pos) const
{
double lat,lon;
UtmToLatLonWGS84(m_xzone,m_yzone,m_easting,m_northing,
lat,lon);
pos.Clear();
pos.SetLatitude(lat);
pos.SetLongitude(lon);
}
 
//=======================================================================
/**
@func turn a position into a UTM point
@parm position
@rdesc true if it was in range
*/
//=======================================================================
 
bool CUtmPoint::FromPosition(const CPosition &pos)
{
Clear();
if(!IsPositionInUtmSpace(pos))
return false;
LatLonToUtmWGS84(m_xzone,m_yzone,m_easting,m_northing,
pos.GetLatitude(),
pos.GetLongitude());
return true;
}