Frames

Untitled

0
1import * as THREE from 'three';
2import julian from 'julian';
3
4import { rescaleArray, rescaleXYZ } from './Scale';
5import { EphemerisTable } from './EphemerisTable';
6
7const sin = Math.sin;
8const cos = Math.cos;
9const sqrt = Math.sqrt;
10
11const DEFAULT_LEAD_TRAIL_YEARS = 10;
12const DEFAULT_SAMPLE_POINTS = 360;
13const DEFAULT_ORBIT_PATH_SETTINGS = {
14 leadDurationYears: DEFAULT_LEAD_TRAIL_YEARS,
15 trailDurationYears: DEFAULT_LEAD_TRAIL_YEARS,
16 numberSamplePoints: DEFAULT_SAMPLE_POINTS,
17};
18
19/**
20 * Special cube root function that assumes input is always positive.
21 */
22function cbrt(x) {
23 return Math.exp(Math.log(x) / 3.0);
24}
25
26/**
27 * Enum of orbital types.
28 */
29export const OrbitType = Object.freeze({
30 UNKNOWN: 0,
31 PARABOLIC: 1,
32 HYPERBOLIC: 2,
33 ELLIPTICAL: 3,
34 TABLE: 4,
35});
36
37/**
38 * Get the type of orbit. Returns one of OrbitType.PARABOLIC, HYPERBOLIC,
39 * ELLIPTICAL, or UNKNOWN.
40 * @return {OrbitType} Name of orbit type
41 */
42export function getOrbitType(ephem) {
43 if (ephem instanceof EphemerisTable) {
44 return OrbitType.TABLE;
45 }
46
47 let e = ephem.get('e');
48 if (e > 0.8 && e < 1.2) {
49 return OrbitType.PARABOLIC;
50 } else if (e > 1.2) {
51 return OrbitType.HYPERBOLIC;
52 } else {
53 return OrbitType.ELLIPTICAL;
54 }
55}
56
57/**
58 * A class that builds a visual representation of a Kepler orbit.
59 * @example
60 * const orbit = new Spacekit.Orbit({
61 * ephem: new Spacekit.Ephem({...}),
62 * options: {
63 * color: 0xFFFFFF,
64 * eclipticLineColor: 0xCCCCCC,
65 * },
66 * });
67 */
68export class Orbit {
69 /**
70 * @param {(Ephem | EphemerisTable)} ephem The ephemerides that define this orbit.
71 * @param {Object} options
72 * @param {Object} options.color The color of the orbital ellipse.
73 * @param {Object} options.orbitPathSettings settings for the path
74 * @param {Object} options.orbitPathSettings.leadDurationYears orbit path lead time in years
75 * @param {Object} options.orbitPathSettings.trailDurationYears orbit path trail time in years
76 * @param {Object} options.orbitPathSettings.numberSamplePoints number of points to use when drawing the orbit line
77 * Only applicable for non-elliptical and ephemeris table orbits.
78 * @param {Object} options.eclipticLineColor The color of lines drawn
79 * perpendicular to the ecliptic in order to illustrate depth (defaults to
80 * 0x333333).
81 */
82 constructor(ephem, options) {
83 /**
84 * Ephem object
85 * @type {Ephem}
86 */
87 this._ephem = ephem;
88
89 /**
90 * Options (see class definition for details)
91 */
92 this._options = options || {};
93
94 /**
95 * configuring orbit path lead/trail data
96 */
97 if (!this._options.orbitPathSettings) {
98 this._options.orbitPathSettings = JSON.parse(
99 JSON.stringify(DEFAULT_ORBIT_PATH_SETTINGS),
100 );
101 }
102
103 if (!this._options.orbitPathSettings.leadDurationYears) {
104 this._options.orbitPathSettings.leadDurationYears = DEFAULT_LEAD_TRAIL_YEARS;
105 }
106
107 if (!this._options.orbitPathSettings.trailDurationYears) {
108 this._options.orbitPathSettings.trailDurationYears = DEFAULT_LEAD_TRAIL_YEARS;
109 }
110
111 if (!this._options.orbitPathSettings.numberSamplePoints) {
112 this._options.orbitPathSettings.numberSamplePoints = DEFAULT_SAMPLE_POINTS;
113 }
114
115 /**
116 * Cached orbital points.
117 * @type {Array.<THREE.Vector3>}
118 */
119 this._orbitPoints = null;
120
121 /**
122 * Cached ecliptic drop lines.
123 * @type {Array.<THREE.Vector3>}
124 */
125 this._eclipticDropLines = null;
126
127 /**
128 * Cached orbit shape.
129 * @type {THREE.Line}
130 */
131 this._orbitShape = null;
132
133 /**
134 * Time span of the drawn orbit line
135 */
136 this._orbitStart = 0;
137 this._orbitStop = 0;
138 }
139
140 /**
141 * Get heliocentric position of object at a given JD.
142 * @param {Number} jd Date value in JD.
143 * @param {boolean} debug Set true for debug output.
144 * @return {Array.<Number>} [X, Y, Z] coordinates
145 */
146 getPositionAtTime(jd, debug) {
147 // Note: logic below must match the vertex shader.
148
149 // This position calculation is used to create orbital ellipses.
150 switch (getOrbitType(this._ephem)) {
151 case OrbitType.PARABOLIC:
152 return this.getPositionAtTimeNearParabolic(jd, debug);
153 case OrbitType.HYPERBOLIC:
154 return this.getPositionAtTimeHyperbolic(jd, debug);
155 case OrbitType.ELLIPTICAL:
156 return this.getPositionAtTimeElliptical(jd, debug);
157 case OrbitType.TABLE:
158 return this.getPositionAtTimeTable(jd, debug);
159 }
160 throw new Error('No handler for this type of orbit');
161 }
162
163 getPositionAtTimeParabolic(jd, debug) {
164 // See https://stjarnhimlen.se/comp/ppcomp.html#17
165 const eph = this._ephem;
166
167 // The Guassian gravitational constant
168 const k = 0.01720209895;
169
170 // Perihelion distance
171 const q = eph.get('q');
172
173 // Compute time since perihelion
174 const d = jd - eph.get('tp');
175
176 const H = (d * (k / sqrt(2))) / sqrt(q * q * q);
177 const h = 1.5 * H;
178 const g = sqrt(1.0 + h * h);
179 const s = cbrt(g + h) - cbrt(g - h);
180
181 // True anomaly
182 const v = 2.0 * Math.atan(s);
183 // Heliocentric distance
184 const r = q * (1.0 + s * s);
185
186 return this.vectorToHeliocentric(v, r);
187 }
188
189 getPositionAtTimeNearParabolic(jd, debug) {
190 // See https://stjarnhimlen.se/comp/ppcomp.html#17
191 const eph = this._ephem;
192
193 // The Guassian gravitational constant
194 const k = 0.01720209895;
195
196 // Eccentricity
197 const e = eph.get('e');
198
199 // Perihelion distance
200 const q = eph.get('q');
201
202 // Compute time since perihelion
203 const d = jd - eph.get('tp');
204
205 const a = 0.75 * d * k * sqrt((1 + e) / (q * q * q));
206 const b = sqrt(1 + a * a);
207 const W = cbrt(b + a) - cbrt(b - a);
208 const f = (1 - e) / (1 + e);
209
210 const a1 = 2 / 3 + (2 / 5) * W * W;
211 const a2 = 7 / 5 + (33 / 35) * W * W + (37 / 175) * W ** 4;
212 const a3 =
213 W * W * (432 / 175 + (956 / 1125) * W * W + (84 / 1575) * W ** 4);
214
215 const C = (W * W) / (1 + W * W);
216 const g = f * C * C;
217 const w = W * (1 + f * C * (a1 + a2 * g + a3 * g * g));
218
219 // True anomaly
220 const v = 2 * Math.atan(w);
221 // Heliocentric distance
222 const r = (q * (1 + w * w)) / (1 + w * w * f);
223
224 return this.vectorToHeliocentric(v, r);
225 }
226
227 getPositionAtTimeHyperbolic(jd, debug) {
228 // See https://stjarnhimlen.se/comp/ppcomp.html#17
229 const eph = this._ephem;
230
231 // Eccentricity
232 const e = eph.get('e');
233
234 // Perihelion distance
235 const q = eph.get('q');
236
237 // Semimajor axis
238 const a = eph.get('a');
239
240 // Mean anomaly
241 const ma = eph.get('ma');
242
243 // Calculate mean anomaly at jd
244 const n = eph.get('n', 'rad');
245 const epoch = eph.get('epoch');
246 const d = jd - epoch;
247
248 const M = ma + n * d;
249
250 let F0 = M;
251 for (let count = 0; count < 100; count++) {
252 const F1 =
253 (M + e * (F0 * Math.cosh(F0) - Math.sinh(F0))) /
254 (e * Math.cosh(F0) - 1);
255 const lastdiff = Math.abs(F1 - F0);
256 F0 = F1;
257
258 if (lastdiff < 0.0000001) {
259 break;
260 }
261 }
262 const F = F0;
263
264 const v = 2 * Math.atan(sqrt((e + 1) / (e - 1))) * Math.tanh(F / 2);
265 const r = (a * (1 - e * e)) / (1 + e * cos(v));
266
267 return this.vectorToHeliocentric(v, r);
268 }
269
270 getPositionAtTimeElliptical(jd, debug) {
271 const eph = this._ephem;
272
273 // Eccentricity
274 const e = eph.get('e');
275
276 // Mean anomaly
277 const ma = eph.get('ma', 'rad');
278
279 // Calculate mean anomaly at jd
280 const n = eph.get('n', 'rad');
281 const epoch = eph.get('epoch');
282 const d = jd - epoch;
283
284 const M = ma + n * d;
285 if (debug) {
286 console.info('period=', eph.get('period'));
287 console.info('n=', n);
288 console.info('ma=', ma);
289 console.info('d=', d);
290 console.info('M=', M);
291 }
292
293 // Estimate eccentric and true anom using iterative approx
294 let E0 = M;
295 for (let count = 0; count < 100; count++) {
296 const E1 = M + e * sin(E0);
297 const lastdiff = Math.abs(E1 - E0);
298 E0 = E1;
299
300 if (lastdiff < 0.0000001) {
301 break;
302 }
303 }
304 const E = E0;
305 const v = 2 * Math.atan(sqrt((1 + e) / (1 - e)) * Math.tan(E / 2));
306
307 // Radius vector, in AU
308 const a = eph.get('a');
309 const r = (a * (1 - e * e)) / (1 + e * cos(v));
310
311 return this.vectorToHeliocentric(v, r);
312 }
313
314 getPositionAtTimeTable(jd, debug) {
315 const point = this._ephem.getPositionAtTime(jd);
316 return rescaleXYZ(point[0], point[1], point[2]);
317 }
318
319 /**
320 * Given true anomaly and heliocentric distance, returns the scaled heliocentric coordinates (X, Y, Z)
321 * @param {Number} v True anomaly
322 * @param {Number} r Heliocentric distance
323 * @return {Array.<Number>} Heliocentric coordinates
324 */
325 vectorToHeliocentric(v, r) {
326 const eph = this._ephem;
327
328 // Inclination, Longitude of ascending node, Longitude of perihelion
329 const i = eph.get('i', 'rad');
330 const o = eph.get('om', 'rad');
331 const p = eph.get('wBar', 'rad');
332
333 // Heliocentric coords
334 const X = r * (cos(o) * cos(v + p - o) - sin(o) * sin(v + p - o) * cos(i));
335 const Y = r * (sin(o) * cos(v + p - o) + cos(o) * sin(v + p - o) * cos(i));
336 const Z = r * (sin(v + p - o) * sin(i));
337
338 return rescaleXYZ(X, Y, Z);
339 }
340
341 /**
342 * Returns whether the requested epoch is within the current orbit's definition
343 * @param jd
344 * @returns {boolean|boolean} true if it is within the orbit span, false if not
345 */
346 timeInRenderedOrbitSpan(jd) {
347 return jd >= this._orbitStart && jd <= this._orbitStop;
348 }
349
350 /**
351 * Calculates, caches, and returns the orbit state for this orbit around this time
352 * @param jd center time of the orbit (only used for ephemeris table ephemeris)
353 * @param forceCompute forces the recomputing of the orbit on this call
354 * @returns {THREE.Line}
355 */
356 getOrbitShape(jd, forceCompute = false) {
357 // For hyperbolic and parabolic orbits, decide on a time range to draw
358 // them.
359 // TODO(ian): Should we compute around current position, not time of perihelion?
360 const orbitType = getOrbitType(this._ephem);
361 const tp = orbitType === OrbitType.TABLE ? jd : this._ephem.get('tp');
362 const centerDate = tp ? tp : julian.toJulianDay(new Date());
363 const startJd =
364 centerDate - this._options.orbitPathSettings.trailDurationYears * 365.0;
365 const endJd =
366 centerDate + this._options.orbitPathSettings.leadDurationYears * 365.0;
367 const step =
368 (endJd - startJd) / this._options.orbitPathSettings.numberSamplePoints;
369
370 this._orbitStart = startJd;
371 this._orbitStop = endJd;
372
373 if (forceCompute) {
374 this._orbitShape = undefined;
375 this._eclipticDropLines = undefined;
376 }
377
378 if (this._orbitShape) {
379 return this._orbitShape;
380 }
381
382 switch (orbitType) {
383 case OrbitType.HYPERBOLIC:
384 return this.getLine(
385 this.getPositionAtTimeHyperbolic.bind(this),
386 startJd,
387 endJd,
388 step,
389 );
390 case OrbitType.PARABOLIC:
391 return this.getLine(
392 this.getPositionAtTimeNearParabolic.bind(this),
393 startJd,
394 endJd,
395 step,
396 );
397 case OrbitType.ELLIPTICAL:
398 return this.getEllipse();
399 case OrbitType.TABLE:
400 return this.getTableOrbit(startJd, endJd, step);
401 }
402 throw new Error('Unknown orbit shape');
403 }
404
405 /**
406 * Compute a line between a given date range.
407 * @private
408 */
409 getLine(orbitFn, startJd, endJd, step) {
410 const points = [];
411 for (let jd = startJd; jd <= endJd; jd += step) {
412 const pos = orbitFn(jd);
413 points.push(new THREE.Vector3(pos[0], pos[1], pos[2]));
414 }
415
416 const pointsGeometry = new THREE.Geometry();
417 pointsGeometry.vertices = points;
418
419 return this.generateAndCacheOrbitShape(pointsGeometry);
420 }
421
422 /**
423 * Returns the orbit for a table lookup orbit definition
424 * @private
425 * @param startJd start of orbit in JDate format
426 * @param stopJd end of orbit in JDate format
427 * @param step step size in days
428 * @returns {THREE.Line}
429 */
430 getTableOrbit(startJd, stopJd, step) {
431 const rawPoints = this._ephem.getPositions(startJd, stopJd, step);
432 const points = rawPoints
433 .map(values => rescaleArray(values))
434 .map(values => new THREE.Vector3(values[0], values[1], values[2]));
435 const pointGeometry = new THREE.Geometry();
436 pointGeometry.vertices = points;
437
438 return this.generateAndCacheOrbitShape(pointGeometry);
439 }
440
441 /**
442 * @private
443 * @return {THREE.Line} The ellipse object that represents this orbit.
444 */
445 getEllipse() {
446 const pointGeometry = this.getEllipsePoints();
447 return this.generateAndCacheOrbitShape(pointGeometry);
448 }
449
450 /**
451 * @private
452 * @return {Array.<THREE.Vector3>} List of points
453 */
454 getEllipsePoints() {
455 const eph = this._ephem;
456
457 const period = eph.get('period');
458 const ecc = eph.get('e');
459 const step = period / this._options.orbitPathSettings.numberSamplePoints;
460
461 const pts = [];
462 let prevPos;
463 for (let time = 0; time < period; time += step) {
464 const pos = this.getPositionAtTime(time);
465 if (isNaN(pos[0]) || isNaN(pos[1]) || isNaN(pos[2])) {
466 console.error(
467 'NaN position value - you may have bad or incomplete data in the following ephemeris:',
468 );
469 console.error(eph);
470 }
471 const vector = new THREE.Vector3(pos[0], pos[1], pos[2]);
472 if (
473 prevPos &&
474 Math.abs(prevPos[0] - pos[0]) +
475 Math.abs(prevPos[1] - pos[1]) +
476 Math.abs(prevPos[2] - pos[2]) >
477 120
478 ) {
479 // Don't render bogus or very large ellipses.
480 points.vertices = [];
481 return points;
482 }
483 prevPos = pos;
484 pts.push(vector);
485 }
486 pts.push(pts[0]);
487
488 const pointGeometry = new THREE.Geometry();
489 pointGeometry.vertices = pts;
490 return pointGeometry;
491 }
492
493 /**
494 * @private
495 */
496 generateAndCacheOrbitShape(pointGeometry) {
497 this._orbitPoints = pointGeometry;
498 this._orbitShape = new THREE.Line(
499 pointGeometry,
500 new THREE.LineBasicMaterial({
501 color: new THREE.Color(this._options.color || 0x444444),
502 }),
503 THREE.LineStrip,
504 );
505 return this._orbitShape;
506 }
507
508 /**
509 * A geometry containing line segments that run between the orbit ellipse and
510 * the ecliptic plane of the solar system. This is a useful visual effect
511 * that makes it easy to tell when an orbit goes below or above the ecliptic
512 * plane.
513 * @return {THREE.Geometry} A geometry with many line segments.
514 */
515 getLinesToEcliptic() {
516 if (this._eclipticDropLines) {
517 return this._eclipticDropLines;
518 }
519
520 if (!this._orbitPoints) {
521 this.getOrbitShape();
522 }
523 const points = this._orbitPoints;
524 const geometry = new THREE.Geometry();
525
526 points.vertices.forEach(vertex => {
527 geometry.vertices.push(vertex);
528 geometry.vertices.push(new THREE.Vector3(vertex.x, vertex.y, 0));
529 });
530
531 this._eclipticDropLines = new THREE.LineSegments(
532 geometry,
533 new THREE.LineBasicMaterial({
534 color: this._options.eclipticLineColor || 0x333333,
535 }),
536 THREE.LineStrip,
537 );
538
539 return this._eclipticDropLines;
540 }
541
542 /**
543 * Get the color of this orbit.
544 * @return {Number} The hexadecimal color of the orbital ellipse.
545 */
546 getHexColor() {
547 return this._orbitShape.material.color.getHex();
548 }
549
550 /**
551 * @param {Number} hexVal The hexadecimal color of the orbital ellipse.
552 */
553 setHexColor(hexVal) {
554 this._orbitShape.material.color = new THREE.Color(hexVal);
555 }
556
557 /**
558 * Get the visibility of this orbit.
559 * @return {boolean} Whether the orbital ellipse is visible. Note that
560 * although the ellipse may not be visible, it is still present in the
561 * underlying Scene and Simultation.
562 */
563 getVisibility() {
564 return this._orbitShape.visible;
565 }
566
567 /**
568 * Change the visibility of this orbit.
569 * @param {boolean} val Whether to show the orbital ellipse.
570 */
571 setVisibility(val) {
572 this._orbitShape.visible = val;
573 }
574}
575