Fix retrograd movements and earth rotation

examples/leo-circular-orbit-tracing-example.cc

@@ -16,7 +16,7 @@ void CourseChange (std::string context, Ptr<const MobilityModel> position)
 {
   Vector pos = position->GetPosition ();
   Ptr<const Node> node = position->GetObject<Node> ();
-  outfile << Simulator::Now () << "," << node->GetId () << "," << pos.x << "," << pos.y << "," << pos.z << "," << position->GetVelocity ().GetLength() <<  std::endl;
+  outfile << Simulator::Now () << ":" << node->GetId () << ":" << pos.x << ":" << pos.y << ":" << pos.z << ":" << position->GetVelocity ().GetLength() << std::endl;
 }
 
 int main(int argc, char *argv[])
@@ -42,7 +42,7 @@ int main(int argc, char *argv[])
 
   outfile << "Time,Satellite,x,y,z,Speed" << std::endl;
 
-  Simulator::Stop (Minutes (60.0));
+  Simulator::Stop (Hours (2.0));
   Simulator::Run ();
   Simulator::Destroy ();
 }

model/leo-circular-orbit-mobility-model.cc

@@ -76,7 +76,6 @@ DotProduct (const Vector3D &l, const Vector3D &r)
 double
 LeoCircularOrbitMobilityModel::GetSpeed () const
 {
-  // force non-retrograd movement
   return sqrt (LEO_EARTH_GM_KM_E10 / m_orbitHeight) * 1e5;
 }
 
@@ -91,39 +90,48 @@ LeoCircularOrbitMobilityModel::DoGetVelocity () const
 Vector3D
 LeoCircularOrbitMobilityModel::PlaneNorm () const
 {
-  int sign = 1;
+  double lat = CalcLatitude ();
-  // ensure correct gradient (not against earth rotation)
+  return Vector3D (sin (-m_inclination) * cos (lat),
-  if (m_inclination < M_PI/4)
+  		   sin (-m_inclination) * sin (lat),
-    {
+  		   cos (m_inclination));
-      sign = -1;
-    }
-  return Vector3D (sign * sin (-m_inclination) * cos (m_latitude),
-  		   sign * sin (-m_inclination) * sin (m_latitude),
-  		   sign * cos (m_inclination));
 }
 
 double
 LeoCircularOrbitMobilityModel::GetProgress (Time t) const
 {
   // TODO use nanos or ms instead? does it give higher precision?
-  return (2 * M_PI * ((GetSpeed () * t.GetSeconds ()) / LEO_EARTH_RAD_M)) + m_offset;
+  int sign = 1;
+  // ensure correct gradient (not against earth rotation)
+  if (m_inclination > M_PI/2)
+    {
+      sign = -1;
+    }
+  return sign * (2 * M_PI * ((GetSpeed () * t.GetSeconds ()) / LEO_EARTH_RAD_M)) + m_offset;
 }
 
 Vector3D
 LeoCircularOrbitMobilityModel::RotatePlane (double a, const Vector3D &x) const
 {
-  Vector3D n = m_plane;
+  Vector3D n = PlaneNorm ();
 
   return Product (DotProduct (n, x), n)
     + Product (cos (a), CrossProduct (CrossProduct (n, x), n))
     + Product (sin (a), CrossProduct (n, x));
 }
 
+double
+LeoCircularOrbitMobilityModel::CalcLatitude () const
+{
+  return m_latitude + ((Simulator::Now ().GetDouble () / Hours (24).GetDouble ()) * 2 * M_PI);
+}
+
 Vector
 LeoCircularOrbitMobilityModel::CalcPosition (Time t) const
 {
-  Vector3D x = Product (m_orbitHeight, Vector3D (cos (m_inclination) * cos (m_latitude),
+  double lat = CalcLatitude ();
-  			       cos (m_inclination) * sin (m_latitude),
+  // account for orbit latitude and earth rotation offset
+  Vector3D x = Product (m_orbitHeight, Vector3D (cos (m_inclination) * cos (lat),
+  			       cos (m_inclination) * sin (lat),
   			       sin (m_inclination)));
 
   return RotatePlane (GetProgress (t), x);
@@ -131,8 +139,6 @@ LeoCircularOrbitMobilityModel::CalcPosition (Time t) const
 
 Vector LeoCircularOrbitMobilityModel::Update ()
 {
-  m_plane = PlaneNorm ();
-
   m_position = CalcPosition (Simulator::Now ());
   NotifyCourseChange ();
 
@@ -180,13 +186,13 @@ void LeoCircularOrbitMobilityModel::SetAltitude (double h)
 
 double LeoCircularOrbitMobilityModel::GetInclination () const
 {
-  return (m_inclination / (M_PI)) * 180.0;
+  return (m_inclination / M_PI) * 180.0;
 }
 
 void LeoCircularOrbitMobilityModel::SetInclination (double incl)
 {
   NS_ASSERT_MSG (incl != 0.0, "Plane must not be orthogonal to axis");
-  m_inclination = (incl / 180) * M_PI;
+  m_inclination = (incl / 180.0) * M_PI;
   Update ();
 }
 

model/leo-circular-orbit-mobility-model.h

@@ -65,11 +65,6 @@ private:
    */
   double m_offset;
 
-  /**
-   * Normal vector of orbital plane.
-   */
-  Vector3D m_plane;
-
   /**
    * Current position
    */
@@ -112,6 +107,12 @@ private:
    */
   Vector CalcPosition (Time t) const;
 
+  /**
+   * Calc the latitude depending on simulation time inside ITRF coordinate
+   * system
+   */
+  double CalcLatitude () const;
+
   Vector Update ();
 };
 

utils/plot-satellites.gnuplot

@@ -1,20 +1,31 @@
-set datafile separator comma
+set datafile separator ",:" 
 set key autotitle columnheader
-set terminal gif animate
+set terminal gif animate size 2000,2000
 set output 'output.gif'
 
 unset xtics
 unset ytics
 unset ztics
 unset border
+set xrange [-8e6:8e6]
+set yrange [-8e6:8e6]
+set zrange [-8e6:8e6]
+set parametric
+set isosamples 100,100 
 unset key
 set hidden3d
-set view equal xyz
 
 # number of nodes per time slot
-n=1200
+EARTH=6.370e6
+sats=ARG1
+ground=ARG2
+numsats=ARG3
+numsamples=ARG4
 
-do for [j=0:100] {
+do for [j=0:numsamples] {
 	set title 'time '.j
-	splot 'leo-circular-orbit-tracing-example.csv' using 3:4:5:2 every ::(j*n)::((j+1)*n)
+	splot [-pi:pi][-pi/2:pi/2] EARTH*cos(u)*cos(v), EARTH*sin(u)*cos(v), EARTH*sin(v), \
+	      ground using 1:2:3 lt rgb "green", \
+	      sats using 3:4:5:2 every ::(j*numsats)::((j+1)*numsats) lt rgb "blue"
+	      
 }