5 changed files with 178 additions and 68 deletions
@ -1,26 +1,44 @@
|
||||
|
||||
OBJDIR = build
|
||||
SRCDIR = src
|
||||
LIBDIR = ext/tangerine
|
||||
LIB = $(LIBDIR)/build/libTangerine.a
|
||||
CXX = g++
|
||||
CXXFLAGS = -std=c++11 -g -ggdb3 -Wall -I$(LIBDIR)/include/
|
||||
LDFLAGS = -lSDL2 -lGLEW -lGL -lassimp
|
||||
BIN = bin/orbital_shipping
|
||||
OBJDIR := build
|
||||
SRCDIR := src
|
||||
LIBDIR := ext/tangerine
|
||||
TESTDIR := tests
|
||||
BINDIR := bin
|
||||
LIB := $(LIBDIR)/build/libTangerine.a
|
||||
|
||||
SOURCES = $(wildcard $(SRCDIR)/*.cpp)
|
||||
OBJECTS = $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(SOURCES))
|
||||
CXX := g++
|
||||
CXXFLAGS := -std=c++11 -g -ggdb3 -Wall -I$(LIBDIR)/include/
|
||||
LDFLAGS := -lSDL2 -lGLEW -lGL -lassimp
|
||||
|
||||
COMMON_SOURCES := $(SRCDIR)/orbits.cpp
|
||||
COMMON_OBJECTS := $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(COMMON_SOURCES))
|
||||
|
||||
all: $(OBJECTS) |
||||
$(CXX) -o $(BIN) $(LDFLAGS) $(OBJDIR)/*.o $(LIB)
|
||||
TEST_SOURCES := $(TESTDIR)/orbit_test.cpp
|
||||
TEST_OBJECTS := $(patsubst $(TESTDIR)/%.cpp, $(OBJDIR)/%.o, $(TEST_SOURCES))
|
||||
TEST_BIN := $(BINDIR)/test_orbit
|
||||
|
||||
SOURCES := $(SRCDIR)/main.cpp
|
||||
OBJECTS := $(patsubst $(SRCDIR)/%.cpp, $(OBJDIR)/%.o, $(SOURCES))
|
||||
BIN := $(BINDIR)/orbital_shipping
|
||||
|
||||
|
||||
all: $(BIN) tests |
||||
.PHONY: all |
||||
|
||||
-include $(OBJDIR)/*.d |
||||
|
||||
$(OBJECTS): $(OBJDIR)/%.o : $(SRCDIR)/%.cpp |
||||
$(CXX) $(CXXFLAGS) -c -MMD $< -o $@
|
||||
$(BIN): $(COMMON_OBJECTS) |
||||
$(CXX) $(CXXFLAGS) -c -MMD $(SOURCES) -o $(OBJECTS)
|
||||
$(CXX) -o $(BIN) $(LDFLAGS) $(COMMON_OBJECTS) $(OBJECTS) $(LIB)
|
||||
|
||||
tests: $(COMMON_OBJECTS) |
||||
$(CXX) $(CXXFLAGS) -c -MMD $(TEST_SOURCES) -o $(TEST_OBJECTS)
|
||||
$(CXX) -o $(TEST_BIN) $(LDFLAGS) $(COMMON_OBJECTS) $(TEST_OBJECTS) $(LIB)
|
||||
.PHONY: tests |
||||
|
||||
$(COMMON_OBJECTS): $(COMMON_SOURCES) |
||||
$(CXX) $(CXXFLAGS) -c -MMD $(SRCDIR)/orbits.cpp -o $@
|
||||
|
||||
clean: |
||||
rm -rf $(OBJDIR)/* $(BIN)
|
||||
rm -rf $(OBJDIR)/* $(BIN) bin/test_orbit
|
||||
.PHONY: clean |
||||
|
||||
@ -0,0 +1,97 @@
|
||||
|
||||
#include "dumbLog.h" |
||||
|
||||
#include "orbits.h" |
||||
|
||||
|
||||
double |
||||
getPropagatedPosition(orbital_elements orbit, |
||||
double initial_anom, |
||||
unsigned int time_step) |
||||
{ |
||||
ellipse_parameters ep = orbit.ep; |
||||
|
||||
// TODO: break out steps with descriptive function names
|
||||
// initial eccentric anomaly
|
||||
double E1 = 2 * atan(sqrt((1 - ep.e) / (1 + ep.e)) * tan(initial_anom / 2)); |
||||
|
||||
// initial mean anomaly
|
||||
double M1 = E1 - ep.e * sin(E1); |
||||
|
||||
// mean motion
|
||||
double n = sqrt(orbit.mu / pow(ep.a, 3)); |
||||
|
||||
// propagated mean anomaly
|
||||
double M2 = M1 + n * (time_step); |
||||
|
||||
// guess starting value for E2
|
||||
double E2_1 = M2 + ep.e * sin(M2) + ((pow(ep.e, 2) / 2) * sin(2 * M2)); |
||||
|
||||
// test if guess is a solution to kepler's equation
|
||||
const double ACCEPTABLE_ERROR = 0.00000001; |
||||
double E2_test = E2_1; |
||||
|
||||
for (uint i = 0; i < 10; i++) { |
||||
double err = E2_test - ep.e * sin(E2_test) - M2; |
||||
|
||||
LOG(Debug) << "-------------------------\n"; |
||||
LOG(Debug) << "index: " << i << "\n"; |
||||
LOG(Debug) << "E2_test: " << E2_test << "\n"; |
||||
LOG(Debug) << "err: " << err << "\n"; |
||||
LOG(Debug) << "err: " << err << "\n"; |
||||
|
||||
if (fabs(err) < ACCEPTABLE_ERROR) |
||||
break; |
||||
|
||||
// compute derivative of the error function
|
||||
double derr = 1 - ep.e * cos(E2_test); |
||||
LOG(Debug) << "derr: " << derr << "\n"; |
||||
|
||||
// use Newton's method to compute next trial value of E2
|
||||
E2_test = E2_test - (err / derr); |
||||
} |
||||
|
||||
return E2_test; |
||||
} |
||||
|
||||
double |
||||
testTrialEccAnomaly(double trial_value, double e, double M2) |
||||
{ |
||||
double err = trial_value - e * sin(trial_value) - M2; |
||||
return err; |
||||
} |
||||
|
||||
ellipse_parameters |
||||
constructEllipse(double a, double b) |
||||
{ |
||||
assert(a > 0 && b > 0 && a >= b); |
||||
ellipse_parameters ep = { a, b }; |
||||
ep.c = sqrt(a * a - b * b); |
||||
ep.e = ep.c / ep.a; |
||||
ep.f1.x = -1 * ep.c; |
||||
ep.f2.x = ep.c; |
||||
return ep; |
||||
} |
||||
|
||||
ellipse_3d |
||||
constructEllipse3D(ellipse_parameters ep, uint vert_count) |
||||
{ |
||||
assert(ep.a > 0 && ep.b > 0 && |
||||
ep.a >= ep.b && |
||||
vert_count > 0); |
||||
|
||||
ellipse_3d e3d = { ep, nullptr, vert_count}; |
||||
e3d.vertices = UTIL_ALLOC(vert_count, glm::vec3); |
||||
double angle = 2 * M_PI / vert_count; |
||||
|
||||
for (uint i = 0; i < vert_count; i++) { |
||||
double a = angle * i; |
||||
// NOTE: solving for distance in polar coordinates
|
||||
double r = ep.b / sqrt(1 - (ep.e * pow(cos(a), 2))); |
||||
// NOTE: converting from polar to rectangular
|
||||
e3d.vertices[i] = glm::vec3(cos(a) * r, sin(a) * r, 0); |
||||
} |
||||
|
||||
return e3d; |
||||
} |
||||
|
||||
@ -0,0 +1,34 @@
|
||||
|
||||
#include "dumbLog.h" |
||||
|
||||
#include "../src/orbits.h" |
||||
|
||||
|
||||
void |
||||
testOrbitPropagation() |
||||
{ |
||||
orbital_elements orbit = {}; |
||||
ellipse_parameters ep = {}; |
||||
ep.a = 26564.5; // meters
|
||||
ep.e = 0.7411; |
||||
orbit.ep = ep; |
||||
|
||||
// NOTE: solve for gravitational parameter from mean motion
|
||||
double mean_motion = 0.00014582; // rad/s
|
||||
orbit.mu = pow(mean_motion, 2) * pow(ep.a, 3); |
||||
|
||||
double initial_anom = 260 * M_PI / 180; // radians
|
||||
unsigned int time_step = 60 * 50; // seconds
|
||||
double next_pos = getPropagatedPosition(orbit, initial_anom, time_step); |
||||
|
||||
LOG(Debug) << "E2_1: " << next_pos << "\n"; |
||||
|
||||
return; |
||||
} |
||||
|
||||
int |
||||
main() |
||||
{ |
||||
testOrbitPropagation(); |
||||
return 0; |
||||
} |
||||
Loading…
Reference in new issue