1 // samporbit.cpp
2 //
3 // Copyright (C) 2002-2008, Celestia Development Team
4 // Original version by Chris Laurel <claurel@gmail.com>
5 //
6 // Trajectories based on unevenly spaced cartesian positions.
7 //
8 // This program is free software; you can redistribute it and/or
9 // modify it under the terms of the GNU General Public License
10 // as published by the Free Software Foundation; either version 2
11 // of the License, or (at your option) any later version.
12
13 #include <cmath>
14 #include <string>
15 #include <algorithm>
16 #include <vector>
17 #include <iostream>
18 #include <fstream>
19 #include <limits>
20 #include <celmath/mathlib.h>
21 #include <celengine/astro.h>
22 #include <celengine/orbit.h>
23 #include <celengine/samporbit.h>
24
25 using namespace std;
26
27 // Trajectories are sampled adaptively for rendering. MaxSampleInterval
28 // is the maximum time (in days) between samples. The threshold angle
29 // is the maximum angle allowed between path segments.
30 static const double MinSampleInterval = 1.0 / 1440.0; // one minute
31 static const double MaxSampleInterval = 50.0;
32 static const double SampleThresholdAngle = 2.0;
33
34 // Position-only sample
35 template <typename T> struct Sample
36 {
37 double t;
38 T x, y, z;
39 };
40
41 // Position + velocity sample
42 template <typename T> struct SampleXYZV
43 {
44 double t;
45 Vector3<T> position;
46 Vector3<T> velocity;
47 };
48
49
operator <(const Sample<T> & a,const Sample<T> & b)50 template <typename T> bool operator<(const Sample<T>& a, const Sample<T>& b)
51 {
52 return a.t < b.t;
53 }
54
operator <(const SampleXYZV<T> & a,const SampleXYZV<T> & b)55 template <typename T> bool operator<(const SampleXYZV<T>& a, const SampleXYZV<T>& b)
56 {
57 return a.t < b.t;
58 }
59
60
61 template <typename T> class SampledOrbit : public CachingOrbit
62 {
63 public:
64 SampledOrbit(TrajectoryInterpolation);
65 virtual ~SampledOrbit();
66
67 void addSample(double t, double x, double y, double z);
68 void setPeriod();
69
70 double getPeriod() const;
71 double getBoundingRadius() const;
72 Point3d computePosition(double jd) const;
73 Vec3d computeVelocity(double jd) const;
74
75 bool isPeriodic() const;
76 void getValidRange(double& begin, double& end) const;
77
78 virtual void sample(double, double, int, OrbitSampleProc& proc) const;
79
80 private:
81 vector<Sample<T> > samples;
82 double boundingRadius;
83 double period;
84 mutable int lastSample;
85
86 TrajectoryInterpolation interpolation;
87 };
88
89
SampledOrbit(TrajectoryInterpolation _interpolation)90 template <typename T> SampledOrbit<T>::SampledOrbit(TrajectoryInterpolation _interpolation) :
91 boundingRadius(0.0),
92 period(1.0),
93 lastSample(0),
94 interpolation(_interpolation)
95 {
96 }
97
98
~SampledOrbit()99 template <typename T> SampledOrbit<T>::~SampledOrbit()
100 {
101 }
102
103
addSample(double t,double x,double y,double z)104 template <typename T> void SampledOrbit<T>::addSample(double t, double x, double y, double z)
105 {
106 double r = sqrt(x * x + y * y + z * z);
107 if (r > boundingRadius)
108 boundingRadius = r;
109
110 Sample<T> samp;
111 samp.x = (T) x;
112 samp.y = (T) y;
113 samp.z = (T) z;
114 samp.t = t;
115 samples.insert(samples.end(), samp);
116 }
117
getPeriod() const118 template <typename T> double SampledOrbit<T>::getPeriod() const
119 {
120 return samples[samples.size() - 1].t - samples[0].t;
121 }
122
123
isPeriodic() const124 template <typename T> bool SampledOrbit<T>::isPeriodic() const
125 {
126 return false;
127 }
128
129
getValidRange(double & begin,double & end) const130 template <typename T> void SampledOrbit<T>::getValidRange(double& begin, double& end) const
131 {
132 begin = samples[0].t;
133 end = samples[samples.size() - 1].t;
134 }
135
136
getBoundingRadius() const137 template <typename T> double SampledOrbit<T>::getBoundingRadius() const
138 {
139 return boundingRadius;
140 }
141
142
cubicInterpolate(const Vec3d & p0,const Vec3d & v0,const Vec3d & p1,const Vec3d & v1,double t)143 static Vec3d cubicInterpolate(const Vec3d& p0, const Vec3d& v0,
144 const Vec3d& p1, const Vec3d& v1,
145 double t)
146 {
147 return p0 + (((2.0 * (p0 - p1) + v1 + v0) * (t * t * t)) +
148 ((3.0 * (p1 - p0) - 2.0 * v0 - v1) * (t * t)) +
149 (v0 * t));
150 }
151
152
cubicInterpolateVelocity(const Vec3d & p0,const Vec3d & v0,const Vec3d & p1,const Vec3d & v1,double t)153 static Vec3d cubicInterpolateVelocity(const Vec3d& p0, const Vec3d& v0,
154 const Vec3d& p1, const Vec3d& v1,
155 double t)
156 {
157 return ((2.0 * (p0 - p1) + v1 + v0) * (3.0 * t * t)) +
158 ((3.0 * (p1 - p0) - 2.0 * v0 - v1) * (2.0 * t)) +
159 v0;
160 }
161
162
computePosition(double jd) const163 template <typename T> Point3d SampledOrbit<T>::computePosition(double jd) const
164 {
165 Vec3d pos;
166 if (samples.size() == 0)
167 {
168 pos = Vec3d(0.0, 0.0, 0.0);
169 }
170 else if (samples.size() == 1)
171 {
172 pos = Vec3d(samples[0].x, samples[1].y, samples[2].z);
173 }
174 else
175 {
176 Sample<T> samp;
177 samp.t = jd;
178 int n = lastSample;
179
180 if (n < 1 || n >= (int) samples.size() || jd < samples[n - 1].t || jd > samples[n].t)
181 {
182 typename vector<Sample<T> >::const_iterator iter = lower_bound(samples.begin(),
183 samples.end(),
184 samp);
185 if (iter == samples.end())
186 n = samples.size();
187 else
188 n = iter - samples.begin();
189
190 lastSample = n;
191 }
192
193 if (n == 0)
194 {
195 pos = Vec3d(samples[n].x, samples[n].y, samples[n].z);
196 }
197 else if (n < (int) samples.size())
198 {
199 if (interpolation == TrajectoryInterpolationLinear)
200 {
201 Sample<T> s0 = samples[n - 1];
202 Sample<T> s1 = samples[n];
203
204 double t = (jd - s0.t) / (s1.t - s0.t);
205 pos = Vec3d(Mathd::lerp(t, (double) s0.x, (double) s1.x),
206 Mathd::lerp(t, (double) s0.y, (double) s1.y),
207 Mathd::lerp(t, (double) s0.z, (double) s1.z));
208 }
209 else if (interpolation == TrajectoryInterpolationCubic)
210 {
211 Sample<T> s0, s1, s2, s3;
212 if (n > 1)
213 s0 = samples[n - 2];
214 else
215 s0 = samples[n - 1];
216 s1 = samples[n - 1];
217 s2 = samples[n];
218 if (n < (int) samples.size() - 1)
219 s3 = samples[n + 1];
220 else
221 s3 = samples[n];
222
223 double h = s2.t - s1.t;
224 double ih = 1.0 / h;
225 double t = (jd - s1.t) * ih;
226 Vec3d p0(s1.x, s1.y, s1.z);
227 Vec3d p1(s2.x, s2.y, s2.z);
228
229 Vec3d v10((double) s1.x - (double) s0.x,
230 (double) s1.y - (double) s0.y,
231 (double) s1.z - (double) s0.z);
232 Vec3d v21((double) s2.x - (double) s1.x,
233 (double) s2.y - (double) s1.y,
234 (double) s2.z - (double) s1.z);
235 Vec3d v32((double) s3.x - (double) s2.x,
236 (double) s3.y - (double) s2.y,
237 (double) s3.z - (double) s2.z);
238
239 // Estimate velocities by averaging the differences at adjacent spans
240 // (except at the end spans, where we just use a single velocity.)
241 Vec3d v0;
242 if (n > 1)
243 {
244 v0 = v10 * (0.5 / (s1.t - s0.t)) + v21 * (0.5 * ih);
245 v0 *= h;
246 }
247 else
248 {
249 v0 = v21;
250 }
251
252 Vec3d v1;
253 if (n < (int) samples.size() - 1)
254 {
255 v1 = v21 * (0.5 * ih) + v32 * (0.5 / (s3.t - s2.t));
256 v1 *= h;
257 }
258 else
259 {
260 v1 = v21;
261 }
262
263 pos = cubicInterpolate(p0, v0, p1, v1, t);
264 }
265 else
266 {
267 // Unknown interpolation type
268 pos = Vec3d(0.0, 0.0, 0.0);
269 }
270 }
271 else
272 {
273 pos = Vec3d(samples[n - 1].x, samples[n - 1].y, samples[n - 1].z);
274 }
275 }
276
277 // Add correction for Celestia's coordinate system
278 return Point3d(pos.x, pos.z, -pos.y);
279 }
280
281
computeVelocity(double jd) const282 template <typename T> Vec3d SampledOrbit<T>::computeVelocity(double jd) const
283 {
284 Vec3d vel;
285 if (samples.size() < 2)
286 {
287 vel = Vec3d(0.0, 0.0, 0.0);
288 }
289 else
290 {
291 Sample<T> samp;
292 samp.t = jd;
293 int n = lastSample;
294
295 if (n < 1 || n >= (int) samples.size() || jd < samples[n - 1].t || jd > samples[n].t)
296 {
297 typename vector<Sample<T> >::const_iterator iter = lower_bound(samples.begin(),
298 samples.end(),
299 samp);
300 if (iter == samples.end())
301 n = samples.size();
302 else
303 n = iter - samples.begin();
304 lastSample = n;
305 }
306
307 if (n == 0)
308 {
309 vel = Vec3d(0.0, 0.0, 0.0);
310 }
311 else if (n < (int) samples.size())
312 {
313 if (interpolation == TrajectoryInterpolationLinear)
314 {
315 Sample<T> s0 = samples[n - 1];
316 Sample<T> s1 = samples[n];
317
318 double dt = (s1.t - s0.t);
319 return (Vec3d(s1.x, s1.y, s1.z) - Vec3d(s0.x, s0.y, s0.z)) * (1.0 / dt);
320 }
321 else if (interpolation == TrajectoryInterpolationCubic)
322 {
323 Sample<T> s0, s1, s2, s3;
324 if (n > 1)
325 s0 = samples[n - 2];
326 else
327 s0 = samples[n - 1];
328 s1 = samples[n - 1];
329 s2 = samples[n];
330 if (n < (int) samples.size() - 1)
331 s3 = samples[n + 1];
332 else
333 s3 = samples[n];
334
335 double h = s2.t - s1.t;
336 double ih = 1.0 / h;
337 double t = (jd - s1.t) * ih;
338 Vec3d p0(s1.x, s1.y, s1.z);
339 Vec3d p1(s2.x, s2.y, s2.z);
340
341 Vec3d v10((double) s1.x - (double) s0.x,
342 (double) s1.y - (double) s0.y,
343 (double) s1.z - (double) s0.z);
344 Vec3d v21((double) s2.x - (double) s1.x,
345 (double) s2.y - (double) s1.y,
346 (double) s2.z - (double) s1.z);
347 Vec3d v32((double) s3.x - (double) s2.x,
348 (double) s3.y - (double) s2.y,
349 (double) s3.z - (double) s2.z);
350
351 // Estimate velocities by averaging the differences at adjacent spans
352 // (except at the end spans, where we just use a single velocity.)
353 Vec3d v0;
354 if (n > 1)
355 {
356 v0 = v10 * (0.5 / (s1.t - s0.t)) + v21 * (0.5 * ih);
357 v0 *= h;
358 }
359 else
360 {
361 v0 = v21;
362 }
363
364 Vec3d v1;
365 if (n < (int) samples.size() - 1)
366 {
367 v1 = v21 * (0.5 * ih) + v32 * (0.5 / (s3.t - s2.t));
368 v1 *= h;
369 }
370 else
371 {
372 v1 = v21;
373 }
374
375 vel = cubicInterpolateVelocity(p0, v0, p1, v1, t);
376 vel *= 1.0 / h;
377 }
378 else
379 {
380 // Unknown interpolation type
381 vel = Vec3d(0.0, 0.0, 0.0);
382 }
383 }
384 else
385 {
386 vel = Vec3d(0.0, 0.0, 0.0);
387 }
388 }
389
390 return Vec3d(vel.x, vel.z, -vel.y);
391 }
392
393
sample(double start,double t,int,OrbitSampleProc & proc) const394 template <typename T> void SampledOrbit<T>::sample(double start, double t, int,
395 OrbitSampleProc& proc) const
396 {
397 double cosThresholdAngle = cos(degToRad(SampleThresholdAngle));
398 double dt = MinSampleInterval;
399 double end = start + t;
400 double current = start;
401
402 proc.sample(current, positionAtTime(current));
403
404 while (current < end)
405 {
406 double dt2 = dt;
407
408 Point3d goodpt;
409 double gooddt = dt;
410 Point3d pos0 = positionAtTime(current);
411 goodpt = positionAtTime(current + dt2);
412 while (1)
413 {
414 Point3d pos1 = positionAtTime(current + dt2);
415 Point3d pos2 = positionAtTime(current + dt2 * 2.0);
416 Vec3d vec1 = pos1 - pos0;
417 Vec3d vec2 = pos2 - pos0;
418
419 vec1.normalize();
420 vec2.normalize();
421 double dot = vec1 * vec2;
422
423 if (dot > cosThresholdAngle && dt2 < MaxSampleInterval)
424 {
425 gooddt = dt2;
426 goodpt = pos1;
427 dt2 *= 2.0;
428 }
429 else
430 {
431 proc.sample(current + gooddt, goodpt);
432 break;
433 }
434 }
435 current += gooddt;
436 }
437 }
438
439
440 // Sampled orbit with positions and velocities
441 template <typename T> class SampledOrbitXYZV : public CachingOrbit
442 {
443 public:
444 SampledOrbitXYZV(TrajectoryInterpolation);
445 virtual ~SampledOrbitXYZV();
446
447 void addSample(double t, const Vec3d& position, const Vec3d& velocity);
448 void setPeriod();
449
450 double getPeriod() const;
451 double getBoundingRadius() const;
452 Point3d computePosition(double jd) const;
453 Vec3d computeVelocity(double jd) const;
454
455 bool isPeriodic() const;
456 void getValidRange(double& begin, double& end) const;
457
458 virtual void sample(double, double, int, OrbitSampleProc& proc) const;
459
460 private:
461 vector<SampleXYZV<T> > samples;
462 double boundingRadius;
463 double period;
464 mutable int lastSample;
465
466 TrajectoryInterpolation interpolation;
467 };
468
469
SampledOrbitXYZV(TrajectoryInterpolation _interpolation)470 template <typename T> SampledOrbitXYZV<T>::SampledOrbitXYZV(TrajectoryInterpolation _interpolation) :
471 boundingRadius(0.0),
472 period(1.0),
473 lastSample(0),
474 interpolation(_interpolation)
475 {
476 }
477
478
~SampledOrbitXYZV()479 template <typename T> SampledOrbitXYZV<T>::~SampledOrbitXYZV()
480 {
481 }
482
483
484 // Add a new sample to the trajectory:
485 // Position in km
486 // Velocity in km/Julian day
addSample(double t,const Vec3d & position,const Vec3d & velocity)487 template <typename T> void SampledOrbitXYZV<T>::addSample(double t, const Vec3d& position, const Vec3d& velocity)
488 {
489 double r = position.length();
490 if (r > boundingRadius)
491 boundingRadius = r;
492
493 SampleXYZV<T> samp;
494 samp.t = t;
495 samp.position = Vector3<T>((T) position.x, (T) position.y, (T) position.z);
496 samp.velocity = Vector3<T>((T) velocity.x, (T) velocity.y, (T) velocity.z);
497 samples.push_back(samp);
498 }
499
getPeriod() const500 template <typename T> double SampledOrbitXYZV<T>::getPeriod() const
501 {
502 return samples[samples.size() - 1].t - samples[0].t;
503 }
504
505
isPeriodic() const506 template <typename T> bool SampledOrbitXYZV<T>::isPeriodic() const
507 {
508 return false;
509 }
510
511
getValidRange(double & begin,double & end) const512 template <typename T> void SampledOrbitXYZV<T>::getValidRange(double& begin, double& end) const
513 {
514 begin = samples[0].t;
515 end = samples[samples.size() - 1].t;
516 }
517
518
getBoundingRadius() const519 template <typename T> double SampledOrbitXYZV<T>::getBoundingRadius() const
520 {
521 return boundingRadius;
522 }
523
524
computePosition(double jd) const525 template <typename T> Point3d SampledOrbitXYZV<T>::computePosition(double jd) const
526 {
527 Vec3d pos;
528 if (samples.size() == 0)
529 {
530 pos = Vec3d(0.0, 0.0, 0.0);
531 }
532 else if (samples.size() == 1)
533 {
534 pos = Vec3d(samples[0].position.x, samples[1].position.y, samples[2].position.z);
535 }
536 else
537 {
538 SampleXYZV<T> samp;
539 samp.t = jd;
540 int n = lastSample;
541
542 if (n < 1 || n >= (int) samples.size() || jd < samples[n - 1].t || jd > samples[n].t)
543 {
544 typename vector<SampleXYZV<T> >::const_iterator iter = lower_bound(samples.begin(),
545 samples.end(),
546 samp);
547 if (iter == samples.end())
548 n = samples.size();
549 else
550 n = iter - samples.begin();
551
552 lastSample = n;
553 }
554
555 if (n == 0)
556 {
557 pos = Vec3d(samples[n].position.x, samples[n].position.y, samples[n].position.z);
558 }
559 else if (n < (int) samples.size())
560 {
561 SampleXYZV<T> s0 = samples[n - 1];
562 SampleXYZV<T> s1 = samples[n];
563
564 if (interpolation == TrajectoryInterpolationLinear)
565 {
566 double t = (jd - s0.t) / (s1.t - s0.t);
567
568 pos = Vec3d(Mathd::lerp(t, (double) s0.position.x, (double) s1.position.x),
569 Mathd::lerp(t, (double) s0.position.y, (double) s1.position.y),
570 Mathd::lerp(t, (double) s0.position.z, (double) s1.position.z));
571 }
572 else if (interpolation == TrajectoryInterpolationCubic)
573 {
574 double h = s1.t - s0.t;
575 double ih = 1.0 / h;
576 double t = (jd - s0.t) * ih;
577
578 Vec3d p0(s0.position.x, s0.position.y, s0.position.z);
579 Vec3d v0(s0.velocity.x, s0.velocity.y, s0.velocity.z);
580 Vec3d p1(s1.position.x, s1.position.y, s1.position.z);
581 Vec3d v1(s1.velocity.x, s1.velocity.y, s1.velocity.z);
582
583 pos = cubicInterpolate(p0, v0 * h, p1, v1 * h, t);
584 }
585 else
586 {
587 // Unknown interpolation type
588 pos = Vec3d(0.0, 0.0, 0.0);
589 }
590 }
591 else
592 {
593 pos = Vec3d(samples[n - 1].position.x, samples[n - 1].position.y, samples[n - 1].position.z);
594 }
595 }
596
597 // Add correction for Celestia's coordinate system
598 return Point3d(pos.x, pos.z, -pos.y);
599 }
600
601
602 // Velocity is computed as the derivative of the interpolating function
603 // for position.
computeVelocity(double jd) const604 template <typename T> Vec3d SampledOrbitXYZV<T>::computeVelocity(double jd) const
605 {
606 Vec3d vel(0.0, 0.0, 0.0);
607
608 if (samples.size() >= 2)
609 {
610 SampleXYZV<T> samp;
611 samp.t = jd;
612 int n = lastSample;
613
614 if (n < 1 || n >= (int) samples.size() || jd < samples[n - 1].t || jd > samples[n].t)
615 {
616 typename vector<SampleXYZV<T> >::const_iterator iter = lower_bound(samples.begin(),
617 samples.end(),
618 samp);
619 if (iter == samples.end())
620 n = samples.size();
621 else
622 n = iter - samples.begin();
623
624 lastSample = n;
625 }
626
627 if (n > 0 && n < (int) samples.size())
628 {
629 SampleXYZV<T> s0 = samples[n - 1];
630 SampleXYZV<T> s1 = samples[n];
631
632 if (interpolation == TrajectoryInterpolationLinear)
633 {
634 double h = s1.t - s0.t;
635 vel = Vec3d(s1.position.x - s0.position.x,
636 s1.position.y - s0.position.y,
637 s1.position.z - s0.position.z) * (1.0 / h) * astro::daysToSecs(1.0);
638 }
639 else if (interpolation == TrajectoryInterpolationCubic)
640 {
641 double h = s1.t - s0.t;
642 double ih = 1.0 / h;
643 double t = (jd - s0.t) * ih;
644
645 Vec3d p0(s0.position.x, s0.position.y, s0.position.z);
646 Vec3d p1(s1.position.x, s1.position.y, s1.position.z);
647 Vec3d v0(s0.velocity.x, s0.velocity.y, s0.velocity.z);
648 Vec3d v1(s1.velocity.x, s1.velocity.y, s1.velocity.z);
649
650 vel = cubicInterpolateVelocity(p0, v0 * h, p1, v1 * h, t) * ih;
651 }
652 else
653 {
654 // Unknown interpolation type
655 vel = Vec3d(0.0, 0.0, 0.0);
656 }
657 }
658 }
659
660 // Add correction for Celestia's coordinate system
661 return Vec3d(vel.x, vel.z, -vel.y);
662 }
663
664
sample(double start,double t,int,OrbitSampleProc & proc) const665 template <typename T> void SampledOrbitXYZV<T>::sample(double start, double t, int,
666 OrbitSampleProc& proc) const
667 {
668 double cosThresholdAngle = cos(degToRad(SampleThresholdAngle));
669 double dt = MinSampleInterval;
670 double end = start + t;
671 double current = start;
672
673 proc.sample(current, positionAtTime(current));
674
675 while (current < end)
676 {
677 double dt2 = dt;
678
679 Point3d goodpt;
680 double gooddt = dt;
681 Point3d pos0 = positionAtTime(current);
682 goodpt = positionAtTime(current + dt2);
683 while (1)
684 {
685 Point3d pos1 = positionAtTime(current + dt2);
686 Point3d pos2 = positionAtTime(current + dt2 * 2.0);
687 Vec3d vec1 = pos1 - pos0;
688 Vec3d vec2 = pos2 - pos0;
689
690 vec1.normalize();
691 vec2.normalize();
692 double dot = vec1 * vec2;
693
694 if (dot > 1.0)
695 dot = 1.0;
696 else if (dot < -1.0)
697 dot = -1.0;
698
699 if (dot > cosThresholdAngle && dt2 < MaxSampleInterval)
700 {
701 gooddt = dt2;
702 goodpt = pos1;
703 dt2 *= 2.0;
704 }
705 else
706 {
707 proc.sample(current + gooddt, goodpt);
708 break;
709 }
710 }
711 current += gooddt;
712 }
713 }
714
715
716 // Scan past comments. A comment begins with the # character and ends
717 // with a newline. Return true if the stream state is good. The stream
718 // position will be at the first non-comment, non-whitespace character.
SkipComments(istream & in)719 static bool SkipComments(istream& in)
720 {
721 bool inComment = false;
722 bool done = false;
723
724 int c = in.get();
725 while (!done)
726 {
727 if (in.eof())
728 {
729 done = true;
730 }
731 else
732 {
733 if (inComment)
734 {
735 if (c == '\n')
736 inComment = false;
737 }
738 else
739 {
740 if (c == '#')
741 {
742 inComment = true;
743 }
744 else if (!isspace(c))
745 {
746 in.unget();
747 done = true;
748 }
749 }
750 }
751
752 if (!done)
753 c = in.get();
754 }
755
756 return in.good();
757 }
758
759
760 // Load an ASCII xyz trajectory file. The file contains records with 4 double
761 // precision values each:
762 //
763 // 1: TDB time
764 // 2: Position x
765 // 3: Position y
766 // 4: Position z
767 //
768 // Positions are in kilometers.
769 //
770 // The numeric data may be preceeded by a comment block. Commented lines begin
771 // with a #; data is read start fromt the first non-whitespace character outside
772 // of a comment.
773
LoadSampledOrbit(const string & filename,TrajectoryInterpolation interpolation,T)774 template <typename T> SampledOrbit<T>* LoadSampledOrbit(const string& filename, TrajectoryInterpolation interpolation, T)
775 {
776 ifstream in(filename.c_str());
777 if (!in.good())
778 return NULL;
779
780 if (!SkipComments(in))
781 return NULL;
782
783 SampledOrbit<T>* orbit = NULL;
784
785 orbit = new SampledOrbit<T>(interpolation);
786
787 double lastSampleTime = -numeric_limits<double>::infinity();
788 while (in.good())
789 {
790 double tdb, x, y, z;
791 in >> tdb;
792 in >> x;
793 in >> y;
794 in >> z;
795
796 if (in.good())
797 {
798 // Skip samples with duplicate times; such trajectories are invalid, but
799 // are unfortunately used in some existing add-ons.
800 if (tdb != lastSampleTime)
801 {
802 orbit->addSample(tdb, x, y, z);
803 lastSampleTime = tdb;
804 }
805 }
806 }
807
808 return orbit;
809 }
810
811
812 // Load an xyzv sampled trajectory file. The file contains records with 7 double
813 // precision values:
814 //
815 // 1: TDB time
816 // 2: Position x
817 // 3: Position y
818 // 4: Position z
819 // 5: Velocity x
820 // 6: Velocity y
821 // 7: Velocity z
822 //
823 // Positions are in kilometers, velocities are kilometers per second.
824 //
825 // The numeric data may be preceeded by a comment block. Commented lines begin
826 // with a #; data is read start fromt the first non-whitespace character outside
827 // of a comment.
828
LoadSampledOrbitXYZV(const string & filename,TrajectoryInterpolation interpolation,T)829 template <typename T> SampledOrbitXYZV<T>* LoadSampledOrbitXYZV(const string& filename, TrajectoryInterpolation interpolation, T)
830 {
831 ifstream in(filename.c_str());
832 if (!in.good())
833 return NULL;
834
835 if (!SkipComments(in))
836 return NULL;
837
838 SampledOrbitXYZV<T>* orbit = NULL;
839
840 orbit = new SampledOrbitXYZV<T>(interpolation);
841
842 double lastSampleTime = -numeric_limits<double>::infinity();
843 while (in.good())
844 {
845 double tdb = 0.0;
846 Vec3d position;
847 Vec3d velocity;
848
849 in >> tdb;
850 in >> position.x;
851 in >> position.y;
852 in >> position.z;
853 in >> velocity.x;
854 in >> velocity.y;
855 in >> velocity.z;
856
857 // Convert velocities from km/sec to km/Julian day
858 velocity = velocity * astro::daysToSecs(1.0);
859
860 if (in.good())
861 {
862 if (tdb != lastSampleTime)
863 {
864 orbit->addSample(tdb, position, velocity);
865 lastSampleTime = tdb;
866 }
867 }
868 }
869
870 return orbit;
871 }
872
873
874 /*! Load a trajectory file containing single precision positions.
875 */
LoadSampledTrajectorySinglePrec(const string & filename,TrajectoryInterpolation interpolation)876 Orbit* LoadSampledTrajectorySinglePrec(const string& filename, TrajectoryInterpolation interpolation)
877 {
878 return LoadSampledOrbit(filename, interpolation, 0.0f);
879 }
880
881
882 /*! Load a trajectory file containing double precision positions.
883 */
LoadSampledTrajectoryDoublePrec(const string & filename,TrajectoryInterpolation interpolation)884 Orbit* LoadSampledTrajectoryDoublePrec(const string& filename, TrajectoryInterpolation interpolation)
885 {
886 return LoadSampledOrbit(filename, interpolation, 0.0);
887 }
888
889
890 /*! Load a trajectory file with single precision positions and velocities.
891 */
LoadXYZVTrajectorySinglePrec(const string & filename,TrajectoryInterpolation interpolation)892 Orbit* LoadXYZVTrajectorySinglePrec(const string& filename, TrajectoryInterpolation interpolation)
893 {
894 return LoadSampledOrbitXYZV(filename, interpolation, 0.0f);
895 }
896
897
898 /*! Load a trajectory file with double precision positions and velocities.
899 */
LoadXYZVTrajectoryDoublePrec(const string & filename,TrajectoryInterpolation interpolation)900 Orbit* LoadXYZVTrajectoryDoublePrec(const string& filename, TrajectoryInterpolation interpolation)
901 {
902 return LoadSampledOrbitXYZV(filename, interpolation, 0.0);
903 }
904