add conversion from equatorial to horizontal coordinates

README.md

@@ -4,7 +4,21 @@ Typescript library for astronomy calculations.
 Based on formulas from:
 Meeus, Jean: Astronomical Algorithms. Second edition. Richmond, Virginia, 1998.
 
-#### Calculating Greenwich Mean Sidereal Time (GMST)
+#### Content
+
+The main purpose of this library is the transformation of equatorial star coordinates to horizontal star coordinates.
+
+This allows to calculate the position of a star in the sky at a given observation time and place (on Earth).
+
+The library is not meant to be a precise tool for professional purposes, only as a sufficiently accurate approximation for amateur astronomy. Test results were in the range of < 0.5° deviation from https://stellarium.org/. 
+
+#### Usage
+
+For a full usage example see the coordinateConversion test in Gmst.test.ts - the following text describes the steps
+of this test in more detail.
+
+
+##### Calculating Greenwich Mean Sidereal Time (GMST)
 
 Get the julian day number (JDN) at 00:00:00 UT (Universal Time) by using:
 
@@ -26,7 +40,7 @@ The result can contain several days and can be negative. Thus:
 
 The result is the number of seconds that has passed on the current day, these need to be converted to hours/minutes/seconds to get the time of the day (which is the GMST).
 
-#### Calculating Local Mean Sidereal Time
+##### Calculating Local Mean Sidereal Time (LMST)
 
 Get the GMST as above.
 
@@ -36,4 +50,20 @@ Then use the geographical longitude of the current location to determine how man
 
 For example, for the longitude of Berlin (13.41°E) at 16:15:27 the local MST is
 
-``let lmst_berlin = gmst_at_16_15_27 + 13.41 / 360 * 86400 ``
+``let lmst_berlin = gmst_at_16_15_27 + 13.41 / 360 * 86400 ``
+
+##### Calculate Horizontal coordinates of a given star
+
+Get the geographical latitude (in radians):
+
+``let berlinLatitudeRad = 52.52 * 2 * Math.PI / 360;``
+
+Look up the Equatorial coordinates of the star (declination and right ascension), for example in Wikipedia. Then create an EquatorialCoordinate object for them, e.g.
+
+``let polaris = EquatorialCoordinate.fromTimeAndDegree(89, 15, 50.8, 2, 31, 49.09);``
+
+Then use the conversion function to obtain a Horizontal coordinate (altitude and azimuth). The LMST must be converted to radians before.
+
+``let polarisHorizontalCoordinates = HorizontalCoordinate.fromEquatorialCoordinate(polaris, lmst / 86400 * 2 * Math.PI, berlinLatitudeRad);``
+
+The result is in polarisHorizontalCoordinates.altitude and polarisHorizontalCoordinates.azimuth - it is in radians.

src/DateUtils.ts

@@ -15,4 +15,5 @@ export function timeToDayFraction(hours : number, minutes : number, seconds : nu
 {
     let result : number = (hours * 60 * 60 + minutes * 60 + seconds) / (60 * 60 * 24);
     return result;
-}
+}
+

src/StarCoordinates.ts

@@ -39,10 +39,47 @@ export class EquatorialCoordinate
      */
     public static fromTimeAndDegree(declDegree, declMinutes, declSeconds, raHours, raMinutes, raSeconds) : EquatorialCoordinate
     {
-        return {
-            declination : (declDegree + declMinutes / 60.0 + declSeconds / (60.0 * 60.0)) / 360.0 * 2.0 * Math.PI,
-            rightAscension : (raHours * 3600.0 + raMinutes * 60.0 + raSeconds) / (60.0 * 60.0 * 24.0) * 2.0 * Math.PI
-        }
+        return new EquatorialCoordinate(
+            (declDegree + declMinutes / 60.0 + declSeconds / (60.0 * 60.0)) / 360.0 * 2.0 * Math.PI,
+            (raHours * 3600.0 + raMinutes * 60.0 + raSeconds) / (60.0 * 60.0 * 24.0) * 2.0 * Math.PI);
+    }
+
+    /**
+     *
+     * @param declination in radians
+     * @param rightAscension in radians
+     */
+    constructor(declination: number, rightAscension: number) {
+        this.declination = declination;
+        this.rightAscension = rightAscension;
+    }
+
+    /**
+     * Calculate Azimuth of this EquatorialCoordinate. Using the formula from
+     * Meeus, Jean: Astronomical Algorithms. Second edition. Richmond, Virginia, 1998. Page 93.
+     * @param lmst Local Mean Sidereal Time, in radians(!): 86400 seconds == 2 * PI
+     * @param phi observer latitude (on Earth), in radians
+     * @return Azimuth in radians
+     */
+    public calculateAzimuth(lmst : number, phi : number) : number {
+        // H hour angle = theta - alpha = local sidereal time - right ascension
+        let H = lmst - this.rightAscension;
+        let result = Math.atan(Math.sin(H) / (Math.sin(phi) * Math.cos(H) - Math.cos(phi) * Math.tan(this.declination)));
+        if (result <0 ) result += 2 * Math.PI;
+        return result;
+    }
+
+    /**
+     * Calculate Altitude of this EquatorialCoordinate. Using the formula from
+     * Meeus, Jean: Astronomical Algorithms. Second edition. Richmond, Virginia, 1998. Page 93.
+     * @param lmst Local Mean Sidereal Time, in radians(!): 86400 seconds == 2 * PI
+     * @param phi observer latitude (on Earth), in radians
+     * @return Altitude in radians
+     */
+    public calculateAltitude(lmst : number, phi : number) : number {
+        // H hour angle = theta - alpha = local sidereal time - right ascension
+        let H = lmst - this.rightAscension;
+        return Math.asin(Math.sin(phi) * Math.sin(this.declination) + Math.cos(phi) * Math.cos(this.declination) * Math.cos(H));
     }
 }
 
