Add AASHyperbolic class

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Changelog
 
+## 2.13.0 (2025-12-07)
+- Addition of a new CAAHyperbolicObjectElements class.
+- Addition of a new CAAHyperbolicObjectDetails class.
+- Addition of a new CAAHyperbolic class.
+
 ## 2.12.1 (2024-12-18)
 - Self documentation (in [PR#19](https://github.com/jsauve/AASharp/pull/19))
 - Extension methods for easy date/time manipulation

Src/AASharp/AASHyperbolic.cs

@@ -0,0 +1,135 @@
+using System;
+
+namespace AASharp {
+    /// <summary>
+    /// This class provides for calculation of the position of an object in an orbit which is hyperbolic i.e. with an eccentricity greater than one. One example being the 2025 interstellar comet 3I/ATLAS.
+    /// </summary>
+    public static class AASHyperbolic {
+        private static double CalculateKeplers(double M, double e, double epsilon = 0.000001) {
+            var bRecalc = true;
+            var Hk = M;
+            while (bRecalc) {
+                var Hk1 = Hk - ((e * Math.Sinh(Hk) - Hk - M) / ((e * Math.Cosh(Hk)) - 1));
+
+                //Prepare for the next loop around
+                bRecalc = (Math.Abs(Hk1 - Hk) > epsilon);
+                Hk = Hk1;
+            }
+
+            return Hk;
+        }
+
+        private static void CalculateTrueAnomalyAndRadius(double JD, ref AASHyperbolicObjectElements elements, ref double v, ref double r, double epsilon = 0.000001) {
+            var a = elements.q / (elements.e - 1);
+            var n = Math.Sqrt(0.0002959122082855911025 / Math.Pow(Math.Abs(a), 3));
+            var M = (JD - elements.T) * n;
+            var H = CalculateKeplers(M, elements.e, epsilon);
+            v = 2 * Math.Atan(Math.Sqrt((elements.e + 1) / (elements.e - 1)) * Math.Tanh(H / 2));
+            r = Math.Abs(a) * (elements.e * Math.Cosh(H) - 1);
+        }
+
+        /// <param name="JD">The date in Dynamical time to calculate for.</param>
+        /// <param name="elements">An instance of the AASHyperbolicObjectElements class with the orbit elements.</param>
+        /// <param name="bHighPrecision">If true then use the full VSOP87 theory instead of the truncated version as provided in Meeus's book.</param>
+        /// <param name="epsilon">The epsilon value passed to the CalculateKeplers method.</param>
+        /// <returns></returns>
+        public static AASHyperbolicObjectDetails Calculate(double JD, ref AASHyperbolicObjectElements elements, bool bHighPrecision, double epsilon = 0.000001) {
+            var Epsilon = AASNutation.MeanObliquityOfEcliptic(elements.JDEquinox);
+            var JD0 = JD;
+
+            //What will be the return value
+            var details = new AASHyperbolicObjectDetails();
+
+            Epsilon = AASCoordinateTransformation.DegreesToRadians(Epsilon);
+            var omega = AASCoordinateTransformation.DegreesToRadians(elements.omega);
+            var w = AASCoordinateTransformation.DegreesToRadians(elements.w);
+            var i = AASCoordinateTransformation.DegreesToRadians(elements.i);
+
+            var sinEpsilon = Math.Sin(Epsilon);
+            var cosEpsilon = Math.Cos(Epsilon);
+            var sinOmega = Math.Sin(omega);
+            var cosOmega = Math.Cos(omega);
+            var cosi = Math.Cos(i);
+            var sini = Math.Sin(i);
+
+            var F = cosOmega;
+            var G =  sinOmega * cosEpsilon;
+            var H = sinOmega * sinEpsilon;
+            var P = -sinOmega * cosi;
+            var Q = (cosOmega * cosi * cosEpsilon) - (sini * sinEpsilon);
+            var R = (cosOmega * cosi * sinEpsilon) + (sini * cosEpsilon);
+            var a = Math.Sqrt((F * F) + (P * P));
+            var b = Math.Sqrt((G * G) + (Q * Q));
+            var c = Math.Sqrt((H * H) + (R * R));
+            var A = Math.Atan2(F, P);
+            var B = Math.Atan2(G, Q);
+            var C = Math.Atan2(H, R);
+
+            var SunCoord = AASSun.EquatorialRectangularCoordinatesAnyEquinox(JD, elements.JDEquinox, bHighPrecision);
+            for (int j = 0; j < 2; j++) {
+                var v = 0.0;
+                var r = 0.0;
+                CalculateTrueAnomalyAndRadius(JD0, ref elements, ref v, ref r, epsilon);
+                var x = r * a * Math.Sin(A + w + v);
+                var y = r * b * Math.Sin(B + w + v);
+                var z = r * c * Math.Sin(C + w + v);
+
+                if (j == 0) {
+                    details.HeliocentricRectangularEquatorial.X = x;
+                    details.HeliocentricRectangularEquatorial.Y = y;
+                    details.HeliocentricRectangularEquatorial.Z = z;
+
+                    //Calculate the heliocentric ecliptic coordinates also
+                    var u = w + v;
+                    var cosu = Math.Cos(u);
+                    var sinu = Math.Sin(u);
+
+                    details.HeliocentricRectangularEcliptical.X = r * ((cosOmega * cosu) - (sinOmega * sinu * cosi));
+                    details.HeliocentricRectangularEcliptical.Y = r * ((sinOmega * cosu) + (cosOmega * sinu * cosi));
+                    details.HeliocentricRectangularEcliptical.Z = r * sini * sinu;
+
+                    details.HeliocentricEclipticLongitude = AASCoordinateTransformation.MapTo0To360Range(AASCoordinateTransformation.RadiansToDegrees(Math.Atan2(details.HeliocentricRectangularEcliptical.Y, details.HeliocentricRectangularEcliptical.X)));
+                    details.HeliocentricEclipticLatitude = AASCoordinateTransformation.RadiansToDegrees(Math.Asin(details.HeliocentricRectangularEcliptical.Z / r));
+                }
+
+                var psi = SunCoord.X + x;
+                var psi2 = psi * psi;
+                var nu = SunCoord.Y + y;
+                var nu2 = nu * nu;
+                var sigma = SunCoord.Z + z;
+
+                var Alpha = Math.Atan2(nu, psi);
+                Alpha = AASCoordinateTransformation.RadiansToDegrees(Alpha);
+                var Delta = Math.Atan2(sigma, Math.Sqrt(psi2 + nu2));
+                Delta = AASCoordinateTransformation.RadiansToDegrees(Delta);
+                var Distance = Math.Sqrt(psi2 + nu2 + (sigma * sigma));
+                var Distance2 = Distance * Distance;
+
+                if (j == 0) {
+                    details.TrueGeocentricRA = AASCoordinateTransformation.MapTo0To24Range(Alpha / 15);
+                    details.TrueGeocentricDeclination = Delta;
+                    details.TrueGeocentricDistance = Distance;
+                    details.TrueGeocentricLightTime = AASElliptical.DistanceToLightTime(Distance);
+                }
+                else {
+                    details.AstrometricGeocentricRA = AASCoordinateTransformation.MapTo0To24Range(Alpha / 15);
+                    details.AstrometricGeocentricDeclination = Delta;
+                    details.AstrometricGeocentricDistance = Distance;
+                    details.AstrometricGeocentricLightTime = AASElliptical.DistanceToLightTime(Distance);
+
+                    var RES = Math.Sqrt((SunCoord.X * SunCoord.X) + (SunCoord.Y * SunCoord.Y) + (SunCoord.Z * SunCoord.Z));
+                    var RES2 = RES * RES;
+                    var r2 = r * r;
+
+                    details.Elongation = AASCoordinateTransformation.RadiansToDegrees(Math.Acos((RES2 + Distance2 - r2) / (2 * RES * Distance)));
+                    details.PhaseAngle = AASCoordinateTransformation.RadiansToDegrees(Math.Acos((r2 + Distance2 - RES2) / (2 * r * Distance)));
+                }
+
+                if (j == 0) //Prepare for the next loop around
+                    JD0 = JD - details.TrueGeocentricLightTime;
+            }
+
+            return details;
+        }
+    }
+}

Src/AASharp/AASHyperbolicObjectDetails.cs

@@ -0,0 +1,81 @@
+namespace AASharp {
+    /// <summary>
+    /// An instance of this class is returned by AASHyperbolic.Calculate method.
+    /// </summary>
+    public class AASHyperbolicObjectDetails {
+        public AASHyperbolicObjectDetails() {
+            HeliocentricRectangularEquatorial = new AAS3DCoordinate();
+            HeliocentricRectangularEcliptical = new AAS3DCoordinate();
+        }
+
+        /// <summary>
+        /// 3D rectangular heliocentric equatorial coordinates of the object.
+        /// </summary>
+        public AAS3DCoordinate HeliocentricRectangularEquatorial { get; set; }
+
+        /// <summary>
+        /// 3D rectangular heliocentric ecliptical coordinates of the object.
+        /// </summary>
+        public AAS3DCoordinate HeliocentricRectangularEcliptical { get; set; }
+
+        /// <summary>
+        /// The heliocentric ecliptical longitude in degrees of the object.
+        /// </summary>
+        public double HeliocentricEclipticLongitude { get; set; }
+
+        /// <summary>
+        /// The heliocentric ecliptical latitude in degrees of the object.
+        /// </summary>
+        public double HeliocentricEclipticLatitude { get; set; }
+
+        /// <summary>
+        /// The geocentric right ascension of the object as an hour angle (i.e. without light time correction applied).
+        /// </summary>
+        public double TrueGeocentricRA { get; set; }
+
+        /// <summary>
+        /// The geocentric declination of the object in degrees (i.e. without light time correction applied).
+        /// </summary>
+        public double TrueGeocentricDeclination { get; set; }
+
+        /// <summary>
+        /// The true distance in astronomical units between the Earth and the object.
+        /// </summary>
+        public double TrueGeocentricDistance { get; set; }
+
+        /// <summary>
+        /// The light travel time in days from the Earth to the object.
+        /// </summary>
+        public double TrueGeocentricLightTime { get; set; }
+
+        /// <summary>
+        /// The geocentric right ascension of the object as an hour angle (i.e. with light time correction applied).
+        /// </summary>
+        public double AstrometricGeocentricRA { get; set; }
+
+        /// <summary>
+        /// The geocentric declination of the object in degrees (i.e. with light time correction applied).
+        /// </summary>
+        public double AstrometricGeocentricDeclination { get; set; }
+
+        /// <summary>
+        /// The observed distance of the object in astronomical units.
+        /// </summary>
+        public double AstrometricGeocentricDistance { get; set; }
+
+        /// <summary>
+        /// The observed light travel time of the object in days.
+        /// </summary>
+        public double AstrometricGeocentricLightTime { get; set; }
+
+        /// <summary>
+        /// The elongation of the object to the Sun in degrees.
+        /// </summary>
+        public double Elongation { get; set; }
+
+        /// <summary>
+        /// The phase angle (the angle Sun - object - Earth) in degrees.
+        /// </summary>
+        public double PhaseAngle { get; set; }
+    }
+}

Src/AASharp/AASHyperbolicObjectElements.cs

@@ -0,0 +1,41 @@
+namespace AASharp {
+    /// <summary>
+    /// An instance of this class is required as a parameter for the AASHyperbolic.Calculate method.
+    /// </summary>
+    public class AASHyperbolicObjectElements {
+        /// <summary>
+        /// The perihelion distance in astronomical units.
+        /// </summary>
+        public double q { get; set; }
+
+        /// <summary>
+        /// The inclination of the plane of the orbit in degrees.
+        /// </summary>
+        public double i { get; set; }
+
+        /// <summary>
+        /// The argument of the perihelion in degrees.
+        /// </summary>
+        public double w { get; set; }
+
+        /// <summary>
+        /// The longitude of the ascending node in degrees.
+        /// </summary>
+        public double omega { get; set; }
+
+        /// <summary>
+        /// The Julian Date for which <b>equatorial coordinates</b> should be calculated for.
+        /// </summary>
+        public double JDEquinox { get; set; }
+
+        /// <summary>
+        /// The Julian date of the time of passage in perihelion.
+        /// </summary>
+        public double T { get; set; }
+
+        /// <summary>
+        /// Eccentricity of the orbit.
+        /// </summary>
+        public double e { get; set; }
+    }
+}

Src/AASharp/AASharp.csproj

@@ -9,9 +9,9 @@ AASharp is a C# port of PJ Naughter's static C++ library, AA+. Naughter's librar
     </Description>
     <GeneratePackageOnBuild>true</GeneratePackageOnBuild>
     <PackageProjectUrl>https://github.com/jsauve/AASharp</PackageProjectUrl>
-    <AssemblyVersion>2.12.1</AssemblyVersion>
-    <FileVersion>2.12.1</FileVersion>
-    <Version>2.12.1</Version>
+    <AssemblyVersion>2.13.0</AssemblyVersion>
+    <FileVersion>2.13.0</FileVersion>
+    <Version>2.13.0</Version>
     <Authors>Joe Sauve, Julien Tschäppät</Authors>
     <RepositoryUrl>https://github.com/jsauve/AASharp</RepositoryUrl>
     <PackageTags>Astronomy,Algorithm</PackageTags>

Src/AASharpExamples/Program.cs

@@ -996,6 +996,26 @@ static void Main(string[] args)
             elements6.JDEquinox = elements6.T;
             AASNearParabolicObjectDetails details4 = AASNearParabolic.Calculate(elements6.T - 63.6954, ref elements6, bHighPrecision);
 
+            // Test out the AASHyperbolic class
+            Console.WriteLine();
+            Console.WriteLine("Information about the comet 3I/Atlas with a hyperbolic orbit:");
+            var tt_now = DateTime.UtcNow.ToTerrestrialDynamical();
+
+            // Orbital elements for 3I/ATLAS comet:
+            // https://minorplanetcenter.net/db_search/show_object?utf8=%E2%9C%93&object_id=3I
+            var _3I_ATLAS_elements = new AASHyperbolicObjectElements() {
+                q = 1.3564534,          // perihelion distance in AU
+                i = 175.11293,          // inclination in degrees
+                w = 128.00774,          // argument of perihelion in degrees
+                omega = 322.15483,      // longitude of ascending node in degrees
+                JDEquinox = tt_now,     // equinox of elements
+                T = AASDynamicalTime.UTC2TT(new AASDate(2025, 10, 29.48307, true).Julian),   // time of perihelion passage
+                e = 6.1394269			// eccentricity
+            };
+            var res = AASHyperbolic.Calculate(tt_now, ref _3I_ATLAS_elements, true);
+            Console.WriteLine($"True geocentric RA (hour angles): {res.TrueGeocentricRA}");
+            Console.WriteLine($"True geocentric DE (degrees):     {res.TrueGeocentricDeclination}");
+
             // test extensions
             Console.WriteLine();