diff --git a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Epoch.scala b/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Epoch.scala index d454ab1f4e..29a97c0d54 100644 --- a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Epoch.scala +++ b/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 } - - - diff --git a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/JulianDate.scala b/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/JulianDate.scala new file mode 100644 index 0000000000..f97c0a9da9 --- /dev/null +++ b/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 +} diff --git a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/ProperMotion.scala b/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/ProperMotion.scala index 95278bfeb6..87a29de68b 100644 --- a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/ProperMotion.scala +++ b/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), diff --git a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Target.scala b/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Target.scala index ef6b9a9830..519efcc946 100644 --- a/bundle/edu.gemini.spModel.core/src/main/scala/edu/gemini/spModel/core/Target.scala +++ b/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 + )) + } /** diff --git a/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/Arbitraries.scala b/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/Arbitraries.scala index 11a5d079cb..53c2488e94 100644 --- a/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/Arbitraries.scala +++ b/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/Arbitraries.scala @@ -1,15 +1,16 @@ package edu.gemini.spModel.core +import java.time.{LocalDate, LocalDateTime, LocalTime, Year, ZoneId} + import edu.gemini.spModel.core.WavelengthConversions._ import org.scalacheck._ import org.scalacheck.Gen._ import org.scalacheck.Arbitrary._ - import squants.motion.VelocityConversions._ import squants.radio.IrradianceConversions._ import squants.radio.SpectralIrradianceConversions._ - -import scalaz.{ ==>>, Order } +import scala.collection.JavaConverters._ +import scalaz.{==>>, Order} import scalaz.std.anyVal._ trait Arbitraries { @@ -101,7 +102,7 @@ trait Arbitraries { } yield Magnitude(value, band, error, system) } - implicit val arbDistribution = Arbitrary[SpectralDistribution] { + implicit val arbDistribution: Arbitrary[SpectralDistribution] = Arbitrary[SpectralDistribution] { Gen.oneOf( BlackBody(8000), BlackBody(10000), @@ -115,7 +116,7 @@ trait Arbitraries { LibraryNonStar.GammaDra ) } - implicit val arbProfile = Arbitrary[SpatialProfile] { + implicit val arbProfile: Arbitrary[SpatialProfile] = Arbitrary[SpatialProfile] { Gen.oneOf( PointSource, UniformSource, @@ -192,3 +193,62 @@ trait Arbitraries { case _ => sys.error("generated empty list!?") }) } + +// Arbitrary but reasonable dates and times. +trait ArbTime { + + implicit val arbZoneId: Arbitrary[ZoneId] = + Arbitrary { + oneOf(ZoneId.getAvailableZoneIds.asScala.toSeq).map(ZoneId.of) + } + + implicit val arbYear: Arbitrary[Year] = + Arbitrary { + choose(2000, 2020).map(Year.of) + } + + implicit val arbLocalDate: Arbitrary[LocalDate] = + Arbitrary { + for { + y <- arbitrary[Year] + d <- choose(1, y.length) + } yield LocalDate.ofYearDay(y.getValue, d) + } + + implicit val arbLocalTime: Arbitrary[LocalTime] = + Arbitrary { + for { + h <- choose(0, 23) + m <- choose(0, 59) + s <- choose(0, 59) + n <- choose(0, 999999999) + } yield LocalTime.of(h, m, s, n) + } + + implicit val arbLocalDateTime: Arbitrary[LocalDateTime] = + Arbitrary { + for { + d <- arbitrary[LocalDate] + t <- arbitrary[LocalTime] + } yield LocalDateTime.of(d, t) + } + + implicit val cogLocalDate: Cogen[LocalDate] = + Cogen[(Int, Int)].contramap(d => (d.getYear, d.getDayOfYear)) +} + +object ArbTime extends ArbTime + +trait ArbJulianDate { + import ArbTime._ + + implicit val arbJulianDate: Arbitrary[JulianDate] = + Arbitrary { + arbitrary[LocalDateTime].map(JulianDate.ofLocalDateTime) + } + + implicit val cogJulianDate: Cogen[JulianDate] = + Cogen[(Int, Long)].contramap(jd => (jd.dayNumber, jd.nanoAdjustment)) +} + +object ArbJulianDate extends ArbJulianDate \ No newline at end of file diff --git a/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/EpochSpec.scala b/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/EpochSpec.scala new file mode 100644 index 0000000000..a86cb4a7d3 --- /dev/null +++ b/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/EpochSpec.scala @@ -0,0 +1,22 @@ +package edu.gemini.spModel.core + +import scalaz._ +import Scalaz._ +import org.scalacheck.Prop._ +import org.specs2.ScalaCheck +import org.specs2.mutable.Specification + +final class EpochSpec extends Specification with ScalaCheck with Arbitraries { + "Epoch" should { + "give an offset of 0 to its own year" ! { + math.abs(Epoch.J2000.untilEpochYear(Epoch.J2000.year)) < 1e-9 + } + } + + "provide sanity in offsets using short" ! { + forAll { (s: Short) => + Epoch.J2000.untilEpochYear(Epoch.J2000.year + s.toDouble) shouldEqual s.toDouble + } + } +} + diff --git a/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/JulianDateSpec.scala b/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/JulianDateSpec.scala new file mode 100644 index 0000000000..71abefb772 --- /dev/null +++ b/bundle/edu.gemini.spModel.core/src/test/scala/edu/gemini/spModel/core/JulianDateSpec.scala @@ -0,0 +1,62 @@ +package edu.gemini.spModel.core + +import scalaz._ +import Scalaz._ +import org.scalacheck.Properties +import org.scalacheck.Prop.forAll +import org.scalacheck.Prop._ +import org.scalacheck.Gen +import org.specs2.ScalaCheck +import org.specs2.mutable.Specification +import AlmostEqual.AlmostEqualOps + +import java.time.LocalDateTime + +final class JulianDateSpec extends Specification with ScalaCheck with Arbitraries with Helpers { + + import ArbTime._ + import ArbJulianDate._ + + // Old Epoch.Scheme.toJulianDay algorithm + private def oldEpochSchemeJulianDate(dt: LocalDateTime): Double = { + import scala.math.floor + + val a = floor((14.0 - dt.getMonthValue) / 12.0) + val y = dt.getYear + 4800.0 - a + val m = dt.getMonthValue + 12 * a - 3.0 + dt.getDayOfMonth + + floor((153.0 * m + 2.0) / 5.0) + + 365 * y + + floor(y / 4.0) - + floor(y / 100.0) + + floor(y / 400.0) - + 32045.0 + } + + "JulianDate" should { + "respect Eq" ! + forAll { (a: JulianDate, b: JulianDate) => + a.equals(b) shouldEqual Equal[JulianDate].equal(a, b) + } + + "have dayNumber match old calculation" ! + forAll { (ldt: LocalDateTime) => + val jd0 = oldEpochSchemeJulianDate(ldt) + val jd1 = JulianDate.ofLocalDateTime(ldt).dayNumber.toDouble + + jd0 shouldEqual jd1 + } + + "Some specific dates compared to USNO calculations" ! { + JulianDate.ofLocalDateTime(LocalDateTime.of(1918, 11, 11, 11, 0, 0)).toDouble ~= 2421908.958333 + JulianDate.ofLocalDateTime(LocalDateTime.of(1969, 7, 21, 2, 56, 15)).toDouble ~= 2440423.622396 + JulianDate.ofLocalDateTime(LocalDateTime.of(2001, 9, 11, 8, 46, 0)).toDouble ~= 2452163.865278 + JulianDate.ofLocalDateTime(LocalDateTime.of(2345, 6, 7, 12, 0, 0)).toDouble ~= 2577711.000000 + } + + "Modified JulianDate should almost equal JulianDate - 2400000.5" ! + forAll { (j: JulianDate) => + j.toModifiedDouble ~= (j.toDouble - 2400000.5) + } + } +} \ No newline at end of file diff --git a/bundle/jsky.app.ot/src/main/scala/jsky/app/ot/gemini/editor/targetComponent/details2/SiderealDetailEditor.scala b/bundle/jsky.app.ot/src/main/scala/jsky/app/ot/gemini/editor/targetComponent/details2/SiderealDetailEditor.scala index 7902599e3f..8ff18cd28e 100644 --- a/bundle/jsky.app.ot/src/main/scala/jsky/app/ot/gemini/editor/targetComponent/details2/SiderealDetailEditor.scala +++ b/bundle/jsky.app.ot/src/main/scala/jsky/app/ot/gemini/editor/targetComponent/details2/SiderealDetailEditor.scala @@ -192,7 +192,7 @@ object RedshiftRepresentations { val repr: Map[RedshiftRepresentations, String] = Map( - RadialVelocity -> "km/sec", + RadialVelocity -> "km/s", RedshiftZ -> "", ApparentRadialVelocity -> "km/s" )