Thrusters

Makefile

@@ -2,4 +2,4 @@ CC=g++
 CFLAGS=--std=c++11 -lm
 
 simtest: simtest.cpp Vessel.cpp CelestialBody.cpp orb.hpp Universe.cpp
-	$(CC) $(CFLAGS) -o simtest simtest.cpp Vessel.cpp CelestialBody.cpp Universe.cpp orb.cpp
+	$(CC) $(CFLAGS) -o simtest simtest.cpp Vessel.cpp CelestialBody.cpp Universe.cpp orb.cpp Thruster.cpp

Thruster.cpp

@@ -0,0 +1,49 @@
+#include "Thruster.hpp"
+
+Thruster::Thruster(Vector3& pos, Vector3& dir, double max_thrust, double isp)
+{
+  this->pos = pos;
+  this->dir = dir;
+  this->max_thrust = max_thrust;
+  this->Isp = isp;
+}
+
+Thruster::~Thruster()
+{
+
+}
+
+Vector3 Thruster::getCurrentThrust()
+{
+  return (level*max_thrust/norm3(dir))*dir;
+}
+
+double Thruster::getCurrentThrustMag()
+{
+  return (level*max_thrust);
+}
+
+void Thruster::setLevel(double lvl)
+{
+  level = lvl;
+}
+
+double Thruster::getLevel()
+{
+  return level;
+}
+
+double Thruster::getIsp()
+{
+  return Isp;
+}
+
+Vector3 Thruster::getPosition()
+{
+  return this->pos;
+}
+
+Vector3 Thruster::getDirection()
+{
+  return this->dir;
+}

Thruster.hpp

@@ -0,0 +1,44 @@
+#ifndef THRUSTER_H
+#define THRUSTER_H
+
+#include "orb.hpp"
+#include <vector>
+
+enum ThrusterGroupType {
+  MAIN,
+  RETRO
+};
+
+struct ThrusterGroup {
+
+  ThrusterGroupType type;
+  std::vector<ThrusterHandle> thrusters;
+};
+
+class Thruster {
+
+  Vector3 pos;
+  Vector3 dir;
+
+  double max_thrust;
+  double Isp;
+  double level; //Thrust level between 0-1
+
+public:
+  Thruster(Vector3& pos, Vector3& dir, double max_thrust, double isp);
+  ~Thruster();
+
+// Thrust Handling
+  Vector3 getCurrentThrust();
+  double  getCurrentThrustMag();
+  double  getLevel();
+  void    setLevel(double lvl);
+  double  getIsp();
+
+// Force Vectors
+  Vector3 getPosition();
+  Vector3 getDirection();
+
+};
+
+#endif

Universe.cpp

@@ -7,7 +7,7 @@
 
 Universe::Universe(double init_mjd)
 {
-  this->mjd = init_mjd;
+  this->mjd  = init_mjd;
   this->simt = 0.0;
 }
 
