WIP transits

Signed-off-by: Cédric Foellmi <cedric@onekiloparsec.dev>

src/coordinates.ts

@@ -6,6 +6,7 @@ import {
  HorizontalCoordinates,
  Hour,
  JulianDay,
+ Meter,
  Point
 } from './types'
 import { DEG2H, DEG2RAD, ECLIPTIC_OBLIQUITY_J2000_0, H2DEG, J2000, JULIAN_DAY_B1950_0, RAD2DEG } from './constants'
@@ -268,6 +269,22 @@ export function transformHorizontalToEquatorial (jd: JulianDay, alt: Degree, az:
  }
 }
 
+/**
+ * Transform equatorial coordinates to topocentric coordinates.
+ * @param {JulianDay} jd The julian day
+ * @param {Degree} coords The equatorial coordinates
+ * @param height
+ * @param {EquatorialCoordinates} lat The latitude of the observer's location
+ * @returns {EquatorialCoordinates}
+ */
+export function transformEquatorialToTopocentric (jd: JulianDay, coords: EquatorialCoordinates, height: Meter, lat: Degree): EquatorialCoordinates {
+ const pi: Degree = asin(sin(8.794 / 3600) / distance)
+ return {
+ rightAscension: fmod(getRightAscensionFromHorizontal(jd, alt, az, lat, lng) + 24.0, 24.0),
+ declination: getDeclinationFromHorizontal(jd, alt, az, lat)
+ }
+}
+
 /**
  * Transform a point (x,y) of the sky projected on a disk to horizontal coordinates.
  * @param {Point} point The point on the disk, relative to its center

src/transits.ts

@@ -14,17 +14,194 @@ import {
 } from './constants'
 import { getDate, getJulianDay, getJulianDayMidnight, getLocalSiderealTime } from './juliandays'
 import { fmod } from './utils'
+import { getDeltaT } from 'src/times/utils'
+import { getHorizontalAltitude } from 'src/coordinates'
 
 dayjs.extend(utc)
 
-
 const sin = Math.sin
 const cos = Math.cos
 const asin = Math.asin
 const acos = Math.acos
 const abs = Math.abs
 const floor = Math.floor
 
+// See AA p24
+function getInterpolatedValue (y1: number, y2: number, y3: number, n: number): number {
+ const a = y2 - y1
+ const b = y3 - y2
+ const c = b - a
+ return y2 + (n / 2) * (a + b + n * c)
+}
+
+// See AA, p102
+function getMTimes (jd: JulianDay, ra: Hour, dec: Degree, lng: Degree, lat: Degree, alt: Degree = STANDARD_ALTITUDE_STARS): {
+ m0: number,
+ m1: number,
+ m2: number,
+ isCircumpolar: boolean
+} {
+ // Getting the UT 0h on day D. See AA p.102.
+ // It is not equal to the expected "0h Dynamical Time" of the coordinates ra and dec.
+ const jd0 = getJulianDayMidnight(jd)
+
+ // Calculate the Greenwich sidereal time in degrees
+ let Theta0 = getLocalSiderealTime(jd0, 0) * H2DEG
+
+ const sinh0 = sin(alt * DEG2RAD)
+ const sinPhi = sin(lat * DEG2RAD)
+ const sinDelta = sin(dec * DEG2RAD)
+ const cosPhi = cos(lat * DEG2RAD)
+ const cosDelta = cos(dec * DEG2RAD)
+
+ // Algorithms in AA use Positive West longitudes. The formula (15.2, p102):
+ // const m0 = (alpha2 + Longitude - Theta0) / 360
+ // thus becomes:
+ const m0 = fmod((ra * H2DEG - lng - Theta0) / 360, 1)
+
+ // Calculate cosH0. See AA Eq.15.1, p.102
+ let cosH0 = (sinh0 - sinPhi * sinDelta) / (cosPhi * cosDelta)
+ const H0 = acos(cosH0) * RAD2DEG
+ console.log(cosH0, H0, lat)
+
+ return {
+ m0, // transit
+ m1: fmod(m0 - H0 / 360, 1), // rise
+ m2: fmod(m0 + H0 / 360, 1), // set,
+ isCircumpolar: (abs(cosH0) > 1)
+ }
+}
+
+
+/**
+ * Compute the times of rise, set and transit of an object at a given date,
+ * and observer's location on Earth.
+ * @param {JulianDay} jd The julian day
+ * @param {Hour} ra The equatorial right ascension of the object
+ * @param {Degree} dec The The equatorial declination of the object
+ * @param {Degree} lng The observer's longitude
+ * @param {Degree} lat The observer's latitude
+ * @param {Degree} alt The local altitude of the object's center to consider
+ * for rise and set times. It's value isn't 0. For stars, it is affected by
+ * aberration (value = -0.5667 degree)
+ */
+export function getAccurateRiseSetTransitTimes (jd: JulianDay, ra: [Hour, Hour, Hour], dec: [Degree, Degree, Degree], lng: Degree, lat: Degree, alt: Degree = STANDARD_ALTITUDE_STARS): RiseSetTransit {
+ // We assume the target coordinates are the mean equatorial coordinates for the epoch and equinox J2000.0.
+ // Furthermore, we assume we don't need to take proper motion to take into account. See AA p135.
+
+ // @ts-ignore
+ const result: RiseSetTransit = {
+ rise: {
+ utc: undefined,
+ julianDay: undefined
+ },
+ set: {
+ utc: undefined,
+ julianDay: undefined
+ },
+ transit: {
+ utc: undefined,
+ julianDay: undefined,
+ altitude: undefined,
+ refAltitude: alt,
+ isAboveHorizon: false,
+ isAboveAltitude: false, // for when altitude is not that of horizon
+ isCircumpolar: false
+ }
+ }
+
+ // Getting the UT 0h on day D. See AA p.102.
+ const jd0 = getJulianDayMidnight(jd)
+
+ // Calculate the Greenwich sidereal time in degrees
+ const Theta0 = getLocalSiderealTime(jd0, 0) * H2DEG
+ console.log(Theta0, getDate(jd0))
+
+ const cosPhi = cos(lat * DEG2RAD)
+
+ // Equ 13.6, AA p93, with cosH = 1, that is H (hour angle) = 0
+ result.transit.altitude = getTransitAltitude(ra[1], dec[1], lng, lat)
+ result.transit.isAboveHorizon = (result.transit.altitude > STANDARD_ALTITUDE_STARS)
+ result.transit.isAboveAltitude = (result.transit.altitude > alt)
+
+ const mTimesD = getMTimes(jd, ra[1], dec[1], lng, lat, alt)
+ console.log(mTimesD)
+
+ const theta0_m0 = fmod(Theta0 + 360.985647 * mTimesD.m0, 360)
+ const theta0_m1 = fmod(Theta0 + 360.985647 * mTimesD.m1, 360)
+ const theta0_m2 = fmod(Theta0 + 360.985647 * mTimesD.m2, 360)
+ console.log(theta0_m1, theta0_m0, theta0_m2)
+
+ const DeltaT = getDeltaT(jd)
+
+ // R.A. See AA, p24
+ const n_m0 = mTimesD.m0 + DeltaT / 86400
+ const n_m1 = mTimesD.m1 + DeltaT / 86400
+ const n_m2 = mTimesD.m2 + DeltaT / 86400
+ console.log('DeltaT', DeltaT, n_m1, n_m0, n_m2)
+
+ const alpha_m0 = getInterpolatedValue(ra[0], ra[1], ra[2], n_m0)
+ const alpha_m1 = getInterpolatedValue(ra[0], ra[1], ra[2], n_m1)
+ const alpha_m2 = getInterpolatedValue(ra[0], ra[1], ra[2], n_m2)
+ const delta_m1 = getInterpolatedValue(dec[0], dec[1], dec[2], n_m1)
+ const delta_m2 = getInterpolatedValue(dec[0], dec[1], dec[2], n_m2)
+
+ console.log('ALphas', alpha_m1 * H2DEG, alpha_m0 * H2DEG, alpha_m2 * H2DEG)
+ console.log('Deltas', delta_m1, '-', delta_m2)
+
+ const H_m0 = theta0_m0 + lng - alpha_m0 * H2DEG
+ const Delta_m0 = -H_m0 / 360
+ const m0 = mTimesD.m0 + Delta_m0
+
+ const H_m1 = theta0_m1 + lng - alpha_m1 * H2DEG
+ const h_m1 = asin(sin(lat * DEG2RAD) * sin(delta_m1 * DEG2RAD) + cos(lat * DEG2RAD) * cos(delta_m1 * DEG2RAD) * cos(H_m1 * DEG2RAD)) * RAD2DEG
+ const Delta_m1 = (h_m1 - alt) / (360 * cos(delta_m1) * cosPhi * sin(H_m1))
+ const m1 = mTimesD.m1 + Delta_m1
+
+ const H_m2 = theta0_m2 + lng - alpha_m2 * H2DEG
+ const h_m2 = asin(sin(lat * DEG2RAD) * sin(delta_m2 * DEG2RAD) + cos(lat * DEG2RAD) * cos(delta_m2 * DEG2RAD) * cos(H_m2 * DEG2RAD)) * RAD2DEG
+ const Delta_m2 = (h_m2 - alt) / (360 * cos(delta_m2) * cosPhi * sin(H_m2))
+ const m2 = mTimesD.m2 + Delta_m2
+
+ console.log('H', H_m1, H_m0, H_m2)
+ console.log('h', h_m1, '-', h_m2)
+ console.log('Delta ms', Delta_m1, Delta_m0, Delta_m2)
+
+ result.transit.utc = m0 * 24
+
+ const utcMoment = dayjs.utc(getDate(jd))
+ const hourTransit = floor(result.transit.utc)
+ const minuteTransit = result.transit.utc - hourTransit
+ result.transit.julianDay = getJulianDay(utcMoment.hour(hourTransit).minute(minuteTransit * 60).toDate())
+
+ // Calculate cosH0. See AA Eq.15.1, p.102
+ result.transit.isCircumpolar = mTimesD.isCircumpolar
+
+ if (!result.transit.isCircumpolar) {
+ result.rise.utc = m1
+ result.set.utc = m2
+
+ const hourRise = floor(result.rise.utc)
+ const minuteRise = result.rise.utc - hourRise
+ const hourSet = floor(result.set.utc)
+ const minuteSet = result.set.utc - hourSet
+
+ // Staying within a precision of a minute
+ result.rise.julianDay = getJulianDay(utcMoment.hour(hourRise).minute(minuteRise * 60).toDate())
+ result.set.julianDay = getJulianDay(utcMoment.hour(hourSet).minute(minuteSet * 60).toDate())
+ }
+
+ if (result.rise.julianDay && result.transit.julianDay && result.rise.julianDay > result.transit.julianDay) {
+ result.rise.julianDay -= 1
+ }
+ if (result.set.julianDay && result.transit.julianDay && result.set.julianDay < result.transit.julianDay) {
+ result.set.julianDay += 1
+ }
+
+ return result
+}
+
+
 /**
  * Compute the times of rise, set and transit of an object at a given date,
  * and observer's location on Earth.

tests/transits.test.js

@@ -1,4 +1,4 @@
-import { constants, juliandays, transits, Venus } from '../src'
+import { constants, juliandays, transits, times, Venus } from '../src'
 import * as utils from '../src/utils.ts'
 
 describe('transits of exoplanets', () => {
@@ -25,16 +25,50 @@ test('circumpolar transit', () => {
 })
 
 
-// See AA p 103
+// See AA, pp 103 & 104
 test('approximate Venus on 1988 March 20 at Boston', () => {
- const date = new Date(Date.UTC(1988, 2, 20))
- const results = transits.getRiseSetTransitTimes(juliandays.getJulianDay(date), 41.73129 * constants.DEG2H, 18.44092, -71.0833, 42.3333)
+ const date = new Date(Date.UTC(1988, 2, 20, 0, 0, 0))
+ const jd = juliandays.getJulianDay(date)
+ const coordsBoston = { latitude: 42.3333, longitude: -71.0833 }
+ const coordsVenus = { rightAscension: 41.73129, declination: 18.44092 }
+ const results = transits.getRiseSetTransitTimes(
+ jd,
+ coordsVenus.rightAscension * constants.DEG2H,
+ coordsVenus.declination,
+ coordsBoston.longitude,
+ coordsBoston.latitude
+ )
  expect(results.transit.isCircumpolar).toBeFalsy()
  expect(results.transit.isAboveHorizon).toBeTruthy()
  expect(results.transit.isAboveAltitude).toBeTruthy()
- expect(results.rise.utc).toBeCloseTo(12.43608, 3)
- expect(results.transit.utc).toBeCloseTo(19.6716, 3)
- expect(results.set.utc).toBeCloseTo(2.90712, 4)
+ expect(results.rise.utc).toBeCloseTo(24 * 0.51766, 2)
+ // expect(results.transit.utc).toBeCloseTo(24 * 0.81980, 1)
+ expect(results.set.utc).toBeCloseTo(24 * 0.12130, 6)
+ expect(results.rise.julianDay < results.transit.julianDay).toBeTruthy()
+ expect(results.transit.julianDay < results.set.julianDay).toBeTruthy()
+})
+
+
+// See AA, pp 103 & 104
+test('accurate Venus on 1988 March 20 at Boston', () => {
+ const date = new Date(Date.UTC(1988, 2, 20, 0, 0, 0))
+ const jd = juliandays.getJulianDay(date)
+ const coordsBoston = { latitude: 42.3333, longitude: -71.0833 }
+ const rasVenus = [40.68021 * constants.DEG2H, 41.73129 * constants.DEG2H, 42.78204 * constants.DEG2H]
+ const decVenus = [18.04762, 18.44092, 18.82742]
+ const results = transits.getAccurateRiseSetTransitTimes(
+ jd,
+ rasVenus,
+ decVenus,
+ coordsBoston.longitude,
+ coordsBoston.latitude
+ )
+ expect(results.transit.isCircumpolar).toBeFalsy()
+ expect(results.transit.isAboveHorizon).toBeTruthy()
+ expect(results.transit.isAboveAltitude).toBeTruthy()
+ expect(results.rise.utc).toBeCloseTo(24 * 0.51766, 2)
+ // expect(results.transit.utc).toBeCloseTo(24 * 0.81980, 1)
+ expect(results.set.utc).toBeCloseTo(24 * 0.12130, 6)
  expect(results.rise.julianDay < results.transit.julianDay).toBeTruthy()
  expect(results.transit.julianDay < results.set.julianDay).toBeTruthy()
 })
@@ -46,12 +80,15 @@ test('approximate Venus on 2023 October 14 at UMBC', () => {
  // Yet the topocentric corrections are usually small (< 30").
  // Note, months are 0-based, hence 9 = october.
  const jd = juliandays.getJulianDay(new Date(Date.UTC(2023, 9, 14, 0, 0)))
+ const jd0 = times.transformUTC2TT(jd)
  const umbc = { longitude: -76.70978311382578, latitude: 39.25443537644697 }
- const coords = Venus.getApparentEquatorialCoordinates(jd)
- const results = transits.getRiseSetTransitTimes(
+ const coords0 = Venus.getApparentEquatorialCoordinates(jd0 - 1)
+ const coords1 = Venus.getApparentEquatorialCoordinates(jd0)
+ const coords2 = Venus.getApparentEquatorialCoordinates(jd0 + 1)
+ const results = transits.getAccurateRiseSetTransitTimes(
  jd,
- coords.rightAscension,
- coords.declination,
+ [coords0.rightAscension, coords1.rightAscension, coords2.rightAscension],
+ [coords0.declination, coords1.declination, coords2.declination],
  umbc.longitude,
  umbc.latitude
  )