Merge pull request #1728 from sraaphorst/ghost/REL-3179

REL-3179: Proper motion calculations for GHOST

bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Epoch.scala

@@ -1,18 +1,54 @@
 package edu.gemini.spModel.core
 
-import scalaz._, Scalaz._
-
-/** Epoch in Gregorian (?) years. */
-final case class Epoch(year: Double)
-
+import java.time._
+
+import scalaz._
+import Scalaz._
+
+/**
+ * An epoch, the astronomer's equivalent of `Instant`, based on a fractional year in some temporal
+ * scheme (Julian, in this case) that determines year zero and the length of a year. The only
+ * meaningful operation for an `Epoch` is to ask the elapsed epoch-years between it and some other
+ * point in time. We need this for proper motion corrections because velocities are measured in
+ * motion per epoch-year.
+ *
+ * @see The Wikipedia [[https://en.wikipedia.org/wiki/Epoch_(astronomy) article]]
+ */
+final case class Epoch(year: Double) {
+  /** Offset in epoch-years from this `Epoch` to the given `Instant`. */
+  def untilInstant(i: Instant): Double =
+    untilLocalDateTime(LocalDateTime.ofInstant(i, ZoneOffset.UTC))
+
+  /** Offset in epoch-years from this `Epoch` to the given `LocalDateTime`. */
+  def untilLocalDateTime(ldt: LocalDateTime): Double =
+    untilJulianDay(Epoch.toJulianDay(ldt))
+
+  /** Offset in epoch-years from this `Epoch` to the given fractional Julian day. */
+  def untilJulianDay(jd: Double): Double =
+    Epoch.JulianYearBasis + (jd - Epoch.JulianBasis) / Epoch.JulianLengthOfYear
+
+  /** Offset in epoch-years from this `Epoch` to the given epoch year under the same scheme. */
+  def untilEpochYear(epochYear: Double): Double =
+    epochYear - year
+}
 
 object Epoch {
+  /**
+   * Standard J2000 epoch and Julian constants.
+   */
+  val JulianYearBasis: Double = 2000.0
+  val JulianBasis: Double = 2451545.0
+  val JulianLengthOfYear: Double = 365.25
+  val J2000: Epoch = Epoch(JulianYearBasis)
+
+  def toJulianDay(dt: LocalDateTime): Double =
+    JulianDate.ofLocalDateTime(dt).dayNumber.toDouble
 
   val year: Epoch @> Double = Lens.lensu((a, b) => a.copy(year = b), _.year)
 
-  val J2000 = Epoch(2000.0)
+  implicit val EpochOrder: Order[Epoch] =
+    Order.orderBy(_.year)
 
+  implicit val EpochShow: Show[Epoch] =
+    Show.showFromToString
 }
-
-
-

bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/JulianDate.scala