@@ -75,4 +112,29 @@ export class HorizontalCoordinate {
      */
     public azimuth : number;
 
+    /**
+     * Calculate HorizontalCoordinate from EquatorialCoordinate. Using the formula from
+     * Meeus, Jean: Astronomical Algorithms. Second edition. Richmond, Virginia, 1998. Page 93.
+     * @param eq Equatorial coordinate of the star
+     * @param lmst Local Mean Sidereal Time, in radians(!): 86400 seconds == 2 * PI
+     * @param phi observer latitude (on Earth), in radians
+     * @return HorizontalCoordinate
+     */
+    public static fromEquatorialCoordinate(eq : EquatorialCoordinate, lmst : number, phi : number) : HorizontalCoordinate
+    {
+        return new HorizontalCoordinate(
+            eq.calculateAltitude(lmst, phi),
+            eq.calculateAzimuth(lmst, phi)
+        )
+    }
+
+    /**
+     *
+     * @param altitude in radians
+     * @param azimuth in radians
+     */
+    constructor(altitude: number, azimuth: number) {
+        this.altitude = altitude;
+        this.azimuth = azimuth;
+    }
 }

tests/Gmst.test.ts

@@ -1,6 +1,9 @@
 
 import {
-    EquatorialCoordinate, greenwichMeanSiderealTimeToLocalMeanSiderealTime, gregorianDateToJulianDayNumber,
+    EquatorialCoordinate,
+    greenwichMeanSiderealTimeToLocalMeanSiderealTime,
+    gregorianDateToJulianDayNumber,
+    HorizontalCoordinate,
     julianDayNumberToGreenwichMeanSiderealTime,
     timeToDayFraction
 } from "../src";
@@ -61,7 +64,8 @@ describe('starCoordinates', () => {
 
 });
 
-/* data for the next test:
+/* Data for the next test:
+
 Current date (local): 2022-11-06 20:00:29
 Current date (UTC): 2022-11-06 19:00:29
 Julian date (current UTC): 2459890.2920023147
@@ -78,10 +82,20 @@ Polaris azimuth: 0.017127232361280796rad = 0.9813181290412728° =0°58'0.8790877
 Polaris altitude: 0.9242232902032961rad = 52.95409385634355°
 Mizar azimuth: 5.906936254323217rad = 338.44251722553537° =338°26'0.5510335321221183''
 Mizar altitude: 0.37716758735416495rad = 21.61011092452545°
+
+Due to rounding errors in intermediate calculation steps of the test, the values are slightly off the expected value above.
+The result was compared with the display in https://stellarium.org/
+The data don't match the Stellarium display exactly as more complex formulas than Meeus's are used there. But the
+result is sufficiently precise for practical purposes in amateur astronomy. The observed discrepancy of the result
+coordinate was in the range of < 0.5°.
+
  */
-describe('localMST', () => {
+describe('coordinateConversion', () => {
+
+    test('polaris', () => {
 
-    test('local mst', () => {
+        // this is an example of the full calculation of a star position, from equatorial coordinates to
+        // horizontal coordinates
 
         let jdnUtc0 = gregorianDateToJulianDayNumber(2022, 11, 6);
         expect(jdnUtc0).toBe(2459889.5);
@@ -95,6 +109,22 @@ describe('localMST', () => {
         let lmst = greenwichMeanSiderealTimeToLocalMeanSiderealTime(gmst, 13.41 /* Berlin longitude*/);
         expect(lmst).toBe(82681.57318261731); // LMST Berlin is 22:58:01 at 19:00:29 UTC
 
+        const berlinLatitudeRad = 52.52 * 2 * Math.PI / 360;
+
+        let polaris = EquatorialCoordinate.fromTimeAndDegree(89, 15, 50.8, 2, 31, 49.09);
+
+        let alt = polaris.calculateAltitude(lmst / 86400 * 2 * Math.PI, berlinLatitudeRad);
+        let az = polaris.calculateAzimuth(lmst / 86400 * 2 * Math.PI, berlinLatitudeRad);
+
+        // this is the final result - we have converted equatorial coordinates to horizontal coordinates.
+        expect(alt).toBe(0.9242257590670654); // in radians - for the polar star, this is by definition approximately the same as berlinLatitudeRad
+        expect(az).toBe(0.017126264029972325); // in radians
+
+        // do the same calculation again, just to test fromEquatorialCoordinate()
+        let polarisHorizontalCoordinates =
+            HorizontalCoordinate.fromEquatorialCoordinate(polaris, lmst / 86400 * 2 * Math.PI, berlinLatitudeRad);
+        expect(polarisHorizontalCoordinates.altitude).toBe(alt);
+        expect(polarisHorizontalCoordinates.azimuth).toBe(az);
     });
 
 });