diff --git a/README.md b/README.md index 156f6e7c..12761a8a 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Examples -
+
Code & Docs
diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 42f37609..6b8568fe 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -44,6 +44,7 @@ def _cbrt(x): KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. +AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. # Jupiter radius data are nominal values obtained from: # https://www.iau.org/static/resolutions/IAU2015_English.pdf @@ -223,6 +224,9 @@ class Vector: def __neg__(self): return Vector(-self.x, -self.y, -self.z, self.t) + def __truediv__(self, scalar): + return Vector(self.x/scalar, self.y/scalar, self.z/scalar, self.t) + def format(self, coord_format): """Returns a custom format string representation of the vector.""" layout = '({:' + coord_format + '}, {:' + coord_format + '}, {:' + coord_format + '}, {})' @@ -290,6 +294,14 @@ class StateVector: self.t ) + def Position(self): + """Extracts a position vector from this state vector.""" + return Vector(self.x, self.y, self.z, self.t) + + def Velocity(self): + """Extracts a velocity vector from this state vector.""" + return Vector(self.vx, self.vy, self.vz, self.t) + @enum.unique class Body(enum.Enum): """The celestial bodies supported by Astronomy Engine calculations. @@ -310,6 +322,14 @@ class Body(enum.Enum): Moon: The Earth's moon. EMB: The Earth/Moon Barycenter. SSB: The Solar System Barycenter. + Star1: User-defined star 1. + Star2: User-defined star 2. + Star3: User-defined star 3. + Star4: User-defined star 4. + Star5: User-defined star 5. + Star6: User-defined star 6. + Star7: User-defined star 7. + Star8: User-defined star 8. """ Invalid = -1 Mercury = 0 @@ -325,6 +345,78 @@ class Body(enum.Enum): Moon = 10 EMB = 11 SSB = 12 + Star1 = 101 + Star2 = 102 + Star3 = 103 + Star4 = 104 + Star5 = 105 + Star6 = 106 + Star7 = 107 + Star8 = 108 + +class _StarDef: + def __init__(self, ra, dec, dist): + self.ra = ra + self.dec = dec + self.dist = dist + +_StarTable = [ + _StarDef(0.0, 0.0, AU_PER_LY), # Star1 + _StarDef(0.0, 0.0, AU_PER_LY), # Star2 + _StarDef(0.0, 0.0, AU_PER_LY), # Star3 + _StarDef(0.0, 0.0, AU_PER_LY), # Star4 + _StarDef(0.0, 0.0, AU_PER_LY), # Star5 + _StarDef(0.0, 0.0, AU_PER_LY), # Star6 + _StarDef(0.0, 0.0, AU_PER_LY), # Star7 + _StarDef(0.0, 0.0, AU_PER_LY), # Star8 +] + +def _IsUserDefinedStar(body): + return Body.Star1.value <= body.value <= Body.Star8.value + +def _GetUserDefinedStar(body): + return _StarTable[body.value - Body.Star1.value] if _IsUserDefinedStar(body) else None + +def DefineStar(body, ra, dec, distanceLightYears): + """Assign equatorial coordinates to a user-defined star. + + Some Astronomy Engine functions allow their `body` parameter to + be a user-defined fixed point in the sky, loosely called a "star". + This function assigns a right ascension, declination, and distance + to one of the eight user-defined stars `Body.Star1`..`Body.Star8`. + + A star that has not been defined through a call to `DefineStar` + defaults to the coordinates RA=0, DEC=0 and a heliocentric distance of 1 light-year. + Once defined, the star keeps the given coordinates until + a subsequent call to `DefineStar` replaces the coordinates with new values. + + Parameters + ---------- + body: Body + One of the eight user-defined star identifiers: `Body.Star1`, `Body.Star2`, ..., `Body.Star8`. + ra: float + The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). + The value is in units of sidereal hours, and must be within the half-open range [0, 24). + dec: float + The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). + The value is in units of degrees north (positive) or south (negative) of the J2000 equator, + and must be within the closed range [-90, +90]. + distanceLightYears: float + The distance between the star and the Sun, expressed in light-years. + This value is used to calculate the tiny parallax shift as seen by an observer on Earth. + If you don't know the distance to the star, using a large value like 1000 will generally work well. + The minimum allowed distance is 1 light-year, which is required to provide certain internal optimizations. + """ + star = _GetUserDefinedStar(body) + if star is None: + raise Error('Body must be a user-defined star, not {}.'.format(body)) + if not (0.0 <= ra < 24.0): + raise Error('Invalid right ascension: {}'.format(ra)) + if not (-90.0 <= dec <= +90.0): + raise Error('Invalid declination: {}'.format(dec)) + star.ra = ra + star.dec = dec + star.dist = distanceLightYears * AU_PER_LY def BodyCode(name): """Finds the Body enumeration value, given the name of a body. @@ -4444,6 +4536,10 @@ def HelioVector(body, time): if body == Body.SSB: return _CalcSolarSystemBarycenter(time) + star = _GetUserDefinedStar(body) + if star is not None: + return VectorFromSphere(Spherical(star.dec, 15.0*star.ra, star.dist), time) + raise InvalidBodyError(body) @@ -4630,6 +4726,24 @@ def BackdatePosition(time, observerBody, targetBody, aberration): Its `t` field holds the time that light left the observed body to arrive at the observer at the observation time. """ + if _IsUserDefinedStar(targetBody): + # This is a user-defined star, which must be treated as a special case. + # First, we assume its heliocentric position does not change with time. + # Second, we assume its heliocentric position has already been corrected + # for light-travel time, its coordinates given as it appears on Earth at the present. + # Therefore, no backdating is applied. + tvec = HelioVector(targetBody, time) + if aberration: + # (Observer velocity) - (light vector) = (Aberration-corrected direction to target body). + # Note that this is an approximation, because technically the light vector should + # be measured in barycentric coordinates, not heliocentric. The error is very small. + ostate = HelioState(observerBody, time) + rvec = tvec - ostate.Position() + s = C_AUDAY / rvec.Length() # conversion factor from relative distance to speed of light + return rvec + ostate.Velocity()/s + # No correction is needed. Simply return the star's current position as seen from the observer. + return tvec - HelioVector(observerBody, time) + if aberration: # With aberration, `BackdatePosition` will calculate `observerPos` at different times. # Therefore, do not waste time calculating it now. @@ -6548,6 +6662,13 @@ def _MaxAltitudeSlope(body, latitude): elif body in [Body.Jupiter, Body.Saturn, Body.Uranus, Body.Neptune, Body.Pluto]: deriv_ra = -0.2 deriv_dec = +0.2 + elif _IsUserDefinedStar(body): + # The minimum allowed heliocentric distance of a user-defined star + # is one light-year. This can cause a tiny amount of parallax (about 0.001 degrees). + # Also, including stellar aberration (22 arcsec = 0.006 degrees), we provide a + # generous safety buffer of 0.008 degrees. + deriv_ra = -0.008 + deriv_dec = +0.008 elif body == Body.Earth: raise EarthNotAllowedError() else: diff --git a/generate/ctest.c b/generate/ctest.c index 83fbebd3..437dbb15 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -7115,11 +7115,11 @@ static int StarRiseSetCulmCase( cdiff = MINUTES_PER_DAY * ABS(expectedCulmTime.ut - culm.time.ut); sdiff = MINUTES_PER_DAY * ABS(expectedSetTime.ut - set.time.ut); - DEBUG("StarRiseSetCulmCase(%s): minutes rdiff = %0.4lf, cdiff = %0.4lf, sdiff = %0.4lf\n", starName, rdiff, cdiff, sdiff); + DEBUG("C StarRiseSetCulmCase(%s): minutes rdiff = %0.4lf, cdiff = %0.4lf, sdiff = %0.4lf\n", starName, rdiff, cdiff, sdiff); - if (rdiff > 0.5) FAIL("StarRiseSetCulmCase(%s): exccessive rise time error = %0.4lf minutes.\n", starName, rdiff); - if (cdiff > 0.5) FAIL("StarRiseSetCulmCase(%s): exccessive culm time error = %0.4lf minutes.\n", starName, cdiff); - if (sdiff > 0.5) FAIL("StarRiseSetCulmCase(%s): exccessive set time error = %0.4lf minutes.\n", starName, sdiff); + if (rdiff > 0.5) FAIL("C StarRiseSetCulmCase(%s): exccessive rise time error = %0.4lf minutes.\n", starName, rdiff); + if (cdiff > 0.5) FAIL("C StarRiseSetCulmCase(%s): exccessive culm time error = %0.4lf minutes.\n", starName, cdiff); + if (sdiff > 0.5) FAIL("C StarRiseSetCulmCase(%s): exccessive set time error = %0.4lf minutes.\n", starName, sdiff); error = 0; fail: @@ -7136,7 +7136,7 @@ static int StarRiseSetCulm(void) CHECK(StarRiseSetCulmCase("Canopus", 6.3992, -52.6956, 310.0, observer, 2022, 11, 21, 4, 17, 7, 44, 11, 11)); CHECK(StarRiseSetCulmCase("Canopus", 6.3992, -52.6956, 310.0, observer, 2022, 11, 25, 4, 1, 7, 28, 10, 56)); - printf("StarRiseSetCulm: PASS\n"); + printf("C StarRiseSetCulm: PASS\n"); fail: return error; } diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs index 0330c63c..40914edb 100644 --- a/generate/dotnet/csharp_test/csharp_test.cs +++ b/generate/dotnet/csharp_test/csharp_test.cs @@ -24,20 +24,19 @@ namespace csharp_test static int Pass(string name) { - Console.WriteLine(name + ": PASS"); + Console.WriteLine("C# " + name + ": PASS"); return 0; } static int Fail(string message) { - Console.WriteLine(message); + Console.WriteLine("C# " + message); return 1; } static bool BoolFail(string message) { - Console.WriteLine(message); - return false; + return 0 == Fail(message); } struct Test diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 8fe7daa4..cd4544e9 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -44,6 +44,7 @@ def _cbrt(x): KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. +AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. # Jupiter radius data are nominal values obtained from: # https://www.iau.org/static/resolutions/IAU2015_English.pdf @@ -223,6 +224,9 @@ class Vector: def __neg__(self): return Vector(-self.x, -self.y, -self.z, self.t) + def __truediv__(self, scalar): + return Vector(self.x/scalar, self.y/scalar, self.z/scalar, self.t) + def format(self, coord_format): """Returns a custom format string representation of the vector.""" layout = '({:' + coord_format + '}, {:' + coord_format + '}, {:' + coord_format + '}, {})' @@ -290,6 +294,14 @@ class StateVector: self.t ) + def Position(self): + """Extracts a position vector from this state vector.""" + return Vector(self.x, self.y, self.z, self.t) + + def Velocity(self): + """Extracts a velocity vector from this state vector.""" + return Vector(self.vx, self.vy, self.vz, self.t) + @enum.unique class Body(enum.Enum): """The celestial bodies supported by Astronomy Engine calculations. @@ -310,6 +322,14 @@ class Body(enum.Enum): Moon: The Earth's moon. EMB: The Earth/Moon Barycenter. SSB: The Solar System Barycenter. + Star1: User-defined star 1. + Star2: User-defined star 2. + Star3: User-defined star 3. + Star4: User-defined star 4. + Star5: User-defined star 5. + Star6: User-defined star 6. + Star7: User-defined star 7. + Star8: User-defined star 8. """ Invalid = -1 Mercury = 0 @@ -325,6 +345,78 @@ class Body(enum.Enum): Moon = 10 EMB = 11 SSB = 12 + Star1 = 101 + Star2 = 102 + Star3 = 103 + Star4 = 104 + Star5 = 105 + Star6 = 106 + Star7 = 107 + Star8 = 108 + +class _StarDef: + def __init__(self, ra, dec, dist): + self.ra = ra + self.dec = dec + self.dist = dist + +_StarTable = [ + _StarDef(0.0, 0.0, AU_PER_LY), # Star1 + _StarDef(0.0, 0.0, AU_PER_LY), # Star2 + _StarDef(0.0, 0.0, AU_PER_LY), # Star3 + _StarDef(0.0, 0.0, AU_PER_LY), # Star4 + _StarDef(0.0, 0.0, AU_PER_LY), # Star5 + _StarDef(0.0, 0.0, AU_PER_LY), # Star6 + _StarDef(0.0, 0.0, AU_PER_LY), # Star7 + _StarDef(0.0, 0.0, AU_PER_LY), # Star8 +] + +def _IsUserDefinedStar(body): + return Body.Star1.value <= body.value <= Body.Star8.value + +def _GetUserDefinedStar(body): + return _StarTable[body.value - Body.Star1.value] if _IsUserDefinedStar(body) else None + +def DefineStar(body, ra, dec, distanceLightYears): + """Assign equatorial coordinates to a user-defined star. + + Some Astronomy Engine functions allow their `body` parameter to + be a user-defined fixed point in the sky, loosely called a "star". + This function assigns a right ascension, declination, and distance + to one of the eight user-defined stars `Body.Star1`..`Body.Star8`. + + A star that has not been defined through a call to `DefineStar` + defaults to the coordinates RA=0, DEC=0 and a heliocentric distance of 1 light-year. + Once defined, the star keeps the given coordinates until + a subsequent call to `DefineStar` replaces the coordinates with new values. + + Parameters + ---------- + body: Body + One of the eight user-defined star identifiers: `Body.Star1`, `Body.Star2`, ..., `Body.Star8`. + ra: float + The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). + The value is in units of sidereal hours, and must be within the half-open range [0, 24). + dec: float + The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). + The value is in units of degrees north (positive) or south (negative) of the J2000 equator, + and must be within the closed range [-90, +90]. + distanceLightYears: float + The distance between the star and the Sun, expressed in light-years. + This value is used to calculate the tiny parallax shift as seen by an observer on Earth. + If you don't know the distance to the star, using a large value like 1000 will generally work well. + The minimum allowed distance is 1 light-year, which is required to provide certain internal optimizations. + """ + star = _GetUserDefinedStar(body) + if star is None: + raise Error('Body must be a user-defined star, not {}.'.format(body)) + if not (0.0 <= ra < 24.0): + raise Error('Invalid right ascension: {}'.format(ra)) + if not (-90.0 <= dec <= +90.0): + raise Error('Invalid declination: {}'.format(dec)) + star.ra = ra + star.dec = dec + star.dist = distanceLightYears * AU_PER_LY def BodyCode(name): """Finds the Body enumeration value, given the name of a body. @@ -2402,6 +2494,10 @@ def HelioVector(body, time): if body == Body.SSB: return _CalcSolarSystemBarycenter(time) + star = _GetUserDefinedStar(body) + if star is not None: + return VectorFromSphere(Spherical(star.dec, 15.0*star.ra, star.dist), time) + raise InvalidBodyError(body) @@ -2588,6 +2684,24 @@ def BackdatePosition(time, observerBody, targetBody, aberration): Its `t` field holds the time that light left the observed body to arrive at the observer at the observation time. """ + if _IsUserDefinedStar(targetBody): + # This is a user-defined star, which must be treated as a special case. + # First, we assume its heliocentric position does not change with time. + # Second, we assume its heliocentric position has already been corrected + # for light-travel time, its coordinates given as it appears on Earth at the present. + # Therefore, no backdating is applied. + tvec = HelioVector(targetBody, time) + if aberration: + # (Observer velocity) - (light vector) = (Aberration-corrected direction to target body). + # Note that this is an approximation, because technically the light vector should + # be measured in barycentric coordinates, not heliocentric. The error is very small. + ostate = HelioState(observerBody, time) + rvec = tvec - ostate.Position() + s = C_AUDAY / rvec.Length() # conversion factor from relative distance to speed of light + return rvec + ostate.Velocity()/s + # No correction is needed. Simply return the star's current position as seen from the observer. + return tvec - HelioVector(observerBody, time) + if aberration: # With aberration, `BackdatePosition` will calculate `observerPos` at different times. # Therefore, do not waste time calculating it now. @@ -4506,6 +4620,13 @@ def _MaxAltitudeSlope(body, latitude): elif body in [Body.Jupiter, Body.Saturn, Body.Uranus, Body.Neptune, Body.Pluto]: deriv_ra = -0.2 deriv_dec = +0.2 + elif _IsUserDefinedStar(body): + # The minimum allowed heliocentric distance of a user-defined star + # is one light-year. This can cause a tiny amount of parallax (about 0.001 degrees). + # Also, including stellar aberration (22 arcsec = 0.006 degrees), we provide a + # generous safety buffer of 0.008 degrees. + deriv_ra = -0.008 + deriv_dec = +0.008 elif body == Body.Earth: raise EarthNotAllowedError() else: diff --git a/generate/test.py b/generate/test.py index 3a4e2b5a..cc124005 100755 --- a/generate/test.py +++ b/generate/test.py @@ -10,6 +10,7 @@ import astronomy Verbose = False SECONDS_PER_DAY = 86400.0 +MINUTES_PER_DAY = 1440.0 def Debug(text): if Verbose: @@ -3030,6 +3031,57 @@ def LunarFraction(): #----------------------------------------------------------------------------------------------------------- +def StarRiseSetCulmCase(starName, ra, dec, distLy, observer, year, month, day, riseHour, riseMinute, culmHour, culmMinute, setHour, setMinute): + func = 'StarRiseSetCulmCase({})'.format(starName) + + # Calculate expected event times. + expectedRiseTime = astronomy.Time.Make(year, month, day, riseHour, riseMinute, 0.0) + expectedCulmTime = astronomy.Time.Make(year, month, day, culmHour, culmMinute, 0.0) + expectedSetTime = astronomy.Time.Make(year, month, day, setHour, setMinute, 0.0) + + # Define a custom star object. + astronomy.DefineStar(astronomy.Body.Star1, ra, dec, distLy) + + # Use Astronomy Engine to search for event times. + searchTime = astronomy.Time.Make(year, month, day, 0, 0, 0.0) + + rise = astronomy.SearchRiseSet(astronomy.Body.Star1, observer, astronomy.Direction.Rise, searchTime, 1.0) + if rise is None: + return Fail(func, 'Star rise search failed.') + + culm = astronomy.SearchHourAngle(astronomy.Body.Star1, observer, 0.0, searchTime, +1) + if culm is None: + return Fail(func, 'Star culmination search failed.') + + set = astronomy.SearchRiseSet(astronomy.Body.Star1, observer, astronomy.Direction.Set, searchTime, 1.0) + if set is None: + return Fail(func, 'Star set search failed.') + + # Compare expected times with calculated times. + rdiff = MINUTES_PER_DAY * vabs(expectedRiseTime.ut - rise.ut) + cdiff = MINUTES_PER_DAY * vabs(expectedCulmTime.ut - culm.time.ut) + sdiff = MINUTES_PER_DAY * vabs(expectedSetTime.ut - set.ut) + + Debug("{}: minutes rdiff = {:0.4f}, cdiff = {:0.4f}, sdiff = {:0.4f}".format(func, rdiff, cdiff, sdiff)) + + if rdiff > 0.5: return Fail(func, "excessive rise time error = {:0.4f} minutes.".format(rdiff)) + if cdiff > 0.5: return Fail(func, "excessive culm time error = {:0.4f} minutes.".format(cdiff)) + if sdiff > 0.5: return Fail(func, "excessive set time error = {:0.4f} minutes.".format(sdiff)) + return 0 + + +def StarRiseSetCulm(): + observer = astronomy.Observer(+25.77, -80.19, 0.0) + return ( + StarRiseSetCulmCase("Sirius", 6.7525, -16.7183, 8.6, observer, 2022, 11, 21, 2, 37, 8, 6, 13, 34) or + StarRiseSetCulmCase("Sirius", 6.7525, -16.7183, 8.6, observer, 2022, 11, 25, 2, 22, 7, 50, 13, 18) or + StarRiseSetCulmCase("Canopus", 6.3992, -52.6956, 310.0, observer, 2022, 11, 21, 4, 17, 7, 44, 11, 11) or + StarRiseSetCulmCase("Canopus", 6.3992, -52.6956, 310.0, observer, 2022, 11, 25, 4, 1, 7, 28, 10, 56) or + Pass("StarRiseSetCulm") + ) + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, 'axis': Axis, @@ -3066,6 +3118,7 @@ UnitTests = { 'seasons187': SeasonsIssue187, 'sidereal': SiderealTime, 'solar_fraction': SolarFraction, + 'star_risesetculm': StarRiseSetCulm, 'time': AstroTime, 'topostate': TopoState, 'transit': Transit, diff --git a/generate/version.txt b/generate/version.txt index a39c0b78..348fc11e 100644 --- a/generate/version.txt +++ b/generate/version.txt @@ -1 +1 @@ -2.1.11 +2.1.12 diff --git a/source/csharp/astronomy.csproj b/source/csharp/astronomy.csproj index 12576986..9bfe71fd 100644 --- a/source/csharp/astronomy.csproj +++ b/source/csharp/astronomy.csproj @@ -4,7 +4,7 @@ true true CosineKitty.AstronomyEngine - 2.1.11 + 2.1.12 https://github.com/cosinekitty/astronomy Don Cross Astronomy Engine diff --git a/source/js/package.json b/source/js/package.json index c926fee3..cbc1273b 100644 --- a/source/js/package.json +++ b/source/js/package.json @@ -1,6 +1,6 @@ { "name": "astronomy-engine", - "version": "2.1.11", + "version": "2.1.12", "description": "Astronomy calculation for Sun, Moon, and planets.", "author": "Donald Cross", "license": "MIT", diff --git a/source/kotlin/README.md b/source/kotlin/README.md index bb4d3e27..c45489f4 100644 --- a/source/kotlin/README.md +++ b/source/kotlin/README.md @@ -25,7 +25,7 @@ allprojects { Now add the dependency: ```kotlin dependencies { - implementation("io.github.cosinekitty:astronomy:2.1.11") + implementation("io.github.cosinekitty:astronomy:2.1.12") } ``` diff --git a/source/kotlin/build.gradle.kts b/source/kotlin/build.gradle.kts index 3b4761df..79a831cc 100644 --- a/source/kotlin/build.gradle.kts +++ b/source/kotlin/build.gradle.kts @@ -6,7 +6,7 @@ plugins { } group = "io.github.cosinekitty" -version = "2.1.11" +version = "2.1.12" repositories { mavenCentral() diff --git a/source/python/README.md b/source/python/README.md index 882d1f0e..5013c609 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -198,6 +198,13 @@ However, because Python defines the [angular conversion functions](https://docs. --- + +### `AU_PER_LY = 63241.07708807546` + +**The number of astronomical units in one light-year.** + +--- + ### `CALLISTO_RADIUS_KM = 2410.3` @@ -881,6 +888,18 @@ The state vector also includes a time stamp. | `float` | `vz` | The z-component of the velocity, measured in AU/day. | | [`Time`](#Time) | `t` | The date and time at which the position and velocity vectors are valid. | +#### member functions + + +### StateVector.Position(self) + +Extracts a position vector from this state vector. + + +### StateVector.Velocity(self) + +Extracts a velocity vector from this state vector. + --- @@ -1086,6 +1105,14 @@ two cases applies to a particular apsis event. | `Moon` | The Earth's moon. | | `EMB` | The Earth/Moon Barycenter. | | `SSB` | The Solar System Barycenter. | +| `Star1` | User-defined star 1. | +| `Star2` | User-defined star 2. | +| `Star3` | User-defined star 3. | +| `Star4` | User-defined star 4. | +| `Star5` | User-defined star 5. | +| `Star6` | User-defined star 6. | +| `Star7` | User-defined star 7. | +| `Star8` | User-defined star 8. | --- @@ -1421,6 +1448,29 @@ body to arrive at the observer at the observation time. --- + +### DefineStar(body, ra, dec, distanceLightYears) + +**Assign equatorial coordinates to a user-defined star.** + +Some Astronomy Engine functions allow their `body` parameter to +be a user-defined fixed point in the sky, loosely called a "star". +This function assigns a right ascension, declination, and distance +to one of the eight user-defined stars `Body.Star1`..`Body.Star8`. +A star that has not been defined through a call to `DefineStar` +defaults to the coordinates RA=0, DEC=0 and a heliocentric distance of 1 light-year. +Once defined, the star keeps the given coordinates until +a subsequent call to `DefineStar` replaces the coordinates with new values. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Body`](#Body) | `body` | One of the eight user-defined star identifiers: `Body.Star1`, `Body.Star2`, ..., `Body.Star8`. | +| `float` | `ra` | The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). The value is in units of sidereal hours, and must be within the half-open range [0, 24). | +| `float` | `dec` | The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). The value is in units of degrees north (positive) or south (negative) of the J2000 equator, and must be within the closed range [-90, +90]. | +| `float` | `distanceLightYears` | The distance between the star and the Sun, expressed in light-years. This value is used to calculate the tiny parallax shift as seen by an observer on Earth. If you don't know the distance to the star, using a large value like 1000 will generally work well. The minimum allowed distance is 1 light-year, which is required to provide certain internal optimizations. | + +--- + ### DeltaT_EspenakMeeus(ut) diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index 42f37609..6b8568fe 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -44,6 +44,7 @@ def _cbrt(x): KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. +AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. # Jupiter radius data are nominal values obtained from: # https://www.iau.org/static/resolutions/IAU2015_English.pdf @@ -223,6 +224,9 @@ class Vector: def __neg__(self): return Vector(-self.x, -self.y, -self.z, self.t) + def __truediv__(self, scalar): + return Vector(self.x/scalar, self.y/scalar, self.z/scalar, self.t) + def format(self, coord_format): """Returns a custom format string representation of the vector.""" layout = '({:' + coord_format + '}, {:' + coord_format + '}, {:' + coord_format + '}, {})' @@ -290,6 +294,14 @@ class StateVector: self.t ) + def Position(self): + """Extracts a position vector from this state vector.""" + return Vector(self.x, self.y, self.z, self.t) + + def Velocity(self): + """Extracts a velocity vector from this state vector.""" + return Vector(self.vx, self.vy, self.vz, self.t) + @enum.unique class Body(enum.Enum): """The celestial bodies supported by Astronomy Engine calculations. @@ -310,6 +322,14 @@ class Body(enum.Enum): Moon: The Earth's moon. EMB: The Earth/Moon Barycenter. SSB: The Solar System Barycenter. + Star1: User-defined star 1. + Star2: User-defined star 2. + Star3: User-defined star 3. + Star4: User-defined star 4. + Star5: User-defined star 5. + Star6: User-defined star 6. + Star7: User-defined star 7. + Star8: User-defined star 8. """ Invalid = -1 Mercury = 0 @@ -325,6 +345,78 @@ class Body(enum.Enum): Moon = 10 EMB = 11 SSB = 12 + Star1 = 101 + Star2 = 102 + Star3 = 103 + Star4 = 104 + Star5 = 105 + Star6 = 106 + Star7 = 107 + Star8 = 108 + +class _StarDef: + def __init__(self, ra, dec, dist): + self.ra = ra + self.dec = dec + self.dist = dist + +_StarTable = [ + _StarDef(0.0, 0.0, AU_PER_LY), # Star1 + _StarDef(0.0, 0.0, AU_PER_LY), # Star2 + _StarDef(0.0, 0.0, AU_PER_LY), # Star3 + _StarDef(0.0, 0.0, AU_PER_LY), # Star4 + _StarDef(0.0, 0.0, AU_PER_LY), # Star5 + _StarDef(0.0, 0.0, AU_PER_LY), # Star6 + _StarDef(0.0, 0.0, AU_PER_LY), # Star7 + _StarDef(0.0, 0.0, AU_PER_LY), # Star8 +] + +def _IsUserDefinedStar(body): + return Body.Star1.value <= body.value <= Body.Star8.value + +def _GetUserDefinedStar(body): + return _StarTable[body.value - Body.Star1.value] if _IsUserDefinedStar(body) else None + +def DefineStar(body, ra, dec, distanceLightYears): + """Assign equatorial coordinates to a user-defined star. + + Some Astronomy Engine functions allow their `body` parameter to + be a user-defined fixed point in the sky, loosely called a "star". + This function assigns a right ascension, declination, and distance + to one of the eight user-defined stars `Body.Star1`..`Body.Star8`. + + A star that has not been defined through a call to `DefineStar` + defaults to the coordinates RA=0, DEC=0 and a heliocentric distance of 1 light-year. + Once defined, the star keeps the given coordinates until + a subsequent call to `DefineStar` replaces the coordinates with new values. + + Parameters + ---------- + body: Body + One of the eight user-defined star identifiers: `Body.Star1`, `Body.Star2`, ..., `Body.Star8`. + ra: float + The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). + The value is in units of sidereal hours, and must be within the half-open range [0, 24). + dec: float + The right ascension to be assigned to the star, expressed in J2000 equatorial coordinates (EQJ). + The value is in units of degrees north (positive) or south (negative) of the J2000 equator, + and must be within the closed range [-90, +90]. + distanceLightYears: float + The distance between the star and the Sun, expressed in light-years. + This value is used to calculate the tiny parallax shift as seen by an observer on Earth. + If you don't know the distance to the star, using a large value like 1000 will generally work well. + The minimum allowed distance is 1 light-year, which is required to provide certain internal optimizations. + """ + star = _GetUserDefinedStar(body) + if star is None: + raise Error('Body must be a user-defined star, not {}.'.format(body)) + if not (0.0 <= ra < 24.0): + raise Error('Invalid right ascension: {}'.format(ra)) + if not (-90.0 <= dec <= +90.0): + raise Error('Invalid declination: {}'.format(dec)) + star.ra = ra + star.dec = dec + star.dist = distanceLightYears * AU_PER_LY def BodyCode(name): """Finds the Body enumeration value, given the name of a body. @@ -4444,6 +4536,10 @@ def HelioVector(body, time): if body == Body.SSB: return _CalcSolarSystemBarycenter(time) + star = _GetUserDefinedStar(body) + if star is not None: + return VectorFromSphere(Spherical(star.dec, 15.0*star.ra, star.dist), time) + raise InvalidBodyError(body) @@ -4630,6 +4726,24 @@ def BackdatePosition(time, observerBody, targetBody, aberration): Its `t` field holds the time that light left the observed body to arrive at the observer at the observation time. """ + if _IsUserDefinedStar(targetBody): + # This is a user-defined star, which must be treated as a special case. + # First, we assume its heliocentric position does not change with time. + # Second, we assume its heliocentric position has already been corrected + # for light-travel time, its coordinates given as it appears on Earth at the present. + # Therefore, no backdating is applied. + tvec = HelioVector(targetBody, time) + if aberration: + # (Observer velocity) - (light vector) = (Aberration-corrected direction to target body). + # Note that this is an approximation, because technically the light vector should + # be measured in barycentric coordinates, not heliocentric. The error is very small. + ostate = HelioState(observerBody, time) + rvec = tvec - ostate.Position() + s = C_AUDAY / rvec.Length() # conversion factor from relative distance to speed of light + return rvec + ostate.Velocity()/s + # No correction is needed. Simply return the star's current position as seen from the observer. + return tvec - HelioVector(observerBody, time) + if aberration: # With aberration, `BackdatePosition` will calculate `observerPos` at different times. # Therefore, do not waste time calculating it now. @@ -6548,6 +6662,13 @@ def _MaxAltitudeSlope(body, latitude): elif body in [Body.Jupiter, Body.Saturn, Body.Uranus, Body.Neptune, Body.Pluto]: deriv_ra = -0.2 deriv_dec = +0.2 + elif _IsUserDefinedStar(body): + # The minimum allowed heliocentric distance of a user-defined star + # is one light-year. This can cause a tiny amount of parallax (about 0.001 degrees). + # Also, including stellar aberration (22 arcsec = 0.006 degrees), we provide a + # generous safety buffer of 0.008 degrees. + deriv_ra = -0.008 + deriv_dec = +0.008 elif body == Body.Earth: raise EarthNotAllowedError() else: diff --git a/source/python/setup.py b/source/python/setup.py index 7b2399d9..083c5fb2 100644 --- a/source/python/setup.py +++ b/source/python/setup.py @@ -6,7 +6,7 @@ def _LoadFile(filename): setup( name='astronomy-engine', - version='2.1.11', + version='2.1.12', description='Astronomy calculation for Sun, Moon, and planets.', long_description=_LoadFile('README.md'), long_description_content_type='text/markdown',