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 const double EARTH_RADIUS = 6371e3;
31 
43 static const double EARTH_SEMIMAJOR_AXIS = 6378137;
45 
47 static const double EARTH_GRS80_ECCENTRICITY = 0.0818191910428158;
48 
50 static const double EARTH_WGS84_ECCENTRICITY = 0.0818191908426215;
51 
52 Vector
54  double longitude,
55  double altitude,
56  EarthSpheroidType sphType)
57 {
59  double latitudeRadians = 0.01745329 * latitude;
60  double longitudeRadians = 0.01745329 * longitude;
61  double a; // semi-major axis of earth
62  double e; // first eccentricity of earth
63  if (sphType == SPHERE)
64  {
65  a = EARTH_RADIUS;
66  e = 0;
67  }
68  else if (sphType == GRS80)
69  {
72  }
73  else // if sphType == WGS84
74  {
77  }
78 
79  double Rn = a / (sqrt (1 - pow (e, 2) * pow (sin (latitudeRadians), 2))); // radius of
80  // curvature
81  double x = (Rn + altitude) * cos (latitudeRadians) * cos (longitudeRadians);
82  double y = (Rn + altitude) * cos (latitudeRadians) * sin (longitudeRadians);
83  double z = ((1 - pow (e, 2)) * Rn + altitude) * sin (latitudeRadians);
84  Vector cartesianCoordinates = Vector (x, y, z);
85  return cartesianCoordinates;
86 }
87 
88 std::list<Vector>
90  double originLongitude,
91  double maxAltitude,
92  int numPoints,
93  double maxDistFromOrigin,
95 {
97  // fixes divide by zero case and limits latitude bounds
98  if (originLatitude >= 90)
99  {
100  NS_LOG_WARN ("origin latitude must be less than 90. setting to 89.999");
101  originLatitude = 89.999;
102  }
103  else if (originLatitude <= -90)
104  {
105  NS_LOG_WARN ("origin latitude must be greater than -90. setting to -89.999");
106  originLatitude = -89.999;
107  }
108 
109  // restricts maximum altitude from being less than zero (below earth's surface).
110  // sets maximum altitude equal to zero if parameter is set to be less than zero.
111  if (maxAltitude < 0)
112  {
113  NS_LOG_WARN ("maximum altitude must be greater than or equal to 0. setting to 0");
114  maxAltitude = 0;
115  }
116 
117  double originLatitudeRadians = originLatitude * (M_PI / 180);
118  double originLongitudeRadians = originLongitude * (M_PI / 180);
119  double originColatitude = (M_PI / 2) - originLatitudeRadians;
120 
121  double a = maxDistFromOrigin / EARTH_RADIUS; // maximum alpha allowed
122  // (arc length formula)
123  if (a > M_PI)
124  {
125  a = M_PI; // pi is largest alpha possible (polar angle from origin that
126  // points can be generated within)
127  }
128 
129  std::list<Vector> generatedPoints;
130  for (int i = 0; i < numPoints; i++)
131  {
132  // random distance from North Pole (towards center of earth)
133  double d = uniRand->GetValue (0, EARTH_RADIUS - EARTH_RADIUS * cos (a));
134  // random angle in latitude slice (wrt Prime Meridian), radians
135  double phi = uniRand->GetValue (0, M_PI * 2);
136  // random angle from Center of Earth (wrt North Pole), radians
137  double alpha = acos((EARTH_RADIUS - d) / EARTH_RADIUS);
138 
139  // shift coordinate system from North Pole referred to origin point referred
140  // reference: http://en.wikibooks.org/wiki/General_Astronomy/Coordinate_Systems
141  double theta = M_PI / 2 - alpha; // angle of elevation of new point wrt
142  // origin point (latitude in coordinate
143  // system referred to origin point)
144  double randPointLatitude = asin(sin(theta)*cos(originColatitude) +
145  cos(theta)*sin(originColatitude)*sin(phi));
146  // declination
147  double intermedLong = asin((sin(randPointLatitude)*cos(originColatitude) -
148  sin(theta)) / (cos(randPointLatitude)*sin(originColatitude)));
149  // right ascension
150  intermedLong = intermedLong + M_PI / 2; // shift to longitude 0
151 
152  //flip / mirror point if it has phi in quadrant II or III (wasn't
153  //resolved correctly by arcsin) across longitude 0
154  if (phi > (M_PI / 2) && phi <= ((3 * M_PI) / 2))
155  intermedLong = -intermedLong;
156 
157  // shift longitude to be referenced to origin
158  double randPointLongitude = intermedLong + originLongitudeRadians;
159 
160  // random altitude above earth's surface
161  double randAltitude = uniRand->GetValue (0, maxAltitude);
162 
164  (randPointLatitude * (180/M_PI),
165  randPointLongitude * (180/M_PI),
166  randAltitude,
167  SPHERE);
168  // convert coordinates
169  // from geographic to cartesian
170 
171  generatedPoints.push_back (pointPosition); //add generated coordinate
172  //points to list
173  }
174  return generatedPoints;
175 }
176 
177 } // namespace ns3
178 
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition: log.h:202
#define NS_LOG_FUNCTION_NOARGS()
Output the name of the function.
static const double EARTH_WGS84_ECCENTRICITY
Earth&#39;s first eccentricity as defined by WGS84.
Every class exported by the ns3 library is enclosed in the ns3 namespace.
double GetValue(double min, double max)
Get the next random value, as a double in the specified range .
static const double EARTH_RADIUS
Earth&#39;s radius in meters if modeled as a perfect sphere.
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)...
#define NS_LOG_WARN(msg)
Use NS_LOG to output a message of level LOG_WARN.
Definition: log.h:262
static const double EARTH_SEMIMAJOR_AXIS
GRS80 and WGS84 sources.
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 const double EARTH_GRS80_ECCENTRICITY
Earth&#39;s first eccentricity as defined by GRS80.