@@ -0,0 +1,99 @@
+package edu.gemini.spModel.core
+
+import java.time.{Instant, LocalDateTime}
+import java.time.ZoneOffset.UTC
+
+import scalaz._
+import Scalaz._
+
+sealed abstract case class JulianDate(dayNumber: Int,
+                                      nanoAdjustment: Long) {
+
+  import JulianDate._
+
+  // Guaranteed by the JulianDate constructors, double checked here.
+  assert(dayNumber >= 0, s"dayNumber >= 0")
+  assert(nanoAdjustment >= MinAdjustment, s"nanoAdjustment >= $MinAdjustment")
+  assert(nanoAdjustment <= MaxAdjustment, s"nanoAdjustment <= $MaxAdjustment")
+
+  /** Julian date value as a Double, including Julian Day Number and fractional
+   * day since the preceding noon.
+   */
+  val toDouble: Double =
+    dayNumber + nanoAdjustment.toDouble / NanoPerDay.toDouble
+
+  /** Modified Julian Date (MJD) double.  This is logically the same as
+   * `toDouble - 2400000.5`. MJD was introduced to preserve a bit of floating
+   * point decimal precision in calculations that use Julian dates.  It also
+   * makes it easier to directly implement algorithms that work with MJD.
+   *
+   * @see http://tycho.usno.navy.mil/mjd.html
+   */
+  def toModifiedDouble: Double = {
+    val h = SecondsPerHalfDay.toLong * Billion
+    val d = dayNumber      - 2400000
+    val n = nanoAdjustment - h
+
+    val (dʹ,nʹ) = if (n >= MinAdjustment) (d, n)
+    else (d - 1, n + SecondsPerDay.toLong * Billion)
+
+    dʹ + nʹ.toDouble / NanoPerDay.toDouble
+  }
+}
+
+object JulianDate {
+
+  /** JulianDate and related constants. */
+  val SecondsPerDay: Int = // 86400
+    24 * 60 * 60
+
+  val SecondsPerHalfDay: Int = // 43200
+    SecondsPerDay / 2
+
+  val Billion: Long = 1000000000
+  val NanoPerDay: Long = SecondsPerDay.toLong * Billion
+
+  val MinAdjustment: Long = -SecondsPerHalfDay.toLong * Billion
+  val MaxAdjustment: Long = SecondsPerHalfDay.toLong * Billion - 1
+
+  /** J2000 reference epoch as Julian Date. */
+  val J2000: JulianDate = // JulianDate(2451545,0)
+    JulianDate.ofLocalDateTime(
+      LocalDateTime.of(2000, 1, 1, 12, 0, 0)
+    )
+
+  /** Convert an `Instant` to a Julian Date.
+   */
+  def ofInstant(i: Instant): JulianDate =
+    ofLocalDateTime(LocalDateTime.ofInstant(i, UTC))
+
+  /** JulianDate from a `LocalDateTime` assumed to represent a time at UTC.
+   */
+  def ofLocalDateTime(ldt: LocalDateTime): JulianDate = {
+    val y = ldt.getYear
+    val m = ldt.getMonthValue
+    val d = ldt.getDayOfMonth
+
+    // Julian Day Number algorithm from:
+    // Fliegel, H.F. and Van Flandern, T.C. (1968). "A Machine Algorithm for
+    // Processing Calendar Dates" Communications of the Association of Computing
+    // Machines ll, 6sT.
+
+    // Yes, integer division.  -1 for Jan and Feb. 0 for Mar - Dec.
+    val t = (m - 14) / 12
+
+    // Julian Day Number (integer division).
+    val jdn = (1461 * (y + 4800 + t)) / 4 +
+      (367 * (m - 2 - 12 * t)) / 12 -
+      (3 * ((y + 4900 + t) / 100)) / 4 +
+      d - 32075
+
+    // Whole seconds since midnight
+    val secs = ldt.getHour * 3600 + ldt.getMinute * 60 + ldt.getSecond
+    val adj = (secs - SecondsPerHalfDay).toLong * Billion + ldt.getNano
+
+    new JulianDate(jdn, adj) {}
+  }
+
+  implicit val JulianDateEquals: Equal[JulianDate] = Equal.equalA
+}

bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/ProperMotion.scala

@@ -1,8 +1,12 @@
 package edu.gemini.spModel.core
 
+import java.time.Instant
+
 import scalaz._
 import Scalaz._
 
