diff --git a/README.md b/README.md index df7dec9..ad1e1c9 100644 --- a/README.md +++ b/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 `` \ No newline at end of file +``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. \ No newline at end of file diff --git a/src/DateUtils.ts b/src/DateUtils.ts index e3ed466..e69e1a2 100644 --- a/src/DateUtils.ts +++ b/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; -} \ No newline at end of file +} + diff --git a/src/StarCoordinates.ts b/src/StarCoordinates.ts index 2106f2c..1716feb 100644 --- a/src/StarCoordinates.ts +++ b/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; + } } \ No newline at end of file diff --git a/tests/Gmst.test.ts b/tests/Gmst.test.ts index 7e10af7..0154ea9 100644 --- a/tests/Gmst.test.ts +++ b/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); }); }); \ No newline at end of file