clean up the interface of the Orbit structure

master
John Bartholomew 2013-03-29 16:10:03 +00:00
parent 3df2beea35
commit f0e5c25745
6 changed files with 169 additions and 141 deletions

View File

@ -6,18 +6,65 @@
#include "Game.h"
#include "Pi.h"
double Orbit::MeanAnomalyAtTime(double time) const {
const double e = eccentricity;
if(e >= 0 && e < 1) { // elliptic orbit
return 2.0*M_PI*time / Period() + orbitalPhaseAtStart; // mean anomaly
static double calc_orbital_period(double semiMajorAxis, double centralMass)
{
return 2.0*M_PI*sqrt((semiMajorAxis*semiMajorAxis*semiMajorAxis)/(G*centralMass));
}
static double calc_orbital_period_gravpoint(double semiMajorAxis, double totalMass, double bodyMass)
{
// variable names according to the formula in:
// http://en.wikipedia.org/wiki/Barycentric_coordinates_(astronomy)#Two-body_problem
//
// We have a 2-body orbital system, represented as a gravpoint (at the barycentre),
// plus two bodies, each orbiting that gravpoint.
// We need to compute the orbital period, given the semi-major axis of one body's orbit
// around the gravpoint, the total mass of the system, and the mass of the body.
//
// According to Kepler, the orbital period P is defined by:
//
// P = 2*pi * sqrt( a**3 / G*(M1 + M2) )
//
// where a is the semi-major axis of the orbit, M1 is the mass of the primary and M2 is
// the mass of the secondary. But we don't have that semi-major axis value, we have the
// the semi-major axis for the orbit of the secondary around the gravpoint, instead.
//
// So, this first computes the semi-major axis of the secondary's orbit around the primary,
// and then uses the above formula to compute the orbital period.
const double r1 = semiMajorAxis;
const double m2 = (totalMass - bodyMass);
const double a = r1 * totalMass / m2;
const double a3 = a*a*a;
return 2.0 * M_PI * sqrt(a3 / (G * totalMass));
}
static double calc_velocity_area_per_sec(double semiMajorAxis, double centralMass, double eccentricity) {
if(eccentricity < 1)
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(1 - eccentricity * eccentricity)/ calc_orbital_period(semiMajorAxis, centralMass);
else
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(eccentricity * eccentricity - 1)/ calc_orbital_period(semiMajorAxis, centralMass);
}
static double calc_velocity_area_per_sec_gravpoint(double semiMajorAxis, double totalMass, double bodyMass, double eccentricity) {
if(eccentricity < 1) {
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(1 - eccentricity * eccentricity)/ calc_orbital_period_gravpoint(semiMajorAxis, totalMass, bodyMass);
} else {
return -2 * time * velocityAreaPerSecond / semiMajorAxis / semiMajorAxis / sqrt(e*e-1) + orbitalPhaseAtStart; // mean anomaly
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(eccentricity * eccentricity - 1)/ calc_orbital_period_gravpoint(semiMajorAxis, totalMass, bodyMass);
}
}
double Orbit::MeanAnomalyAtTime(double time) const {
const double e = m_eccentricity;
if(e >= 0 && e < 1) { // elliptic orbit
return 2.0*M_PI*time / Period() + m_orbitalPhaseAtStart; // mean anomaly
} else {
return -2 * time * m_velocityAreaPerSecond / m_semiMajorAxis / m_semiMajorAxis / sqrt(e*e-1) + m_orbitalPhaseAtStart; // mean anomaly
}
}
vector3d Orbit::OrbitalPosAtTime(double t) const
{
const double e = eccentricity;
const double e = m_eccentricity;
double r = 0, cos_v = 0, sin_v = 0;
const double M = MeanAnomalyAtTime(t); // mean anomaly
@ -29,7 +76,7 @@ vector3d Orbit::OrbitalPosAtTime(double t) const
E = E - (E-e*(sin(E))-M) / (1.0 - e*cos(E));
}
// heliocentric distance
r = semiMajorAxis * (1.0 - e*cos(E));
r = m_semiMajorAxis * (1.0 - e*cos(E));
// true anomaly (angle of orbit position)
cos_v = (cos(E) - e) / (1.0 - e*cos(E));
sin_v = (sqrt(1.0-e*e)*sin(E))/ (1.0 - e*cos(E));
@ -46,7 +93,7 @@ vector3d Orbit::OrbitalPosAtTime(double t) const
ch = sqrt(1 + sh*sh);
// heliocentric distance
r = semiMajorAxis * (e*ch - 1.0);
r = m_semiMajorAxis * (e*ch - 1.0);
// true anomaly (angle of orbit position)
cos_v = (ch - e) / (1.0 - e*ch);
@ -54,7 +101,7 @@ vector3d Orbit::OrbitalPosAtTime(double t) const
}
vector3d pos = vector3d(-cos_v*r, sin_v*r, 0);
pos = rotMatrix * pos;
pos = m_orient * pos;
return pos;
}
@ -63,16 +110,16 @@ vector3d Orbit::OrbitalPosAtTime(double t) const
// to be taken into account
vector3d Orbit::EvenSpacedPosTrajectory(double angle) const
{
const double e = eccentricity;
const double e = m_eccentricity;
vector3d pos = vector3d(0.0f,0.0f,0.0f);
if(e < 1) {
double v = 2*M_PI*angle +TrueAnomaly(MeanAnomalyAtTime(Pi::game->GetTime()));
const double r = semiMajorAxis * (1 - e*e) / (1 + e*cos(v));
const double r = m_semiMajorAxis * (1 - e*e) / (1 + e*cos(v));
pos = vector3d(-cos(v)*r, sin(v)*r, 0);
} else {
double v = 2*M_PI*angle +TrueAnomaly(MeanAnomalyAtTime(Pi::game->GetTime()));
double r = semiMajorAxis * (e*e - 1) / (1 + e*cos(v));
double r = m_semiMajorAxis * (e*e - 1) / (1 + e*cos(v));
// planet is in infinity
if(v <= - acos(-1/e)) {
@ -86,13 +133,13 @@ vector3d Orbit::EvenSpacedPosTrajectory(double angle) const
pos = vector3d(-cos(v)*r, sin(v)*r, 0);
}
pos = rotMatrix * pos;
pos = m_orient * pos;
return pos;
}
double Orbit::Period() const {
if(eccentricity < 1 && eccentricity >= 0) {
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(1 - eccentricity * eccentricity)/ velocityAreaPerSecond;
if(m_eccentricity < 1 && m_eccentricity >= 0) {
return M_PI * m_semiMajorAxis * m_semiMajorAxis * sqrt(1 - m_eccentricity * m_eccentricity)/ m_velocityAreaPerSecond;
} else { // hyperbola.. period makes no sense, should not be used
assert(0);
return 0;
@ -100,7 +147,7 @@ double Orbit::Period() const {
}
double Orbit::TrueAnomaly(double MeanAnomaly) const {
const double e = eccentricity, M = MeanAnomaly;
const double e = m_eccentricity, M = MeanAnomaly;
double cos_v, sin_v, v;
if(e >= 0 && e < 1) { // elliptic orbit
// eccentric anomaly
@ -136,7 +183,7 @@ double Orbit::TrueAnomaly(double MeanAnomaly) const {
double Orbit::MeanAnomalyFromTrueAnomaly(double trueAnomaly) const {
double M_t0;
const double e = eccentricity;
const double e = m_eccentricity;
if(e > 0 && e < 1) {
M_t0 = 2*atan(tan(trueAnomaly/2)*sqrt((1-e)/(1+e)));
M_t0 = M_t0 - e*sin(M_t0);
@ -146,73 +193,40 @@ double Orbit::MeanAnomalyFromTrueAnomaly(double trueAnomaly) const {
// in time decrease their true anomaly. Yes, it is confusing.
M_t0 = 2*atanh(tan(trueAnomaly/2)*sqrt((e-1)/(1+e)));
M_t0 = M_t0 - e*sinh(M_t0);
M_t0 += Pi::game->GetTime() * 2 * velocityAreaPerSecond / semiMajorAxis / semiMajorAxis / sqrt(e*e-1);
M_t0 += Pi::game->GetTime() * 2 * m_velocityAreaPerSecond / m_semiMajorAxis / m_semiMajorAxis / sqrt(e*e-1);
}
return M_t0;
}
vector3d Orbit::Apogeum() const {
if(eccentricity < 1) {
return semiMajorAxis * (1 + eccentricity) * (rotMatrix * vector3d(1,0,0));
if(m_eccentricity < 1) {
return m_semiMajorAxis * (1 + m_eccentricity) * (m_orient * vector3d(1,0,0));
} else {
return vector3d(0,0,0);
}
}
vector3d Orbit::Perigeum() const {
if(eccentricity < 1) {
return semiMajorAxis * (1 - eccentricity) * (rotMatrix * vector3d(-1,0,0));
if(m_eccentricity < 1) {
return m_semiMajorAxis * (1 - m_eccentricity) * (m_orient * vector3d(-1,0,0));
} else {
return semiMajorAxis * (eccentricity - 1) * (rotMatrix * vector3d(-1,0,0));
return m_semiMajorAxis * (m_eccentricity - 1) * (m_orient * vector3d(-1,0,0));
}
}
double Orbit::calc_velocity_area_per_sec(double semiMajorAxis, double centralMass, double eccentricity) {
if(eccentricity < 1)
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(1 - eccentricity * eccentricity)/ calc_orbital_period(semiMajorAxis, centralMass);
else
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(eccentricity * eccentricity - 1)/ calc_orbital_period(semiMajorAxis, centralMass);
}
double Orbit::calc_velocity_area_per_sec_gravpoint(double semiMajorAxis, double totalMass, double bodyMass, double eccentricity) {
if(eccentricity < 1) {
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(1 - eccentricity * eccentricity)/ calc_orbital_period_gravpoint(semiMajorAxis, totalMass, bodyMass);
} else {
return M_PI * semiMajorAxis * semiMajorAxis * sqrt(eccentricity * eccentricity - 1)/ calc_orbital_period_gravpoint(semiMajorAxis, totalMass, bodyMass);
}
}
double Orbit::calc_orbital_period(double semiMajorAxis, double centralMass)
void Orbit::SetShapeAroundBarycentre(double semiMajorAxis, double totalMass, double bodyMass, double eccentricity)
{
return 2.0*M_PI*sqrt((semiMajorAxis*semiMajorAxis*semiMajorAxis)/(G*centralMass));
m_semiMajorAxis = semiMajorAxis;
m_eccentricity = eccentricity;
m_velocityAreaPerSecond = calc_velocity_area_per_sec_gravpoint(semiMajorAxis, totalMass, bodyMass, eccentricity);
}
double Orbit::calc_orbital_period_gravpoint(double semiMajorAxis, double totalMass, double bodyMass)
void Orbit::SetShapeAroundPrimary(double semiMajorAxis, double centralMass, double eccentricity)
{
// variable names according to the formula in:
// http://en.wikipedia.org/wiki/Barycentric_coordinates_(astronomy)#Two-body_problem
//
// We have a 2-body orbital system, represented as a gravpoint (at the barycentre),
// plus two bodies, each orbiting that gravpoint.
// We need to compute the orbital period, given the semi-major axis of one body's orbit
// around the gravpoint, the total mass of the system, and the mass of the body.
//
// According to Kepler, the orbital period P is defined by:
//
// P = 2*pi * sqrt( a**3 / G*(M1 + M2) )
//
// where a is the semi-major axis of the orbit, M1 is the mass of the primary and M2 is
// the mass of the secondary. But we don't have that semi-major axis value, we have the
// the semi-major axis for the orbit of the secondary around the gravpoint, instead.
//
// So, this first computes the semi-major axis of the secondary's orbit around the primary,
// and then uses the above formula to compute the orbital period.
const double r1 = semiMajorAxis;
const double m2 = (totalMass - bodyMass);
const double a = r1 * totalMass / m2;
const double a3 = a*a*a;
return 2.0 * M_PI * sqrt(a3 / (G * totalMass));
m_semiMajorAxis = semiMajorAxis;
m_eccentricity = eccentricity;
m_velocityAreaPerSecond = calc_velocity_area_per_sec(semiMajorAxis, centralMass, eccentricity);
}
Orbit Orbit::FromBodyState(const vector3d &pos, const vector3d &vel, double mass)
@ -230,29 +244,29 @@ Orbit Orbit::FromBodyState(const vector3d &pos, const vector3d &vel, double mass
double EE = vel.LengthSqr()/2 - mass*G/r_now;
if (mass <= 1e-3 || r_now <= 1e-3 || v_now <= 1e-3 || fabs(EE) <= 1e-12 || 1.0-ang.z*ang.z/LL/LL < 0) {
ret.eccentricity = 0;
ret.semiMajorAxis = 0;
ret.velocityAreaPerSecond = 0;
ret.orbitalPhaseAtStart = 0;
ret.rotMatrix = matrix3x3d::Identity();
ret.m_eccentricity = 0;
ret.m_semiMajorAxis = 0;
ret.m_velocityAreaPerSecond = 0;
ret.m_orbitalPhaseAtStart = 0;
ret.m_orient = matrix3x3d::Identity();
return ret;
}
// http://en.wikipedia.org/wiki/Orbital_eccentricity
ret.eccentricity = 1 + 2*EE*LL*LL/pow(mass*G, 2);
if (ret.eccentricity < 1e-12) ret.eccentricity = 1e-12;
if (ret.eccentricity == 1.0) ret.eccentricity = 1-1e-6;
ret.eccentricity = sqrt(ret.eccentricity);
ret.m_eccentricity = 1 + 2*EE*LL*LL/pow(mass*G, 2);
if (ret.m_eccentricity < 1e-12) ret.m_eccentricity = 1e-12;
if (ret.m_eccentricity == 1.0) ret.m_eccentricity = 1-1e-6;
ret.m_eccentricity = sqrt(ret.m_eccentricity);
// lines represent these quantities:
// (e M G)^2
// M G (e - 1) / 2 EE, always positive (EE and (e-1) change sign
// M G / 2 EE,
// which is a (http://en.wikipedia.org/wiki/Semi-major_axis); a of hyperbola is taken as positive here
ret.semiMajorAxis = 2*EE*LL*LL + pow(mass*G, 2);
if (ret.semiMajorAxis < 0) ret.semiMajorAxis = 0;
ret.semiMajorAxis = (sqrt(ret.semiMajorAxis ) - mass*G)/2/EE;
ret.semiMajorAxis = ret.semiMajorAxis/fabs(1.0-ret.eccentricity);
ret.m_semiMajorAxis = 2*EE*LL*LL + pow(mass*G, 2);
if (ret.m_semiMajorAxis < 0) ret.m_semiMajorAxis = 0;
ret.m_semiMajorAxis = (sqrt(ret.m_semiMajorAxis ) - mass*G)/2/EE;
ret.m_semiMajorAxis = ret.m_semiMajorAxis/fabs(1.0-ret.m_eccentricity);
// The formulas for rotation matrix were derived based on following assumptions:
// 1. Trajectory follows Kepler's law and vector {-r cos(v), -r sin(v), 0}, r(t) and v(t) are parameters.
@ -269,14 +283,14 @@ Orbit Orbit::FromBodyState(const vector3d &pos, const vector3d &vel, double mass
double off = 0, ccc = 0;
matrix3x3d mat;
if (ret.eccentricity < 1) {
off = ret.semiMajorAxis*(1 - ret.eccentricity*ret.eccentricity) - r_now;
if (ret.m_eccentricity < 1) {
off = ret.m_semiMajorAxis*(1 - ret.m_eccentricity*ret.m_eccentricity) - r_now;
} else {
off = ret.semiMajorAxis*(-1 + ret.eccentricity*ret.eccentricity) - r_now;
off = ret.m_semiMajorAxis*(-1 + ret.m_eccentricity*ret.m_eccentricity) - r_now;
}
// correct sign of offset is given by sign pos.Dot(vel) (heading towards apohelion or perihelion?]
off = Clamp(off/(r_now * ret.eccentricity), -1 + 1e-6,1 - 1e-6);
off = Clamp(off/(r_now * ret.m_eccentricity), -1 + 1e-6,1 - 1e-6);
off = -pos.Dot(vel)/fabs(pos.Dot(vel))*acos(off);
ccc = acos(-pos.z/r_now/sin(angle1)) * i;
@ -291,10 +305,10 @@ Orbit Orbit::FromBodyState(const vector3d &pos, const vector3d &vel, double mass
// matrix3x3d::RotateX(M_PI) and minus sign before offset changes solution above, derived for orbits {-r cos(v), -r sin(v), 0}
// to {-r cos(v), -r sin(v), 0}
ret.rotMatrix = matrix3x3d::RotateZ(angle2) * matrix3x3d::RotateY(angle1) * matrix3x3d::RotateZ(cc - offset) * matrix3x3d::RotateX(M_PI);
ret.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec(ret.semiMajorAxis, mass,ret.eccentricity);
ret.m_orient = matrix3x3d::RotateZ(angle2) * matrix3x3d::RotateY(angle1) * matrix3x3d::RotateZ(cc - offset) * matrix3x3d::RotateX(M_PI);
ret.m_velocityAreaPerSecond = calc_velocity_area_per_sec(ret.m_semiMajorAxis, mass,ret.m_eccentricity);
ret.orbitalPhaseAtStart = ret.MeanAnomalyFromTrueAnomaly(-offset);
ret.m_orbitalPhaseAtStart = ret.MeanAnomalyFromTrueAnomaly(-offset);
return ret;
}

View File

@ -7,18 +7,28 @@
#include "vector3.h"
#include "matrix3x3.h"
struct Orbit {
class Orbit {
public:
// note: the resulting Orbit is at the given position at t=0
static Orbit FromBodyState(const vector3d &position, const vector3d &velocity, double central_mass);
Orbit(): orbitalPhaseAtStart(0.0) {};
Orbit():
m_eccentricity(0.0),
m_semiMajorAxis(0.0),
m_orbitalPhaseAtStart(0.0),
m_velocityAreaPerSecond(0.0),
m_orient(matrix3x3d::Identity())
{}
void SetShapeAroundBarycentre(double semiMajorAxis, double totalMass, double bodyMass, double eccentricity);
void SetShapeAroundPrimary(double semiMajorAxis, double totalMass, double eccentricity);
void SetPlane(const matrix3x3d &orient) { m_orient = orient; }
void SetPhase(double orbitalPhaseAtStart) { m_orbitalPhaseAtStart = orbitalPhaseAtStart; }
vector3d OrbitalPosAtTime(double t) const;
// 0.0 <= t <= 1.0. Not for finding orbital pos
vector3d EvenSpacedPosTrajectory(double angle) const;
/* duplicated from SystemBody... should remove probably */
static double calc_orbital_period(double semiMajorAxis, double centralMass);
static double calc_orbital_period_gravpoint(double semiMajorAxis, double totalMass, double bodyMass);
static double calc_velocity_area_per_sec(double semiMajorAxis, double centralMass, double eccentricity);
static double calc_velocity_area_per_sec_gravpoint(double semiMajorAxis, double totalMass, double bodyMass, double eccentricity);
vector3d EvenSpacedPosTrajectory(double t) const;
double Period() const;
double TrueAnomaly(double MeanAnomaly) const;
@ -26,12 +36,20 @@ struct Orbit {
double MeanAnomalyAtTime(double time) const;
vector3d Apogeum() const;
vector3d Perigeum() const;
double eccentricity;
double semiMajorAxis;
double orbitalPhaseAtStart; // 0 to 2 pi radians
// basic accessors
double GetEccentricity() const { return m_eccentricity; }
double GetSemiMajorAxis() const { return m_semiMajorAxis; }
double GetOrbitalPhaseAtStart() const { return m_orbitalPhaseAtStart; }
const matrix3x3d &GetPlane() const { return m_orient; }
private:
double m_eccentricity;
double m_semiMajorAxis;
double m_orbitalPhaseAtStart; // 0 to 2 pi radians
/* dup " " --------------------------------------- */
double velocityAreaPerSecond; // seconds
matrix3x3d rotMatrix;
double m_velocityAreaPerSecond; // seconds
matrix3x3d m_orient;
};
#endif

View File

@ -371,7 +371,7 @@ static void RelocateStarportIfUnderwaterOrBuried(SystemBody *sbody, Frame *frame
const double radius = planet->GetSystemBody()->GetRadius();
// suggested position
rot = sbody->orbit.rotMatrix;
rot = sbody->orbit.GetPlane();
pos = rot * vector3d(0,1,0);
// Check if height varies too much around the starport center
@ -567,7 +567,7 @@ static Frame *MakeFrameFor(SystemBody *sbody, Body *b, Frame *f)
vector3d pos;
Planet *planet = static_cast<Planet*>(rotFrame->GetBody());
RelocateStarportIfUnderwaterOrBuried(sbody, rotFrame, planet, pos, rot);
sbody->orbit.rotMatrix = rot;
sbody->orbit.SetPlane(rot);
b->SetPosition(pos * planet->GetTerrainHeight(pos));
b->SetOrient(rot);
return rotFrame;

View File

@ -88,7 +88,7 @@ void SystemInfoView::OnBodyViewed(SystemBody *b)
_add_label_and_value(Lang::ORBITAL_PERIOD, data);
_add_label_and_value(Lang::PERIAPSIS_DISTANCE, format_distance(b->orbMin.ToDouble()*AU, 3));
_add_label_and_value(Lang::APOAPSIS_DISTANCE, format_distance(b->orbMax.ToDouble()*AU, 3));
_add_label_and_value(Lang::ECCENTRICITY, stringf("%0{f.2}", b->orbit.eccentricity));
_add_label_and_value(Lang::ECCENTRICITY, stringf("%0{f.2}", b->orbit.GetEccentricity()));
if (b->type != SystemBody::TYPE_STARPORT_ORBITAL) {
_add_label_and_value(Lang::AXIAL_TILT, stringf(Lang::N_DEGREES, formatarg("angle", b->axialTilt.ToDouble() * (180.0/M_PI))));
if (b->rotationPeriod != 0) {

View File

@ -146,10 +146,10 @@ void SystemView::PutOrbit(const Orbit *orbit, const vector3d &offset, const Colo
if (num_vertices > 1) {
// don't close the loop for hyperbolas and parabolas and crashed ellipses
if ((orbit->eccentricity <= 1.0) || (num_vertices < int(COUNTOF(vts))))
m_renderer->DrawLines(num_vertices, vts, color, LINE_LOOP);
else
if ((orbit->GetEccentricity() > 1.0) || (num_vertices < int(COUNTOF(vts))))
m_renderer->DrawLines(num_vertices, vts, color, LINE_STRIP);
else
m_renderer->DrawLines(num_vertices, vts, color, LINE_LOOP);
}
}
@ -175,7 +175,7 @@ void SystemView::OnClickObject(const SystemBody *b)
if (b->parent) {
desc += std::string(Lang::SEMI_MAJOR_AXIS);
desc += ":\n";
data += format_distance(b->orbit.semiMajorAxis)+"\n";
data += format_distance(b->orbit.GetSemiMajorAxis())+"\n";
desc += std::string(Lang::ORBITAL_PERIOD);
desc += ":\n";
@ -247,8 +247,8 @@ void SystemView::PutBody(const SystemBody *b, const vector3d &offset, const matr
if (b->children.size()) {
for(std::vector<SystemBody*>::const_iterator kid = b->children.begin(); kid != b->children.end(); ++kid) {
if (is_zero_general((*kid)->orbit.semiMajorAxis)) continue;
if ((*kid)->orbit.semiMajorAxis * m_zoom < ROUGH_SIZE_OF_TURD) {
if (is_zero_general((*kid)->orbit.GetSemiMajorAxis())) continue;
if ((*kid)->orbit.GetSemiMajorAxis() * m_zoom < ROUGH_SIZE_OF_TURD) {
PutOrbit(&((*kid)->orbit), offset, Color(0.f, 1.f, 0.f, 1.f));
}

View File

@ -783,8 +783,7 @@ static void position_settlement_on_planet(SystemBody *b)
// used for orientation on planet surface
double r2 = r.Double(); // function parameter evaluation order is implementation-dependent
double r1 = r.Double(); // can't put two rands in the same expression
b->orbit.rotMatrix = matrix3x3d::RotateZ(2*M_PI*r1) *
matrix3x3d::RotateY(2*M_PI*r2);
b->orbit.SetPlane(matrix3x3d::RotateZ(2*M_PI*r1) * matrix3x3d::RotateY(2*M_PI*r2));
}
double SystemBody::GetMaxChildOrbitalDistance() const
@ -929,26 +928,27 @@ void StarSystem::CustomGetKidsOf(SystemBody *parent, const std::vector<CustomSys
kid->orbitalPhaseAtStart = csbody->orbitalPhaseAtStart;
kid->axialTilt = csbody->axialTilt;
kid->semiMajorAxis = csbody->semiMajorAxis;
kid->orbit.eccentricity = csbody->eccentricity.ToDouble();
kid->orbit.semiMajorAxis = csbody->semiMajorAxis.ToDouble() * AU;
if(parent->type == SystemBody::TYPE_GRAVPOINT) // generalize Kepler's law to multiple stars
kid->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec_gravpoint(kid->orbit.semiMajorAxis, parent->GetMass(), kid->GetMass(), kid->orbit.eccentricity);
else
kid->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec(kid->orbit.semiMajorAxis, parent->GetMass(), kid->orbit.eccentricity);
kid->orbit.orbitalPhaseAtStart = csbody->orbitalPhaseAtStart.ToDouble();
if (csbody->heightMapFilename.length() > 0) {
kid->heightMapFilename = csbody->heightMapFilename.c_str();
kid->heightMapFractal = csbody->heightMapFractal;
}
if(parent->type == SystemBody::TYPE_GRAVPOINT) // generalize Kepler's law to multiple stars
kid->orbit.SetShapeAroundBarycentre(csbody->semiMajorAxis.ToDouble() * AU, parent->GetMass(), kid->GetMass(), csbody->eccentricity.ToDouble());
else
kid->orbit.SetShapeAroundPrimary(csbody->semiMajorAxis.ToDouble() * AU, parent->GetMass(), csbody->eccentricity.ToDouble());
kid->orbit.SetPhase(csbody->orbitalPhaseAtStart.ToDouble());
if (kid->type == SystemBody::TYPE_STARPORT_SURFACE) {
kid->orbit.rotMatrix = matrix3x3d::RotateY(csbody->longitude) *
matrix3x3d::RotateX(-0.5*M_PI + csbody->latitude);
kid->orbit.SetPlane(matrix3x3d::RotateY(csbody->longitude) * matrix3x3d::RotateX(-0.5*M_PI + csbody->latitude));
} else {
if (kid->orbit.semiMajorAxis < 1.2 * parent->GetRadius()) {
if (kid->orbit.GetSemiMajorAxis() < 1.2 * parent->GetRadius()) {
Error("%s's orbit is too close to its parent", csbody->name.c_str());
}
double offset = csbody->want_rand_offset ? rand.Double(2*M_PI) : (csbody->orbitalOffset.ToDouble()*M_PI);
kid->orbit.rotMatrix = matrix3x3d::RotateY(offset) * matrix3x3d::RotateX(-0.5*M_PI + csbody->latitude);
kid->orbit.SetPlane(matrix3x3d::RotateY(offset) * matrix3x3d::RotateX(-0.5*M_PI + csbody->latitude));
}
if (kid->GetSuperType() == SystemBody::SUPERTYPE_STARPORT) {
(*outHumanInfestedness)++;
@ -1055,19 +1055,16 @@ void StarSystem::MakeBinaryPair(SystemBody *a, SystemBody *b, fixed minDist, Ran
mul *= 2;
} while (a->semiMajorAxis < minDist);
a->orbit.eccentricity = a->eccentricity.ToDouble();
a->orbit.semiMajorAxis = AU * (a->semiMajorAxis * a0).ToDouble();
a->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec_gravpoint(a->orbit.semiMajorAxis, a->GetMass() + b->GetMass(), a->GetMass(), a->orbit.eccentricity);
const double total_mass = a->GetMass() + b->GetMass();
const double e = a->eccentricity.ToDouble();
a->orbit.SetShapeAroundBarycentre(AU * (a->semiMajorAxis * a0).ToDouble(), total_mass, a->GetMass(), e);
b->orbit.SetShapeAroundBarycentre(AU * (a->semiMajorAxis * a1).ToDouble(), total_mass, b->GetMass(), e);
const float rotX = -0.5f*float(M_PI);//(float)(rand.Double()*M_PI/2.0);
const float rotY = static_cast<float>(rand.Double(M_PI));
a->orbit.rotMatrix = matrix3x3d::RotateY(rotY) * matrix3x3d::RotateX(rotX);
b->orbit.rotMatrix = matrix3x3d::RotateY(rotY-M_PI) * matrix3x3d::RotateX(rotX);
b->orbit.eccentricity = a->eccentricity.ToDouble();
b->orbit.semiMajorAxis = AU * (a->semiMajorAxis * a1).ToDouble();
b->orbit.velocityAreaPerSecond = a->orbit.velocityAreaPerSecond;
b->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec_gravpoint(b->orbit.semiMajorAxis, a->GetMass() + b->GetMass(), b->GetMass(), b->orbit.eccentricity);
a->orbit.SetPlane(matrix3x3d::RotateY(rotY) * matrix3x3d::RotateX(rotX));
b->orbit.SetPlane(matrix3x3d::RotateY(rotY-M_PI) * matrix3x3d::RotateX(rotX));
fixed orbMin = a->semiMajorAxis - a->eccentricity*a->semiMajorAxis;
fixed orbMax = 2*a->semiMajorAxis - orbMin;
@ -1658,17 +1655,16 @@ void StarSystem::MakePlanetsAround(SystemBody *primary, Random &rand)
planet->mass = mass;
planet->rotationPeriod = fixed(rand.Int32(1,200), 24);
planet->orbit.eccentricity = ecc.ToDouble();
planet->orbit.semiMajorAxis = semiMajorAxis.ToDouble() * AU;
if(primary->type == SystemBody::TYPE_GRAVPOINT) // generalize Kepler's law to multiple stars
planet->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec_gravpoint(planet->orbit.semiMajorAxis, primary->GetMass(), planet->GetMass(), planet->orbit.eccentricity);
const double e = ecc.ToDouble();
if(primary->type == SystemBody::TYPE_GRAVPOINT)
planet->orbit.SetShapeAroundBarycentre(semiMajorAxis.ToDouble() * AU, primary->GetMass(), planet->GetMass(), e);
else
planet->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec(planet->orbit.semiMajorAxis, primary->GetMass(), planet->orbit.eccentricity);
planet->orbit.SetShapeAroundPrimary(semiMajorAxis.ToDouble() * AU, primary->GetMass(), e);
double r1 = rand.Double(2*M_PI); // function parameter evaluation order is implementation-dependent
double r2 = rand.NDouble(5); // can't put two rands in the same expression
planet->orbit.rotMatrix = matrix3x3d::RotateY(r1) *
matrix3x3d::RotateX(-0.5*M_PI + r2*M_PI/2.0);
planet->orbit.SetPlane(matrix3x3d::RotateY(r1) * matrix3x3d::RotateX(-0.5*M_PI + r2*M_PI/2.0));
planet->orbMin = periapsis;
planet->orbMax = apoapsis;
@ -2122,10 +2118,10 @@ void SystemBody::PopulateAddStations(StarSystem *system)
sp->semiMajorAxis = orbMinS;
sp->eccentricity = fixed(0);
sp->axialTilt = fixed(0);
sp->orbit.eccentricity = 0;
sp->orbit.semiMajorAxis = sp->semiMajorAxis.ToDouble()*AU;
sp->orbit.velocityAreaPerSecond = Orbit::calc_velocity_area_per_sec(sp->orbit.semiMajorAxis, this->mass.ToDouble() * EARTH_MASS, sp->orbit.eccentricity);
sp->orbit.rotMatrix = matrix3x3d::Identity();
sp->orbit.SetShapeAroundPrimary(sp->semiMajorAxis.ToDouble()*AU, this->mass.ToDouble() * EARTH_MASS, 0.0);
sp->orbit.SetPlane(matrix3x3d::Identity());
children.insert(children.begin(), sp);
system->m_spaceStations.push_back(sp);
sp->orbMin = sp->semiMajorAxis;
@ -2139,7 +2135,7 @@ void SystemBody::PopulateAddStations(StarSystem *system)
SystemPath path2 = sp2->path;
*sp2 = *sp;
sp2->path = path2;
sp2->orbit.rotMatrix = matrix3x3d::RotateZ(M_PI);
sp2->orbit.SetPlane(matrix3x3d::RotateZ(M_PI));
sp2->name = gen_unique_station_name(sp, system, namerand);
children.insert(children.begin(), sp2);
system->m_spaceStations.push_back(sp2);