diff --git a/Makefile b/Makefile index bab3d1f..2aa80a2 100644 --- a/Makefile +++ b/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 diff --git a/Thruster.cpp b/Thruster.cpp new file mode 100644 index 0000000..4f89006 --- /dev/null +++ b/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; +} diff --git a/Thruster.hpp b/Thruster.hpp new file mode 100644 index 0000000..3107cb1 --- /dev/null +++ b/Thruster.hpp @@ -0,0 +1,44 @@ +#ifndef THRUSTER_H +#define THRUSTER_H + +#include "orb.hpp" +#include + +enum ThrusterGroupType { + MAIN, + RETRO +}; + +struct ThrusterGroup { + + ThrusterGroupType type; + std::vector 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 diff --git a/Universe.cpp b/Universe.cpp index 8e1d33b..57426ae 100644 --- a/Universe.cpp +++ b/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; diff --git a/Vessel.cpp b/Vessel.cpp index 585d22c..ce8c785 100644 --- a/Vessel.cpp +++ b/Vessel.cpp @@ -1,4 +1,6 @@ #include "Vessel.hpp" +#include "Thruster.hpp" +#include 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); +} diff --git a/Vessel.hpp b/Vessel.hpp index 72f0c2b..2a4b759 100644 --- a/Vessel.hpp +++ b/Vessel.hpp @@ -2,10 +2,14 @@ #define VESSEL_H #include "orb.hpp" +#include "Thruster.hpp" #include +#define g0 9.81 + class CelestialBody; + class Vessel { public: @@ -13,8 +17,10 @@ class Vessel { LinearStateVector lin_state; AngularStateVector rot_state; -// Forces +// Forces + Thrusters std::vector force_stack; + std::vector thrusters; + std::vector 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 diff --git a/orb.hpp b/orb.hpp index 70b698e..8b3aca7 100644 --- a/orb.hpp +++ b/orb.hpp @@ -2,7 +2,8 @@ #define ORB_H #include -typedef void *ObjHandle; +typedef void * ObjHandle; +typedef void * ThrusterHandle; const double pi = 3.14159; typedef struct Vector3 { diff --git a/simtest.cpp b/simtest.cpp index 313a011..dd02f53 100644 --- a/simtest.cpp +++ b/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;