From b742aac138055afd3e995b306d88a3fb8e5d83d1 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Wed, 1 Nov 2023 16:38:30 -0400 Subject: [PATCH] Go: added function Horizon --- generate/template/astronomy.go | 134 ++++++++++++++++++++++++++++---- source/golang/README.md | 50 +++++++----- source/golang/astronomy.go | 134 ++++++++++++++++++++++++++++---- source/golang/astronomy_test.go | 43 ++++++++++ 4 files changed, 315 insertions(+), 46 deletions(-) diff --git a/generate/template/astronomy.go b/generate/template/astronomy.go index 036de886..beefee64 100644 --- a/generate/template/astronomy.go +++ b/generate/template/astronomy.go @@ -310,7 +310,7 @@ type CalendarDateTime struct { // CalendarFromDays converts a J2000 day value to a Gregorian calendar date and time. func CalendarFromDays(ut float64) (*CalendarDateTime, error) { if !isfinite(ut) { - return nil, errors.New("Non-finite value of ut") + return nil, errors.New("non-finite value of ut") } // Adapted from the NOVAS C 3.1 function cal_date(). @@ -352,11 +352,11 @@ func CalendarFromDays(ut float64) (*CalendarDateTime, error) { year := int(100*(n-49) + m + k - 400*c) if year < -999999 || year > +999999 { - return nil, errors.New("The supplied time is too far from the year 2000 to be represented.") + return nil, errors.New("the supplied time is too far from the year 2000 to be represented.") } if month < 1 || month > 12 || day < 1 || day > 31 { - return nil, errors.New("Internal error: invalid calendar date calculated.") + return nil, errors.New("internal error: invalid calendar date calculated.") } return &CalendarDateTime{year, month, day, hour, minute, second}, nil @@ -396,7 +396,7 @@ func TimeFromCalendar(year, month, day, hour, minute int, second float64) AstroT // https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf func Atmosphere(elevationMeters float64) (AtmosphereInfo, error) { if !isfinite(elevationMeters) || elevationMeters < -500.0 || elevationMeters > 100000.0 { - return AtmosphereInfo{}, errors.New("Invalid elevation") + return AtmosphereInfo{}, errors.New("invalid elevation") } const P0 = 101325.0 // pressure at sea level [pascals] const T0 = 288.15 // temperature at sea level [kelvins] @@ -1034,7 +1034,7 @@ func CorrectLightTravel(posFunc PositionFunction, time AstroTime) (*AstroVector, // values for stellar aberration angles. lt := pos.Length() / SpeedOfLightAuPerDay if lt > 1.0 { - return nil, errors.New("Object is too distant for light-travel solver.") + return nil, errors.New("object is too distant for light-travel solver.") } ltime2 := time.AddDays(-lt) dt := math.Abs(ltime2.Tt - ltime.Tt) @@ -1043,7 +1043,115 @@ func CorrectLightTravel(posFunc PositionFunction, time AstroTime) (*AstroVector, } ltime = ltime2 } - return nil, errors.New("Light travel time correction did not converge.") + return nil, errors.New("light travel time correction did not converge.") +} + +func Horizon(time AstroTime, observer Observer, ra, dec float64, refraction Refraction) (*Topocentric, error) { + sinlat := dsin(observer.Latitude) + coslat := dcos(observer.Latitude) + sinlon := dsin(observer.Longitude) + coslon := dcos(observer.Longitude) + sindc := dsin(dec) + cosdc := dcos(dec) + sinra := dsin(ra * 15.0) + cosra := dcos(ra * 15.0) + + // Calculate three mutually perpendicular unit vectors + // in equatorial coordinates: uze, une, uwe. + // + // uze = The direction of the observer's local zenith (straight up). + // une = The direction toward due north on the observer's horizon. + // uwe = The direction toward due west on the observer's horizon. + // + // HOWEVER, these are uncorrected for the Earth's rotation due to the time of day. + // + // The components of these 3 vectors are as follows: + // x = direction from center of Earth toward 0 degrees longitude (the prime meridian) on equator. + // y = direction from center of Earth toward 90 degrees west longitude on equator. + // z = direction from center of Earth toward the north pole. + uze := AstroVector{coslat * coslon, coslat * sinlon, sinlat, time} + une := AstroVector{-sinlat * coslon, -sinlat * sinlon, coslat, time} + uwe := AstroVector{sinlon, -coslon, 0.0, time} + + // Correct the vectors uze, une, uwe for the Earth's rotation by calculating + // sidereal time. Call spin() for each uncorrected vector to rotate about + // the Earth's axis to yield corrected unit vectors uz, un, uw. + // Multiply sidereal hours by -15 to convert to degrees and flip eastward + // rotation of the Earth to westward apparent movement of objects with time. + angle := -15.0 * SiderealTime(&time) + uz := spin(angle, uze) + un := spin(angle, une) + uw := spin(angle, uwe) + + // Convert angular equatorial coordinates (RA, DEC) to + // cartesian equatorial coordinates in 'p', using the + // same orientation system as uze, une, uwe. + p := AstroVector{cosdc * cosra, cosdc * sinra, sindc, time} + + // Use dot products of p with the zenith, north, and west + // vectors to obtain the cartesian coordinates of the body in + // the observer's horizontal orientation system. + // pz = zenith component [-1, +1] + // pn = north component [-1, +1] + // pw = west component [-1, +1] + pz := p.X*uz.X + p.Y*uz.Y + p.Z*uz.Z + pn := p.X*un.X + p.Y*un.Y + p.Z*un.Z + pw := p.X*uw.X + p.Y*uw.Y + p.Z*uw.Z + + // proj is the "shadow" of the body vector along the observer's flat ground. + proj := math.Hypot(pn, pw) + + // Calculate az = azimuth (compass direction clockwise from East.) + var az float64 + if proj > 0.0 { + // If the body is not exactly straight up/down, it has an azimuth. + // Invert the angle to produce degrees eastward from north. + az = -datan2(pw, pn) + if az < 0.0 { + az += 360.0 + } + } else { + // The body is straight up/down, so it does not have an azimuth. + // Report an arbitrary but reasonable value. + az = 0.0 + } + + // zd = the angle of the body away from the observer's zenith, in degrees. + zd := datan2(proj, pz) + horRa := ra + horDec := dec + + if refraction == NormalRefraction || refraction == JplHorizonsRefraction { + zd0 := zd + refr := RefractionAngle(refraction, 90.0-zd) + zd -= refr + + if refr > 0.0 && zd > 3.0e-4 { + sinzd := dsin(zd) + coszd := dcos(zd) + sinzd0 := dsin(zd0) + coszd0 := dcos(zd0) + + prx := ((p.X-coszd0*uz.X)/sinzd0)*sinzd + uz.X*coszd + pry := ((p.Y-coszd0*uz.Y)/sinzd0)*sinzd + uz.Y*coszd + prz := ((p.Z-coszd0*uz.Z)/sinzd0)*sinzd + uz.Z*coszd + + proj = math.Hypot(prx, pry) + if proj > 0.0 { + horRa = datan2(pry, prx) / 15.0 + if horRa < 0.0 { + horRa += 24.0 + } + } else { + horRa = 0.0 + } + horDec = datan2(prz, proj) + } + } else if refraction != NoRefraction { + return nil, errors.New("unsupported refraction option.") + } + + return &Topocentric{az, 90.0 - zd, horRa, horDec}, nil } type jupiterMoon struct { @@ -1298,16 +1406,16 @@ func userDefinedStar(body Body) *starDef { func DefineStar(body Body, ra, dec, distanceLightYears float64) error { star := getStar(body) if star == nil { - return errors.New("Invalid body value for a user-defined star.") + return errors.New("invalid body value for a user-defined star.") } if !isfinite(ra) || ra < 0.0 || ra >= 24.0 { - return errors.New("Invalid right ascension for user-defined star.") + return errors.New("invalid right ascension for user-defined star.") } if !isfinite(dec) || dec < -90.0 || dec > +90.0 { - return errors.New("Invalid declination for user-defined star.") + return errors.New("invalid declination for user-defined star.") } if !isfinite(distanceLightYears) || distanceLightYears < 1.0 { - return errors.New("Invalid heliocentric distance for user-defined star. Must be at least 1 light year.") + return errors.New("invalid heliocentric distance for user-defined star. Must be at least 1 light year.") } star.ra = ra star.dec = dec @@ -1709,10 +1817,10 @@ func CombineRotation(a, b RotationMatrix) RotationMatrix { // in the range [-360, +360]. func Pivot(rotation RotationMatrix, axis int, angle float64) (*RotationMatrix, error) { if axis < 0 || axis > 2 { - return nil, errors.New("Invalid coordinate axis for Pivot. Must be 0..2.") + return nil, errors.New("invalid coordinate axis for Pivot. Must be 0..2.") } if !isfinite(angle) { - return nil, errors.New("Angle is not a finite number. Not valid for Pivot.") + return nil, errors.New("angle is not a finite number. Not valid for Pivot.") } c := dcos(angle) s := dsin(angle) @@ -2459,7 +2567,7 @@ func visualMagnitude(body Body, phase, helioDist, geoDist float64) (float64, err c1 = +4.00 break default: - return 0.0, errors.New("Invalid body for visual magnitude") + return 0.0, errors.New("invalid body for visual magnitude") } x := phase / 100.0 mag := c0 + x*(c1+x*(c2+x*c3)) diff --git a/source/golang/README.md b/source/golang/README.md index 148f9e49..bb367a09 100644 --- a/source/golang/README.md +++ b/source/golang/README.md @@ -76,6 +76,7 @@ It provides a suite of well\-tested functions for calculating positions of the S - [func \(state StateVector\) Velocity\(\) AstroVector](<#StateVector.Velocity>) - [type TimeFormat](<#TimeFormat>) - [type Topocentric](<#Topocentric>) + - [func Horizon\(time AstroTime, observer Observer, ra, dec float64, refraction Refraction\) \(\*Topocentric, error\)](<#Horizon>) ## Constants @@ -172,7 +173,7 @@ const ( ``` -## func [AngleBetween]() +## func [AngleBetween]() ```go func AngleBetween(avec AstroVector, bvec AstroVector) float64 @@ -190,7 +191,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 @@ -199,7 +200,7 @@ func DefineStar(body Body, ra, dec, distanceLightYears float64) error -## func [DegreesFromRadians]() +## func [DegreesFromRadians]() ```go func DegreesFromRadians(radians float64) float64 @@ -217,7 +218,7 @@ func Dot(a, b AstroVector) float64 Returns the scalar dot product of two vectors. -## func [MassProduct]() +## func [MassProduct]() ```go func MassProduct(body Body) float64 @@ -226,7 +227,7 @@ func MassProduct(body Body) float64 Returns the product of mass and universal gravitational constant of a Solar System body. For problems involving the gravitational interactions of Solar System bodies, it is helpful to know the product GM, where G = the universal gravitational constant and M = the mass of the body. In practice, GM is known to a higher precision than either G or M alone, and thus using the product results in the most accurate results. This function returns the product GM in the units au^3/day^2. The values come from page 10 of a JPL memorandum regarding the DE405/LE405 ephemeris: https://web.archive.org/web/20120220062549/http://iau-comm4.jpl.nasa.gov/de405iom/de405iom.pdf -## func [ObserverGravity]() +## func [ObserverGravity]() ```go func ObserverGravity(latitude, height float64) float64 @@ -235,7 +236,7 @@ func ObserverGravity(latitude, height float64) float64 Calculates the gravitational acceleration experienced by an observer on the Earth. This function implements the WGS 84 Ellipsoidal Gravity Formula. The result is a combination of inward gravitational acceleration with outward centrifugal acceleration, as experienced by an observer in the Earth's rotating frame of reference. The resulting value increases toward the Earth's poles and decreases toward the equator, consistent with changes of the weight measured by a spring scale of a fixed mass moved to different latitudes and heights on the Earth. The latitude is of the observer in degrees north or south of the equator. By formula symmetry, positive latitudes give the same answer as negative latitudes, so the sign does not matter. The height is specified above the sea level geoid in meters. No range checking is done; however, accuracy is only valid in the range 0 to 100000 meters. The return value is the gravitational acceleration expressed in meters per second squared. -## func [PlanetOrbitalPeriod]() +## func [PlanetOrbitalPeriod]() ```go func PlanetOrbitalPeriod(body Body) float64 @@ -244,7 +245,7 @@ func PlanetOrbitalPeriod(body Body) float64 PlanetOrbitalPeriod returns the average number of days it takes for a planet to orbit the Sun. -## func [RadiansFromDegrees]() +## func [RadiansFromDegrees]() ```go func RadiansFromDegrees(degrees float64) float64 @@ -262,7 +263,7 @@ func RefractionAngle(refraction Refraction, altitude float64) float64 RefractionAngle calculates the amount of "lift" to an altitude angle caused by atmospheric refraction. Given an altitude angle and a refraction option, calculates the amount of "lift" caused by atmospheric refraction. This is the number of degrees higher in the sky an object appears due to the lensing of the Earth's atmosphere. This function works best near sea level. To correct for higher elevations, call Atmosphere for that elevation and multiply the refraction angle by the resulting relative density. The refraction parameter specifies which refraction correction to use. If set to NormalRefraction, uses a well\-behaved refraction model that works well for all valid values \(\-90 to \+90\) of altitude. If set to JplHorizonsRefraction, this function returns a value compatible with the JPL Horizons tool. This is provided for internal unit tests that compare against JPL Horizons data. Any other value, including NoRefraction, causes this function to return 0.0. The return value is a non\-negative value expressed in degrees of refraction above the horizontal. -## func [SiderealTime]() +## func [SiderealTime]() ```go func SiderealTime(time *AstroTime) float64 @@ -407,7 +408,7 @@ For common use cases, it is simpler to use BackdatePosition for calculating the For geocentric calculations, GeoVector also backdates the returned position vector for light travel time, only it returns the observation time in the returned vector's \`t\` field rather than the backdated time. -### func [GeoMoon]() +### func [GeoMoon]() ```go func GeoMoon(time AstroTime) AstroVector @@ -416,7 +417,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 @@ -581,7 +582,7 @@ type JupiterMoonsInfo struct { ``` -### func [JupiterMoons]() +### func [JupiterMoons]() ```go func JupiterMoons(time AstroTime) JupiterMoonsInfo @@ -681,7 +682,7 @@ type RotationMatrix struct { ``` -### func [CombineRotation]() +### func [CombineRotation]() ```go func CombineRotation(a, b RotationMatrix) RotationMatrix @@ -708,7 +709,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 [Pivot]() +### func [Pivot]() ```go func Pivot(rotation RotationMatrix, axis int, angle float64) (*RotationMatrix, error) @@ -717,7 +718,7 @@ func Pivot(rotation RotationMatrix, axis int, angle float64) (*RotationMatrix, e Pivot re\-orients a rotation matrix by pivoting it by an angle around one of its axes. Given a rotation matrix, a selected coordinate axis, and an angle in degrees, this function pivots the rotation matrix by that angle around that coordinate axis. For example, if you have rotation matrix that converts ecliptic coordinates \(ECL\) to horizontal coordinates \(HOR\), but you really want to convert ECL to the orientation of a telescope camera pointed at a given body, you can use Pivot twice: \(1\) pivot around the zenith axis by the body's azimuth, then \(2\) pivot around the western axis by the body's altitude angle. The resulting rotation matrix will then reorient ECL coordinates to the orientation of your telescope camera. The axis parameter is an integer that selects which axis to pivot about: 0=x, 1=y, 2=z. The angle parameter is an angle in degrees indicating the amount of rotation around the specified axis. Positive angles indicate rotation counterclockwise as seen from the positive direction along that axis, looking towards the origin point of the orientation system. Any finite number of degrees is allowed, but best precision will result from keeping angle in the range \[\-360, \+360\]. -### func [RotationEclEqj]() +### func [RotationEclEqj]() ```go func RotationEclEqj() RotationMatrix @@ -726,7 +727,7 @@ func RotationEclEqj() RotationMatrix Calculates a rotation matrix from J2000 mean ecliptic \(ECL\) to J2000 mean equator \(EQJ\). -### func [RotationEqdEqj]() +### func [RotationEqdEqj]() ```go func RotationEqdEqj(time *AstroTime) RotationMatrix @@ -735,7 +736,7 @@ func RotationEqdEqj(time *AstroTime) RotationMatrix Calculates a rotation matrix that converts equator\-of\-date \(EQD\) to J2000 mean equator \(EQJ\). -### func [RotationEqjEcl]() +### func [RotationEqjEcl]() ```go func RotationEqjEcl() RotationMatrix @@ -744,7 +745,7 @@ func RotationEqjEcl() RotationMatrix Calculates a rotation matrix from J2000 mean equator \(EQJ\) to J2000 mean ecliptic \(ECL\). -### func [RotationEqjGal]() +### func [RotationEqjGal]() ```go func RotationEqjGal() RotationMatrix @@ -753,7 +754,7 @@ func RotationEqjGal() RotationMatrix Calculates a rotation matrix from J2000 mean equator \(EQJ\) to galactic \(GAL\). -### func [RotationGalEqj]() +### func [RotationGalEqj]() ```go func RotationGalEqj() RotationMatrix @@ -762,7 +763,7 @@ func RotationGalEqj() RotationMatrix -## type [SearchContext]() +## type [SearchContext]() @@ -826,7 +827,7 @@ type StateVector struct { ``` -### func [RotateState]() +### func [RotateState]() ```go func RotateState(rotation RotationMatrix, state StateVector) StateVector @@ -886,4 +887,13 @@ type Topocentric struct { } ``` + +### func [Horizon]() + +```go +func Horizon(time AstroTime, observer Observer, ra, dec float64, refraction Refraction) (*Topocentric, error) +``` + + + Generated by [gomarkdoc]() diff --git a/source/golang/astronomy.go b/source/golang/astronomy.go index eb1ebefc..88be6633 100644 --- a/source/golang/astronomy.go +++ b/source/golang/astronomy.go @@ -310,7 +310,7 @@ type CalendarDateTime struct { // CalendarFromDays converts a J2000 day value to a Gregorian calendar date and time. func CalendarFromDays(ut float64) (*CalendarDateTime, error) { if !isfinite(ut) { - return nil, errors.New("Non-finite value of ut") + return nil, errors.New("non-finite value of ut") } // Adapted from the NOVAS C 3.1 function cal_date(). @@ -352,11 +352,11 @@ func CalendarFromDays(ut float64) (*CalendarDateTime, error) { year := int(100*(n-49) + m + k - 400*c) if year < -999999 || year > +999999 { - return nil, errors.New("The supplied time is too far from the year 2000 to be represented.") + return nil, errors.New("the supplied time is too far from the year 2000 to be represented.") } if month < 1 || month > 12 || day < 1 || day > 31 { - return nil, errors.New("Internal error: invalid calendar date calculated.") + return nil, errors.New("internal error: invalid calendar date calculated.") } return &CalendarDateTime{year, month, day, hour, minute, second}, nil @@ -396,7 +396,7 @@ func TimeFromCalendar(year, month, day, hour, minute int, second float64) AstroT // https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf func Atmosphere(elevationMeters float64) (AtmosphereInfo, error) { if !isfinite(elevationMeters) || elevationMeters < -500.0 || elevationMeters > 100000.0 { - return AtmosphereInfo{}, errors.New("Invalid elevation") + return AtmosphereInfo{}, errors.New("invalid elevation") } const P0 = 101325.0 // pressure at sea level [pascals] const T0 = 288.15 // temperature at sea level [kelvins] @@ -1034,7 +1034,7 @@ func CorrectLightTravel(posFunc PositionFunction, time AstroTime) (*AstroVector, // values for stellar aberration angles. lt := pos.Length() / SpeedOfLightAuPerDay if lt > 1.0 { - return nil, errors.New("Object is too distant for light-travel solver.") + return nil, errors.New("object is too distant for light-travel solver.") } ltime2 := time.AddDays(-lt) dt := math.Abs(ltime2.Tt - ltime.Tt) @@ -1043,7 +1043,115 @@ func CorrectLightTravel(posFunc PositionFunction, time AstroTime) (*AstroVector, } ltime = ltime2 } - return nil, errors.New("Light travel time correction did not converge.") + return nil, errors.New("light travel time correction did not converge.") +} + +func Horizon(time AstroTime, observer Observer, ra, dec float64, refraction Refraction) (*Topocentric, error) { + sinlat := dsin(observer.Latitude) + coslat := dcos(observer.Latitude) + sinlon := dsin(observer.Longitude) + coslon := dcos(observer.Longitude) + sindc := dsin(dec) + cosdc := dcos(dec) + sinra := dsin(ra * 15.0) + cosra := dcos(ra * 15.0) + + // Calculate three mutually perpendicular unit vectors + // in equatorial coordinates: uze, une, uwe. + // + // uze = The direction of the observer's local zenith (straight up). + // une = The direction toward due north on the observer's horizon. + // uwe = The direction toward due west on the observer's horizon. + // + // HOWEVER, these are uncorrected for the Earth's rotation due to the time of day. + // + // The components of these 3 vectors are as follows: + // x = direction from center of Earth toward 0 degrees longitude (the prime meridian) on equator. + // y = direction from center of Earth toward 90 degrees west longitude on equator. + // z = direction from center of Earth toward the north pole. + uze := AstroVector{coslat * coslon, coslat * sinlon, sinlat, time} + une := AstroVector{-sinlat * coslon, -sinlat * sinlon, coslat, time} + uwe := AstroVector{sinlon, -coslon, 0.0, time} + + // Correct the vectors uze, une, uwe for the Earth's rotation by calculating + // sidereal time. Call spin() for each uncorrected vector to rotate about + // the Earth's axis to yield corrected unit vectors uz, un, uw. + // Multiply sidereal hours by -15 to convert to degrees and flip eastward + // rotation of the Earth to westward apparent movement of objects with time. + angle := -15.0 * SiderealTime(&time) + uz := spin(angle, uze) + un := spin(angle, une) + uw := spin(angle, uwe) + + // Convert angular equatorial coordinates (RA, DEC) to + // cartesian equatorial coordinates in 'p', using the + // same orientation system as uze, une, uwe. + p := AstroVector{cosdc * cosra, cosdc * sinra, sindc, time} + + // Use dot products of p with the zenith, north, and west + // vectors to obtain the cartesian coordinates of the body in + // the observer's horizontal orientation system. + // pz = zenith component [-1, +1] + // pn = north component [-1, +1] + // pw = west component [-1, +1] + pz := p.X*uz.X + p.Y*uz.Y + p.Z*uz.Z + pn := p.X*un.X + p.Y*un.Y + p.Z*un.Z + pw := p.X*uw.X + p.Y*uw.Y + p.Z*uw.Z + + // proj is the "shadow" of the body vector along the observer's flat ground. + proj := math.Hypot(pn, pw) + + // Calculate az = azimuth (compass direction clockwise from East.) + var az float64 + if proj > 0.0 { + // If the body is not exactly straight up/down, it has an azimuth. + // Invert the angle to produce degrees eastward from north. + az = -datan2(pw, pn) + if az < 0.0 { + az += 360.0 + } + } else { + // The body is straight up/down, so it does not have an azimuth. + // Report an arbitrary but reasonable value. + az = 0.0 + } + + // zd = the angle of the body away from the observer's zenith, in degrees. + zd := datan2(proj, pz) + horRa := ra + horDec := dec + + if refraction == NormalRefraction || refraction == JplHorizonsRefraction { + zd0 := zd + refr := RefractionAngle(refraction, 90.0-zd) + zd -= refr + + if refr > 0.0 && zd > 3.0e-4 { + sinzd := dsin(zd) + coszd := dcos(zd) + sinzd0 := dsin(zd0) + coszd0 := dcos(zd0) + + prx := ((p.X-coszd0*uz.X)/sinzd0)*sinzd + uz.X*coszd + pry := ((p.Y-coszd0*uz.Y)/sinzd0)*sinzd + uz.Y*coszd + prz := ((p.Z-coszd0*uz.Z)/sinzd0)*sinzd + uz.Z*coszd + + proj = math.Hypot(prx, pry) + if proj > 0.0 { + horRa = datan2(pry, prx) / 15.0 + if horRa < 0.0 { + horRa += 24.0 + } + } else { + horRa = 0.0 + } + horDec = datan2(prz, proj) + } + } else if refraction != NoRefraction { + return nil, errors.New("unsupported refraction option.") + } + + return &Topocentric{az, 90.0 - zd, horRa, horDec}, nil } type jupiterMoon struct { @@ -1298,16 +1406,16 @@ func userDefinedStar(body Body) *starDef { func DefineStar(body Body, ra, dec, distanceLightYears float64) error { star := getStar(body) if star == nil { - return errors.New("Invalid body value for a user-defined star.") + return errors.New("invalid body value for a user-defined star.") } if !isfinite(ra) || ra < 0.0 || ra >= 24.0 { - return errors.New("Invalid right ascension for user-defined star.") + return errors.New("invalid right ascension for user-defined star.") } if !isfinite(dec) || dec < -90.0 || dec > +90.0 { - return errors.New("Invalid declination for user-defined star.") + return errors.New("invalid declination for user-defined star.") } if !isfinite(distanceLightYears) || distanceLightYears < 1.0 { - return errors.New("Invalid heliocentric distance for user-defined star. Must be at least 1 light year.") + return errors.New("invalid heliocentric distance for user-defined star. Must be at least 1 light year.") } star.ra = ra star.dec = dec @@ -1709,10 +1817,10 @@ func CombineRotation(a, b RotationMatrix) RotationMatrix { // in the range [-360, +360]. func Pivot(rotation RotationMatrix, axis int, angle float64) (*RotationMatrix, error) { if axis < 0 || axis > 2 { - return nil, errors.New("Invalid coordinate axis for Pivot. Must be 0..2.") + return nil, errors.New("invalid coordinate axis for Pivot. Must be 0..2.") } if !isfinite(angle) { - return nil, errors.New("Angle is not a finite number. Not valid for Pivot.") + return nil, errors.New("angle is not a finite number. Not valid for Pivot.") } c := dcos(angle) s := dsin(angle) @@ -2459,7 +2567,7 @@ func visualMagnitude(body Body, phase, helioDist, geoDist float64) (float64, err c1 = +4.00 break default: - return 0.0, errors.New("Invalid body for visual magnitude") + return 0.0, errors.New("invalid body for visual magnitude") } x := phase / 100.0 mag := c0 + x*(c1+x*(c2+x*c3)) diff --git a/source/golang/astronomy_test.go b/source/golang/astronomy_test.go index 5f7555e6..c766ee27 100644 --- a/source/golang/astronomy_test.go +++ b/source/golang/astronomy_test.go @@ -394,3 +394,46 @@ func TestRefraction(t *testing.T) { refractionCheck(t, NormalRefraction, -80.0, 0.0726495079641601) refractionCheck(t, NormalRefraction, -90.0, 0.0) } + +func TestHorizon(t *testing.T) { + // Test case generated using the Python version of Astronomy Engine. + // 2023-11-01T18:07:13.348Z + time := TimeFromCalendar(2023, 11, 1, 18, 7, 13.348) + ttExpected := 8705.255869201963 + if diff := math.Abs(ttExpected - time.Tt); diff > 1.0e-16 { + t.Errorf("excessive tt diff = %g", diff) + } + utExpected := 8705.255015601851 + if diff := math.Abs(utExpected - time.Ut); diff > 1.0e-16 { + t.Errorf("excessive ut diff = %g", diff) + } + + observer := Observer{Latitude: 35.0, Longitude: -87.0, Height: 0.0} + ra := 12.0 + dec := 5.0 + topo, err := Horizon(time, observer, ra, dec, NormalRefraction) + if err != nil { + t.Fatal(err) + } else { + azimuthExpected := 245.2118420327564 + altitudeExpected := 38.42409530728063 + raExpected := 12.001062402159066 + decExpected := 5.014149039382782 + + if diff := math.Abs(topo.Azimuth - azimuthExpected); diff > 1.0e-13 { + t.Errorf("excessive azimuth diff = %g, calculated = %g", diff, topo.Azimuth) + } + + if diff := math.Abs(topo.Altitude - altitudeExpected); diff > 3.0e-14 { + t.Errorf("excessive altitude diff = %g, calculated = %g", diff, topo.Altitude) + } + + if diff := math.Abs(topo.Ra - raExpected); diff > 1.0e-16 { + t.Errorf("excessive RA diff = %g, calculated = %g", diff, topo.Ra) + } + + if diff := math.Abs(topo.Dec - decExpected); diff > 1.0e-16 { + t.Errorf("excessive DEC diff = %g, calculated = %g", diff, topo.Dec) + } + } +}