{\r
public static class Constants\r
{\r
- #region fields\r
- // based on WGS84\r
- public const double SemiMajorAxisA = 6378137.00000;\r
- public const double SemiMajorAxisB = 6356752.31425;\r
- public const double FirstEccentricitySquared = 6.69437999013e-3;\r
- public const double SecondEccentricitySquared = 6.73949674226e-3;\r
- public const double Pi = 3.1415926535898;\r
+ #region type definitions\r
+ public static class Wgs84\r
+ {\r
+ public const double SemiMajorAxisA = 6378137.00000;\r
+ public const double SemiMajorAxisB = 6356752.31425;\r
+ public const double FirstEccentricitySquared = 6.69437999013e-3;\r
+ public const double SecondEccentricitySquared = 6.73949674226e-3;\r
+ public const double InverseFlattening = 298.257223563;\r
+ public const double Flattening = 3.3528106647474807e-3;\r
+ public const double Pi = 3.1415926535898;\r
+ public const double SpeedOfLight = 299792458; \r
+ }\r
+\r
+ public static class Grs80\r
+ {\r
+ public const double SemiMajorAxisA = 6378137.00000;\r
+ public const double SemiMajorAxisB = 6356752.31414;\r
+ public const double FirstEccentricitySquared = 6.69438002290e-3;\r
+ public const double SecondEccentricitySquared = 6.73949677548e-3;\r
+ public const double InverseFlattening = 298.257222101;\r
+ public const double Flattening = 3.3528106811823189e-3;\r
+ }\r
+\r
+ public static class Bessel\r
+ {\r
+ }\r
#endregion\r
}\r
}\r
using System;\r
using System.Collections.Generic;\r
using System.Text;\r
+using Const = Yubeshi.Constants.Wgs84;\r
\r
namespace Yubeshi\r
{\r
{\r
#region fields\r
private double degree;\r
- private static readonly double degToRad = Constants.Pi / 180.0;\r
- private static readonly double radToDeg = 180.0 / Constants.Pi;\r
+ private static readonly double degToRad = Const.Pi / 180.0;\r
+ private static readonly double radToDeg = 180.0 / Const.Pi;\r
#endregion\r
\r
#region constructors\r
using System;\r
using System.Collections.Generic;\r
using System.Text;\r
+using Const = Yubeshi.Constants.Wgs84;\r
\r
namespace Yubeshi\r
{\r
public class EcefCoordinate\r
{\r
#region fields\r
- private const double a = Constants.SemiMajorAxisA;\r
- private const double b = Constants.SemiMajorAxisB;\r
- private const double e1sq = Constants.FirstEccentricitySquared;\r
- private const double e2sq = Constants.SecondEccentricitySquared;\r
- private const double pi = Constants.Pi;\r
+ private const double a = Const.SemiMajorAxisA;\r
+ private const double b = Const.SemiMajorAxisB;\r
+ private const double e1sq = Const.FirstEccentricitySquared;\r
+ private const double e2sq = Const.SecondEccentricitySquared;\r
+ private const double pi = Const.Pi;\r
#endregion\r
\r
#region constructors\r
double phi = Math.Atan2(Z + e2sq * b * sin * sin * sin,\r
p - e1sq * a * cos * cos * cos);\r
double sinPhi = Math.Sin(phi);\r
+ double cosPhi = Math.Cos(phi);\r
double khi = Math.Sqrt(1.0 - e1sq * sinPhi * sinPhi);\r
double n = a / khi;\r
- double h = p / Math.Cos(phi) - n;\r
- return new GeodeticCoordinate(phi * 180 / pi, lambda * 180 / pi, h);\r
+ double ht;\r
+ if (cosPhi < 0.5)\r
+ {\r
+ ht = Z / sinPhi - n * (1.0 - e1sq);\r
+ }\r
+ else\r
+ { \r
+ ht = p / Math.Cos(phi) - n;\r
+ }\r
+ Height h = new Height(ht, Height.Base.Wgs84Ellipsoid);\r
+ return new GeodeticCoordinate(\r
+ phi * 180 / pi, lambda * 180 / pi, h);\r
}\r
#endregion\r
}\r
\r
#region constructors\r
\r
- public EnuCoordinate(double e, double n, double u)\r
+ public EnuCoordinate(double e, double n, Height u)\r
: this(e, n, u, null)\r
{\r
}\r
\r
- public EnuCoordinate(double e, double n, double u, \r
+ public EnuCoordinate(double e, double n, Height u, \r
EcefCoordinate origin)\r
{\r
E = e;\r
public class GeodeticCoordinate\r
{\r
#region fields\r
- private const double a = Constants.SemiMajorAxisA;\r
- private const double b = Constants.SemiMajorAxisB;\r
- private const double e1sq = Constants.FirstEccentricitySquared;\r
- private const double pi = Constants.Pi;\r
+ private const double a = Constants.Wgs84.SemiMajorAxisA;\r
+ private const double b = Constants.Wgs84.SemiMajorAxisB;\r
+ private const double e1sq = Constants.Wgs84.FirstEccentricitySquared;\r
+ private const double pi = Constants.Wgs84.Pi;\r
\r
#endregion\r
#region constructors\r
}\r
\r
public GeodeticCoordinate(\r
- Degree latitude, Degree longitude, double altitude)\r
+ Degree latitude, Degree longitude, Height altitude)\r
{\r
Latitude = latitude;\r
Longitude = longitude;\r
#region properties\r
\r
/// <summary>\r
- /// in Degree\r
+ /// geographic latitude in degree\r
/// </summary>\r
public Degree Latitude\r
{\r
get;\r
set;\r
}\r
+\r
+ /// <summary>\r
+ /// geocentric latitude in degree\r
+ /// </summary>\r
+ public Degree GeocentricLatitude\r
+ {\r
+ get\r
+ {\r
+ return Math.Atan((1 - e1sq) * Math.Tan(Latitude.Radian));\r
+ }\r
+ }\r
+\r
/// <summary>\r
- /// in Degree\r
+ /// longitude in degree\r
/// </summary>\r
public Degree Longitude\r
{\r
}\r
\r
/// <summary>\r
- /// in metre\r
+ /// altitude in metre\r
/// </summary>\r
- public double Altitude\r
+ public Height Altitude\r
{\r
get;\r
set;\r
[Test]\r
public void ToGeodeticCoordinate()\r
{\r
- GeodeticCoordinate geo = C.SkyTreeTop;\r
- EcefCoordinate ecef = new EcefCoordinate(-3961181.340,\r
- 3346187.629,\r
- 3702487.226);\r
+ GeodeticCoordinate geo = C.SkytreeTop;\r
+ EcefCoordinate ecef = new EcefCoordinate(-3961187.702,\r
+ 3346179.551,\r
+ 3702497.996);\r
\r
GeodeticCoordinate converted = ecef.ToGeodeticCoordinate();\r
\r
Assert.AreEqual(geo.Longitude, converted.Longitude, 1e-6);\r
Assert.AreEqual(geo.Altitude, converted.Altitude, 2e-4);\r
\r
+ // North Pole\r
+ geo = new GeodeticCoordinate(90.0, 0.0, \r
+ new Height(0, Height.Base.Wgs84Ellipsoid));\r
+ ecef = geo.ToEcefCoordinate();\r
+ ecef = new EcefCoordinate(0, 0, Constants.Wgs84.SemiMajorAxisB);\r
+\r
+ converted = ecef.ToGeodeticCoordinate();\r
+\r
+ Assert.AreEqual(geo.Latitude, converted.Latitude, 1e-6);\r
+ Assert.AreEqual(geo.Longitude, converted.Longitude, 1e-6);\r
+ Assert.AreEqual(geo.Altitude, converted.Altitude, 2e-4);\r
+\r
}\r
}\r
}\r
[Test]\r
public void ZeroDistance()\r
{\r
- EnuCoordinate enu = new EnuCoordinate(C.SkyTreeTop, C.SkyTreeTop);\r
+ EnuCoordinate enu = new EnuCoordinate(C.SkytreeTop, C.SkytreeTop);\r
Assert.AreEqual(0.0, enu.E);\r
Assert.AreEqual(0.0, enu.N);\r
Assert.AreEqual(0.0, enu.U);\r
public void UpOnly()\r
{\r
EnuCoordinate enu = new EnuCoordinate(\r
- C.SkyTreeTop, C.SkyTreeBottom);\r
+ C.SkytreeTop, C.SkytreeBottom);\r
\r
Assert.AreEqual(0.0, enu.E, metreError);\r
Assert.AreEqual(0.0, enu.N, metreError);\r
[Test]\r
public void TwoPoints()\r
{\r
- GeodeticCoordinate c1 = C.SkyTreeTop;\r
+ GeodeticCoordinate c1 = C.SkytreeTop;\r
GeodeticCoordinate c2 = C.TokyoTowerTop;\r
GeodeticCoordinate c0 = new GeodeticCoordinate(\r
0.5 * (c1.Latitude + c2.Latitude),\r
// c2->c1 45.965 + 180 deg\r
// 45.9455 = 0.5 * (45.926 + 45.965)\r
Assert.AreEqual(45.946, direction, 1.0 / 60.0);\r
- Assert.AreEqual(283.0, du, 1e-3);\r
+ Assert.AreEqual(289.0, du, 1e-3);\r
}\r
\r
[Test]\r
public void ToEcefCoordinate()\r
{\r
EnuCoordinate enu = new EnuCoordinate(\r
- C.SkyTreeTop, C.TokyoTowerTop);\r
+ C.SkytreeTop, C.TokyoTowerTop);\r
EcefCoordinate ecef1 = enu.ToEcefCoordinate();\r
- EcefCoordinate ecef2 = C.SkyTreeTop.ToEcefCoordinate();\r
+ EcefCoordinate ecef2 = C.SkytreeTop.ToEcefCoordinate();\r
\r
Assert.AreEqual(ecef2.X, ecef1.X, metreError);\r
Assert.AreEqual(ecef2.Y, ecef1.Y, metreError);\r
Assert.AreEqual(ecef.Y, converted.Y, 1e-3);\r
Assert.AreEqual(ecef.Z, converted.Z, 1e-3);\r
\r
-\r
-\r
}\r
}\r
}\r
{\r
public class SampleCoordinates\r
{\r
- public static readonly GeodeticCoordinate SkyTreeTop = \r
- new GeodeticCoordinate(35.710058, 139.810719, 634.0);\r
+ public static readonly GeodeticCoordinate SkytreeTop =\r
+ new GeodeticCoordinate(35.710139, 139.810833, \r
+ new Height(640.0, Height.Base.Geoid));\r
\r
- public static readonly GeodeticCoordinate SkyTreeBottom =\r
- new GeodeticCoordinate(35.710058, 139.810719, 0.0);\r
+ public static readonly GeodeticCoordinate SkytreeBottom =\r
+ new GeodeticCoordinate(35.710139, 139.810833, \r
+ new Height(6.0, Height.Base.Geoid));\r
\r
public static readonly GeodeticCoordinate TokyoTowerTop =\r
- new GeodeticCoordinate(35.658611, 139.745556, 351.0);\r
+ new GeodeticCoordinate(35.658611, 139.745556,\r
+ new Height(351.0, Height.Base.Geoid));\r
}\r
}\r