diff --git a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts index 86c45c6..07537fe 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.spec.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.spec.ts @@ -4,7 +4,7 @@ import { describe, expect, it } from 'vitest'; import { DEFAULT_EPOCH_JD, GM_SUN_AU3_PER_DAY2 } from '../../shared/astro/constants'; import { keplerRates } from '../../shared/astro/kepler'; import { eclipticToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; -import { BodyRecord } from '../../shared/models/body.model'; +import { BodyRecord, RotationalElements } from '../../shared/models/body.model'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; import { SystemOrbitsRenderer } from './system-orbits-renderer'; @@ -377,8 +377,8 @@ describe('SystemOrbitsRenderer exoplanet propagation', () => { }); }); -describe('rotation', () => { - /** Earth, near enough: a day of 23.934 h, tipped 23.44 degrees off its orbit. */ +describe('rotation without IAU elements', () => { + /** A body with a day of 23.934 h and no pole: Eris, Haumea and Makemake are drawn this way. */ function spinning(overrides: Partial = {}): BodyRecord { return { id: 'earth', @@ -389,7 +389,6 @@ describe('rotation', () => { orbit: { semiMajorAxisAu: 1, eccentricity: 0.0167, inclinationDeg: 0, longitudeOfAscendingNodeDeg: 0, argumentOfPeriapsisDeg: 0, meanAnomalyAtEpochDeg: 0, epochJd: DEFAULT_EPOCH_JD }, rates: keplerRates(1, GM_SUN_AU3_PER_DAY2), orbitSource: 'test', rotationPeriodHours: 23.934, - obliquityDeg: 23.4392911, ...overrides }; } @@ -428,18 +427,9 @@ describe('rotation', () => { return axis.normalize().dot(new THREE.Vector3(0, 0, 1).applyQuaternion(renderer.referenceFrame)); } - it('turns Venus backwards, as Horizons gives it: a negative rate and an obliquity past 90', () => { - // Both say retrograde, in two conventions. Applied together they cancelled into a forward - // turn, which is how Venus and Uranus used to be drawn. - const venus = spinning({ id: 'venus', rotationPeriodHours: -5832.54, obliquityDeg: 177.3 }); - - expect(spinSense(spinning())).toBeGreaterThan(0.9); - expect(spinSense(venus)).toBeLessThan(-0.9); - }); - - it('reads the sign of the period only where no obliquity says which way the pole points', () => { - expect(spinSense(spinning({ rotationPeriodHours: -23.934, obliquityDeg: undefined }))).toBeLessThan(-0.9); - expect(spinSense(spinning({ rotationPeriodHours: 23.934, obliquityDeg: undefined }))).toBeGreaterThan(0.9); + it('turns it about its orbit’s normal, backwards for a negative period', () => { + expect(spinSense(spinning({ rotationPeriodHours: -23.934 }))).toBeLessThan(-0.99); + expect(spinSense(spinning({ rotationPeriodHours: 23.934 }))).toBeGreaterThan(0.99); }); it('leaves a body with no published rotation still', () => { @@ -514,6 +504,18 @@ describe('solar-system bodies against Horizons', () => { uranus: {kind: 'planet', orbit: {semiMajorAxisAu: 19.18797948, eccentricity: 0.0468574, inclinationDeg: 0.77298127, longitudeOfAscendingNodeDeg: 73.96250215, argumentOfPeriapsisDeg: 98.47154226, meanAnomalyAtEpochDeg: 141.76872184, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.011731557178644764, longitudeOfAscendingNodeDegPerDay: 0.0000015714439425051334, argumentOfPeriapsisDegPerDay: 9.65718275154004e-7, semiMajorAxisAuPerDay: -5.600273785078713e-9, eccentricityPerDay: -4.2436687200547574e-10, inclinationDegPerDay: -4.932375085557837e-8, meanAnomalyTerms: {b: 0.00058331, c: -0.97731848, s: 0.17689245, f: 7.67025}}}, titania: {kind: 'moon', orbit: {semiMajorAxisAu: 0.002916485361445723, eccentricity: 0.0011, inclinationDeg: 0.079, longitudeOfAscendingNodeDeg: 279.771, argumentOfPeriapsisDeg: 284.4, meanAnomalyAtEpochDeg: 24.614, epochJd: 2444239.5}, rates: {meanMotionDegPerDay: 41.3514246, longitudeOfAscendingNodeDegPerDay: -0.005044947168524978, argumentOfPeriapsisDegPerDay: 0.006102004540272753}, laplacePole: {raDeg: 77.311, decDeg: 15.175}, parentBodyId: 'uranus'}, charon: {kind: 'moon', orbit: {semiMajorAxisAu: 0.00013095774631236113, eccentricity: 0.0002, inclinationDeg: 0.08, longitudeOfAscendingNodeDeg: 26.928, argumentOfPeriapsisDeg: 146.106, meanAnomalyAtEpochDeg: 131.07, epochJd: 2451545}, rates: {meanMotionDegPerDay: 56.362521, longitudeOfAscendingNodeDegPerDay: -0.00010926638529337138, argumentOfPeriapsisDegPerDay: 0.00009683851540842405}, laplacePole: {raDeg: 132.993, decDeg: -6.163}, parentBodyId: 'pluto', massRatio: 0.1220485755631374}, + venus: {kind: 'planet', orbit: {semiMajorAxisAu: 0.72332102, eccentricity: 0.00676399, inclinationDeg: 3.39777545, longitudeOfAscendingNodeDeg: 76.67261496, argumentOfPeriapsisDeg: 55.094942169999996, meanAnomalyAtEpochDeg: 50.21215136999999, epochJd: 2451545}, rates: {meanMotionDegPerDay: 1.6021304750882956, longitudeOfAscendingNodeDegPerDay: -0.000007467261875427789, argumentOfPeriapsisDegPerDay: 0.000009022264750171116, semiMajorAxisAuPerDay: -7.118412046543463e-12, eccentricityPerDay: -1.398220396988364e-9, inclinationDegPerDay: 1.1908008213552361e-8}}, + mars: {kind: 'planet', orbit: {semiMajorAxisAu: 1.52371243, eccentricity: 0.09336511, inclinationDeg: 1.85181869, longitudeOfAscendingNodeDeg: 49.71320984, argumentOfPeriapsisDeg: -73.63065768, meanAnomalyAtEpochDeg: 19.3493162, epochJd: 2451545}, rates: {meanMotionDegPerDay: 0.5240328362061601, longitudeOfAscendingNodeDegPerDay: -0.000007351794934976044, argumentOfPeriapsisDegPerDay: 0.00001973334866529774, semiMajorAxisAuPerDay: 2.6557152635181385e-11, eccentricityPerDay: 2.5048596851471597e-9, inclinationDegPerDay: -1.9842765229295004e-7}}, + }; + // The IAU WGCCRE 2015 rotational elements bodies.json carries for them, from pck00011.tpc. + const ROTATION: Record = { + venus: {poleRaDeg: [272.76, 0, 0], poleDecDeg: [67.16, 0, 0], primeMeridianDeg: [160.2, -1.4813688, 0]}, + earth: {poleRaDeg: [0, -0.641, 0], poleDecDeg: [90, -0.557, 0], primeMeridianDeg: [190.147, 360.9856235, 0]}, + mars: {poleRaDeg: [317.269202, -0.10927547, 0], poleDecDeg: [54.432516, -0.05827105, 0], primeMeridianDeg: [176.049863, 350.891982443297, 0], terms: [{angleDeg: [79.398797, 0.5042615, 0], ra: 0.419057, dec: 0, pm: 0}, {angleDeg: [166.325722, 0.5042615, 0], ra: 0, dec: 1.591274, pm: 0}, {angleDeg: [95.391654, 0.5042615, 0], ra: 0, dec: 0, pm: 0.584542}]}, + jupiter: {poleRaDeg: [268.056595, -0.006499, 0], poleDecDeg: [64.495303, 0.002413, 0], primeMeridianDeg: [284.95, 870.536, 0]}, + uranus: {poleRaDeg: [257.311, 0, 0], poleDecDeg: [-15.175, 0, 0], primeMeridianDeg: [203.81, -501.1600928, 0]}, + pluto: {poleRaDeg: [132.993, 0, 0], poleDecDeg: [-6.163, 0, 0], primeMeridianDeg: [302.695, 56.3625225, 0]}, + moon: {poleRaDeg: [269.9949, 0.0031, 0], poleDecDeg: [66.5392, 0.013, 0], primeMeridianDeg: [38.3213, 13.17635815, -1.4e-12], terms: [{angleDeg: [125.045, -1935.5364525], ra: -3.8787, dec: 1.5419, pm: 3.561}, {angleDeg: [250.089, -3871.072905], ra: -0.1204, dec: 0.0239, pm: 0.1208}, {angleDeg: [260.008, 475263.3328725], ra: 0.07, dec: -0.0278, pm: -0.0642}, {angleDeg: [176.625, 487269.629985], ra: -0.0172, dec: 0.0068, pm: 0.0158}, {angleDeg: [357.529, 35999.0509575], ra: 0, dec: 0, pm: 0.0252}]}, }; // Each ceiling sits just above what these elements measure on that date: Earth 0.003 degrees, // Jupiter 0.063, Saturn 0.164, Pluto 0.054, the Moon 0.72 (no mean ellipse has its evection or @@ -534,7 +536,7 @@ describe('solar-system bodies against Horizons', () => { ]; function record(id: string): BodyRecord { - return { id, systemStarId: 0, name: id, radiusKm: 1000, orbitSource: 'test', ...RECORDS[id] }; + return { id, systemStarId: 0, name: id, radiusKm: 1000, orbitSource: 'test', ...RECORDS[id], rotationalElements: ROTATION[id] }; } const renderer = new SystemOrbitsRenderer(Object.keys(RECORDS).map(record), []); @@ -589,4 +591,95 @@ describe('solar-system bodies against Horizons', () => { expect(Math.abs(moon.position.clone().normalize().dot(normal))).toBeLessThan(1e-9); } }); + + const JUNE_1_2025_NOON_UTC = 2460828.0; + + /** + * The tilt of a body's drawn spin from the orbit it is drawn going round, in degrees: its angular + * velocity, read off the sphere a quarter of an hour apart, against its orbit line's normal. Past + * 90 is a body turning backwards against its orbit. + */ + function drawnObliquity(id: string): number { + const marker = renderer.members.find((member) => member.id === id)!.marker; + const line = renderer.object.children[renderer.object.children.indexOf(marker) - 1]; + expect(line.name).toBe('orbit-line'); + renderer.update(JUNE_1_2025_NOON_UTC); + const start = marker.quaternion.clone(); + renderer.update(JUNE_1_2025_NOON_UTC + 0.01); + const turn = marker.quaternion.clone().multiply(start.invert()); + const spin = new THREE.Vector3(turn.x, turn.y, turn.z).multiplyScalar(Math.sign(turn.w)); + return (spin.angleTo(new THREE.Vector3(0, 0, 1).applyQuaternion(line.quaternion)) * 180) / Math.PI; + } + + it('turns Venus, Uranus and Pluto backwards against their orbits, at the tilts Horizons gives', () => { + // The IAU names a planet's north pole by the side of the solar system it lies on, so Venus's W + // and Uranus's run backwards; Pluto's pole follows the right-hand rule instead, and points + // south. Either way the spin read off the drawn sphere is past 90 degrees from the orbit's pole. + expect(drawnObliquity('venus')).toBeCloseTo(177.3, 0); + expect(drawnObliquity('uranus')).toBeCloseTo(97.77, 0); + expect(drawnObliquity('pluto')).toBeCloseTo(119.6, 0); + expect(drawnObliquity('earth')).toBeCloseTo(23.44, 0); + }); + + /** + * Where on its drawn sphere a body faces a point, as east longitude and latitude on its map: read + * from the texture coordinates where a ray from that point meets the sphere, so the map's own + * convention is part of what is measured. + */ + function facing(id: string, point: THREE.Vector3): { eastDeg: number; latDeg: number } { + const marker = renderer.members.find((member) => member.id === id)!.marker as THREE.Mesh; + const centre = worldPosition(id); + const towards = point.clone().sub(centre).normalize(); + const radius = (marker.geometry as THREE.SphereGeometry).parameters.radius; + const hit = new THREE.Raycaster(centre.clone().addScaledVector(towards, radius * 4), towards.clone().negate()).intersectObject(marker)[0]; + return { eastDeg: (hit.uv!.x - 0.5) * 360, latDeg: (hit.uv!.y - 0.5) * 180 }; + } + + function worldPosition(id: string): THREE.Vector3 { + const marker = renderer.members.find((member) => member.id === id)!.marker; + marker.updateWorldMatrix(true, false); + return marker.getWorldPosition(new THREE.Vector3()); + } + + /** Degrees between two longitudes, the short way round. */ + const apart = (a: number, b: number): number => Math.abs(((((a - b) % 360) + 540) % 360) - 180); + + it('lights Earth where the Sun really stands: within 4 degrees of Greenwich at noon UTC', () => { + // The equation of time is all that separates them: on 1 June 2025 it puts the Sun over 0.53 W, + // and the drawn sphere has it over 0.43 W. + renderer.update(JUNE_1_2025_NOON_UTC); + expect(Math.abs(facing('earth', new THREE.Vector3()).eastDeg)).toBeLessThan(4); + }); + + // Horizons' observer quantities 14 and 15 at 2025-06-01 12:00 UTC, from Earth's centre (from the + // Sun's, for Earth): the sub-observer and sub-solar longitude and latitude, east-positive for + // Earth and the Moon and west-positive for Mars and Jupiter, as each is printed. Horizons gives + // each body as it was when the light now arriving left it, so it is drawn that much earlier. Its + // latitudes are planetodetic, on the body's flattened figure, which a sphere does not have, so the + // drawn latitude is put on that figure before they are compared: without it they differ by what + // the flattening makes of them, 0.14 degrees on Earth, 0.26 on Mars and 0.33 on Jupiter. + // + // Measured: every longitude within 0.09 degrees and every latitude within 0.03, but for the + // Moon's face towards Earth, 0.70 and 0.09 out because its mean orbit is (its evection alone is + // 1.27 degrees); its face towards the Sun is within 0.002. + const SUB_POINTS: Array<[id: string, observer: string | undefined, lightMinutes: number, west: boolean, flattening: number, observerLon: number, observerLat: number, sunLon: number, sunLat: number, maxObserverDeg: number]> = [ + ['earth', undefined, 8.43351424, false, 1 / 298.257, 1.5855, 22.261204, 1.579501, 22.260426, 0.1], + ['mars', 'earth', 14.13295841, true, 1 - 3376.2 / 3396.19, 307.365389, 21.27653, 269.287887, 25.451264, 0.1], + ['moon', 'earth', 0.02150549, false, 0, 7.256763, -3.462104, 116.285934, 1.503004, 0.8], + ['jupiter', 'earth', 50.70337676, true, 1 - 66854 / 71492, 251.139846, 2.58787, 247.855871, 2.572658, 0.1] + ]; + + for (const [id, observer, lightMinutes, west, flattening, observerLon, observerLat, sunLon, sunLat, maxObserverDeg] of SUB_POINTS) { + it(`faces ${observer ?? 'the Sun'} and the Sun with the points Horizons gives on ${id}`, () => { + renderer.update(JUNE_1_2025_NOON_UTC - lightMinutes / 1440); + const seen = facing(id, observer ? worldPosition(observer) : new THREE.Vector3()); + const lit = facing(id, new THREE.Vector3()); + const east = (longitude: number): number => (west ? -longitude : longitude); + const planetodetic = (latDeg: number): number => (Math.atan(Math.tan((latDeg * Math.PI) / 180) / (1 - flattening) ** 2) * 180) / Math.PI; + expect(apart(seen.eastDeg, east(observerLon))).toBeLessThan(maxObserverDeg); + expect(Math.abs(planetodetic(seen.latDeg) - observerLat)).toBeLessThan(maxObserverDeg); + expect(apart(lit.eastDeg, east(sunLon))).toBeLessThan(0.1); + expect(Math.abs(planetodetic(lit.latDeg) - sunLat)).toBeLessThan(0.05); + }); + } }); diff --git a/src/app/features/galaxy-system/system-orbits-renderer.ts b/src/app/features/galaxy-system/system-orbits-renderer.ts index bd1996e..2c55b24 100644 --- a/src/app/features/galaxy-system/system-orbits-renderer.ts +++ b/src/app/features/galaxy-system/system-orbits-renderer.ts @@ -5,8 +5,9 @@ import { PlanetAppearance } from '../../shared/astro/planet-appearance'; import { planetTexture } from '../../shared/rendering/procedural-planet-texture'; import { bodyTexturePath, loadCachedTexture } from '../../shared/rendering/texture-catalog'; import { isPropagatableOrbit, keplerRates, meanElementsAt, orbitEllipsePoints, positionAtEpoch, resolveGravitationalParameter, resolveOrbitalElements } from '../../shared/astro/kepler'; -import { CartesianCoordinates, laplacePlaneToEquatorial, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; -import { BodyRecord, MeanElementRates, OrbitalElements } from '../../shared/models/body.model'; +import { CartesianCoordinates, OBLIQUITY_J2000_DEG } from '../../shared/astro/coordinates'; +import { BodyRecord, MeanElementRates, OrbitalElements, RotationalElements } from '../../shared/models/body.model'; +import { bodyOrientation, poleFrame } from '../../shared/rendering/body-orientation'; import { bodyMarkerRadiusAu, systemGridRingsAu } from './system-framing'; import { PolarGridPlane, TetherField } from './grid-plane'; import { ExoplanetRecord } from '../../shared/models/exoplanet.model'; @@ -51,19 +52,12 @@ const ECLIPTIC_FRAME = new THREE.Quaternion().setFromAxisAngle(new THREE.Vector3 /** * Rotation carrying a moon's element frame into the scene: its local Laplace plane where JPL - * gives one, the ecliptic otherwise. Built from the three axes {@link laplacePlaneToEquatorial} - * sends, so the scene and the ETL's check against Horizons share the one conversion. + * gives one, the ecliptic otherwise. {@link poleFrame} builds it from the axes + * `laplacePlaneToEquatorial` sends, so the scene and the ETL's check against Horizons share the + * one conversion. */ function moonFrame(body: BodyRecord): THREE.Quaternion { - const pole = body.laplacePole; - if (!pole) { - return ECLIPTIC_FRAME.clone(); - } - const axis = (x: number, y: number, z: number): THREE.Vector3 => { - const turned = laplacePlaneToEquatorial({ x, y, z }, pole); - return new THREE.Vector3(turned.x, turned.y, turned.z); - }; - return new THREE.Quaternion().setFromRotationMatrix(new THREE.Matrix4().makeBasis(axis(1, 0, 0), axis(0, 1, 0), axis(0, 0, 1))); + return body.laplacePole ? poleFrame(body.laplacePole) : ECLIPTIC_FRAME.clone(); } const X_AXIS = new THREE.Vector3(1, 0, 0); @@ -253,32 +247,20 @@ const SPIN_AXIS = new THREE.Vector3(0, 1, 0); const HOURS_PER_DAY = 24; /** - * How a body is turned at a given date: its own sidereal rotation, about its own axis. + * How a body the IAU gives no rotational elements for is turned at a given date — Eris, Haumea + * and Makemake, whose periods are known and whose poles are not: at its own sidereal rate, about + * its orbit's normal, backwards for a negative period. None of them has an obliquity, so none is + * applied. The phase is arbitrary: each body starts at its elements' epoch in the shortest + * rotation of +Y onto its axis, and turns from there. Exoplanets have no published rotation at + * all, and are left still. * - * The obliquity fixes how far the pole leans from the orbit normal, and nothing more: which way - * it leans needs the pole's right ascension, which the Horizons pages this reads do not carry. The - * lean is taken about the orbit's ascending node because that is the one line the elements name, - * not because the data says so — so the tilt is real and its azimuth is not. It is the node on the - * date drawn, so the lean follows the orbit as the node turns. Likewise the phase: each body - * starts at its elements' epoch (J2000 for the planets) in an arbitrary orientation, the - * shortest rotation of +Y onto its axis, and turns from there. The rate and the sense are real; - * the face towards the camera is not. - * - * Horizons states a retrograde spin twice over, in two conventions: an obliquity past 90 degrees - * (Venus 177.3, Uranus 97.8) and a negative rate. Either one alone turns the body backwards, and - * both together cancel into a forward turn — which is how Venus and Uranus were drawn. Where an - * obliquity is given it carries the sense, and the period is taken as a magnitude; the sign of the - * period is only read for a body with no obliquity at all. + * Every other body is turned by {@link bodyOrientation}. */ -function spinFor(elements: OrbitalElements, frame: THREE.Quaternion, rotationPeriodHours: number, obliquityDeg: number | undefined, daysSinceEpoch: number): THREE.Quaternion { +function spinFor(elements: OrbitalElements, frame: THREE.Quaternion, rotationPeriodHours: number, daysSinceEpoch: number): THREE.Quaternion { const node = elements.longitudeOfAscendingNodeDeg * DEG_TO_RAD; const inclination = elements.inclinationDeg * DEG_TO_RAD; - const nodeDirection = new THREE.Vector3(Math.cos(node), Math.sin(node), 0); - const axis = new THREE.Vector3(Math.sin(inclination) * Math.sin(node), -Math.sin(inclination) * Math.cos(node), Math.cos(inclination)) - .applyAxisAngle(nodeDirection, (obliquityDeg ?? 0) * DEG_TO_RAD) - .applyQuaternion(frame); - const period = obliquityDeg === undefined ? rotationPeriodHours : Math.abs(rotationPeriodHours); - const turns = (daysSinceEpoch * HOURS_PER_DAY) / period; + const axis = new THREE.Vector3(Math.sin(inclination) * Math.sin(node), -Math.sin(inclination) * Math.cos(node), Math.cos(inclination)).applyQuaternion(frame); + const turns = (daysSinceEpoch * HOURS_PER_DAY) / rotationPeriodHours; return new THREE.Quaternion() .setFromUnitVectors(SPIN_AXIS, axis) .multiply(new THREE.Quaternion().setFromAxisAngle(SPIN_AXIS, turns * 2 * Math.PI)); @@ -297,7 +279,7 @@ interface TrackedTopLevelBody { position: THREE.Vector3; /** Sidereal rotation, where the catalogue publishes one; negative is retrograde. */ rotationPeriodHours?: number; - obliquityDeg?: number; + rotationalElements?: RotationalElements; } interface TrackedMoon { @@ -310,7 +292,7 @@ interface TrackedMoon { pivot: THREE.Group; parentId: string; rotationPeriodHours?: number; - obliquityDeg?: number; + rotationalElements?: RotationalElements; /** * Where the moon and its planet go round a barycentre outside the planet (Charon): the moon's * mass over the planet's, and the planet's own small orbit round that point. @@ -398,7 +380,7 @@ export class SystemOrbitsRenderer { } // A body reaches here only when it has no parentBodyId, so `kind` is 'planet' or 'dwarf'. const kind: SystemMemberKind = body.kind; - const tracked = this.addTopLevelBody(body.id, kind, body.orbit, body.rates, body.radiusKm, ECLIPTIC_FRAME, appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, obliquityDeg: body.obliquityDeg }); + const tracked = this.addTopLevelBody(body.id, kind, body.orbit, body.rates, body.radiusKm, ECLIPTIC_FRAME, appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, elements: body.rotationalElements }); members.push({ id: body.id, kind, marker: tracked.marker }); } @@ -411,7 +393,7 @@ export class SystemOrbitsRenderer { if (!parentTracked) { continue; // orphaned moon reference; skip rather than crash. } - const moon = this.addMoon(body.id, body.orbit, body.rates, body.radiusKm, parentTracked, moonFrame(body), appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, obliquityDeg: body.obliquityDeg }, body.massRatio); + const moon = this.addMoon(body.id, body.orbit, body.rates, body.radiusKm, parentTracked, moonFrame(body), appearanceForBody(body, bodies, hostLuminositySolar), { periodHours: body.rotationPeriodHours, elements: body.rotationalElements }, body.massRatio); members.push({ id: body.id, kind: 'moon', marker: moon.marker, parentId: parent.id }); } @@ -484,8 +466,10 @@ export class SystemOrbitsRenderer { body.position.set(orbital.x, orbital.y, orbital.z).applyQuaternion(body.frame); body.marker.position.copy(body.position); orientOrbit(body.orbitLine.quaternion, current, body.frame); - if (body.rotationPeriodHours) { - body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, body.obliquityDeg, epochJd - body.elements.epochJd)); + if (body.rotationalElements) { + bodyOrientation(body.rotationalElements, epochJd, body.marker.quaternion); + } else if (body.rotationPeriodHours) { + body.marker.quaternion.copy(spinFor(current, body.frame, body.rotationPeriodHours, epochJd - body.elements.epochJd)); } } @@ -507,8 +491,10 @@ export class SystemOrbitsRenderer { moon.marker.position.multiplyScalar(1 / (1 + massRatio)); parentOrbitLine.quaternion.copy(moon.orbitLine.quaternion); } - if (moon.rotationPeriodHours) { - moon.marker.quaternion.copy(spinFor(current, moon.frame, moon.rotationPeriodHours, moon.obliquityDeg, epochJd - moon.elements.epochJd)); + if (moon.rotationalElements) { + bodyOrientation(moon.rotationalElements, epochJd, moon.marker.quaternion); + } else if (moon.rotationPeriodHours) { + moon.marker.quaternion.copy(spinFor(current, moon.frame, moon.rotationPeriodHours, epochJd - moon.elements.epochJd)); } } @@ -567,7 +553,7 @@ export class SystemOrbitsRenderer { radiusKm: number | undefined, frame: THREE.Quaternion, appearance?: PlanetAppearance, - rotation?: { periodHours?: number; obliquityDeg?: number } + rotation?: { periodHours?: number; elements?: RotationalElements } ): TrackedTopLevelBody { const orbitLine = buildOrbitLine(elements, kind, frame); const marker = buildMarker(id, kind, radiusKm, appearance, this.deferSurface); @@ -575,7 +561,7 @@ export class SystemOrbitsRenderer { this.trackDisposable(orbitLine.geometry, orbitLine.material as THREE.Material); this.trackDisposable(marker.geometry, marker.material as THREE.Material); - const tracked: TrackedTopLevelBody = { id, kind, elements, rates, marker, orbitLine, frame, position: new THREE.Vector3(), rotationPeriodHours: rotation?.periodHours, obliquityDeg: rotation?.obliquityDeg }; + const tracked: TrackedTopLevelBody = { id, kind, elements, rates, marker, orbitLine, frame, position: new THREE.Vector3(), rotationPeriodHours: rotation?.periodHours, rotationalElements: rotation?.elements }; this.topLevelBodies.push(tracked); return tracked; } @@ -588,7 +574,7 @@ export class SystemOrbitsRenderer { parent: TrackedTopLevelBody, frame: THREE.Quaternion, appearance?: PlanetAppearance, - rotation?: { periodHours?: number; obliquityDeg?: number }, + rotation?: { periodHours?: number; elements?: RotationalElements }, massRatio?: number ): TrackedMoon { const pivot = new THREE.Group(); @@ -612,7 +598,7 @@ export class SystemOrbitsRenderer { barycentre = { massRatio, parentOrbitLine }; } - const moon: TrackedMoon = { id, elements, rates, marker, orbitLine, frame, pivot, parentId: parent.id, rotationPeriodHours: rotation?.periodHours, obliquityDeg: rotation?.obliquityDeg, barycentre }; + const moon: TrackedMoon = { id, elements, rates, marker, orbitLine, frame, pivot, parentId: parent.id, rotationPeriodHours: rotation?.periodHours, rotationalElements: rotation?.elements, barycentre }; this.moons.push(moon); return moon; } diff --git a/src/app/shared/astro/constants.ts b/src/app/shared/astro/constants.ts index 9ccb835..41e2569 100644 --- a/src/app/shared/astro/constants.ts +++ b/src/app/shared/astro/constants.ts @@ -19,6 +19,14 @@ export const DEFAULT_EPOCH_JD = 2451545.0; */ export const GM_SUN_AU3_PER_DAY2 = 0.01720209895 * 0.01720209895; +/** + * TT - UTC, in days: 32.184 s plus the 37 leap seconds UTC has taken since 1972, the last at the + * end of 2016. TDB, which ephemerides run on, stays within 2 ms of TT. Held constant, as Horizons + * holds it for dates past the last announced leap second; before 2017 it was smaller, about 29 s + * in 1950. + */ +export const TT_MINUS_UTC_DAYS = 69.184 / 86400; + /** Converts a JS `Date` into a Julian date (days), for driving the Kepler propagator "now". */ export function dateToJulianDate(date: Date = new Date()): number { return date.getTime() / 86400000 + 2440587.5; diff --git a/src/app/shared/rendering/body-orientation.ts b/src/app/shared/rendering/body-orientation.ts new file mode 100644 index 0000000..763460f --- /dev/null +++ b/src/app/shared/rendering/body-orientation.ts @@ -0,0 +1,60 @@ +import * as THREE from 'three/webgpu'; + +import { TT_MINUS_UTC_DAYS } from '../astro/constants'; +import { laplacePlaneToEquatorial } from '../astro/coordinates'; +import { orientationAt } from '../astro/rotational-elements'; +import { RotationalElements } from '../models/body.model'; + +const DEG_TO_RAD = Math.PI / 180; +const Z_AXIS = new THREE.Vector3(0, 0, 1); + +/** + * How a surface map sits on a sphere, settled once for every map the app wraps. + * + * `THREE.SphereGeometry` is built round +Y and runs its u coordinate eastward, anticlockwise seen + * from +Y, from a seam on -X: u = 0.5 faces +X and u = 0.75 faces -Z. Every photograph in + * `texture-catalog.ts` is an equirectangular map centred on longitude 0 with east to the right — + * Greenwich is in the middle of Earth's; on Mars's, Olympus Mons (226.2 E, which is -133.8) sits a + * little over a third of the width left of centre; on the Moon's, Mare Crisium (59 E) is right of + * centre and Mare Orientale (95 W) left of it; on Mercury's, the rayed crater Kuiper (31.5 W, 11 S) + * is just left of centre and below the equator. A map labelled in west longitude, as most planets' + * are, is still drawn with east to the right, as any map of a sphere seen from outside is; only its + * numbers run the other way. So longitude 0 is +X and 90 E is -Z, and a quarter turn about X + * carries that onto the IAU's body-fixed frame: pole +Z, prime meridian +X, 90 E +Y. + * + * The derived surfaces have no meridian of their own, and take the same convention. + */ +export const MAP_TO_BODY = new THREE.Quaternion().setFromAxisAngle(new THREE.Vector3(1, 0, 0), Math.PI / 2); + +const scratchMatrix = new THREE.Matrix4(); +const scratchAxes = [new THREE.Vector3(), new THREE.Vector3(), new THREE.Vector3()]; +const scratchTurn = new THREE.Quaternion(); + +/** + * Sets `target` to the rotation carrying a frame whose +Z is `pole` and whose +X is where its + * equator rises through the ICRF equator into the ICRF: the frame the IAU counts W in, and the + * one JPL refers a moon's Laplace plane to — so both go through {@link laplacePlaneToEquatorial}. + */ +export function poleFrame(pole: { raDeg: number; decDeg: number }, target = new THREE.Quaternion()): THREE.Quaternion { + const [x, y, z] = scratchAxes.map((axis, index) => { + const turned = laplacePlaneToEquatorial({ x: index === 0 ? 1 : 0, y: index === 1 ? 1 : 0, z: index === 2 ? 1 : 0 }, pole); + return axis.set(turned.x, turned.y, turned.z); + }); + return target.setFromRotationMatrix(scratchMatrix.makeBasis(x, y, z)); +} + +/** + * Sets `target` to the rotation carrying a sphere, wrapped in its map as `SphereGeometry` wraps it, + * into the ICRF at the map's own clock: the map onto the body's frame, turned by W about the pole, + * and on to where the pole points. + * + * The clock is UTC and the IAU's elements run on TDB, 69.184 s ahead; in that time Earth turns + * 0.29 degrees, Jupiter 0.70 and Phobos 0.90, so the difference is added here. + */ +export function bodyOrientation(elements: RotationalElements, jdUtc: number, target = new THREE.Quaternion()): THREE.Quaternion { + const { poleRaDeg, poleDecDeg, primeMeridianDeg } = orientationAt(elements, jdUtc + TT_MINUS_UTC_DAYS); + return poleFrame({ raDeg: poleRaDeg, decDeg: poleDecDeg }, target) + .multiply(scratchTurn.setFromAxisAngle(Z_AXIS, primeMeridianDeg * DEG_TO_RAD)) + .multiply(MAP_TO_BODY); +} +