Index: thirdparty/curveplot/include/curveplot.h =================================================================== --- thirdparty/curveplot/include/curveplot.h (revision 4936) +++ thirdparty/curveplot/include/curveplot.h (working copy) @@ -79,6 +79,16 @@ double subdivisionThreshold, double startTime, double endTime) const; + void renderFaded(const Eigen::Transform3d& modelview, + double nearZ, + double farZ, + const Eigen::Vector3d viewFrustumPlaneNormals[], + double subdivisionThreshold, + double startTime, + double endTime, + const Eigen::Vector4f& color, + double fadeStartTime, + double fadeEndTime) const; unsigned int lastUsed() const { return m_lastUsed; } void setLastUsed(unsigned int lastUsed) { m_lastUsed = lastUsed; } @@ -89,6 +99,8 @@ bool empty() const { return m_samples.empty(); } + unsigned int sampleCount() const { return m_samples.size(); } + private: std::deque m_samples; Index: thirdparty/curveplot/src/curveplot.cpp =================================================================== --- thirdparty/curveplot/src/curveplot.cpp (revision 4936) +++ thirdparty/curveplot/src/curveplot.cpp (working copy) @@ -235,6 +235,25 @@ #endif } + inline void vertex(const Vector4d& v, const Vector4f& color) + { +#if USE_VERTEX_BUFFER + data[currentPosition++] = v.cast(); + ++currentStripLength; + if (currentPosition == capacity) + { + flush(); + + data[0] = v.cast(); + currentPosition = 1; + currentStripLength = 1; + } +#else + glColor4fv(color.data()); + glVertex3dv(v.data()); +#endif + } + inline void begin() { #if !USE_VERTEX_BUFFER @@ -423,6 +442,86 @@ return restartCurve; } + // Return the GL restart status: true if the last segment of the + // curve was culled and we need to start a new primitive sequence + // with glBegin(). + bool renderCubicFaded(bool restartCurve, + const Matrix4d& coeff, + double t0, double t1, + const Vector4f& color, + double fadeStart, double fadeRate, + double curveBoundingRadius, + int depth) const + { + const double dt = (t1 - t0) * InvSubdivisionFactor; + double segmentBoundingRadius = curveBoundingRadius * InvSubdivisionFactor; + +#if DEBUG_ADAPTIVE_SPLINE + { + int c = depth % 10; + glColor4f(SplineColors[c][0], SplineColors[c][1], SplineColors[c][2], 1.0f); + ++SegmentCounts[depth]; + } +#endif + + Vector4d lastP = coeff * Vector4d(1.0, t0, t0 * t0, t0 * t0 * t0); + double lastOpacity = (t0 - fadeStart) * fadeRate; + lastOpacity = max(0.0, min(1.0, lastOpacity)); // clamp + + for (unsigned int i = 1; i <= SubdivisionFactor; i++) + { + double t = t0 + dt * i; + Vector4d p = coeff * Vector4d(1.0, t, t * t, t * t * t); + double opacity = (t - fadeStart) * fadeRate; + opacity = max(0.0, min(1.0, opacity)); // clamp + + double minDistance = max(-m_viewFrustum.nearZ(), abs(p.z()) - segmentBoundingRadius); + + if (segmentBoundingRadius >= m_subdivisionThreshold * minDistance) + { + if (m_viewFrustum.cullSphere(p, segmentBoundingRadius)) + { + if (!restartCurve) + { + m_vbuf.end(); + restartCurve = true; + } + } + else + { + restartCurve = renderCubicFaded(restartCurve, + coeff, t - dt, t, + color, + fadeStart, fadeRate, + segmentBoundingRadius, + depth + 1); + } + } + else + { +#if DEBUG_ADAPTIVE_SPLINE + { + int c = depth % 10; + glColor4f(SplineColors[c][0], SplineColors[c][1], SplineColors[c][2], i % 2 ? 0.25f : 1.0f); + } +#endif + + if (restartCurve) + { + m_vbuf.begin(); + m_vbuf.vertex(lastP, Vector4f(color.x(), color.y(), color.z(), color.w() * float(lastOpacity))); + restartCurve = false; + } + + m_vbuf.vertex(p, Vector4f(color.x(), color.y(), color.z(), color.w() * float(opacity))); + } + lastP = p; + lastOpacity = opacity; + } + + return restartCurve; + } + private: HighPrec_VertexBuffer& m_vbuf; HighPrec_Frustum& m_viewFrustum; @@ -554,7 +653,6 @@ vbuf.createVertexBuffer(); vbuf.setup(); - clog << "size: " << m_samples.size() << endl; for (unsigned int i = 1; i < m_samples.size(); i++) { // Transform the points into camera space. @@ -805,3 +903,170 @@ vbuf.flush(); vbuf.finish(); } + + +/** Draw a curve with a fade effect. The curve is at full opacity at fadeStartTime + * and completely transparent at fadeEndTime. fadeStartTime may be greater than + * fadeEndTime--this just means that the fade direction will be reversed. + */ +void +CurvePlot::renderFaded(const Eigen::Transform3d& modelview, + double nearZ, + double farZ, + const Eigen::Vector3d viewFrustumPlaneNormals[], + double subdivisionThreshold, + double startTime, + double endTime, + const Vector4f& color, + double fadeStartTime, + double fadeEndTime) const +{ + // Flag to indicate whether we need to issue a glBegin() + bool restartCurve = true; + + if (m_samples.empty() || endTime <= m_samples.front().t || startTime >= m_samples.back().t) + return; + + // Linear search for the first sample + unsigned int startSample = 0; + while (startSample < m_samples.size() - 1 && startTime < m_samples[startSample].t) + startSample++; + + // Start at the first sample with time <= startTime + if (startSample > 0) + startSample--; + + double fadeDuration = fadeEndTime - fadeStartTime; + double fadeRate = 1.0 / fadeDuration; + + const Vector3d& p0_ = m_samples[startSample].position; + const Vector3d& v0_ = m_samples[startSample].velocity; + Vector4d p0 = modelview * Vector4d(p0_.x(), p0_.y(), p0_.z(), 1.0); + Vector4d v0 = modelview * Vector4d(v0_.x(), v0_.y(), v0_.z(), 0.0); + double opacity0 = (m_samples[startSample].t - fadeStartTime) * fadeRate; + opacity0 = max(0.0, min(1.0, opacity0)); + + HighPrec_Frustum viewFrustum(nearZ, farZ, viewFrustumPlaneNormals); + HighPrec_RenderContext rc(vbuf, viewFrustum, subdivisionThreshold); + + vbuf.createVertexBuffer(); + vbuf.setup(); + + bool firstSegment = true; + bool lastSegment = false; + + for (unsigned int i = startSample + 1; i < m_samples.size() && !lastSegment; i++) + { + // Transform the points into camera space. + const Vector3d& p1_ = m_samples[i].position; + const Vector3d& v1_ = m_samples[i].velocity; + Vector4d p1 = modelview * Vector4d(p1_.x(), p1_.y(), p1_.z(), 1.0); + Vector4d v1 = modelview * Vector4d(v1_.x(), v1_.y(), v1_.z(), 0.0); + double opacity1 = (m_samples[i].t - fadeStartTime) * fadeRate; + opacity1 = max(0.0, min(1.0, opacity1)); + + if (endTime <= m_samples[i].t) + { + lastSegment = true; + } + + // O(t) is an approximating function for this segment of + // the orbit, with 0 <= t <= 1 + // C is the viewer position + // d(t) = |O(t) - C|, the distance from viewer to the + // orbit segment. + + double curveBoundingRadius = m_samples[i].boundingRadius; + + // Estimate the minimum possible distance from the + // curve to the z=0 plane. If the curve is far enough + // away to be approximated as a straight line, we'll just + // render it. Otherwise, it should be a performance win + // to do a sphere-frustum cull test before subdividing and + // rendering segment. + double minDistance = abs(p0.z()) - curveBoundingRadius; + + // Render close segments as splines with adaptive subdivision. The + // subdivisions eliminates kinks between line segments and also + // prevents clipping precision problems that occur when a + // very long line is rendered with a relatively small view + // volume. + if (curveBoundingRadius >= subdivisionThreshold * minDistance || lastSegment || firstSegment) + { + // Skip rendering this section if it lies outside the view + // frustum. + if (viewFrustum.cullSphere(p0, curveBoundingRadius)) + { + if (!restartCurve) + { + vbuf.end(); + restartCurve = true; + } + } + else + { + double dt = m_samples[i].t - m_samples[i - 1].t; + double t0 = 0.0; + double t1 = 1.0; + + if (firstSegment) + { + t0 = (startTime - m_samples[i - 1].t) / dt; + t0 = std::max(0.0, std::min(1.0, t0)); + firstSegment = false; + } + + if (lastSegment) + { + t1 = (endTime - m_samples[i - 1].t) / dt; + } + + Matrix4d coeff = cubicHermiteCoefficients(p0, p1, v0 * dt, v1 * dt); + restartCurve = rc.renderCubicFaded(restartCurve, coeff, + t0, t1, + color, + (fadeStartTime - m_samples[i - 1].t) / dt, fadeRate * dt, + curveBoundingRadius, 1); + } + } + else + { + // Apparent size of curve is small enough that we can approximate + // it as a line. + + // Simple cull test--just check the far plane. This is required because + // apparent clipping precision limitations can cause a GPU to draw lines + // that lie completely beyond the far plane. + if (p0.z() + curveBoundingRadius < farZ) + { + if (!restartCurve) + { + vbuf.end(); + restartCurve = true; + } + } + else + { + if (restartCurve) + { + vbuf.begin(); + vbuf.vertex(p0, Vector4f(color.x(), color.y(), color.z(), color.w() * float(opacity0))); + restartCurve = false; + } + vbuf.vertex(p1, Vector4f(color.x(), color.y(), color.z(), color.w() * float(opacity1))); + } + } + + p0 = p1; + v0 = v1; + opacity0 = opacity1; + } + + if (!restartCurve) + { + vbuf.end(); + } + + vbuf.flush(); + vbuf.finish(); +} Index: src/celephem/samporbit.cpp =================================================================== --- src/celephem/samporbit.cpp (revision 4936) +++ src/celephem/samporbit.cpp (working copy) @@ -21,6 +21,7 @@ #include #include #include +#include using namespace Eigen; using namespace std; @@ -76,7 +77,7 @@ 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 > samples; @@ -392,7 +393,7 @@ } -template void SampledOrbit::sample(double /* start */, double /* t */, int, +template void SampledOrbit::sample(double /* startTime */, double /* endTime */, OrbitSampleProc& proc) const { for (unsigned int i = 0; i < samples.size(); i++) @@ -446,7 +447,7 @@ 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 > samples; @@ -492,7 +493,10 @@ template double SampledOrbitXYZV::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 void SampledOrbitXYZV::sample(double /* start */, double /* t */, int, +template void SampledOrbitXYZV::sample(double /* startTime */, double /* endTime */, OrbitSampleProc& proc) const { for (typename vector >::const_iterator iter = samples.begin(); Index: src/celephem/vsop87.cpp =================================================================== --- src/celephem/vsop87.cpp (revision 4936) +++ src/celephem/vsop87.cpp (working copy) @@ -11084,13 +11084,7 @@ 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 @@ 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); + } + }; Index: src/celephem/orbit.h =================================================================== --- src/celephem/orbit.h (revision 4936) +++ src/celephem/orbit.h (working copy) @@ -34,14 +34,25 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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; Index: src/celephem/orbit.cpp =================================================================== --- src/celephem/orbit.cpp (revision 4936) +++ src/celephem/orbit.cpp (working copy) @@ -27,6 +27,152 @@ 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,68 +432,8 @@ } -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() : lastTime(-1.0e30), positionCacheValid(false), @@ -418,15 +504,6 @@ } -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 @@ } -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 @@ 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 @@ } -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. } Index: src/celengine/render.cpp =================================================================== --- src/celengine/render.cpp (revision 4936) +++ src/celengine/render.cpp (working copy) @@ -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 @@ //#define USE_BLOOM_LISTS #endif +// #define ENABLE_SELF_SHADOW + #ifndef _WIN32 #ifndef TARGET_OS_MAC #include @@ -239,6 +242,11 @@ 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 @@ 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 @@ class OrbitSampler : public OrbitSampleProc { public: - CurvePlot* m_orbitPath; + vector samples; - OrbitSampler(CurvePlot* orbitPath) : m_orbitPath(orbitPath) {}; + OrbitSampler() + { + } + 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::const_iterator iter = samples.begin(); iter != samples.end(); ++iter) + { + plot->addSample(*iter); + } + } + + void insertBackward(CurvePlot* plot) + { + for (vector::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 @@ } #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 @@ } 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) @@ -1914,37 +1953,46 @@ // 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; + const double windowSlack = 0.2; + double startTime = t - orbit->getPeriod(); + double endTime = t; - if (startTime < cachedOrbit->startTime()) + double currentWindowStart = cachedOrbit->startTime(); + double currentWindowEnd = cachedOrbit->endTime(); + double newWindowStart = startTime - orbit->getPeriod() * windowSlack; + double newWindowEnd = endTime + orbit->getPeriod() * 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 +2018,8 @@ 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 +2036,11 @@ if (orbit->isPeriodic()) { - cachedOrbit->render(modelview, - nearZ, farZ, viewFrustumPlaneNormals, - subdivisionThreshold, - t - orbit->getPeriod() / 2.0, t + orbit->getPeriod() / 2.0); + cachedOrbit->renderFaded(modelview, + nearZ, farZ, viewFrustumPlaneNormals, + subdivisionThreshold, + t - orbit->getPeriod(), t, + orbitColor, t - orbit->getPeriod(), t - orbit->getPeriod() * 0.5); } else {