+import scala.math.{atan2, cos, hypot, sin}
+
 /**
  * Specification of proper motion.
  * @param deltaRA angular velocity in right ascension per year
@@ -11,16 +15,125 @@ import Scalaz._
 case class ProperMotion(
   deltaRA: RightAscensionAngularVelocity,
   deltaDec: DeclinationAngularVelocity,
-  epoch: Epoch = Epoch.J2000)
+  epoch: Epoch = Epoch.J2000) {
+
+  /** Coordinates at the current time. */
+  def calculateAt(t: Target, i: Instant): Coordinates =
+    plusYears(t, epoch.untilInstant(i))
+
+  /** Coordinates `elapsedYears` fractional epoch-years after `epoch`. */
+  def plusYears(t: Target, elapsedYears: Double): Coordinates = {
+    val baseCoordinates = t.coords(None).getOrElse(Coordinates.zero)
+
+    val properVelocity: Offset = Offset(OffsetP(Angle.fromDegrees(deltaRA.velocity.toDegreesPerYear)),
+                                        OffsetQ(Angle.fromDegrees(deltaDec.velocity.toDegreesPerYear)))
+
+    val radialVelocity: Double = {
+      val redshift: Option[Redshift] = Target.redshift.get(t).flatten
+      redshift.getOrElse(Redshift.zero).toRadialVelocity.toKilometersPerSecond
+    }
+    val parallax = Parallax.mas.get(Target.parallax.get(t).flatten.getOrElse(Parallax.zero)) / ProperMotion.Million
+    ProperMotion.properMotionCalculator(baseCoordinates,
+                                        Epoch.JulianLengthOfYear,
+                                        properVelocity,
+                                        radialVelocity,
+                                        parallax,
+                                        elapsedYears)
+  }
+}
 
 object ProperMotion {
+  // Some constants we need
+  val AstronomicalUnit:    Long   = 149597870660L
+  private val secsPerDay:  Double = 86400.0
+  private val auPerKm:     Double = 1000.0 / AstronomicalUnit.toDouble
+  private val radsPerAsec: Double = Angle.fromArcsecs(1).toRadians
+  private val TwoPi:       Double = 6.283185307179586476925286766559
+  private val Million:     Double = 1000000.0
+
+  private type Vec2 = (Double, Double)
+  private type Vec3 = (Double, Double, Double)
+
+  // Scalar multiplication and addition for Vec3.
+  private implicit class Vec3Ops(a: Vec3) {
+    def *(d: Double): Vec3 =
+      (a._1 * d, a._2 * d, a._3 * d)
+
+    def |+|(d: Double): Vec3 =
+      (a._1 + d, a._2 + d, a._3 + d)
+
+    def |+|(other: Vec3): Vec3 =
+      (a._1 + other._1, a._2 + other._2, a._3 + other._3)
+  }
+
+  /**
+   * Calculator to determine the new coordinates due to proper motion from epoch after elapsed years.
+   * @param baseCoordinates base coordinates of the target
+   * @param daysPerYear     length of epoch year in fractional days
+   * @param properVelocity  proper velocity per epoch year
+   * @param radialVelocity  radial velocity (km / s, positive if receding)
+   * @param parallax        parallax
+   * @param elapsedYears    elapsed time in epoch years
+   * @return
+   */
+  def properMotionCalculator(baseCoordinates: Coordinates,
+                             daysPerYear: Double,
+                             properVelocity: Offset,
+                             radialVelocity: Double,
+                             parallax: Double,
+                             elapsedYears: Double): Coordinates = {
+    // We want the baseCoordinates in radians.
+    val (ra,  dec)  = (baseCoordinates.ra.toAngle.toSignedDegrees.toRadians, baseCoordinates.dec.toAngle.toSignedDegrees.toRadians)
+    val (dRa, dDec) = (properVelocity.p.toAngle.toSignedDegrees.toRadians,   properVelocity.q.toAngle.toSignedDegrees.toRadians)
+
+    val pos: Vec3 = {
+      val cd = cos(dec)
+      (cos(ra) * cd, sin(ra) * cd, sin(dec))
+    }
+
+    // Change per year due to radial velocity and parallax. The units work out to asec/y.
+    val dPos1: Vec3 =
+        pos            *
+        daysPerYear    *
+        secsPerDay     *
+        radsPerAsec    *
+        auPerKm        *
+        radialVelocity *
+        parallax
+
+    // Change per year due to proper velocity
+    val dPos2 = (
+      -dRa * pos._2 - dDec * cos(ra) * sin(dec),
+       dRa * pos._1 - dDec * sin(ra) * sin(dec),
+                      dDec *           cos(dec)
+    )
+
+    // Our new position (still in polar coordinates). `|+|` here is scalar addition.
+    val pp = pos |+| ((dPos1 |+| dPos2) * elapsedYears)
+
+    // Back to spherical
+    val (x, y, z) = pp
+    val r    = hypot(x, y)
+    val rap  = if (r === 0.0) 0.0 else atan2(y, x)
+    val decp = if (z === 0.0) 0.0 else atan2(z, r)
+    val rapp = {
+      // Normalize to [0 .. 2π)
+      val rem = rap % TwoPi
+      if (rem < 0.0) rem + TwoPi else rem
+    }
+
+    Coordinates(RA.fromAngle(Angle.fromRadians(rapp)),
+                Declination.fromAngle(Angle.fromRadians(decp))
+                  .getOrElse(sys.error(s"Invalid declination: $decp radians")))
+  }
+
   val zero = ProperMotion(RightAscensionAngularVelocity.Zero, DeclinationAngularVelocity.Zero, Epoch.J2000)
 
   val deltaRA : ProperMotion @> RightAscensionAngularVelocity = Lens(t => Store(s => t.copy(deltaRA = s), t.deltaRA))
   val deltaDec: ProperMotion @> DeclinationAngularVelocity    = Lens(t => Store(s => t.copy(deltaDec = s), t.deltaDec))
   val epoch:    ProperMotion @> Epoch                         = Lens(t => Store(s => t.copy(epoch = s), t.epoch))
 
-  // Warning: This assumes the same epoch across proper motion.
+  // Warning: This assumes the same epoch across proper motion, which we take to be J2000.
   // We use it in GhostAsterism to simplify, where we assume same epoch.
   implicit val monoid: Monoid[ProperMotion] = Monoid.instance(
     (pm1,pm2) => ProperMotion(pm1.deltaRA |+| pm2.deltaRA, pm1.deltaDec |+| pm2.deltaDec),

bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Target.scala

@@ -111,6 +111,13 @@ trait TargetLenses {
       PLens.nil.run
     ))
 
+  val parallax: Target @?> Option[Parallax] =
+    PLens(_.fold(
+      PLens.nil.run,
+      pRunTarget(SiderealTarget.parallax.partial),
+      PLens.nil.run
+    ))
+
 }
 
 /**