diff --git a/generate/template/astronomy.go b/generate/template/astronomy.go index 54d789b1..48b37eba 100644 --- a/generate/template/astronomy.go +++ b/generate/template/astronomy.go @@ -14,6 +14,7 @@ import ( const ( DaysPerTropicalYear = 365.24217 // the number of days in one tropical year + SecondsPerDay = 86400.0 // the number of seconds in one day SpeedOfLightAuPerDay = 173.1446326846693 // the speed of light in vacuum expressed in astronomical units per day KmPerAu = 1.4959787069098932e+8 // the number of kilometers in one astronomical unit AuPerLightYear = 63241.07708807546 // the number of astronomical units in one light year @@ -1523,6 +1524,87 @@ func clampIndex(frac float64, nsteps int) int { //--- Pluto calc ends +type ascentInfo struct { + valid bool + tx AstroTime + ty AstroTime + ax float64 + ay float64 +} + +type SearchContext interface { + Eval(time AstroTime) (float64, error) +} + +func findAscent(depth int, context SearchContext, maxDerivAlt float64, t1, t2 AstroTime, a1, a2 float64) ascentInfo { + // See if we can find any time interval where the altitude-diff function + // rises from non-positive to positive. + + if a1 < 0.0 && a2 >= 0.0 { + // Trivial success case: the endpoints already rise through zero. + return ascentInfo{valid: true, tx: t1, ty: t2, ax: a1, ay: a2} + } + + if a1 >= 0.0 && a2 < 0.0 { + // Trivial failure case: Assume Nyquist condition prevents an ascent. + return ascentInfo{valid: false} + } + + if depth > 17 { + // Safety valve: do not allow unlimited recursion. + // This should never happen if the rest of the logic is working correctly, + // so fail the whole search if it does happen. It's a bug! + panic("Excessive recursion in rise/set ascent search.") + } + + // Both altitudes are on the same side of zero: both are negative, or both are non-negative. + // There could be a convex "hill" or a concave "valley" that passes through zero. + // In polar regions sometimes there is a rise/set or set/rise pair within minutes of each other. + // For example, the Moon can be below the horizon, then the very top of it becomes + // visible (moonrise) for a few minutes, then it moves sideways and down below + // the horizon again (moonset). We want to catch these cases. + // However, for efficiency and practicality concerns, because the rise/set search itself + // has a 0.1 second threshold, we do not worry about rise/set pairs that are less than + // one second apart. These are marginal cases that are rendered highly uncertain + // anyway, due to unpredictable atmospheric refraction conditions (air temperature and pressure). + + dt := t2.Ut - t1.Ut + if dt*SecondsPerDay < 1.0 { + return ascentInfo{valid: false} + } + + // Is it possible to reach zero from the altitude that is closer to zero? + da := math.Min(math.Abs(a1), math.Abs(a2)) + + // Without loss of generality, assume |a1| <= |a2|. + // (Reverse the argument in the case |a2| < |a1|.) + // Imagine you have to "drive" from a1 to 0, then back to a2. + // You can't go faster than max_deriv_alt. If you can't reach 0 in half the time, + // you certainly don't have time to reach 0, turn around, and still make your way + // back up to a2 (which is at least as far from 0 than a1 is) in the time interval dt. + // Therefore, the time threshold is half the time interval, or dt/2. + if da > maxDerivAlt*(dt/2) { + // Prune: the altitude cannot change fast enough to reach zero. + return ascentInfo{valid: false} + } + + // Bisect the time interval and evaluate the altitude at the midpoint. + tmid := TimeFromUniversalDays((t1.Ut + t2.Ut) / 2.0) + amid, err := context.Eval(tmid) + if err != nil { + panic("Altitude context should not have failed to evaluate.") + } + + // Recurse to the left interval. + ascent := findAscent(1+depth, context, maxDerivAlt, t1, tmid, a1, amid) + if !ascent.valid { + // Recurse to the right interval. + ascent = findAscent(1+depth, context, maxDerivAlt, tmid, t2, amid, a2) + } + + return ascent +} + //--- Generated code begins here ------------------------------------------------------------------ //$ASTRO_CONSTEL() diff --git a/source/golang/README.md b/source/golang/README.md index df356e0b..4a135b81 100644 --- a/source/golang/README.md +++ b/source/golang/README.md @@ -53,6 +53,7 @@ It provides a suite of well\-tested functions for calculating positions of the S - [func IdentityMatrix\(\) RotationMatrix](<#IdentityMatrix>) - [func InverseRotation\(rotation RotationMatrix\) RotationMatrix](<#InverseRotation>) - [func RotationEqdEqj\(time \*AstroTime\) RotationMatrix](<#RotationEqdEqj>) +- [type SearchContext](<#SearchContext>) - [type SeasonsInfo](<#SeasonsInfo>) - [type Spherical](<#Spherical>) - [type StateVector](<#StateVector>) @@ -70,6 +71,7 @@ It provides a suite of well\-tested functions for calculating positions of the S ```go const ( DaysPerTropicalYear = 365.24217 // the number of days in one tropical year + SecondsPerDay = 86400.0 // the number of seconds in one day SpeedOfLightAuPerDay = 173.1446326846693 // the speed of light in vacuum expressed in astronomical units per day KmPerAu = 1.4959787069098932e+8 // the number of kilometers in one astronomical unit AuPerLightYear = 63241.07708807546 // the number of astronomical units in one light year @@ -155,7 +157,7 @@ const ( ``` -## func [AngleBetween]() +## func [AngleBetween]() ```go func AngleBetween(avec AstroVector, bvec AstroVector) float64 @@ -164,7 +166,7 @@ func AngleBetween(avec AstroVector, bvec AstroVector) float64 AngleBetween calculates the angle in degrees between two vectors. Given a pair of vectors avec and bvec, this function returns the angle in degrees between the vectors in 3D space. The angle is measured in the plane that contains both vectors. The returned value is in the closed range \[0, 180\]. -## func [DaysFromCalendar]() +## func [DaysFromCalendar]() ```go func DaysFromCalendar(year, month, day, hour, minute int, second float64) float64 @@ -173,7 +175,7 @@ func DaysFromCalendar(year, month, day, hour, minute int, second float64) float6 -## func [DefineStar]() +## func [DefineStar]() ```go func DefineStar(body Body, ra, dec, distanceLightYears float64) error @@ -182,7 +184,7 @@ func DefineStar(body Body, ra, dec, distanceLightYears float64) error -## func [DegreesFromRadians]() +## func [DegreesFromRadians]() ```go func DegreesFromRadians(radians float64) float64 @@ -191,7 +193,7 @@ func DegreesFromRadians(radians float64) float64 DegreesFromRadians converts an angle expressed in radians to an angle expressed in degrees. -## func [Dot]() +## func [Dot]() ```go func Dot(a, b AstroVector) float64 @@ -200,7 +202,7 @@ func Dot(a, b AstroVector) float64 Returns the scalar dot product of two vectors. -## func [RadiansFromDegrees]() +## func [RadiansFromDegrees]() ```go func RadiansFromDegrees(degrees float64) float64 @@ -209,7 +211,7 @@ func RadiansFromDegrees(degrees float64) float64 RadiansFromDegrees converts an angle expressed in degrees to an angle expressed in radians. -## func [SiderealTime]() +## func [SiderealTime]() ```go func SiderealTime(time *AstroTime) float64 @@ -218,7 +220,7 @@ func SiderealTime(time *AstroTime) float64 Given a date and time, SiderealTime calculates the rotation of the Earth, represented by the equatorial angle of the Greenwich prime meridian with respect to distant stars \(not the Sun, which moves relative to background stars by almost one degree per day\). This angle is called Greenwich Apparent Sidereal Time \(GAST\). GAST is measured in sidereal hours in the half\-open range \[0, 24\). When GAST = 0, it means the prime meridian is aligned with the of\-date equinox, corrected at that time for precession and nutation of the Earth's axis. In this context, the "equinox" is the direction in space where the Earth's orbital plane \(the ecliptic\) intersects with the plane of the Earth's equator, at the location on the Earth's orbit of the \(seasonal\) March equinox. As the Earth rotates, GAST increases from 0 up to 24 sidereal hours, then starts over at 0. To convert to degrees, multiply the return value by 15. As an optimization, this function caches the sidereal time value in the time parameter. The value is reused later as needed, to avoid redundant calculations. -## type [AstroMoonQuarter]() +## type [AstroMoonQuarter]() @@ -230,7 +232,7 @@ type AstroMoonQuarter struct { ``` -## type [AstroSearchFunc]() +## type [AstroSearchFunc]() @@ -239,7 +241,7 @@ type AstroSearchFunc func(context interface{}, time AstroTime) float64 ``` -## type [AstroTime]() +## type [AstroTime]() @@ -287,7 +289,7 @@ type AstroTime struct { ``` -### func [TimeFromCalendar]() +### func [TimeFromCalendar]() ```go func TimeFromCalendar(year, month, day, hour, minute int, second float64) AstroTime @@ -296,7 +298,7 @@ func TimeFromCalendar(year, month, day, hour, minute int, second float64) AstroT TimeFromCalendar returns an AstroTime value for a date and time expressed in civil UTC. -### func [TimeFromTerrestrialDays]() +### func [TimeFromTerrestrialDays]() ```go func TimeFromTerrestrialDays(tt float64) AstroTime @@ -305,7 +307,7 @@ func TimeFromTerrestrialDays(tt float64) AstroTime TimeFromTerrestrialDays converts a Terrestrial Time \(TT\) day value, also known as ephemeris days, to an AstroTime value that can be used for astronomy calculations. -### func [TimeFromUniversalDays]() +### func [TimeFromUniversalDays]() ```go func TimeFromUniversalDays(ut float64) AstroTime @@ -314,7 +316,7 @@ func TimeFromUniversalDays(ut float64) AstroTime TimeFromUniversalDays converts a UTC number of days since January 1, 2000 into an AstroTime value that can be used for astronomy calculations. -### func \(AstroTime\) [AddDays]() +### func \(AstroTime\) [AddDays]() ```go func (time AstroTime) AddDays(days float64) AstroTime @@ -323,7 +325,7 @@ func (time AstroTime) AddDays(days float64) AstroTime Given an AstroTime value, creates a new AstroTime value that is the specified number of UT days in the future \(positive\) or past \(negative\). -## type [AstroVector]() +## type [AstroVector]() AstroVector represents a position in 3D space at a given time. Usually the distance components are expressed in astronomical units \(AU\). The origin and orientation system depends on context. Occasionally AstroVector is used to represent a velocity vector, in which case the component units are astronomical units per day. @@ -337,7 +339,7 @@ type AstroVector struct { ``` -### func [GeoMoon]() +### func [GeoMoon]() ```go func GeoMoon(time AstroTime) AstroVector @@ -346,7 +348,7 @@ func GeoMoon(time AstroTime) AstroVector GeoMoon calculates the equatorial geocentric position of the Moon at a given time. The returned vector indicates the Moon's center relative to the Earth's center. The vector components are expressed in AU \(astronomical units\). The coordinates are oriented with respect to the Earth's equator at the J2000 epoch. In Astronomy Engine, this orientation is called EQJ. -### func [RotateVector]() +### func [RotateVector]() ```go func RotateVector(rotation RotationMatrix, vector AstroVector) AstroVector @@ -355,7 +357,7 @@ func RotateVector(rotation RotationMatrix, vector AstroVector) AstroVector RotateVector applies a rotation to a vector, yielding a vector in another orientation system. -### func \(AstroVector\) [Length]() +### func \(AstroVector\) [Length]() ```go func (vec AstroVector) Length() float64 @@ -364,7 +366,7 @@ func (vec AstroVector) Length() float64 Returns the length of vec expressed in the same distance units as vec's components. -## type [AtmosphereInfo]() +## type [AtmosphereInfo]() AtmosphereInfo contains information about idealized atmospheric variables at a given elevation. @@ -377,7 +379,7 @@ type AtmosphereInfo struct { ``` -### func [Atmosphere]() +### func [Atmosphere]() ```go func Atmosphere(elevationMeters float64) (AtmosphereInfo, error) @@ -386,7 +388,7 @@ func Atmosphere(elevationMeters float64) (AtmosphereInfo, error) Atmosphere calculates U.S. Standard Atmosphere \(1976\) variables as a function of elevation. elevationMeters is the elevation above sea level at which to calculate atmospheric variables. It must be in the range \-500 to \+100000 or an error will occur. 1. COESA, U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, DC, 1976. 2. Jursa, A. S., Ed., Handbook of Geophysics and the Space Environment, Air Force Geophysics Laboratory, 1985. See: https://hbcp.chemnetbase.com/faces/documents/14_12/14_12_0001.xhtml https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf -## type [AxisInfo]() +## type [AxisInfo]() @@ -400,7 +402,7 @@ type AxisInfo struct { ``` -## type [Body]() +## type [Body]() @@ -409,7 +411,7 @@ type Body int ``` -## type [CalendarDateTime]() +## type [CalendarDateTime]() CalendarDateTime represents a Gregorian calendar date and time within plus or minus 1 million years from the year 0. @@ -425,7 +427,7 @@ type CalendarDateTime struct { ``` -### func [CalendarFromDays]() +### func [CalendarFromDays]() ```go func CalendarFromDays(ut float64) (*CalendarDateTime, error) @@ -434,7 +436,7 @@ func CalendarFromDays(ut float64) (*CalendarDateTime, error) CalendarFromDays converts a J2000 day value to a Gregorian calendar date and time. -## type [DeltaTimeFunc]() +## type [DeltaTimeFunc]() @@ -443,7 +445,7 @@ type DeltaTimeFunc func(ut float64) float64 ``` -## type [EclipseKind]() +## type [EclipseKind]() @@ -452,7 +454,7 @@ type EclipseKind int ``` -## type [Ecliptic]() +## type [Ecliptic]() A location of a body expressed in angular coordinates relative to the plane of the Earth's orbit around the Sun @@ -465,7 +467,7 @@ type Ecliptic struct { ``` -## type [Equatorial]() +## type [Equatorial]() A location of a body expressed in angular coordinates relative to the Earth's equator @@ -479,7 +481,7 @@ type Equatorial struct { ``` -## type [JupiterMoonsInfo]() +## type [JupiterMoonsInfo]() @@ -493,7 +495,7 @@ type JupiterMoonsInfo struct { ``` -### func [JupiterMoons]() +### func [JupiterMoons]() ```go func JupiterMoons(time AstroTime) JupiterMoonsInfo @@ -502,7 +504,7 @@ func JupiterMoons(time AstroTime) JupiterMoonsInfo Calculates Jovicentric positoins and velocities of Jupiter's largest 4 moons. Calculates position and velocity vectors for Jupiter's moons Io, Europa, Ganymede, and Callisto, at the given date and time. The vectors are jovicentric \(relative to the center of Jupiter\). Their orientation is the Earth's equatorial system at the J2000 epoch \(EQJ\). The position components are expressed in astronomical units \(AU\), and the velocity components are in AU/day. To convert to heliocentric position vectors, call HelioVector with Jupiter as the body to get Jupiter's heliocentric position, then add the jovicentric moon positions. Likewise, you can call \#Astronomy.GeoVector to convert to geocentric positions; however, you will have to manually correct for light travel time from the Jupiter system to Earth to figure out what time to pass to \`JupiterMoons\` to get an accurate picture of how Jupiter and its moons look from Earth. -## type [LibrationInfo]() +## type [LibrationInfo]() @@ -518,7 +520,7 @@ type LibrationInfo struct { ``` -## type [NodeEventInfo]() +## type [NodeEventInfo]() @@ -530,7 +532,7 @@ type NodeEventInfo struct { ``` -## type [NodeEventKind]() +## type [NodeEventKind]() @@ -539,7 +541,7 @@ type NodeEventKind int ``` -## type [Observer]() +## type [Observer]() The location of a point on or near the surface of the Earth @@ -552,7 +554,7 @@ type Observer struct { ``` -## type [Refraction]() +## type [Refraction]() @@ -571,7 +573,7 @@ const ( ``` -## type [RotationMatrix]() +## type [RotationMatrix]() RotationMatrix is a 3x3 matrix used to convert a vector from one orientation system to another. @@ -582,7 +584,7 @@ type RotationMatrix struct { ``` -### func [CombineRotation]() +### func [CombineRotation]() ```go func CombineRotation(a, b RotationMatrix) RotationMatrix @@ -591,7 +593,7 @@ func CombineRotation(a, b RotationMatrix) RotationMatrix CombineRotation combines the effects of two consecutive rotation matrices into a single rotation matrix. -### func [IdentityMatrix]() +### func [IdentityMatrix]() ```go func IdentityMatrix() RotationMatrix @@ -600,7 +602,7 @@ func IdentityMatrix() RotationMatrix Creates a rotation matrix that represents no rotation at all. -### func [InverseRotation]() +### func [InverseRotation]() ```go func InverseRotation(rotation RotationMatrix) RotationMatrix @@ -609,7 +611,7 @@ func InverseRotation(rotation RotationMatrix) RotationMatrix Calculates the inverse of a rotation matrix. Given a rotation matrix that performs some coordinate transform, this function returns the matrix that reverses that transform. -### func [RotationEqdEqj]() +### func [RotationEqdEqj]() ```go func RotationEqdEqj(time *AstroTime) RotationMatrix @@ -617,8 +619,19 @@ func RotationEqdEqj(time *AstroTime) RotationMatrix Calculates a rotation matrix that converts equator\-of\-date \(EQD\) to J2000 mean equator \(EQJ\). + +## type [SearchContext]() + + + +```go +type SearchContext interface { + Eval(time AstroTime) (float64, error) +} +``` + -## type [SeasonsInfo]() +## type [SeasonsInfo]() @@ -632,7 +645,7 @@ type SeasonsInfo struct { ``` -## type [Spherical]() +## type [Spherical]() Spherical coordinates for a body in space @@ -645,7 +658,7 @@ type Spherical struct { ``` -## type [StateVector]() +## type [StateVector]() StateVector represents the combined position and velocity of a body at a given moment of time. @@ -662,7 +675,7 @@ type StateVector struct { ``` -### func [RotateState]() +### func [RotateState]() ```go func RotateState(rotation RotationMatrix, state StateVector) StateVector @@ -671,7 +684,7 @@ func RotateState(rotation RotationMatrix, state StateVector) StateVector -### func \(StateVector\) [Position]() +### func \(StateVector\) [Position]() ```go func (state StateVector) Position() AstroVector @@ -680,7 +693,7 @@ func (state StateVector) Position() AstroVector Position returns the position vector inside a state vector. -### func \(StateVector\) [Velocity]() +### func \(StateVector\) [Velocity]() ```go func (state StateVector) Velocity() AstroVector @@ -689,7 +702,7 @@ func (state StateVector) Velocity() AstroVector Position returns the velocity vector inside a state vector. -## type [TimeFormat]() +## type [TimeFormat]() @@ -709,7 +722,7 @@ const ( ``` -## type [Topocentric]() +## type [Topocentric]() A location of a body as seen from an observer's point of view on or near the surface of the Earth. The topocentric position can optionally be corrected for atmospheric refraction. diff --git a/source/golang/astronomy.go b/source/golang/astronomy.go index ec0ba268..a51ef3f5 100644 --- a/source/golang/astronomy.go +++ b/source/golang/astronomy.go @@ -14,6 +14,7 @@ import ( const ( DaysPerTropicalYear = 365.24217 // the number of days in one tropical year + SecondsPerDay = 86400.0 // the number of seconds in one day SpeedOfLightAuPerDay = 173.1446326846693 // the speed of light in vacuum expressed in astronomical units per day KmPerAu = 1.4959787069098932e+8 // the number of kilometers in one astronomical unit AuPerLightYear = 63241.07708807546 // the number of astronomical units in one light year @@ -1523,6 +1524,87 @@ func clampIndex(frac float64, nsteps int) int { //--- Pluto calc ends +type ascentInfo struct { + valid bool + tx AstroTime + ty AstroTime + ax float64 + ay float64 +} + +type SearchContext interface { + Eval(time AstroTime) (float64, error) +} + +func findAscent(depth int, context SearchContext, maxDerivAlt float64, t1, t2 AstroTime, a1, a2 float64) ascentInfo { + // See if we can find any time interval where the altitude-diff function + // rises from non-positive to positive. + + if a1 < 0.0 && a2 >= 0.0 { + // Trivial success case: the endpoints already rise through zero. + return ascentInfo{valid: true, tx: t1, ty: t2, ax: a1, ay: a2} + } + + if a1 >= 0.0 && a2 < 0.0 { + // Trivial failure case: Assume Nyquist condition prevents an ascent. + return ascentInfo{valid: false} + } + + if depth > 17 { + // Safety valve: do not allow unlimited recursion. + // This should never happen if the rest of the logic is working correctly, + // so fail the whole search if it does happen. It's a bug! + panic("Excessive recursion in rise/set ascent search.") + } + + // Both altitudes are on the same side of zero: both are negative, or both are non-negative. + // There could be a convex "hill" or a concave "valley" that passes through zero. + // In polar regions sometimes there is a rise/set or set/rise pair within minutes of each other. + // For example, the Moon can be below the horizon, then the very top of it becomes + // visible (moonrise) for a few minutes, then it moves sideways and down below + // the horizon again (moonset). We want to catch these cases. + // However, for efficiency and practicality concerns, because the rise/set search itself + // has a 0.1 second threshold, we do not worry about rise/set pairs that are less than + // one second apart. These are marginal cases that are rendered highly uncertain + // anyway, due to unpredictable atmospheric refraction conditions (air temperature and pressure). + + dt := t2.Ut - t1.Ut + if dt*SecondsPerDay < 1.0 { + return ascentInfo{valid: false} + } + + // Is it possible to reach zero from the altitude that is closer to zero? + da := math.Min(math.Abs(a1), math.Abs(a2)) + + // Without loss of generality, assume |a1| <= |a2|. + // (Reverse the argument in the case |a2| < |a1|.) + // Imagine you have to "drive" from a1 to 0, then back to a2. + // You can't go faster than max_deriv_alt. If you can't reach 0 in half the time, + // you certainly don't have time to reach 0, turn around, and still make your way + // back up to a2 (which is at least as far from 0 than a1 is) in the time interval dt. + // Therefore, the time threshold is half the time interval, or dt/2. + if da > maxDerivAlt*(dt/2) { + // Prune: the altitude cannot change fast enough to reach zero. + return ascentInfo{valid: false} + } + + // Bisect the time interval and evaluate the altitude at the midpoint. + tmid := TimeFromUniversalDays((t1.Ut + t2.Ut) / 2.0) + amid, err := context.Eval(tmid) + if err != nil { + panic("Altitude context should not have failed to evaluate.") + } + + // Recurse to the left interval. + ascent := findAscent(1+depth, context, maxDerivAlt, t1, tmid, a1, amid) + if !ascent.valid { + // Recurse to the right interval. + ascent = findAscent(1+depth, context, maxDerivAlt, tmid, t2, amid, a2) + } + + return ascent +} + //--- Generated code begins here ------------------------------------------------------------------ var constelNames = [...]constelInfo{