@@ -20,7 +20,7 @@ void Universe::propagate(double step)
 {
   for(Vessel *v : vessel_list)
   {
-    v->preStep(mjd, simt, step);
+    v->update(mjd, simt, step);
     propagate_linear(v, step);
     integrate_rk4_angular(v, step);
   }
@@ -80,6 +80,11 @@ LinearStateVector Universe::eom_linear(double mjd, LinearStateVector state, Vect
 
 ObjHandle Universe::compute_gravity(double mjd, Vessel *ves, Vector3& G)
 {
+  /*
+   * Computes gravitational acceleration vector for Vessel ves 
+   * and places it in Vector3 G.  Returns handle to celestialBody
+   * providing greatest contribution to gravity vector
+   */
   Vector3 r, r_hat, delta_r, accel = {0.0, 0.0, 0.0}; /* object position */
   double r_mag = 0;
   double g_mag = 0;
@@ -91,9 +96,10 @@ ObjHandle Universe::compute_gravity(double mjd, Vessel *ves, Vector3& G)
   for(CelestialBody *cb : celbody_list)
   {
     delta_r = r + (-1.0)*cb->get_ephemeris(mjd);
-    r_mag = norm3(delta_r);
-    r_hat = (1.0/r_mag)*delta_r;
-    g_mag = ((cb->mu) / (r_mag*r_mag));
+    r_mag   = norm3(delta_r);
+    r_hat   = (1.0/r_mag)*delta_r;
+    g_mag   = ((cb->mu) / (r_mag*r_mag));
+
     if(g_mag > g_max)
     {
       g_max_body = cb;

Vessel.cpp

@@ -1,4 +1,6 @@
 #include "Vessel.hpp"
+#include "Thruster.hpp"
+#include <cmath>
 
 Vessel::Vessel()
 {
@@ -15,9 +17,18 @@ int Vessel::preStep(double mjd, double simt, double dt)
   return 0;
 }
 
-void update(double mjd, double simt, double dt)
+void Vessel::update(double mjd, double simt, double dt)
 {
-
+  // Decrease mass according to current thruster usage
+  double deltaM = 0;
+  for(Thruster *th : thrusters)
+  {
+    deltaM += dt*(th->getCurrentThrustMag()/(g0*th->getIsp()));
+  }
+  if(mass-empty_mass < deltaM) {
+    killThrust();
+  }
+  mass = fmax(mass-(empty_mass+deltaM), empty_mass);
 }
 
 void Vessel::addForce(const Vector3& F, const Vector3& pos)
@@ -36,6 +47,11 @@ void Vessel::getTotalForce(Vector3& F)
   {
     F = F + fd->F;
   }
+
+  for(Thruster *th : thrusters)
+  {
+    F = F + th->getCurrentThrust();
+  }
 }
 
 void Vessel::getTotalTorque(Vector3& T)
@@ -45,6 +61,11 @@ void Vessel::getTotalTorque(Vector3& T)
   {
     T = T + cross(fd->pos, fd->F);
   }
+
+  for(Thruster *th : thrusters)
+  {
+    T = T + cross(th->getPosition(), th->getCurrentThrust());
+  }
 }
 
 void Vessel::Local2Global(const Vector3& loc, Vector3& glob)
@@ -58,10 +79,61 @@ void Vessel::Local2Global(const Vector3& loc, Vector3& glob)
 
 double Vessel::interpTotalMass(double step)
 {
+  /* Interpolates the mass of the vehicle during timesteps
+  */
   return mass;
 }
 
 void Vessel::setReference(ObjHandle hRef)
 {
+  /* Sets the vessel reference body */
   ref = (CelestialBody *)hRef;
 }
+
+void Vessel::getLVLH(Vector3& X, Vector3& Y, Vector3& Z)
+{
+  /* Return East, North, Up unit vectors in (x, y, z)
+  representing the Local vertical local horizontal frame
+  Vectors written in global bases */
+}
+
+/* Thruster Management */
+ThrusterHandle Vessel::addThruster(Vector3& pos, Vector3& dir,
+                                   double max_thrust, double Isp)
+{
+  Thruster *th = new Thruster(pos, dir, max_thrust, Isp);
+  thrusters.push_back(th);
+  return (ThrusterHandle)th;
+}
+
+bool Vessel::delThruster(ThrusterHandle th)
+{
+  delete (Thruster *)th;
+  return 1;
+}
+
+void Vessel::killThrust()
+{
+  for(Thruster *th : thrusters)
+  {
+    th->setLevel(0.0);
+  }
+}
+
+// Thruster Group Management
+ThrusterGroup * Vessel::createThrusterGroup(ThrusterGroupType type)
+{
+  ThrusterGroup *th_grp = new ThrusterGroup;
+  th_grp->type = type;
+  return th_grp;
+}
+
+void Vessel::destroyThrusterGroup(ThrusterGroup *th_grp)
+{
+  delete th_grp;
+}
+
+void Vessel::addThrusterToGroup(ThrusterGroup *th_grp, ThrusterHandle th)
+{
+  th_grp->thrusters.push_back(th);
+}

Vessel.hpp

@@ -2,19 +2,25 @@
 #define VESSEL_H
 
 #include "orb.hpp"
+#include "Thruster.hpp"
 #include <vector>
 
+#define g0 9.81
+
 class CelestialBody;
 
+
 class Vessel {
 
 public:
 // Physical State Vectors
   LinearStateVector lin_state;
   AngularStateVector rot_state;
 
-// Forces
+// Forces + Thrusters
   std::vector<ForceData*> force_stack;
+  std::vector<Thruster*> thrusters;
+  std::vector<ThrusterGroup*> thruster_groups;
 
 // Internal handle and current reference body
   ObjHandle h;
@@ -44,7 +50,23 @@ class Vessel {
   void getTotalTorque(Vector3& T);
 
 // Coordinate Transforms
+
+  /* Body Coordinates to Global Frame */
   void Local2Global(const Vector3& loc, Vector3& glob);
+
+  /* Get LVLH Frame */
+  void getLVLH(Vector3& X, Vector3& Y, Vector3& Z);
+
+// Thruster Management
+  ThrusterHandle addThruster(Vector3& pos, Vector3& dir,
+                             double max_thrust, double Isp);
+  bool delThruster(ThrusterHandle th);
+  void killThrust(); /* Kills all currently active thrusters */
+
+// Thruster Groups
+  ThrusterGroup* createThrusterGroup(ThrusterGroupType type);
+  void destroyThrusterGroup(ThrusterGroup *th_grp);
+  void addThrusterToGroup(ThrusterGroup *th_grp, ThrusterHandle th);
 };
 
 #endif

orb.hpp

@@ -2,7 +2,8 @@
 #define ORB_H
 
 #include <cmath>
-typedef void *ObjHandle;
+typedef void * ObjHandle;
+typedef void * ThrusterHandle;
 
 const double pi = 3.14159;
 typedef struct Vector3 {

simtest.cpp

@@ -19,6 +19,7 @@ int main()
   Vessel *Probe = new Vessel();
 
   Probe->mass = 1000.0;
+  Probe->empty_mass = 10.0;
   LinearStateVector ProbeInitState;
   ProbeInitState.r.x = 6871.0e3 + 400e3;
   ProbeInitState.r.y = 0.0;