- Use adaptive sampling to calculate cubic splines for drawing orbits and trajectories
- Added capability to show fading orbit 'trails' (disabled by default) - Added more parameters to control the appearance of orbitssensor-dev
Chris Laurel
parent
6b6bdc4316
commit
bfa2f5fb09
|
|
@ -10,6 +10,7 @@
|
|||
|
||||
#define DEBUG_COALESCE 0
|
||||
#define DEBUG_SECONDARY_ILLUMINATION 0
|
||||
#define DEBUG_ORBIT_CACHE 0
|
||||
|
||||
//#define DEBUG_HDR
|
||||
#ifdef DEBUG_HDR
|
||||
|
|
@ -33,6 +34,8 @@ std::ofstream hdrlog;
|
|||
//#define USE_BLOOM_LISTS
|
||||
#endif
|
||||
|
||||
// #define ENABLE_SELF_SHADOW
|
||||
|
||||
#ifndef _WIN32
|
||||
#ifndef TARGET_OS_MAC
|
||||
#include <config.h>
|
||||
|
|
@ -239,6 +242,11 @@ Color Renderer::EclipticColor (0.5f, 0.1f, 0.1f);
|
|||
Color Renderer::SelectionCursorColor (1.0f, 0.0f, 0.0f);
|
||||
|
||||
|
||||
#if ENABLE_SELF_SHADOW
|
||||
static FramebufferObject* shadowFbo = NULL;
|
||||
#endif
|
||||
|
||||
|
||||
// Some useful unit conversions
|
||||
inline float mmToInches(float mm)
|
||||
{
|
||||
|
|
@ -1133,6 +1141,18 @@ bool Renderer::init(GLContext* _context,
|
|||
genSceneTexture();
|
||||
genBlurTextures();
|
||||
#endif
|
||||
|
||||
#if ENABLE_SELF_SHADOW
|
||||
if (GLEW_EXT_framebuffer_object)
|
||||
{
|
||||
shadowFbo = new FramebufferObject(1024, 1024, FramebufferObject::DepthAttachment);
|
||||
if (!shadowFbo->isValid())
|
||||
{
|
||||
clog << "Error creating shadow FBO.\n";
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
commonDataInitialized = true;
|
||||
}
|
||||
|
||||
|
|
@ -1737,21 +1757,40 @@ static void disableSmoothLines()
|
|||
class OrbitSampler : public OrbitSampleProc
|
||||
{
|
||||
public:
|
||||
CurvePlot* m_orbitPath;
|
||||
vector<CurvePlotSample> samples;
|
||||
|
||||
OrbitSampler()
|
||||
{
|
||||
}
|
||||
|
||||
OrbitSampler(CurvePlot* orbitPath) : m_orbitPath(orbitPath) {};
|
||||
void sample(double t, const Vector3d& position, const Vector3d& velocity)
|
||||
{
|
||||
CurvePlotSample samp;
|
||||
samp.t = t;
|
||||
samp.position = position;
|
||||
samp.velocity = velocity;
|
||||
samp.t = t;
|
||||
m_orbitPath->addSample(samp);
|
||||
};
|
||||
samples.push_back(samp);
|
||||
}
|
||||
|
||||
void insertForward(CurvePlot* plot)
|
||||
{
|
||||
for (vector<CurvePlotSample>::const_iterator iter = samples.begin(); iter != samples.end(); ++iter)
|
||||
{
|
||||
plot->addSample(*iter);
|
||||
}
|
||||
}
|
||||
|
||||
void insertBackward(CurvePlot* plot)
|
||||
{
|
||||
for (vector<CurvePlotSample>::const_reverse_iterator iter = samples.rbegin(); iter != samples.rend(); ++iter)
|
||||
{
|
||||
plot->addSample(*iter);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
void renderOrbitColor(const Body *body, bool selected, float opacity)
|
||||
Vector4f renderOrbitColor(const Body *body, bool selected, float opacity)
|
||||
{
|
||||
Color orbitColor;
|
||||
|
||||
|
|
@ -1803,9 +1842,9 @@ void renderOrbitColor(const Body *body, bool selected, float opacity)
|
|||
}
|
||||
|
||||
#ifdef USE_HDR
|
||||
glColor(orbitColor, 1.f - opacity * orbitColor.alpha());
|
||||
return Vector4f(orbitColor.red(), orbitColor.green(), orbitColor.blue(), 1.0f - opacity * orbitColor.alpha());
|
||||
#else
|
||||
glColor(orbitColor, opacity * orbitColor.alpha());
|
||||
return Vector4f(orbitColor.red(), orbitColor.green(), orbitColor.blue(), opacity * orbitColor.alpha());
|
||||
#endif
|
||||
}
|
||||
|
||||
|
|
@ -1872,17 +1911,17 @@ void Renderer::renderOrbit(const OrbitPathListEntry& orbitPath,
|
|||
}
|
||||
else
|
||||
{
|
||||
startTime = t - orbit->getPeriod() / 2.0;
|
||||
startTime = t - orbit->getPeriod();
|
||||
}
|
||||
|
||||
cachedOrbit = new CurvePlot();
|
||||
cachedOrbit->setLastUsed(frameCount);
|
||||
|
||||
OrbitSampler sampler(cachedOrbit);
|
||||
OrbitSampler sampler;
|
||||
orbit->sample(startTime,
|
||||
orbit->getPeriod(),
|
||||
nSamples,
|
||||
startTime + orbit->getPeriod(),
|
||||
sampler);
|
||||
sampler.insertForward(cachedOrbit);
|
||||
|
||||
// If the orbit cache is full, first try and eliminate some old orbits
|
||||
if (orbitCache.size() > OrbitCacheCullThreshold)
|
||||
|
|
@ -1909,42 +1948,73 @@ void Renderer::renderOrbit(const OrbitPathListEntry& orbitPath,
|
|||
if (cachedOrbit->empty())
|
||||
return;
|
||||
|
||||
//*** Orbit rendering parameters
|
||||
|
||||
// The 'window' is the interval of time for which the orbit will be drawn.
|
||||
|
||||
// End of the orbit window relative to the current simulation time. Units
|
||||
// are orbital periods.
|
||||
const double OrbitWindowEnd = 0.5;
|
||||
|
||||
// Number of orbit periods shown. The orbit window is:
|
||||
// [ t + (OrbitWindowEnd - OrbitPeriodsShown) * T, t + OrbitWindowEnd * T ]
|
||||
// where t is the current simulation time and T is the orbital period.
|
||||
const double OrbitPeriodsShown = 1.0;
|
||||
|
||||
// Fraction of the window over which the orbit fades from opaque to transparent.
|
||||
// Fading is disabled when this value is zero.
|
||||
const double LinearFadeFraction = 0.0;
|
||||
|
||||
// Extra size of the internal sample cache.
|
||||
const double WindowSlack = 0.2;
|
||||
|
||||
//***
|
||||
|
||||
// 'Periodic' orbits are generally not strictly periodic because of perturbations
|
||||
// from other bodies. Here we update the trajectory samples to make sure that the
|
||||
// orbit covers a time range centered at the current time and covering a full revolution.
|
||||
if (orbit->isPeriodic())
|
||||
{
|
||||
double startTime = t - orbit->getPeriod() / 2.0;
|
||||
double endTime = t + orbit->getPeriod() / 2.0;
|
||||
double dt = orbit->getPeriod() / detailOptions.orbitPathSamplePoints;
|
||||
double period = orbit->getPeriod();
|
||||
double endTime = t + period * OrbitWindowEnd;
|
||||
double startTime = endTime - period * OrbitPeriodsShown;
|
||||
|
||||
if (startTime < cachedOrbit->startTime())
|
||||
double currentWindowStart = cachedOrbit->startTime();
|
||||
double currentWindowEnd = cachedOrbit->endTime();
|
||||
double newWindowStart = startTime - period * WindowSlack;
|
||||
double newWindowEnd = endTime + period * WindowSlack;
|
||||
|
||||
if (startTime < currentWindowStart)
|
||||
{
|
||||
cachedOrbit->removeSamplesAfter(endTime + dt);
|
||||
// Remove samples at the end of the time window
|
||||
cachedOrbit->removeSamplesAfter(newWindowEnd);
|
||||
|
||||
double orbitStartTime = cachedOrbit->empty() ? endTime : cachedOrbit->startTime() - dt;
|
||||
do {
|
||||
CurvePlotSample sample;
|
||||
sample.t = orbitStartTime;
|
||||
sample.position = orbit->positionAtTime(orbitStartTime);
|
||||
sample.velocity = orbit->velocityAtTime(orbitStartTime);
|
||||
cachedOrbit->addSample(sample);
|
||||
orbitStartTime -= dt;
|
||||
} while (orbitStartTime > startTime);
|
||||
// Trim the first sample (because it will be duplicated when we sample the orbit.)
|
||||
cachedOrbit->removeSamplesBefore(cachedOrbit->startTime() * (1.0 + 1.0e-15));
|
||||
|
||||
// Add the new samples
|
||||
OrbitSampler sampler;
|
||||
orbit->sample(newWindowStart, min(currentWindowStart, newWindowEnd), sampler);
|
||||
sampler.insertBackward(cachedOrbit);
|
||||
#if DEBUG_ORBIT_CACHE
|
||||
clog << "new sample count: " << cachedOrbit->sampleCount() << endl;
|
||||
#endif
|
||||
}
|
||||
else if (endTime > cachedOrbit->endTime())
|
||||
else if (endTime > currentWindowEnd)
|
||||
{
|
||||
cachedOrbit->removeSamplesBefore(startTime - dt);
|
||||
// Remove samples at the beginning of the time window
|
||||
cachedOrbit->removeSamplesBefore(newWindowStart);
|
||||
|
||||
double orbitEndTime = cachedOrbit->empty() ? startTime : cachedOrbit->endTime() + dt;
|
||||
do {
|
||||
CurvePlotSample sample;
|
||||
sample.t = orbitEndTime;
|
||||
sample.position = orbit->positionAtTime(orbitEndTime);
|
||||
sample.velocity = orbit->velocityAtTime(orbitEndTime);
|
||||
cachedOrbit->addSample(sample);
|
||||
orbitEndTime += dt;
|
||||
} while (orbitEndTime < endTime);
|
||||
// Trim the last sample (because it will be duplicated when we sample the orbit.)
|
||||
cachedOrbit->removeSamplesAfter(cachedOrbit->endTime() * (1.0 - 1.0e-15));
|
||||
|
||||
// Add the new samples
|
||||
OrbitSampler sampler;
|
||||
orbit->sample(max(currentWindowEnd, newWindowStart), newWindowEnd, sampler);
|
||||
sampler.insertForward(cachedOrbit);
|
||||
#if DEBUG_ORBIT_CACHE
|
||||
clog << "new sample count: " << cachedOrbit->sampleCount() << endl;
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1970,7 +2040,8 @@ void Renderer::renderOrbit(const OrbitPathListEntry& orbitPath,
|
|||
highlight = highlightObject.body() == body;
|
||||
else
|
||||
highlight = highlightObject.star() == orbitPath.star;
|
||||
renderOrbitColor(body, highlight, orbitPath.opacity);
|
||||
Vector4f orbitColor = renderOrbitColor(body, highlight, orbitPath.opacity);
|
||||
glColor4fv(orbitColor.data());
|
||||
|
||||
#ifdef STIPPLED_LINES
|
||||
glLineStipple(3, 0x5555);
|
||||
|
|
@ -1987,10 +2058,28 @@ void Renderer::renderOrbit(const OrbitPathListEntry& orbitPath,
|
|||
|
||||
if (orbit->isPeriodic())
|
||||
{
|
||||
cachedOrbit->render(modelview,
|
||||
nearZ, farZ, viewFrustumPlaneNormals,
|
||||
subdivisionThreshold,
|
||||
t - orbit->getPeriod() / 2.0, t + orbit->getPeriod() / 2.0);
|
||||
double period = orbit->getPeriod();
|
||||
double windowEnd = t + period * OrbitWindowEnd;
|
||||
double windowStart = windowEnd - period * OrbitPeriodsShown;
|
||||
double windowDuration = windowEnd - windowStart;
|
||||
|
||||
if (LinearFadeFraction == 0.0f)
|
||||
{
|
||||
cachedOrbit->render(modelview,
|
||||
nearZ, farZ, viewFrustumPlaneNormals,
|
||||
subdivisionThreshold,
|
||||
windowStart, windowEnd);
|
||||
}
|
||||
else
|
||||
{
|
||||
cachedOrbit->renderFaded(modelview,
|
||||
nearZ, farZ, viewFrustumPlaneNormals,
|
||||
subdivisionThreshold,
|
||||
windowStart, windowEnd,
|
||||
orbitColor,
|
||||
windowStart,
|
||||
windowEnd - windowDuration * (1.0 - LinearFadeFraction));
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
|
|||
|
|
@ -27,6 +27,152 @@ using namespace std;
|
|||
static const double ORBITAL_VELOCITY_DIFF_DELTA = 1.0 / 1440.0;
|
||||
|
||||
|
||||
static Vector3d cubicInterpolate(const Vector3d& p0, const Vector3d& v0,
|
||||
const Vector3d& p1, const Vector3d& v1,
|
||||
double t)
|
||||
{
|
||||
return p0 + (((2.0 * (p0 - p1) + v1 + v0) * (t * t * t)) +
|
||||
((3.0 * (p1 - p0) - 2.0 * v0 - v1) * (t * t)) +
|
||||
(v0 * t));
|
||||
}
|
||||
|
||||
|
||||
/** Sample the orbit over the time range [ startTime, endTime ] using the
|
||||
* default sampling parameters for the orbit type.
|
||||
*
|
||||
* Subclasses of orbit should override this method as necessary. The default
|
||||
* implementation uses an adaptive sampling scheme with the following defaults:
|
||||
* tolerance: 1 km
|
||||
* start step: T / 1e5
|
||||
* min step: T / 1e7
|
||||
* max step: T / 100
|
||||
*
|
||||
* Where T is either the mean orbital period for periodic orbits or the valid
|
||||
* time span for aperiodic trajectories.
|
||||
*/
|
||||
void Orbit::sample(double startTime, double endTime, OrbitSampleProc& proc) const
|
||||
{
|
||||
double span = 0.0;
|
||||
if (isPeriodic())
|
||||
{
|
||||
span = getPeriod();
|
||||
}
|
||||
else
|
||||
{
|
||||
double startValidInterval = 0.0;
|
||||
double endValidInterval = 0.0;
|
||||
getValidRange(startValidInterval, endValidInterval);
|
||||
if (startValidInterval == endValidInterval)
|
||||
{
|
||||
span = endValidInterval - startValidInterval;
|
||||
}
|
||||
else
|
||||
{
|
||||
span = endTime - startTime;
|
||||
}
|
||||
}
|
||||
|
||||
AdaptiveSamplingParameters samplingParams;
|
||||
samplingParams.tolerance = 1.0; // kilometers
|
||||
samplingParams.maxStep = span / 100.0;
|
||||
samplingParams.minStep = span / 1.0e7;
|
||||
samplingParams.startStep = span / 1.0e5;
|
||||
|
||||
adaptiveSample(startTime, endTime, proc, samplingParams);
|
||||
}
|
||||
|
||||
|
||||
/** Adaptively sample the orbit over the range [ startTime, endTime ].
|
||||
*/
|
||||
void Orbit::adaptiveSample(double startTime, double endTime, OrbitSampleProc& proc, const AdaptiveSamplingParameters& samplingParams) const
|
||||
{
|
||||
double startStepSize = samplingParams.startStep;
|
||||
double maxStepSize = samplingParams.maxStep;
|
||||
double minStepSize = samplingParams.minStep;
|
||||
double tolerance = samplingParams.tolerance;
|
||||
double t = startTime;
|
||||
const double stepFactor = 1.25;
|
||||
|
||||
Vector3d lastP = positionAtTime(t);
|
||||
Vector3d lastV = velocityAtTime(t);
|
||||
proc.sample(t, lastP, lastV);
|
||||
int sampCount = 0;
|
||||
int nTests = 0;
|
||||
|
||||
while (t < endTime)
|
||||
{
|
||||
// Make sure that we don't go past the end of the sample interval
|
||||
maxStepSize = min(maxStepSize, endTime - t);
|
||||
double dt = min(maxStepSize, startStepSize * 2.0);
|
||||
|
||||
Vector3d p1 = positionAtTime(t + dt);
|
||||
Vector3d v1 = velocityAtTime(t + dt);
|
||||
|
||||
double tmid = t + dt / 2.0;
|
||||
Vector3d pTest = positionAtTime(tmid);
|
||||
Vector3d pInterp = cubicInterpolate(lastP, lastV * dt,
|
||||
p1, v1 * dt,
|
||||
0.5);
|
||||
nTests++;
|
||||
|
||||
double positionError = (pInterp - pTest).norm();
|
||||
|
||||
// Error is greater than tolerance; decrease the step until the
|
||||
// error is within the tolerance.
|
||||
if (positionError > tolerance)
|
||||
{
|
||||
while (positionError > tolerance && dt > minStepSize)
|
||||
{
|
||||
dt /= stepFactor;
|
||||
|
||||
p1 = positionAtTime(t + dt);
|
||||
v1 = velocityAtTime(t + dt);
|
||||
|
||||
tmid = t + dt / 2.0;
|
||||
pTest = positionAtTime(tmid);
|
||||
pInterp = cubicInterpolate(lastP, lastV * dt,
|
||||
p1, v1 * dt,
|
||||
0.5);
|
||||
nTests++;
|
||||
|
||||
positionError = (pInterp - pTest).norm();
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Error is less than the tolerance; increase the step size until the
|
||||
// tolerance is just exceeded.
|
||||
while (positionError < tolerance && dt < maxStepSize)
|
||||
{
|
||||
dt *= stepFactor;
|
||||
|
||||
p1 = positionAtTime(t + dt);
|
||||
v1 = velocityAtTime(t + dt);
|
||||
|
||||
tmid = t + dt / 2.0;
|
||||
pTest = positionAtTime(tmid);
|
||||
pInterp = cubicInterpolate(lastP, lastV * dt,
|
||||
p1, v1 * dt,
|
||||
0.5);
|
||||
nTests++;
|
||||
|
||||
positionError = (pInterp - pTest).norm();
|
||||
}
|
||||
}
|
||||
|
||||
t = t + dt;
|
||||
lastP = p1;
|
||||
lastV = v1;
|
||||
|
||||
proc.sample(t, lastP, lastV);
|
||||
sampCount++;
|
||||
}
|
||||
|
||||
// Statistics for debugging
|
||||
// clog << "Orbit samples: " << sampCount << ", nTests: " << nTests << endl;
|
||||
}
|
||||
|
||||
|
||||
EllipticalOrbit::EllipticalOrbit(double _pericenterDistance,
|
||||
double _eccentricity,
|
||||
double _inclination,
|
||||
|
|
@ -286,66 +432,6 @@ double EllipticalOrbit::getBoundingRadius() const
|
|||
}
|
||||
|
||||
|
||||
void EllipticalOrbit::sample(double startTime, double duration, int nSamples,
|
||||
OrbitSampleProc& proc) const
|
||||
{
|
||||
#if 0
|
||||
{
|
||||
// Sample uniformly in time
|
||||
for (int i = 0; i < nSamples; i++)
|
||||
{
|
||||
double tsamp = startTime + duration / nSamples;
|
||||
proc.sample(tsamp, positionAtTime(tsamp), velocityAtTime(tsamp));
|
||||
}
|
||||
}
|
||||
#endif
|
||||
if (eccentricity >= 1.0 || 1)
|
||||
{
|
||||
// Sample uniformly in eccentric anomaly
|
||||
double t = startTime - epoch;
|
||||
double meanMotion = 2.0 * PI / period;
|
||||
double M0 = meanAnomalyAtEpoch + t * meanMotion;
|
||||
double E0 = eccentricAnomaly(M0);
|
||||
double dE = 2 * PI / (double) nSamples;
|
||||
|
||||
for (int i = 0; i < nSamples; i++)
|
||||
{
|
||||
// Compute the time tag for this sample
|
||||
double E = E0 + dE * i;
|
||||
double M = E - eccentricity * sin(E); // Mean anomaly from ecc anomaly
|
||||
double tsamp = startTime + (M - M0) * period / (2 * PI); // Time from mean anomaly
|
||||
proc.sample(tsamp, positionAtE(E), velocityAtE(E));
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Adaptive sampling of the orbit; more samples in regions of high curvature.
|
||||
// nSamples is the number of samples that will be used for a perfectly circular
|
||||
// orbit. Elliptical orbits will have regions of higher curvature thar require
|
||||
// additional sample points.
|
||||
double E = 0.0;
|
||||
double dE = 2 * PI / (double) nSamples;
|
||||
double w = (1 - square(eccentricity));
|
||||
double M0 = E - eccentricity * sin(E);
|
||||
|
||||
while (E < 2 * PI)
|
||||
{
|
||||
// Compute the time tag for this sample
|
||||
double M = E - eccentricity * sin(E); // Mean anomaly from ecc anomaly
|
||||
double tsamp = startTime + (M - M0) * period / (2 * PI); // Time from mean anomaly
|
||||
|
||||
proc.sample(tsamp, positionAtE(E), velocityAtE(E));
|
||||
|
||||
// Compute the curvature
|
||||
double k = w * pow(square(sin(E)) + w * w * square(cos(E)), -1.5);
|
||||
|
||||
// Step amount based on curvature--constrain it so that we don't end up
|
||||
// taking too many samples anywhere. Clamping the curvature to 20 effectively
|
||||
// limits the numbers of samples to 3*nSamples
|
||||
E += dE / max(min(k, 20.0), 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
CachingOrbit::CachingOrbit() :
|
||||
|
|
@ -418,15 +504,6 @@ Vector3d CachingOrbit::computeVelocity(double jd) const
|
|||
}
|
||||
|
||||
|
||||
void CachingOrbit::sample(double start, double t, int nSamples,
|
||||
OrbitSampleProc& proc) const
|
||||
{
|
||||
double dt = t / (double) nSamples;
|
||||
for (int i = 0; i < nSamples; i++)
|
||||
proc.sample(start + dt * i, positionAtTime(start + dt * i), velocityAtTime(start + dt * i));
|
||||
}
|
||||
|
||||
|
||||
static EllipticalOrbit* StateVectorToOrbit(const Vector3d& position,
|
||||
const Vector3d& v,
|
||||
double mass,
|
||||
|
|
@ -556,17 +633,16 @@ double MixedOrbit::getBoundingRadius() const
|
|||
}
|
||||
|
||||
|
||||
void MixedOrbit::sample(double t0, double t1, int nSamples,
|
||||
OrbitSampleProc& proc) const
|
||||
void MixedOrbit::sample(double startTime, double endTime, OrbitSampleProc& proc) const
|
||||
{
|
||||
Orbit* o;
|
||||
if (t0 < begin)
|
||||
if (startTime < begin)
|
||||
o = beforeApprox;
|
||||
else if (t0 < end)
|
||||
else if (startTime < end)
|
||||
o = primary;
|
||||
else
|
||||
o = afterApprox;
|
||||
o->sample(t0, t1, nSamples, proc);
|
||||
o->sample(startTime, endTime, proc);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -612,8 +688,10 @@ FixedOrbit::getBoundingRadius() const
|
|||
|
||||
|
||||
void
|
||||
FixedOrbit::sample(double, double, int, OrbitSampleProc&) const
|
||||
FixedOrbit::sample(double /* startTime */, double /* endTime */, OrbitSampleProc&) const
|
||||
{
|
||||
// Don't add any samples. This will prevent a fixed trajectory from
|
||||
// every being drawn when orbit visualization is enabled.
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -650,7 +728,7 @@ double SynchronousOrbit::getBoundingRadius() const
|
|||
}
|
||||
|
||||
|
||||
void SynchronousOrbit::sample(double, double, int, OrbitSampleProc&) const
|
||||
void SynchronousOrbit::sample(double /* startTime */, double /* endTime */, OrbitSampleProc&) const
|
||||
{
|
||||
// Empty method--we never want to show a synchronous orbit.
|
||||
}
|
||||
|
|
|
|||
|
|
@ -34,14 +34,25 @@ class Orbit
|
|||
|
||||
virtual double getPeriod() const = 0;
|
||||
virtual double getBoundingRadius() const = 0;
|
||||
virtual void sample(double start, double t,
|
||||
int nSamples, OrbitSampleProc& proc) const = 0;
|
||||
|
||||
virtual void sample(double startTime, double endTime, OrbitSampleProc& proc) const;
|
||||
|
||||
virtual bool isPeriodic() const { return true; };
|
||||
|
||||
// Return the time range over which the orbit is valid; if the orbit
|
||||
// is always valid, begin and end should be equal.
|
||||
virtual void getValidRange(double& begin, double& end) const
|
||||
{ begin = 0.0; end = 0.0; };
|
||||
|
||||
struct AdaptiveSamplingParameters
|
||||
{
|
||||
double tolerance;
|
||||
double startStep;
|
||||
double minStep;
|
||||
double maxStep;
|
||||
};
|
||||
|
||||
void adaptiveSample(double startTime, double endTime, OrbitSampleProc& proc, const AdaptiveSamplingParameters& samplingParameters) const;
|
||||
};
|
||||
|
||||
|
||||
|
|
@ -57,7 +68,6 @@ class EllipticalOrbit : public Orbit
|
|||
virtual Eigen::Vector3d velocityAtTime(double) const;
|
||||
double getPeriod() const;
|
||||
double getBoundingRadius() const;
|
||||
virtual void sample(double, double, int, OrbitSampleProc&) const;
|
||||
|
||||
private:
|
||||
double eccentricAnomaly(double) const;
|
||||
|
|
@ -108,8 +118,6 @@ class CachingOrbit : public Orbit
|
|||
Eigen::Vector3d positionAtTime(double jd) const;
|
||||
Eigen::Vector3d velocityAtTime(double jd) const;
|
||||
|
||||
virtual void sample(double, double, int, OrbitSampleProc& proc) const;
|
||||
|
||||
private:
|
||||
mutable Eigen::Vector3d lastPosition;
|
||||
mutable Eigen::Vector3d lastVelocity;
|
||||
|
|
@ -135,8 +143,7 @@ class MixedOrbit : public Orbit
|
|||
virtual Eigen::Vector3d velocityAtTime(double jd) const;
|
||||
virtual double getPeriod() const;
|
||||
virtual double getBoundingRadius() const;
|
||||
virtual void sample(double t0, double t1,
|
||||
int nSamples, OrbitSampleProc& proc) const;
|
||||
virtual void sample(double startTime, double endTime, OrbitSampleProc& proc) const;
|
||||
|
||||
private:
|
||||
Orbit* primary;
|
||||
|
|
@ -165,7 +172,7 @@ class SynchronousOrbit : public Orbit
|
|||
virtual Eigen::Vector3d positionAtTime(double jd) const;
|
||||
virtual double getPeriod() const;
|
||||
virtual double getBoundingRadius() const;
|
||||
virtual void sample(double, double, int, OrbitSampleProc& proc) const;
|
||||
virtual void sample(double, double, OrbitSampleProc& proc) const;
|
||||
|
||||
private:
|
||||
const Body& body;
|
||||
|
|
@ -187,7 +194,7 @@ class FixedOrbit : public Orbit
|
|||
virtual double getPeriod() const;
|
||||
virtual bool isPeriodic() const;
|
||||
virtual double getBoundingRadius() const;
|
||||
virtual void sample(double, double, int, OrbitSampleProc&) const;
|
||||
virtual void sample(double, double, OrbitSampleProc&) const;
|
||||
|
||||
private:
|
||||
Eigen::Vector3d position;
|
||||
|
|
|
|||
|
|
@ -21,6 +21,7 @@
|
|||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <limits>
|
||||
#include <iomanip>
|
||||
|
||||
using namespace Eigen;
|
||||
using namespace std;
|
||||
|
|
@ -76,7 +77,7 @@ public:
|
|||
bool isPeriodic() const;
|
||||
void getValidRange(double& begin, double& end) const;
|
||||
|
||||
virtual void sample(double, double, int, OrbitSampleProc& proc) const;
|
||||
virtual void sample(double startTime, double endTime, OrbitSampleProc& proc) const;
|
||||
|
||||
private:
|
||||
vector<Sample<T> > samples;
|
||||
|
|
@ -392,7 +393,7 @@ template <typename T> Vector3d SampledOrbit<T>::computeVelocity(double jd) const
|
|||
}
|
||||
|
||||
|
||||
template <typename T> void SampledOrbit<T>::sample(double /* start */, double /* t */, int,
|
||||
template <typename T> void SampledOrbit<T>::sample(double /* startTime */, double /* endTime */,
|
||||
OrbitSampleProc& proc) const
|
||||
{
|
||||
for (unsigned int i = 0; i < samples.size(); i++)
|
||||
|
|
@ -446,7 +447,7 @@ public:
|
|||
bool isPeriodic() const;
|
||||
void getValidRange(double& begin, double& end) const;
|
||||
|
||||
virtual void sample(double, double, int, OrbitSampleProc& proc) const;
|
||||
virtual void sample(double startTime, double endTime, OrbitSampleProc& proc) const;
|
||||
|
||||
private:
|
||||
vector<SampleXYZV<T> > samples;
|
||||
|
|
@ -492,7 +493,10 @@ template <typename T> void SampledOrbitXYZV<T>::addSample(double t, const Vector
|
|||
|
||||
template <typename T> double SampledOrbitXYZV<T>::getPeriod() const
|
||||
{
|
||||
return samples[samples.size() - 1].t - samples[0].t;
|
||||
if (samples.empty())
|
||||
return 0.0;
|
||||
else
|
||||
return samples[samples.size() - 1].t - samples[0].t;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -654,7 +658,7 @@ template <typename T> Vector3d SampledOrbitXYZV<T>::computeVelocity(double jd) c
|
|||
}
|
||||
|
||||
|
||||
template <typename T> void SampledOrbitXYZV<T>::sample(double /* start */, double /* t */, int,
|
||||
template <typename T> void SampledOrbitXYZV<T>::sample(double /* startTime */, double /* endTime */,
|
||||
OrbitSampleProc& proc) const
|
||||
{
|
||||
for (typename vector<SampleXYZV<T> >::const_iterator iter = samples.begin();
|
||||
|
|
|
|||
|
|
@ -11084,13 +11084,7 @@ class VSOP87Orbit : public CachingOrbit
|
|||
r += SumSeries(vsR[i], t) * T;
|
||||
T = t * T;
|
||||
}
|
||||
#if 0
|
||||
// l and b are in units of 1e+8 radians
|
||||
l *= 1.0e-8;
|
||||
b *= 1.0e-8;
|
||||
// r is in units of 1e-8 AU
|
||||
r *= (KM_PER_AU / 100000000.0);
|
||||
#endif
|
||||
|
||||
r *= KM_PER_AU;
|
||||
|
||||
// Corrections for internal coordinate system
|
||||
|
|
@ -11101,6 +11095,28 @@ class VSOP87Orbit : public CachingOrbit
|
|||
cos(b) * r,
|
||||
-sin(l) * sin(b) * r);
|
||||
}
|
||||
|
||||
|
||||
/** Custom implementation of sample() for VSOP87 orbits. The default
|
||||
* implementation runs too slowly and produces too many samples.
|
||||
*/
|
||||
void sample(double startTime, double endTime, OrbitSampleProc& proc) const
|
||||
{
|
||||
double span = getPeriod();
|
||||
|
||||
AdaptiveSamplingParameters samplingParams;
|
||||
samplingParams.tolerance = 1.0; // kilometers
|
||||
|
||||
// startStep, minStep, and maxStep are all set to identical values,
|
||||
// so the adaptive sampling function actuall behaves like uniform
|
||||
// sampling.
|
||||
samplingParams.startStep = span / 150.0;
|
||||
samplingParams.minStep = samplingParams.startStep;
|
||||
samplingParams.maxStep = samplingParams.startStep;
|
||||
|
||||
adaptiveSample(startTime, endTime, proc, samplingParams);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue