refactor: ♻️ refactor common celestial equations out of sun and moon

thkruz · Dec 31, 2023 · 7f65165 · 7f65165
1 parent 7f4146e
commit 7f65165
Show file tree

Hide file tree

Showing 5 changed files with 106 additions and 104 deletions.
diff --git a/src/body/Celestial.ts b/src/body/Celestial.ts
@@ -0,0 +1,78 @@
+import { AzEl, Degrees, RaDec, Radians } from '@src/ootk';
+import { Sun } from './Sun';
+
+export class Celestial {
+  private constructor() {
+    // disable constructor
+  }
+
+  static getStarAzEl(date: Date, lat: Degrees, lon: Degrees, ra: Radians, dec: Radians): AzEl<Radians> {
+    const c: RaDec = {
+      ra,
+      dec,
+    };
+    const azEl = Sun.azEl(date, lat, lon, c);
+
+    const el = <Radians>(azEl.el + Celestial.astroRefraction(azEl.el)); // elevation correction for refraction
+
+    return {
+      az: azEl.az,
+      el,
+    };
+  }
+
+  /**
+   * get astro refraction
+   * @param {Radians} h - elevation
+   * @returns {number} refraction
+   */
+  static astroRefraction(h: Radians): Radians {
+    if (h < 0) {
+      h = <Radians>0;
+    }
+
+    return <Radians>(0.0002967 / Math.tan(h + 0.00312536 / (h + 0.08901179)));
+  }
+
+  /**
+   * get declination
+   * @param {number} l - ecliptic longitude
+   * @param {number} b - ecliptic latitude
+   * @returns {number} declination
+   */
+  static declination(l: number, b: number): Radians {
+    return <Radians>Math.asin(Math.sin(b) * Math.cos(Sun.e) + Math.cos(b) * Math.sin(Sun.e) * Math.sin(l));
+  }
+
+  /**
+   * get right ascension
+   * @param {number} l - ecliptic longitude
+   * @param {number} b - ecliptic latitude
+   * @returns {number} right ascension
+   */
+  static rightAscension(l: number, b: number): Radians {
+    return <Radians>Math.atan2(Math.sin(l) * Math.cos(Sun.e) - Math.tan(b) * Math.sin(Sun.e), Math.cos(l));
+  }
+
+  /**
+   * get elevation (sometimes called altitude)
+   * @param {number} H - siderealTime
+   * @param {Radians} phi - latitude
+   * @param {Radians} dec - The declination of the sun
+   * @returns {Radians} elevation
+   */
+  static elevation(H: number, phi: Radians, dec: Radians): Radians {
+    return <Radians>Math.asin(Math.sin(phi) * Math.sin(dec) + Math.cos(phi) * Math.cos(dec) * Math.cos(H));
+  }
+
+  /**
+   * get azimuth
+   * @param {number} H - siderealTime
+   * @param {Radians} phi - latitude
+   * @param {Radians} dec - The declination of the sun
+   * @returns {Radians} azimuth in rad
+   */
+  static azimuth(H: number, phi: Radians, dec: Radians): Radians {
+    return <Radians>(Math.PI + Math.atan2(Math.sin(H), Math.cos(H) * Math.sin(phi) - Math.tan(dec) * Math.cos(phi)));
+  }
+}
diff --git a/src/body/Star.ts b/src/body/Star.ts
diff --git a/src/body/Sun.ts b/src/body/Sun.ts
@@ -5,6 +5,7 @@ import { angularDiameter, AngularDiameterMethod } from '../operations/functions'
 import { Vector3D } from '../operations/Vector3D';
 import { EpochUTC } from '../time/EpochUTC';
 import { DEG2RAD, MS_PER_DAY } from '../utils/constants';
+import { Celestial } from './Celestial';
 import { Earth } from './Earth';
 
 /**
@@ -14,7 +15,7 @@ export class Sun {
   private static readonly J0_ = 0.0009;
   private static readonly J1970_ = 2440587.5;
   private static readonly J2000_ = 2451545;
-  private static readonly e_ = DEG2RAD * 23.4397;
+  static readonly e = DEG2RAD * 23.4397;
 
   /**
    * Array representing the times for different phases of the sun.
@@ -83,8 +84,8 @@ export class Sun {
     const H = Sun.siderealTime(d, lw) - c.ra;
 
     return {
-      az: Sun.azimuth_(H, phi, c.dec),
-      el: Sun.elevation_(H, phi, c.dec),
+      az: Celestial.azimuth(H, phi, c.dec),
+      el: Celestial.elevation(H, phi, c.dec),
     };
   }
 
@@ -241,7 +242,7 @@ export class Sun {
       const d = Sun.date2jSince2000(newDate);
       const c = Sun.raDec(newDate);
       const H = Sun.siderealTime(d, lw) - c.ra;
-      const newAz = Sun.azimuth_(H, phi, c.dec);
+      const newAz = Celestial.azimuth(H, phi, c.dec);
 
       addval /= 2;
       if (newAz < az) {
@@ -432,8 +433,8 @@ export class Sun {
     const L = Sun.eclipticLongitude(M);
 
     return {
-      dec: Sun.declination_(L, 0),
-      ra: Sun.rightAscension_(L, 0),
+      dec: Celestial.declination(L, 0),
+      ra: Celestial.rightAscension(L, 0),
     };
   }
 
@@ -491,17 +492,6 @@ export class Sun {
     return Sun.J0_ + (Ht + lw) / tau + n;
   }
 
-  /**
-   * get azimuth
-   * @param {number} H - siderealTime
-   * @param {Radians} phi - latitude
-   * @param {Radians} dec - The declination of the sun
-   * @returns {Radians} azimuth in rad
-   */
-  private static azimuth_(H: number, phi: Radians, dec: Radians): Radians {
-    return <Radians>(Math.PI + Math.atan2(Math.sin(H), Math.cos(H) * Math.sin(phi) - Math.tan(dec) * Math.cos(phi)));
-  }
-
   private static calculateJnoon_(lon: Degrees, lat: Degrees, alt: Meters, date: Date) {
     const lw = <Radians>(DEG2RAD * -lon);
     const phi = <Radians>(DEG2RAD * lat);
@@ -512,33 +502,12 @@ export class Sun {
     const ds = Sun.approxTransit_(0, lw, n);
     const M = Sun.solarMeanAnomaly_(ds);
     const L = Sun.eclipticLongitude(M);
-    const dec = Sun.declination_(L, 0);
+    const dec = Celestial.declination(L, 0);
     const Jnoon = Sun.solarTransitJulian(ds, M, L);
 
     return { Jnoon, dh, lw, phi, dec, n, M, L };
   }
 
-  /**
-   * get declination
-   * @param {number} l - ecliptic longitude
-   * @param {number} b - ecliptic latitude
-   * @returns {number} declination
-   */
-  private static declination_(l: number, b: number): Radians {
-    return <Radians>Math.asin(Math.sin(b) * Math.cos(Sun.e_) + Math.cos(b) * Math.sin(Sun.e_) * Math.sin(l));
-  }
-
-  /**
-   * get elevation (sometimes called altitude)
-   * @param {number} H - siderealTime
-   * @param {Radians} phi - latitude
-   * @param {Radians} dec - The declination of the sun
-   * @returns {Radians} elevation
-   */
-  private static elevation_(H: number, phi: Radians, dec: Radians): Radians {
-    return <Radians>Math.asin(Math.sin(phi) * Math.sin(dec) + Math.cos(phi) * Math.cos(dec) * Math.cos(H));
-  }
-
   /**
    * returns set time for the given sun altitude
    * @param {Meters} alt - altitude at 0
@@ -580,16 +549,6 @@ export class Sun {
     return <Degrees>((-2.076 * Math.sqrt(alt)) / 60);
   }
 
-  /**
-   * get right ascension
-   * @param {number} l - ecliptic longitude
-   * @param {number} b - ecliptic latitude
-   * @returns {number} right ascension
-   */
-  private static rightAscension_(l: number, b: number): Radians {
-    return <Radians>Math.atan2(Math.sin(l) * Math.cos(Sun.e_) - Math.tan(b) * Math.sin(Sun.e_), Math.cos(l));
-  }
-
   /**
    * get solar mean anomaly
    * @param {number} d - julian day

diff --git a/src/types/types.ts b/src/types/types.ts
@@ -589,6 +589,7 @@ export type AzEl<Units = Radians> = {
 export type RaDec = {
   dec: Radians;
   ra: Radians;
+  dist?: Kilometers;
 };
 
 export interface RadarSensor extends Sensor {

diff --git a/src/utils/moon-math.ts b/src/utils/moon-math.ts
@@ -3,7 +3,7 @@
  * @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 SunMath for calculating the position of the moon
+ * @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.
@@ -34,10 +34,11 @@
  * Orbital Object ToolKit. If not, see <http://www.gnu.org/licenses/>.
  */
 
-import { Degrees, Radians } from '../ootk';
+import { Celestial } from '@src/body/Celestial';
+import { Sun } from '@src/body/Sun';
+import { Degrees, Kilometers, RaDec, Radians } from '../ootk';
 import { MS_PER_DAY } from '../utils/constants';
 import { DEG2RAD } from './constants';
-import { SunMath } from './sun-math';
 
 type MoonIlluminationData = {
   fraction: number;
@@ -183,8 +184,9 @@ export class MoonMath {
 
     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 d = SunMath.date2jSince2000(new Date(dateValue));
-    const s = SunMath.getSunRaDec(d);
+    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(
@@ -231,11 +233,9 @@ export class MoonMath {
     // eslint-disable-next-line init-declarations
     let phase;
 
-    for (let index = 0; index < MoonMath.moonCycles_.length; index++) {
-      const element = MoonMath.moonCycles_[index];
-
-      if (phaseValue >= element.from && phaseValue <= element.to) {
-        phase = element;
+    for (const moonCycle of MoonMath.moonCycles_) {
+      if (phaseValue >= moonCycle.from && phaseValue <= moonCycle.to) {
+        phase = moonCycle;
         break;
       }
     }
@@ -284,17 +284,17 @@ export class MoonMath {
   static getMoonPosition(date: Date, lat: Degrees, lon: Degrees) {
     const lw = <Radians>(DEG2RAD * -lon);
     const phi = <Radians>(DEG2RAD * lat);
-    const d = SunMath.date2jSince2000(date);
+    const d = Sun.date2jSince2000(date);
     const c = MoonMath.moonCoords(d);
-    const H = SunMath.siderealTime(d, lw) - c.ra;
-    let h = SunMath.elevation(H, phi, c.dec);
+    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 = <Radians>(h + SunMath.astroRefraction(h)); // altitude correction for refraction
+    h = <Radians>(h + Celestial.astroRefraction(h)); // altitude correction for refraction
 
     return {
-      az: SunMath.azimuth(H, phi, c.dec),
+      az: Celestial.azimuth(H, phi, c.dec),
       el: h,
       rng: c.dist,
       parallacticAngle: pa,
@@ -365,7 +365,7 @@ export class MoonMath {
     return new Date(date.getTime() + (h * MS_PER_DAY) / 24);
   }
 
-  static moonCoords(d) {
+  static moonCoords(d): RaDec {
     // geocentric ecliptic coordinates of the moon
 
     const L = DEG2RAD * (218.316 + 13.176396 * d); // ecliptic longitude
@@ -376,9 +376,9 @@ export class MoonMath {
     const dt = 385001 - 20905 * Math.cos(M); // distance to the moon in km
 
     return {
-      ra: SunMath.rightAscension(l, b),
-      dec: SunMath.declination(l, b),
-      dist: dt,
+      ra: Celestial.rightAscension(l, b),
+      dec: Celestial.declination(l, b),
+      dist: dt as Kilometers,
     };
   }