A Discrete-Event Network Simulator
API
geographic-positions.cc
Go to the documentation of this file.
1 /* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
2 /*
3  * Copyright (c) 2014 University of Washington
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License version 2 as
7  * published by the Free Software Foundation;
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17  *
18  * Author: Benjamin Cizdziel <ben.cizdziel@gmail.com>
19  */
20 
21 #include <ns3/log.h>
22 #include <cmath>
23 #include "geographic-positions.h"
24 
25 NS_LOG_COMPONENT_DEFINE ("GeographicPositions");
26 
27 namespace ns3 {
28 
30 static constexpr double EARTH_RADIUS = 6371e3;
31 
44 static constexpr double EARTH_SEMIMAJOR_AXIS = 6378137;
45 
47 static constexpr double EARTH_GRS80_ECCENTRICITY = 0.0818191910428158;
48 
50 static constexpr double EARTH_WGS84_ECCENTRICITY = 0.0818191908426215;
51 
53 static constexpr double DEG2RAD = M_PI / 180.0;
54 
56 static constexpr double RAD2DEG = 180.0 * M_1_PI;
57 
58 Vector
60  double longitude,
61  double altitude,
62  EarthSpheroidType sphType)
63 {
65  double latitudeRadians = DEG2RAD * latitude;
66  double longitudeRadians = DEG2RAD * longitude;
67  double a; // semi-major axis of earth
68  double e; // first eccentricity of earth
69  if (sphType == SPHERE)
70  {
71  a = EARTH_RADIUS;
72  e = 0;
73  }
74  else if (sphType == GRS80)
75  {
78  }
79  else // if sphType == WGS84
80  {
83  }
84 
85  double Rn = a / (sqrt (1 - pow (e, 2) * pow (sin (latitudeRadians), 2))); // radius of
86  // curvature
87  double x = (Rn + altitude) * cos (latitudeRadians) * cos (longitudeRadians);
88  double y = (Rn + altitude) * cos (latitudeRadians) * sin (longitudeRadians);
89  double z = ((1 - pow (e, 2)) * Rn + altitude) * sin (latitudeRadians);
90  Vector cartesianCoordinates = Vector (x, y, z);
91  return cartesianCoordinates;
92 }
93 
94 Vector
96 {
97  NS_LOG_FUNCTION (pos << sphType);
98 
99  double a; // semi-major axis of earth
100  double e; // first eccentricity of earth
101  if (sphType == SPHERE)
102  {
103  a = EARTH_RADIUS;
104  e = 0;
105  }
106  else if (sphType == GRS80)
107  {
110  }
111  else // if sphType == WGS84
112  {
115  }
116 
117  Vector lla, tmp;
118  lla.y = atan2 (pos.y, pos.x); // longitude (rad), in +/- pi
119 
120  double e2 = e * e;
121  // sqrt (pos.x^2 + pos.y^2)
122  double p = CalculateDistance (pos, {0, 0, pos.z});
123  lla.x = atan2 (pos.z, p * (1 - e2)); // init latitude (rad), in +/- pi
124 
125  do
126  {
127  tmp = lla;
128  double N = a / sqrt (1 - e2 * sin (tmp.x) * sin (tmp.x));
129  double v = p / cos (tmp.x);
130  lla.z = v - N; // altitude
131  lla.x = atan2 (pos.z, p * (1 - e2 * N / v));
132  }
133  // 1 m difference is approx 1 / 30 arc seconds = 9.26e-6 deg
134  while (fabs (lla.x - tmp.x) > 0.00000926 * DEG2RAD);
135 
136  lla.x *= RAD2DEG;
137  lla.y *= RAD2DEG;
138 
139  // canonicalize (latitude) x in [-90, 90] and (longitude) y in [-180, 180)
140  if (lla.x > 90.0)
141  {
142  lla.x = 180 - lla.x;
143  lla.y += lla.y < 0 ? 180 : -180;
144  }
145  else if (lla.x < -90.0)
146  {
147  lla.x = -180 - lla.x;
148  lla.y += lla.y < 0 ? 180 : -180;
149  }
150  if (lla.y == 180.0) lla.y = -180;
151 
152  // make sure lat/lon in the right range to double check canonicalization
153  // and conversion routine
154  NS_ASSERT_MSG (-180.0 <= lla.y, "Conversion error: longitude too negative");
155  NS_ASSERT_MSG (180.0 > lla.y, "Conversion error: longitude too positive");
156  NS_ASSERT_MSG (-90.0 <= lla.x, "Conversion error: latitude too negative");
157  NS_ASSERT_MSG (90.0 >= lla.x, "Conversion error: latitude too positive");
158 
159  return lla;
160 }
161 
162 std::list<Vector>
164  double originLongitude,
165  double maxAltitude,
166  int numPoints,
167  double maxDistFromOrigin,
169 {
171  // fixes divide by zero case and limits latitude bounds
172  if (originLatitude >= 90)
173  {
174  NS_LOG_WARN ("origin latitude must be less than 90. setting to 89.999");
175  originLatitude = 89.999;
176  }
177  else if (originLatitude <= -90)
178  {
179  NS_LOG_WARN ("origin latitude must be greater than -90. setting to -89.999");
180  originLatitude = -89.999;
181  }
182 
183  // restricts maximum altitude from being less than zero (below earth's surface).
184  // sets maximum altitude equal to zero if parameter is set to be less than zero.
185  if (maxAltitude < 0)
186  {
187  NS_LOG_WARN ("maximum altitude must be greater than or equal to 0. setting to 0");
188  maxAltitude = 0;
189  }
190 
191  double originLatitudeRadians = originLatitude * DEG2RAD;
192  double originLongitudeRadians = originLongitude * DEG2RAD;
193  double originColatitude = (M_PI_2) - originLatitudeRadians;
194 
195  double a = maxDistFromOrigin / EARTH_RADIUS; // maximum alpha allowed
196  // (arc length formula)
197  if (a > M_PI)
198  {
199  a = M_PI; // pi is largest alpha possible (polar angle from origin that
200  // points can be generated within)
201  }
202 
203  std::list<Vector> generatedPoints;
204  for (int i = 0; i < numPoints; i++)
205  {
206  // random distance from North Pole (towards center of earth)
207  double d = uniRand->GetValue (0, EARTH_RADIUS - EARTH_RADIUS * cos (a));
208  // random angle in latitude slice (wrt Prime Meridian), radians
209  double phi = uniRand->GetValue (0, M_PI * 2);
210  // random angle from Center of Earth (wrt North Pole), radians
211  double alpha = acos((EARTH_RADIUS - d) / EARTH_RADIUS);
212 
213  // shift coordinate system from North Pole referred to origin point referred
214  // reference: http://en.wikibooks.org/wiki/General_Astronomy/Coordinate_Systems
215  double theta = M_PI_2 - alpha; // angle of elevation of new point wrt
216  // origin point (latitude in coordinate
217  // system referred to origin point)
218  double randPointLatitude = asin(sin(theta)*cos(originColatitude) +
219  cos(theta)*sin(originColatitude)*sin(phi));
220  // declination
221  double intermedLong = asin((sin(randPointLatitude)*cos(originColatitude) -
222  sin(theta)) / (cos(randPointLatitude)*sin(originColatitude)));
223  // right ascension
224  intermedLong = intermedLong + M_PI_2; // shift to longitude 0
225 
226  //flip / mirror point if it has phi in quadrant II or III (wasn't
227  //resolved correctly by arcsin) across longitude 0
228  if (phi > (M_PI_2) && phi <= (3 * M_PI_2))
229  intermedLong = -intermedLong;
230 
231  // shift longitude to be referenced to origin
232  double randPointLongitude = intermedLong + originLongitudeRadians;
233 
234  // random altitude above earth's surface
235  double randAltitude = uniRand->GetValue (0, maxAltitude);
236 
238  (randPointLatitude * RAD2DEG,
239  randPointLongitude * RAD2DEG,
240  randAltitude,
241  SPHERE);
242  // convert coordinates
243  // from geographic to cartesian
244 
245  generatedPoints.push_back (pointPosition); //add generated coordinate
246  //points to list
247  }
248  return generatedPoints;
249 }
250 
251 } // namespace ns3
252 
static std::list< Vector > RandCartesianPointsAroundGeographicPoint(double originLatitude, double originLongitude, double maxAltitude, int numPoints, double maxDistFromOrigin, Ptr< UniformRandomVariable > uniRand)
Generates uniformly distributed random points (in ECEF Cartesian coordinates) within a given altitude...
EarthSpheroidType
Spheroid model to use for earth: perfect sphere (SPHERE), Geodetic Reference System 1980 (GRS80),...
static Vector GeographicToCartesianCoordinates(double latitude, double longitude, double altitude, EarthSpheroidType sphType)
Converts earth geographic/geodetic coordinates (latitude and longitude in degrees) with a given altit...
static Vector CartesianToGeographicCoordinates(Vector pos, EarthSpheroidType sphType)
Inverse of GeographicToCartesianCoordinates using [1].
double GetValue(double min, double max)
Get the next random value, as a double in the specified range .
#define NS_ASSERT_MSG(condition, message)
At runtime, in debugging builds, if this condition is not true, the program prints the message to out...
Definition: assert.h:88
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition: log.h:205
#define NS_LOG_FUNCTION_NOARGS()
Output the name of the function.
#define NS_LOG_FUNCTION(parameters)
If log level LOG_FUNCTION is enabled, this macro will output all input parameters separated by ",...
#define NS_LOG_WARN(msg)
Use NS_LOG to output a message of level LOG_WARN.
Definition: log.h:265
Every class exported by the ns3 library is enclosed in the ns3 namespace.
static constexpr double DEG2RAD
Conversion factor: degrees to radians.
static constexpr double EARTH_GRS80_ECCENTRICITY
Earth's first eccentricity as defined by GRS80.
static constexpr double EARTH_RADIUS
Earth's radius in meters if modeled as a perfect sphere.
static constexpr double EARTH_SEMIMAJOR_AXIS
GRS80 and WGS84 sources.
static constexpr double RAD2DEG
Conversion factor: radians to degrees.
double CalculateDistance(const Vector3D &a, const Vector3D &b)
Definition: vector.cc:105
static constexpr double EARTH_WGS84_ECCENTRICITY
Earth's first eccentricity as defined by WGS84.
float alpha
Plot alpha value (transparency)
list x
Random number samples.