diff --git a/src/body/Moon.ts b/src/body/Moon.ts index fed7684..7823fa9 100644 --- a/src/body/Moon.ts +++ b/src/body/Moon.ts @@ -1,11 +1,88 @@ -import { AngularDiameterMethod } from '../ootk-core'; +/** + * @author Theodore Kruczek. + * @description Orbital Object ToolKit (OOTK) is a collection of tools for working + * with satellites and other orbital objects. + * + * Some of the math in this file was originally created by Vladimir Agafonkin. + * Robert Gester's update was referenced for documentation. There were a couple of + * bugs in both versions so there will be some differences if you are migrating from + * either to this library. + * + * @Copyright (c) 2011-2015, Vladimir Agafonkin + * SunCalc is a JavaScript library for calculating sun/moon position and light phases. + * https://github.com/mourner/suncalc + * + * Reworked and enhanced by Robert Gester + * @Copyright (c) 2022 Robert Gester + * https://github.com/hypnos3/suncalc3 + * + * moon calculations, based on http://aa.quae.nl/en/reken/hemelpositie.html formulas + * + * @license MIT License + * + * @Copyright (c) 2022-2024 Theodore Kruczek + * Permission is hereby granted, free of charge, to any person obtaining a copy of this + * software and associated documentation files (the "Software"), to deal in the Software + * without restriction, including without limitation the rights to use, copy, modify, merge, + * publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons + * to whom the Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all copies or + * substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, + * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR + * PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE + * FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR + * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + */ + +import { AngularDiameterMethod, Celestial, Degrees, Kilometers, RaDec, Radians } from '../ootk-core'; import { Vector3D } from '../operations/Vector3D'; import { EpochUTC } from '../time/EpochUTC'; -import { DEG2RAD } from '../utils/constants'; +import { DEG2RAD, MS_PER_DAY } from '../utils/constants'; import { angularDiameter } from '../utils/functions'; import { Earth } from './Earth'; import { Sun } from './Sun'; +type MoonIlluminationData = { + fraction: number; + phase: { + from: number; + to: number; + id: string; + emoji: string; + code: string; + name: string; + weight: number; + css: string; + }; + phaseValue: number; + angle: number; + next: { + value: number; + date: string; + type: string; + newMoon: { + value: number; + date: string; + }; + fullMoon: { + value: number; + date: string; + }; + firstQuarter: { + value: number; + date: string; + }; + thirdQuarter: { + value: number; + date: string; + }; + }; +}; + // / Moon metrics and operations. export class Moon { private constructor() { @@ -19,7 +96,7 @@ export class Moon { static radiusEquator = 1738.0; // / Calculate the Moon's ECI position _(km)_ for a given UTC [epoch]. - static position(epoch: EpochUTC): Vector3D { + static eci(epoch: EpochUTC): Vector3D { const jc = epoch.toJulianCenturies(); const dtr = DEG2RAD; const lamEcl = @@ -69,7 +146,7 @@ export class Moon { static illumination(epoch: EpochUTC, origin?: Vector3D): number { const orig = origin ?? Vector3D.origin; const sunPos = Sun.position(epoch).subtract(orig); - const moonPos = this.position(epoch).subtract(orig); + const moonPos = this.eci(epoch).subtract(orig); const phaseAngle = sunPos.angle(moonPos); return 0.5 * (1 - Math.cos(phaseAngle)); @@ -85,4 +162,379 @@ export class Moon { static diameter(obsPos: Vector3D, moonPos: Vector3D): number { return angularDiameter(this.radiusEquator * 2, obsPos.subtract(moonPos).magnitude(), AngularDiameterMethod.Sphere); } + + /** + * calculations for illumination parameters of the moon, + * based on http://idlastro.gsfc.nasa.gov/ftp/pro/astro/mphase.pro formulas and + * Chapter 48 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. + * @param {number | Date} date Date object or timestamp for calculating moon-illumination + * @return {MoonIlluminationData} result object of moon-illumination + */ + // eslint-disable-next-line max-statements + static getMoonIllumination(date: number | Date): MoonIlluminationData { + const dateValue = date instanceof Date ? date.getTime() : date; + + const lunarDaysMs = 2551442778; // The duration in days of a lunar cycle is 29.53058770576 days. + const firstNewMoon2000 = 947178840000; // first newMoon in the year 2000 2000-01-06 18:14 + const dateObj = new Date(dateValue); + const d = Sun.date2jSince2000(dateObj); + const s = Sun.raDec(dateObj); + const m = Moon.moonCoords(d); + const sdist = 149598000; // distance from Earth to Sun in km + const phi = Math.acos( + Math.sin(s.dec) * Math.sin(m.dec) + Math.cos(s.dec) * Math.cos(m.dec) * Math.cos(s.ra - m.ra), + ); + const inc = Math.atan2(sdist * Math.sin(phi), m.dist - sdist * Math.cos(phi)); + const angle = Math.atan2( + Math.cos(s.dec) * Math.sin(s.ra - m.ra), + Math.sin(s.dec) * Math.cos(m.dec) - Math.cos(s.dec) * Math.sin(m.dec) * Math.cos(s.ra - m.ra), + ); + const phaseValue = 0.5 + (0.5 * inc * (angle < 0 ? -1 : 1)) / Math.PI; + + // calculates the difference in ms between the sirst fullMoon 2000 and given Date + const diffBase = dateValue - firstNewMoon2000; + // Calculate modulus to drop completed cycles + let cycleModMs = diffBase % lunarDaysMs; + // If negative number (date before new moon 2000) add lunarDaysMs + + if (cycleModMs < 0) { + cycleModMs += lunarDaysMs; + } + const nextNewMoon = lunarDaysMs - cycleModMs + dateValue; + let nextFullMoon = lunarDaysMs / 2 - cycleModMs + dateValue; + + if (nextFullMoon < dateValue) { + nextFullMoon += lunarDaysMs; + } + const quater = lunarDaysMs / 4; + let nextFirstQuarter = quater - cycleModMs + dateValue; + + if (nextFirstQuarter < dateValue) { + nextFirstQuarter += lunarDaysMs; + } + let nextThirdQuarter = lunarDaysMs - quater - cycleModMs + dateValue; + + if (nextThirdQuarter < dateValue) { + nextThirdQuarter += lunarDaysMs; + } + /* + * Calculate the fraction of the moon cycle + * const currentfrac = cycleModMs / lunarDaysMs; + */ + const next = Math.min(nextNewMoon, nextFirstQuarter, nextFullMoon, nextThirdQuarter); + // eslint-disable-next-line init-declarations + let phase; + + for (const moonCycle of Moon.moonCycles_) { + if (phaseValue >= moonCycle.from && phaseValue <= moonCycle.to) { + phase = moonCycle; + break; + } + } + + let type = ''; + + if (next === nextNewMoon) { + type = 'newMoon'; + } else if (next === nextFirstQuarter) { + type = 'firstQuarter'; + } else if (next === nextFullMoon) { + type = 'fullMoon'; + } else { + type = 'thirdQuarter'; + } + + return { + fraction: (1 + Math.cos(inc)) / 2, + phase, + phaseValue, + angle, + next: { + value: next, + date: new Date(next).toISOString(), + type, + newMoon: { + value: nextNewMoon, + date: new Date(nextNewMoon).toISOString(), + }, + fullMoon: { + value: nextFullMoon, + date: new Date(nextFullMoon).toISOString(), + }, + firstQuarter: { + value: nextFirstQuarter, + date: new Date(nextFirstQuarter).toISOString(), + }, + thirdQuarter: { + value: nextThirdQuarter, + date: new Date(nextThirdQuarter).toISOString(), + }, + }, + }; + } + + static rae( + date: Date, + lat: Degrees, + lon: Degrees, + ): { + az: Radians; + el: Radians; + rng: Kilometers; + parallacticAngle: Radians; + } { + const lw = (DEG2RAD * -lon); + const phi = (DEG2RAD * lat); + const d = Sun.date2jSince2000(date); + const c = Moon.moonCoords(d); + const H = Sun.siderealTime(d, lw) - c.ra; + let h = Celestial.elevation(H, phi, c.dec); + // formula 14.1 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. + const pa = Math.atan2(Math.sin(H), Math.tan(phi) * Math.cos(c.dec) - Math.sin(c.dec) * Math.cos(H)); + + h = (h + Celestial.astroRefraction(h)); // altitude correction for refraction + + return { + az: Celestial.azimuth(H, phi, c.dec), + el: h, + rng: c.dist, + parallacticAngle: pa as Radians, + }; + } + + /** + * calculations for moon rise/set times are based on http://www.stargazing.net/kepler/moonrise.html article + */ + static getMoonTimes(date: Date, lat: Degrees, lon: Degrees, isUtc = false) { + if (isUtc) { + date.setUTCHours(0, 0, 0, 0); + } else { + date.setHours(0, 0, 0, 0); + } + + const { rise, set, ye } = Moon.calculateRiseSetTimes_(date, lat, lon); + + const result = { + rise: null, + set: null, + ye: null, + alwaysUp: null, + alwaysDown: null, + highest: null, + }; + + if (rise) { + result.rise = new Date(Moon.hoursLater_(date, rise)); + } else { + result.rise = NaN; + } + + if (set) { + result.set = new Date(Moon.hoursLater_(date, set)); + } else { + result.set = NaN; + } + + if (!rise && !set) { + if (ye > 0) { + result.alwaysUp = true; + result.alwaysDown = false; + } else { + result.alwaysUp = false; + result.alwaysDown = true; + } + } else if (rise && set) { + result.alwaysUp = false; + result.alwaysDown = false; + result.highest = new Date(Moon.hoursLater_(date, Math.min(rise, set) + Math.abs(set - rise) / 2)); + } else { + result.alwaysUp = false; + result.alwaysDown = false; + } + + return result; + } + + private static hoursLater_(date: Date, h: number) { + return new Date(date.getTime() + (h * MS_PER_DAY) / 24); + } + + /** + * Calculates the geocentric ecliptic coordinates of the moon. + * + * @param d - The number of days since year 2000. + * @returns An object containing the right ascension, declination, and distance to the moon. + */ + static moonCoords(d: number): RaDec { + const L = DEG2RAD * (218.316 + 13.176396 * d); // ecliptic longitude + const M = DEG2RAD * (134.963 + 13.064993 * d); // mean anomaly + const F = DEG2RAD * (93.272 + 13.22935 * d); // mean distance + const l = L + DEG2RAD * 6.289 * Math.sin(M); // longitude + const b = DEG2RAD * 5.128 * Math.sin(F); // latitude + const dt = 385001 - 20905 * Math.cos(M); // distance to the moon in km + + return { + ra: Celestial.rightAscension(l, b), + dec: Celestial.declination(l, b), + dist: dt as Kilometers, + }; + } + + private static calculateRiseSetTimes_(t: Date, lat: Degrees, lon: Degrees) { + const hc = 0.133 * DEG2RAD; + let h0 = Moon.rae(t, lat, lon).el - hc; + let h1 = 0; + let h2 = 0; + let rise = 0; + let set = 0; + let a = 0; + let b = 0; + let xe = 0; + let ye = 0; + let d = 0; + let roots = 0; + let x1 = 0; + let x2 = 0; + let dx = 0; + + // go in 2-hour chunks, each time seeing if a 3-point quadratic curve crosses zero (which means rise or set) + for (let i = 1; i <= 24; i += 2) { + h1 = Moon.rae(Moon.hoursLater_(t, i), lat, lon).el - hc; + h2 = Moon.rae(Moon.hoursLater_(t, i + 1), lat, lon).el - hc; + + a = (h0 + h2) / 2 - h1; + b = (h2 - h0) / 2; + xe = -b / (2 * a); + ye = (a * xe + b) * xe + h1; + d = b * b - 4 * a * h1; + roots = 0; + + if (d >= 0) { + dx = Math.sqrt(d) / (Math.abs(a) * 2); + x1 = xe - dx; + x2 = xe + dx; + if (Math.abs(x1) <= 1) { + roots++; + } + if (Math.abs(x2) <= 1) { + roots++; + } + if (x1 < -1) { + x1 = x2; + } + } + + if (roots === 1) { + if (h0 < 0) { + rise = i + x1; + } else { + set = i + x1; + } + } else if (roots === 2) { + rise = i + (ye < 0 ? x2 : x1); + set = i + (ye < 0 ? x1 : x2); + } + + if (rise && set) { + break; + } + + h0 = h2; + } + + return { rise, set, ye }; + } + + private static moonCycles_ = [ + { + from: 0, + to: 0.033863193308711, + id: 'newMoon', + emoji: '🌚', + code: ':new_moon_with_face:', + name: 'New Moon', + weight: 1, + css: 'wi-moon-new', + }, + { + from: 0.033863193308711, + to: 0.216136806691289, + id: 'waxingCrescentMoon', + emoji: '🌒', + code: ':waxing_crescent_moon:', + name: 'Waxing Crescent', + weight: 6.3825, + css: 'wi-moon-wax-cres', + }, + { + from: 0.216136806691289, + to: 0.283863193308711, + id: 'firstQuarterMoon', + emoji: '🌓', + code: ':first_quarter_moon:', + name: 'First Quarter', + weight: 1, + css: 'wi-moon-first-quart', + }, + { + from: 0.283863193308711, + to: 0.466136806691289, + id: 'waxingGibbousMoon', + emoji: '🌔', + code: ':waxing_gibbous_moon:', + name: 'Waxing Gibbous', + weight: 6.3825, + css: 'wi-moon-wax-gibb', + }, + { + from: 0.466136806691289, + to: 0.533863193308711, + id: 'fullMoon', + emoji: '🌝', + code: ':full_moon_with_face:', + name: 'Full Moon', + weight: 1, + css: 'wi-moon-full', + }, + { + from: 0.533863193308711, + to: 0.716136806691289, + id: 'waningGibbousMoon', + emoji: '🌖', + code: ':waning_gibbous_moon:', + name: 'Waning Gibbous', + weight: 6.3825, + css: 'wi-moon-wan-gibb', + }, + { + from: 0.716136806691289, + to: 0.783863193308711, + id: 'thirdQuarterMoon', + emoji: '🌗', + code: ':last_quarter_moon:', + name: 'third Quarter', + weight: 1, + css: 'wi-moon-third-quart', + }, + { + from: 0.783863193308711, + to: 0.966136806691289, + id: 'waningCrescentMoon', + emoji: '🌘', + code: ':waning_crescent_moon:', + name: 'Waning Crescent', + weight: 6.3825, + css: 'wi-moon-wan-cres', + }, + { + from: 0.966136806691289, + to: 1, + id: 'newMoon', + emoji: '🌚', + code: ':new_moon_with_face:', + name: 'New Moon', + weight: 1, + css: 'wi-moon-new', + }, + ]; } diff --git a/src/utils/functions.ts b/src/utils/functions.ts index 7ee3d81..0e67c47 100644 --- a/src/utils/functions.ts +++ b/src/utils/functions.ts @@ -6,9 +6,11 @@ import { AngularDiameterMethod, AngularDistanceMethod, DifferentiableFunction, + EciVec3, JacobianFunction, Radians, SpaceObjectType, + Vec3, } from '../types/types'; /** @@ -301,6 +303,16 @@ export function angularDistance( } } +/** + * Calculates the linear distance between two points in three-dimensional space. + * @param pos1 The first position. + * @param pos2 The second position. + * @returns The linear distance between the two positions in kilometers. + */ +export function linearDistance(pos1: Vec3, pos2: Vec3): D { + return Math.sqrt((pos1.x - pos2.x) ** 2 + (pos1.y - pos2.y) ** 2 + (pos1.z - pos2.z) ** 2); +} + /** * Calculates the angular diameter of an object. * @@ -654,3 +666,34 @@ const spaceObjTypeStrMap = { export const spaceObjType2Str = (spaceObjType: SpaceObjectType): string => spaceObjTypeStrMap[spaceObjType] || 'Unknown'; + +export const dopplerFactor = (location: EciVec3, position: EciVec3, velocity: EciVec3): number => { + const mfactor = 7.292115e-5; + const c = 299792.458; // Speed of light in km/s + + const range = { + x: position.x - location.x, + y: position.y - location.y, + z: position.z - location.z, + }; + const distance = Math.sqrt(range.x ** 2 + range.y ** 2 + range.z ** 2); + + const rangeVel = { + x: velocity.x + mfactor * location.y, + y: velocity.y - mfactor * location.x, + z: velocity.z, + }; + + const rangeRate = (range.x * rangeVel.x + range.y * rangeVel.y + range.z * rangeVel.z) / distance; + let dopplerFactor = 0; + + if (rangeRate < 0) { + dopplerFactor = 1 + (rangeRate / c) * sign(rangeRate); + } + + if (rangeRate >= 0) { + dopplerFactor = 1 - (rangeRate / c) * sign(rangeRate); + } + + return dopplerFactor; +}; diff --git a/src/utils/moon-math.ts b/src/utils/moon-math.ts deleted file mode 100644 index 34f8e29..0000000 --- a/src/utils/moon-math.ts +++ /dev/null @@ -1,454 +0,0 @@ -/** - * @author Theodore Kruczek. - * @description Orbital Object ToolKit (OOTK) is a collection of tools for working - * with satellites and other orbital objects. - * - * @file MoonMath is a an extension to Sun for calculating the position of the moon - * and its phases. This was originally created by Vladimir Agafonkin. Robert Gester's - * update was referenced for documentation. There were a couple of bugs in both versions - * so there will be some differences if you are migrating from either to this library. - * - * @Copyright (c) 2011-2015, Vladimir Agafonkin - * SunCalc is a JavaScript library for calculating sun/moon position and light phases. - * https://github.com/mourner/suncalc - * - * Reworked and enhanced by Robert Gester - * @Copyright (c) 2022 Robert Gester - * https://github.com/hypnos3/suncalc3 - * - * moon calculations, based on http://aa.quae.nl/en/reken/hemelpositie.html formulas - * - * @license MIT License - * - * @Copyright (c) 2024 Theodore Kruczek - * Permission is hereby granted, free of charge, to any person obtaining a copy of this - * software and associated documentation files (the "Software"), to deal in the Software - * without restriction, including without limitation the rights to use, copy, modify, merge, - * publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons - * to whom the Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included in all copies or - * substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, - * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR - * PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE - * FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR - * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER - * DEALINGS IN THE SOFTWARE. - */ - -import { Celestial } from '../body/Celestial'; -import { Sun } from '../body/Sun'; -import { Degrees, Kilometers, RaDec, Radians } from '../ootk-core'; -import { MS_PER_DAY } from '../utils/constants'; -import { DEG2RAD } from './constants'; - -type MoonIlluminationData = { - fraction: number; - phase: { - from: number; - to: number; - id: string; - emoji: string; - code: string; - name: string; - weight: number; - css: string; - }; - phaseValue: number; - angle: number; - next: { - value: number; - date: string; - type: string; - newMoon: { - value: number; - date: string; - }; - fullMoon: { - value: number; - date: string; - }; - firstQuarter: { - value: number; - date: string; - }; - thirdQuarter: { - value: number; - date: string; - }; - }; -}; - -export class MoonMath { - private static moonCycles_ = [ - { - from: 0, - to: 0.033863193308711, - id: 'newMoon', - emoji: '🌚', - code: ':new_moon_with_face:', - name: 'New Moon', - weight: 1, - css: 'wi-moon-new', - }, - { - from: 0.033863193308711, - to: 0.216136806691289, - id: 'waxingCrescentMoon', - emoji: '🌒', - code: ':waxing_crescent_moon:', - name: 'Waxing Crescent', - weight: 6.3825, - css: 'wi-moon-wax-cres', - }, - { - from: 0.216136806691289, - to: 0.283863193308711, - id: 'firstQuarterMoon', - emoji: '🌓', - code: ':first_quarter_moon:', - name: 'First Quarter', - weight: 1, - css: 'wi-moon-first-quart', - }, - { - from: 0.283863193308711, - to: 0.466136806691289, - id: 'waxingGibbousMoon', - emoji: '🌔', - code: ':waxing_gibbous_moon:', - name: 'Waxing Gibbous', - weight: 6.3825, - css: 'wi-moon-wax-gibb', - }, - { - from: 0.466136806691289, - to: 0.533863193308711, - id: 'fullMoon', - emoji: '🌝', - code: ':full_moon_with_face:', - name: 'Full Moon', - weight: 1, - css: 'wi-moon-full', - }, - { - from: 0.533863193308711, - to: 0.716136806691289, - id: 'waningGibbousMoon', - emoji: '🌖', - code: ':waning_gibbous_moon:', - name: 'Waning Gibbous', - weight: 6.3825, - css: 'wi-moon-wan-gibb', - }, - { - from: 0.716136806691289, - to: 0.783863193308711, - id: 'thirdQuarterMoon', - emoji: '🌗', - code: ':last_quarter_moon:', - name: 'third Quarter', - weight: 1, - css: 'wi-moon-third-quart', - }, - { - from: 0.783863193308711, - to: 0.966136806691289, - id: 'waningCrescentMoon', - emoji: '🌘', - code: ':waning_crescent_moon:', - name: 'Waning Crescent', - weight: 6.3825, - css: 'wi-moon-wan-cres', - }, - { - from: 0.966136806691289, - to: 1, - id: 'newMoon', - emoji: '🌚', - code: ':new_moon_with_face:', - name: 'New Moon', - weight: 1, - css: 'wi-moon-new', - }, - ]; - - /** - * calculations for illumination parameters of the moon, - * based on http://idlastro.gsfc.nasa.gov/ftp/pro/astro/mphase.pro formulas and - * Chapter 48 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. - * @param {number | Date} date Date object or timestamp for calculating moon-illumination - * @return {MoonIlluminationData} result object of moon-illumination - */ - // eslint-disable-next-line max-statements - static getMoonIllumination(date: number | Date): MoonIlluminationData { - const dateValue = date instanceof Date ? date.getTime() : date; - - const lunarDaysMs = 2551442778; // The duration in days of a lunar cycle is 29.53058770576 days. - const firstNewMoon2000 = 947178840000; // first newMoon in the year 2000 2000-01-06 18:14 - const dateObj = new Date(dateValue); - const d = Sun.date2jSince2000(dateObj); - const s = Sun.raDec(dateObj); - const m = MoonMath.moonCoords(d); - const sdist = 149598000; // distance from Earth to Sun in km - const phi = Math.acos( - Math.sin(s.dec) * Math.sin(m.dec) + Math.cos(s.dec) * Math.cos(m.dec) * Math.cos(s.ra - m.ra), - ); - const inc = Math.atan2(sdist * Math.sin(phi), m.dist - sdist * Math.cos(phi)); - const angle = Math.atan2( - Math.cos(s.dec) * Math.sin(s.ra - m.ra), - Math.sin(s.dec) * Math.cos(m.dec) - Math.cos(s.dec) * Math.sin(m.dec) * Math.cos(s.ra - m.ra), - ); - const phaseValue = 0.5 + (0.5 * inc * (angle < 0 ? -1 : 1)) / Math.PI; - - // calculates the difference in ms between the sirst fullMoon 2000 and given Date - const diffBase = dateValue - firstNewMoon2000; - // Calculate modulus to drop completed cycles - let cycleModMs = diffBase % lunarDaysMs; - // If negative number (date before new moon 2000) add lunarDaysMs - - if (cycleModMs < 0) { - cycleModMs += lunarDaysMs; - } - const nextNewMoon = lunarDaysMs - cycleModMs + dateValue; - let nextFullMoon = lunarDaysMs / 2 - cycleModMs + dateValue; - - if (nextFullMoon < dateValue) { - nextFullMoon += lunarDaysMs; - } - const quater = lunarDaysMs / 4; - let nextFirstQuarter = quater - cycleModMs + dateValue; - - if (nextFirstQuarter < dateValue) { - nextFirstQuarter += lunarDaysMs; - } - let nextThirdQuarter = lunarDaysMs - quater - cycleModMs + dateValue; - - if (nextThirdQuarter < dateValue) { - nextThirdQuarter += lunarDaysMs; - } - /* - * Calculate the fraction of the moon cycle - * const currentfrac = cycleModMs / lunarDaysMs; - */ - const next = Math.min(nextNewMoon, nextFirstQuarter, nextFullMoon, nextThirdQuarter); - // eslint-disable-next-line init-declarations - let phase; - - for (const moonCycle of MoonMath.moonCycles_) { - if (phaseValue >= moonCycle.from && phaseValue <= moonCycle.to) { - phase = moonCycle; - break; - } - } - - let type = ''; - - if (next === nextNewMoon) { - type = 'newMoon'; - } else if (next === nextFirstQuarter) { - type = 'firstQuarter'; - } else if (next === nextFullMoon) { - type = 'fullMoon'; - } else { - type = 'thirdQuarter'; - } - - return { - fraction: (1 + Math.cos(inc)) / 2, - phase, - phaseValue, - angle, - next: { - value: next, - date: new Date(next).toISOString(), - type, - newMoon: { - value: nextNewMoon, - date: new Date(nextNewMoon).toISOString(), - }, - fullMoon: { - value: nextFullMoon, - date: new Date(nextFullMoon).toISOString(), - }, - firstQuarter: { - value: nextFirstQuarter, - date: new Date(nextFirstQuarter).toISOString(), - }, - thirdQuarter: { - value: nextThirdQuarter, - date: new Date(nextThirdQuarter).toISOString(), - }, - }, - }; - } - - static getMoonPosition(date: Date, lat: Degrees, lon: Degrees) { - const lw = (DEG2RAD * -lon); - const phi = (DEG2RAD * lat); - const d = Sun.date2jSince2000(date); - const c = MoonMath.moonCoords(d); - const H = Sun.siderealTime(d, lw) - c.ra; - let h = Celestial.elevation(H, phi, c.dec); - // formula 14.1 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998. - const pa = Math.atan2(Math.sin(H), Math.tan(phi) * Math.cos(c.dec) - Math.sin(c.dec) * Math.cos(H)); - - h = (h + Celestial.astroRefraction(h)); // altitude correction for refraction - - return { - az: Celestial.azimuth(H, phi, c.dec), - el: h, - rng: c.dist, - parallacticAngle: pa, - }; - } - - /** - * calculations for moon rise/set times are based on http://www.stargazing.net/kepler/moonrise.html article - * @param {number|Date} dateValue Date object or timestamp for calculating moon-times - * @param {Degrees} lat latitude for calculating moon-times - * @param {Degrees} lon longitude for calculating moon-times - * @param {boolean} [isUtc] defines if the calculation should be in utc or local time (default is local) - * @return {IMoonTimes} result object of sunTime - */ - static getMoonTimes(dateValue: number | Date, lat: Degrees, lon: Degrees, isUtc = false) { - const t = new Date(dateValue); - - if (isUtc) { - t.setUTCHours(0, 0, 0, 0); - } else { - t.setHours(0, 0, 0, 0); - } - - const { rise, set, ye } = MoonMath.calculateRiseSetTimes_(t, lat, lon); - - const result = { - rise: null, - set: null, - ye: null, - alwaysUp: null, - alwaysDown: null, - highest: null, - }; - - if (rise) { - result.rise = new Date(MoonMath.hoursLater(t, rise)); - } else { - result.rise = NaN; - } - - if (set) { - result.set = new Date(MoonMath.hoursLater(t, set)); - } else { - result.set = NaN; - } - - if (!rise && !set) { - if (ye > 0) { - result.alwaysUp = true; - result.alwaysDown = false; - } else { - result.alwaysUp = false; - result.alwaysDown = true; - } - } else if (rise && set) { - result.alwaysUp = false; - result.alwaysDown = false; - result.highest = new Date(MoonMath.hoursLater(t, Math.min(rise, set) + Math.abs(set - rise) / 2)); - } else { - result.alwaysUp = false; - result.alwaysDown = false; - } - - return result; - } - - static hoursLater(date: Date, h: number) { - return new Date(date.getTime() + (h * MS_PER_DAY) / 24); - } - - static moonCoords(d): RaDec { - // geocentric ecliptic coordinates of the moon - - const L = DEG2RAD * (218.316 + 13.176396 * d); // ecliptic longitude - const M = DEG2RAD * (134.963 + 13.064993 * d); // mean anomaly - const F = DEG2RAD * (93.272 + 13.22935 * d); // mean distance - const l = L + DEG2RAD * 6.289 * Math.sin(M); // longitude - const b = DEG2RAD * 5.128 * Math.sin(F); // latitude - const dt = 385001 - 20905 * Math.cos(M); // distance to the moon in km - - return { - ra: Celestial.rightAscension(l, b), - dec: Celestial.declination(l, b), - dist: dt as Kilometers, - }; - } - - // eslint-disable-next-line max-statements - private static calculateRiseSetTimes_(t: Date, lat: Degrees, lon: Degrees) { - const hc = 0.133 * DEG2RAD; - let h0 = MoonMath.getMoonPosition(t, lat, lon).el - hc; - let h1 = 0; - let h2 = 0; - let rise = 0; - let set = 0; - let a = 0; - let b = 0; - let xe = 0; - let ye = 0; - let d = 0; - let roots = 0; - let x1 = 0; - let x2 = 0; - let dx = 0; - - // go in 2-hour chunks, each time seeing if a 3-point quadratic curve crosses zero (which means rise or set) - for (let i = 1; i <= 24; i += 2) { - h1 = MoonMath.getMoonPosition(MoonMath.hoursLater(t, i), lat, lon).el - hc; - h2 = MoonMath.getMoonPosition(MoonMath.hoursLater(t, i + 1), lat, lon).el - hc; - - a = (h0 + h2) / 2 - h1; - b = (h2 - h0) / 2; - xe = -b / (2 * a); - ye = (a * xe + b) * xe + h1; - d = b * b - 4 * a * h1; - roots = 0; - - if (d >= 0) { - dx = Math.sqrt(d) / (Math.abs(a) * 2); - x1 = xe - dx; - x2 = xe + dx; - if (Math.abs(x1) <= 1) { - roots++; - } - if (Math.abs(x2) <= 1) { - roots++; - } - if (x1 < -1) { - x1 = x2; - } - } - - if (roots === 1) { - if (h0 < 0) { - rise = i + x1; - } else { - set = i + x1; - } - } else if (roots === 2) { - rise = i + (ye < 0 ? x2 : x1); - set = i + (ye < 0 ? x1 : x2); - } - - if (rise && set) { - break; - } - - h0 = h2; - } - - return { rise, set, ye }; - } -} diff --git a/src/utils/utils.ts b/src/utils/utils.ts index 658e858..0a3ce30 100644 --- a/src/utils/utils.ts +++ b/src/utils/utils.ts @@ -1,33 +1,8 @@ import * as Types from '../types/types'; -import { EciVec3, Kilometers } from '../types/types'; -import { MILLISECONDS_TO_DAYS } from './constants'; -import { sign } from './functions'; -import { MoonMath } from './moon-math'; class Utils { - static MoonMath = MoonMath; static Types = Types; - // eslint-disable-next-line max-params - static jday = (year?: number, mon?: number, day?: number, hr?: number, minute?: number, sec?: number) => { - if (!year) { - const now = new Date(); - const jDayStart = new Date(now.getUTCFullYear(), 0, 0); - const jDayDiff = now.getDate() - jDayStart.getDate(); - - return Math.floor(jDayDiff / MILLISECONDS_TO_DAYS); - } - - return ( - 367.0 * year - - Math.floor(7 * (year + Math.floor((mon + 9) / 12.0)) * 0.25) + - Math.floor((275 * mon) / 9.0) + - day + - 1721013.5 + - ((sec / 60.0 + minute) / 60.0 + hr) / 24.0 - ); - }; - static createVec(start: number, stop: number, step: number): number[] { const array = []; @@ -37,41 +12,6 @@ class Utils { return array; } - - static distance(pos1: EciVec3, pos2: EciVec3): Kilometers { - return Math.sqrt((pos1.x - pos2.x) ** 2 + (pos1.y - pos2.y) ** 2 + (pos1.z - pos2.z) ** 2); - } - - static dopplerFactor(location: EciVec3, position: EciVec3, velocity: EciVec3): number { - const mfactor = 7.292115e-5; - const c = 299792.458; // Speed of light in km/s - - const range = { - x: position.x - location.x, - y: position.y - location.y, - z: position.z - location.z, - }; - const distance = Math.sqrt(range.x ** 2 + range.y ** 2 + range.z ** 2); - - const rangeVel = { - x: velocity.x + mfactor * location.y, - y: velocity.y - mfactor * location.x, - z: velocity.z, - }; - - const rangeRate = (range.x * rangeVel.x + range.y * rangeVel.y + range.z * rangeVel.z) / distance; - let dopplerFactor = 0; - - if (rangeRate < 0) { - dopplerFactor = 1 + (rangeRate / c) * sign(rangeRate); - } - - if (rangeRate >= 0) { - dopplerFactor = 1 - (rangeRate / c) * sign(rangeRate); - } - - return dopplerFactor; - } } export { Utils };