|
| 1 | +/** |
| 2 | + * Solar ephemeris — Meeus, "Astronomical Algorithms" Ch. 25. |
| 3 | + * Computes apparent ecliptic longitude, apparent obliquity, and distance. |
| 4 | + * Accurate to ~0.01°. |
| 5 | + * No THREE.js dependency — safe for Web Workers. |
| 6 | + */ |
| 7 | +import { DEG2RAD } from '../constants'; |
| 8 | +import { epochToJulianDateTT, normalizeEpoch } from './epoch'; |
| 9 | + |
| 10 | +export interface SunEcliptic { |
| 11 | + /** Apparent ecliptic longitude (radians), includes nutation + aberration */ |
| 12 | + lambdaApp: number; |
| 13 | + /** Apparent obliquity of the ecliptic (radians), includes nutation */ |
| 14 | + epsilonApp: number; |
| 15 | + /** Distance from Earth (AU) */ |
| 16 | + rAU: number; |
| 17 | +} |
| 18 | + |
| 19 | +export function computeSunEcliptic(epoch: number): SunEcliptic { |
| 20 | + epoch = normalizeEpoch(epoch); |
| 21 | + const jd = epochToJulianDateTT(epoch); |
| 22 | + const T = (jd - 2451545.0) / 36525.0; // Julian centuries from J2000.0 |
| 23 | + const T2 = T * T; |
| 24 | + const T3 = T2 * T; |
| 25 | + |
| 26 | + // Geometric mean longitude of the Sun (degrees) — Meeus eq. 25.2 |
| 27 | + let L0 = 280.46646 + 36000.76983 * T + 0.0003032 * T2; |
| 28 | + L0 = ((L0 % 360) + 360) % 360; |
| 29 | + |
| 30 | + // Mean anomaly of the Sun (degrees) — Meeus eq. 25.3 |
| 31 | + let M = 357.52911 + 35999.05029 * T - 0.0001537 * T2; |
| 32 | + M = ((M % 360) + 360) % 360; |
| 33 | + const Mrad = M * DEG2RAD; |
| 34 | + |
| 35 | + // Eccentricity of Earth's orbit — Meeus eq. 25.4 |
| 36 | + const e = 0.016708634 - 0.000042037 * T - 0.0000001267 * T2; |
| 37 | + |
| 38 | + // Sun's equation of center (degrees) — Meeus p. 164 |
| 39 | + const C = (1.9146 - 0.004817 * T - 0.000014 * T2) * Math.sin(Mrad) |
| 40 | + + (0.019993 - 0.000101 * T) * Math.sin(2 * Mrad) |
| 41 | + + 0.000289 * Math.sin(3 * Mrad); |
| 42 | + |
| 43 | + // Sun's true longitude and true anomaly (degrees) |
| 44 | + const theta = L0 + C; |
| 45 | + const nu = M + C; |
| 46 | + const nuRad = nu * DEG2RAD; |
| 47 | + |
| 48 | + // Sun's radius vector in AU — Meeus eq. 25.5 |
| 49 | + const rAU = (1.000001018 * (1 - e * e)) / (1 + e * Math.cos(nuRad)); |
| 50 | + |
| 51 | + // Longitude of the ascending node of the Moon's mean orbit (degrees) |
| 52 | + // Used for nutation correction — Meeus eq. 25.8 (low-precision nutation) |
| 53 | + const omega = 125.04 - 1934.136 * T; |
| 54 | + const omegaRad = omega * DEG2RAD; |
| 55 | + |
| 56 | + // Apparent longitude (degrees) — corrected for nutation and aberration |
| 57 | + // Meeus eq. 25.8: -0.00569° for aberration, -0.00478°·sin(Ω) for nutation |
| 58 | + const lambdaApp = (theta - 0.00569 - 0.00478 * Math.sin(omegaRad)) * DEG2RAD; |
| 59 | + |
| 60 | + // Mean obliquity of the ecliptic — Meeus eq. 22.2 (in degrees) |
| 61 | + // 23°26'21.448" = 23.4392911° |
| 62 | + const eps0 = 23.4392911 - (46.8150 / 3600) * T - (0.00059 / 3600) * T2 + (0.001813 / 3600) * T3; |
| 63 | + |
| 64 | + // Apparent obliquity (corrected for nutation) — Meeus eq. 25.8 |
| 65 | + const epsilonApp = (eps0 + 0.00256 * Math.cos(omegaRad)) * DEG2RAD; |
| 66 | + |
| 67 | + return { lambdaApp, epsilonApp, rAU }; |
| 68 | +} |
| 69 | + |
| 70 | +/** |
| 71 | + * Compute unit sun direction in **standard ECI** (Earth-Centered Inertial) coordinates. |
| 72 | + * NOT render coords — use calculateSunPosition() from sun.ts for render coords. |
| 73 | + * No THREE.js dependency — safe for Web Workers. |
| 74 | + */ |
| 75 | +export function sunDirectionECI(epoch: number): { x: number; y: number; z: number } { |
| 76 | + const { lambdaApp, epsilonApp } = computeSunEcliptic(epoch); |
| 77 | + |
| 78 | + // Ecliptic to equatorial (ECI). Sun's ecliptic latitude β ≈ 0. |
| 79 | + const cosLam = Math.cos(lambdaApp); |
| 80 | + const sinLam = Math.sin(lambdaApp); |
| 81 | + const cosEps = Math.cos(epsilonApp); |
| 82 | + const sinEps = Math.sin(epsilonApp); |
| 83 | + |
| 84 | + const x = cosLam; |
| 85 | + const y = sinLam * cosEps; |
| 86 | + const z = sinLam * sinEps; |
| 87 | + |
| 88 | + const len = Math.sqrt(x * x + y * y + z * z); |
| 89 | + return { x: x / len, y: y / len, z: z / len }; |
| 90 | +} |
0 commit comments