diff --git a/DataConverter/DataConverter.csproj b/DataConverter/DataConverter.csproj index 2f3dd48..8fd630f 100644 --- a/DataConverter/DataConverter.csproj +++ b/DataConverter/DataConverter.csproj @@ -2,7 +2,7 @@ Exe - net6.0 + net6.0;net7.0 enable enable x64 @@ -26,11 +26,11 @@ - + - + diff --git a/DataConverter/DataReader.cs b/DataConverter/DataReader.cs index 58fcf7a..7d3c80a 100644 --- a/DataConverter/DataReader.cs +++ b/DataConverter/DataReader.cs @@ -54,12 +54,14 @@ public static class DataReader public static List ReadData() { - List VSOP2013DATA = new List(); + List VSOP2013DATA = new(); for (int ip = 0; ip < 9; ip++) { - PlanetTable planet = new PlanetTable(); - planet.body = (VSOPBody)ip; - planet.variables = new VariableTable[6]; + PlanetTable planet = new() + { + body = (VSOPBody)ip, + variables = new VariableTable[6] + }; for (int iv = 0; iv < 6; iv++) { @@ -115,7 +117,6 @@ private static void ReadPlanet(PlanetTable Planet, int ip) ReadTerm(line, ref buffer[i]); } - Planet.variables[H.iv].PowerTables[H.it].header = H; Planet.variables[H.iv].PowerTables[H.it].iv = H.iv; Planet.variables[H.iv].PowerTables[H.it].it = H.it; Planet.variables[H.iv].PowerTables[H.it].Body = (VSOPBody)H.ip; @@ -127,7 +128,7 @@ private static void ReadPlanet(PlanetTable Planet, int ip) private static Header ReadHeader(string line) { - Header H = new Header(); + Header H = new(); int lineptr = 9; H.ip = Convert.ToInt32(line.Substring(lineptr, 3).Trim()) - 1; lineptr += 3; @@ -148,7 +149,7 @@ private static Term ReadTerm(string line, ref Term T) // lineptr = 5; - T.rank = Convert.ToInt32(line.Substring(0, lineptr)); + T.rank = Convert.ToInt32(line[..lineptr]); // lineptr++; // @@ -204,7 +205,7 @@ private static Term ReadTerm(string line, ref Term T) ci = Convert.ToDouble(line.Substring(lineptr, 20).Trim()); lineptr += 20; lineptr++; - ci = ci * Math.Pow(10, Convert.ToDouble(line.Substring(lineptr, 3).Trim())); + ci *= Math.Pow(10, Convert.ToDouble(line.Substring(lineptr, 3).Trim())); T.cc = ci; @@ -216,7 +217,6 @@ private static Term ReadTerm(string line, ref Term T) T.aa += Bufferiphi[j] * ci0[j]; T.bb += Bufferiphi[j] * ci1[j]; } - return T; } diff --git a/DataConverter/Program.cs b/DataConverter/Program.cs index d7116ac..5d89daa 100644 --- a/DataConverter/Program.cs +++ b/DataConverter/Program.cs @@ -1,6 +1,6 @@ -using MessagePack; -using System.Diagnostics; -using System.Runtime.Serialization.Formatters.Binary; +using System.Diagnostics; +using System.IO.Compression; +using MessagePack; namespace VSOP2013.DataConverter { @@ -12,18 +12,17 @@ private static void Main(string[] args) #region Read Original Data - Stopwatch sw = new Stopwatch(); + Stopwatch sw = new(); sw.Start(); List VSOP2013DATA = DataReader.ReadData(); + VSOP2013DATA = VSOP2013DATA.OrderBy(x => x.body).ToList(); sw.Stop(); double ticks = sw.ElapsedTicks; double Freq = Stopwatch.Frequency; double milliseconds = (ticks / Freq) * 1000; - Console.WriteLine($"Data Read OK...Elapsed milliseconds: {milliseconds} ms"); - Console.WriteLine("Press Enter to dump data..."); - Console.ReadLine(); + Console.WriteLine($"Data Read & Convert OK...Elapsed milliseconds: {milliseconds} ms"); #endregion Read Original Data @@ -34,7 +33,7 @@ private static void Main(string[] args) { Directory.CreateDirectory(OutputDirPath); } - DirectoryInfo OutputDir = new DirectoryInfo(Directory.GetCurrentDirectory() + @"\Data"); + DirectoryInfo OutputDir = new(Directory.GetCurrentDirectory() + @"\Data"); sw.Restart(); @@ -45,40 +44,35 @@ private static void Main(string[] args) { File.Delete(filename); } - using (FileStream fs = new FileStream(filename, FileMode.OpenOrCreate)) - { - BinaryFormatter bf = new BinaryFormatter(); - MessagePackSerializer.Serialize(fs, VSOP2013DATA[ip], lz4Options); - } + using FileStream fs = new(filename, FileMode.OpenOrCreate); + using BrotliStream bs = new(fs, CompressionLevel.SmallestSize); + MessagePackSerializer.Serialize(bs, VSOP2013DATA[ip]); }); sw.Stop(); ticks = sw.ElapsedTicks; milliseconds = (ticks / Freq) * 1000; Console.WriteLine($"Data Dumped OK. Elapsed: {milliseconds}ms"); - Console.WriteLine("Press Enter to test dumped data..."); - Console.ReadLine(); #endregion Dump Data #region Test sw.Restart(); - + VSOP2013DATA.Clear(); result = Parallel.For(0, 9, ip => { string filename = Path.Combine(OutputDir.FullName, string.Format("VSOP2013_{0}.BIN", ((VSOPBody)ip).ToString())); - using (FileStream fs = new FileStream(filename, FileMode.OpenOrCreate)) - { - - VSOP2013DATA.Add(MessagePackSerializer.Deserialize(fs, lz4Options)); - } + using FileStream fs = new(filename, FileMode.OpenOrCreate); + using BrotliStream bs = new(fs, CompressionMode.Decompress); + VSOP2013DATA.Add(MessagePackSerializer.Deserialize(bs)); }); sw.Stop(); ticks = sw.ElapsedTicks; milliseconds = (ticks / Freq) * 1000; - Console.WriteLine($"Data Reload Test OK. Elapsed: {milliseconds}ms"); + Console.WriteLine($"Dump Data Reload Test OK. Elapsed: {milliseconds}ms"); + Console.WriteLine("Press Enter to exit..."); Console.ReadLine(); #endregion Test diff --git a/DataConverter/Resources/VSOP2013.ctl b/DataConverter/Resources/VSOP2013.ctl index 130216a..c6826c7 100644 --- a/DataConverter/Resources/VSOP2013.ctl +++ b/DataConverter/Resources/VSOP2013.ctl @@ -1,5 +1,4 @@ - - PLANETARY EPHEMERIS VSOP2013 MERCURY + PLANETARY EPHEMERIS VSOP2013 MERCURY 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -50,7 +49,7 @@ Julian Date JD 2411545.0 (TDB) -0.1300935038 -0.4472876448 -0.0245983783 0.0213663969 -0.0064479843 -0.0024878666 -0.1300936046 -0.4005937206 -0.2004893069 0.0213663956 -0.0049262997 -0.0048474330 - PLANETARY EPHEMERIS VSOP2013 VENUS + PLANETARY EPHEMERIS VSOP2013 VENUS 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -101,7 +100,7 @@ Julian Date JD 2411545.0 (TDB) -0.7183022848 -0.0326544811 0.0410142477 0.0007981222 -0.0202952185 -0.0003234552 -0.7183022964 -0.0462742464 0.0246406381 0.0007981175 -0.0184918375 -0.0083697353 - PLANETARY EPHEMERIS VSOP2013 EMB + PLANETARY EPHEMERIS VSOP2013 EMB 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -152,7 +151,7 @@ Julian Date JD 2411545.0 (TDB) -0.1771590071 0.9672193117 -0.0000009492 -0.0172031083 -0.0031639158 0.0000000254 -0.1771587839 0.8874068590 0.3847367185 -0.0172031091 -0.0029028420 -0.0012585096 - PLANETARY EPHEMERIS VSOP2013 MARS + PLANETARY EPHEMERIS VSOP2013 MARS 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -203,7 +202,7 @@ Julian Date JD 2411545.0 (TDB) 1.3907159210 -0.0134159938 -0.0344678037 0.0006714961 0.0151872482 0.0003016544 1.3907159214 0.0014012149 -0.0369601677 0.0006714996 0.0138140375 0.0063179004 - PLANETARY EPHEMERIS VSOP2013 JUPITER + PLANETARY EPHEMERIS VSOP2013 JUPITER 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -254,7 +253,7 @@ Julian Date JD 2411545.0 (TDB) 4.0011764936 2.9385770224 -0.1017848771 -0.0045683150 0.0064432050 0.0000755810 4.0011771819 2.7365785897 1.0755125254 -0.0045683135 0.0058814622 0.0026323029 - PLANETARY EPHEMERIS VSOP2013 SATURN + PLANETARY EPHEMERIS VSOP2013 SATURN 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -305,7 +304,7 @@ Julian Date JD 2411545.0 (TDB) 6.4064073173 6.5699911538 -0.3690759502 -0.0042923528 0.0038903147 0.0001029490 6.4064088704 6.1746578061 2.2747707349 -0.0042923519 0.0035283446 0.0016419315 - PLANETARY EPHEMERIS VSOP2013 URANUS + PLANETARY EPHEMERIS VSOP2013 URANUS 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -356,7 +355,7 @@ Julian Date JD 2411545.0 (TDB) 14.4318597263 -13.7343125162 -0.2381402777 0.0026781039 0.0026726980 -0.0000247713 14.4318565807 -12.5062632452 -5.6816829828 0.0026781045 0.0024620054 0.0010404106 - PLANETARY EPHEMERIS VSOP2013 NEPTUNE + PLANETARY EPHEMERIS VSOP2013 NEPTUNE 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -407,7 +406,7 @@ Julian Date JD 2411545.0 (TDB) 16.8120537367 -24.9917592776 0.1272247241 0.0025792738 0.0017769013 -0.0000959080 16.8120479567 -22.9801038994 -9.8244204429 0.0025792742 0.0016684245 0.0006188166 - PLANETARY EPHEMERIS VSOP2013 PLUTO + PLANETARY EPHEMERIS VSOP2013 PLUTO 1/ Elliptic Elements: a(au), lambda (radian), k, h, q, p - Dynamical Frame J2000 2/ Ecliptic Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - Dynamical Frame J2000 @@ -456,4 +455,4 @@ Julian Date JD 2411545.0 (TDB) Julian Date JD 2451545.0 (TDB) 39.2648542648 4.1726045776 -0.1758641167 -0.1701234143 -0.0517015914 0.1398654514 -9.8753625435 -27.9588613710 5.8504463318 0.0030287536 -0.0015378008 -0.0007122001 - -9.8753695808 -27.9789262247 -5.7537118247 0.0030287533 -0.0011276087 -0.0012651327 + -9.8753695808 -27.9789262247 -5.7537118247 0.0030287533 -0.0011276087 -0.0012651327 \ No newline at end of file diff --git a/Demo/Demo.csproj b/Demo/Demo.csproj index 1ad34c1..45f5d72 100644 --- a/Demo/Demo.csproj +++ b/Demo/Demo.csproj @@ -2,12 +2,16 @@ Exe - net6.0 + net6.0;net7.0 x64 - + + + + + diff --git a/Demo/PerfTest.cs b/Demo/PerfTest.cs new file mode 100644 index 0000000..2fb3f18 --- /dev/null +++ b/Demo/PerfTest.cs @@ -0,0 +1,39 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Jobs; +using VSOP2013; + +namespace Demo +{ + [SimpleJob(RuntimeMoniker.Net70)] + [MemoryDiagnoser] + public class PerfTest + { + private Calculator vsop; + + private DateTime dt; + VSOPTime vTime; + public PerfTest() { vsop = new Calculator(); } + + [GlobalSetup] + public void init() + { + dt = DateTime.Now.ToUniversalTime(); + dt.ToUniversalTime(); + dt = dt.AddSeconds(-69.184); + vTime =new VSOPTime(dt); + } + + [Benchmark] + public VSOPResult Go() + { + var ell = vsop.GetPlanetPosition(VSOPBody.JUPITER, vTime); + return ell; + } + + } +} diff --git a/Demo/Program.cs b/Demo/Program.cs index 502d51c..428541c 100644 --- a/Demo/Program.cs +++ b/Demo/Program.cs @@ -1,120 +1,156 @@ using System; -using System.Diagnostics; using System.Globalization; -using System.IO; -using System.Threading.Tasks; using VSOP2013; -using VSOP2013.VSOPResult; +using BenchmarkDotNet.Running; +using BenchmarkDotNet.Configs; namespace Demo { - class Program + internal class Program { - static Calculator vsop; - static VSOPTime vTime; - static void Main(string[] args) - { - DateTime Tinput; + private static Calculator vsop; + private static void Main(string[] args) + { vsop = new Calculator(); - - Tinput = DateTime.Now.ToUniversalTime(); - //Console.WriteLine("Parse UTC string that conforms to ISO 8601: 2018-08-18T07:22:16.0000000Z"); - + //Parse Time //while (true) //{ // Console.Write("Input Time As UTC:"); // //string inputT = Console.ReadLine(); - // string inputT = "2000-01-01T12:00:00.0000000Z"; - // CultureInfo culture = CultureInfo.CreateSpecificCulture("en-US"); - // DateTimeStyles style = DateTimeStyles.AdjustToUniversal; - // if (DateTime.TryParse(inputT, culture, style, out Tinput)) break; - // else Console.WriteLine("Invalid Entry..."); - //} - //Tinput = VSOPTime.TDBtoTT(Tinput); - //Tinput = VSOPTime.TTtoTAI(Tinput); - //Tinput = VSOPTime.TAItoUTC(Tinput); - - Console.WriteLine(Tinput.ToString()); - - //Convert UTC to TDB (Barycentric Dynamical Time) - vTime = new VSOPTime(Tinput); - + DateTime dt = DateTime.Now.ToUniversalTime(); + string inputT = "2000-01-01T12:00:00.0000000Z"; + CultureInfo culture = CultureInfo.CreateSpecificCulture("en-US"); + DateTimeStyles style = DateTimeStyles.AdjustToUniversal; + DateTime.TryParse(inputT, culture, style, out dt); + dt.ToUniversalTime(); + dt = dt.AddSeconds(-69.184); + VSOPTime vTime = new VSOPTime(dt); Console.WriteLine(); + Console.WriteLine("Start Substitution..."); - Console.WriteLine("Press Enter To Start Substitution..."); - Console.ReadLine(); - - DynamicalELL[] ResultAll = vsop.CalcAllPlanet(vTime); + VSOPResult_ELL ell; + VSOPResult_XYZ xyz; + VSOPResult_LBR lbr; - foreach(DynamicalELL result in ResultAll) + foreach(VSOPBody body in Enum.GetValues(typeof(VSOPBody))) { - FormattedPrint(result, vTime); + ell = vsop.GetPlanetPosition(body, vTime); + xyz = (VSOPResult_XYZ)ell; + lbr = (VSOPResult_LBR)ell; + FormattedPrint(ell, vTime); + FormattedPrint(xyz, vTime); + FormattedPrint(lbr, vTime); } - Console.WriteLine("Press Enter to Start PerformanceTest."); - Console.ReadLine(); + + - - PerformanceTest(1000); + Console.WriteLine("Press Enter to Start Performance Test..."); Console.ReadLine(); +#if DEBUG + var summary = BenchmarkRunner.Run(new DebugBuildConfig()); +#else + var summary = BenchmarkRunner.Run(); +#endif } - public static void FormattedPrint(DynamicalELL Result, VSOPTime vtime) + public static void FormattedPrint(VSOPResult Result, VSOPTime vtime) { Console.WriteLine("==============================================================="); - Console.WriteLine("PLANETARY EPHEMERIS VSOP2013"); + WriteColorLine(ConsoleColor.Cyan, "PLANETARY EPHEMERIS VSOP2013"); Console.WriteLine("==============================================================="); - Console.WriteLine(Enum.GetName(typeof(VSOPBody), Result.Body) + " at UTC:" + vtime.UTC.ToString()); - Console.WriteLine(); - Console.WriteLine(" Elliptic Elements - Dynamical Frame J2000"); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "semi-major axis (au)", (Result).A)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "mean longitude (rd)", (Result).L)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "k = e*cos(pi) (rd)", (Result).K)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "h = e*sin(pi) (rd)", (Result).H)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "q = sin(i/2)*cos(omega) (rd)", (Result).Q)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "p = sin(i/2)*sin(omega) (rd)", (Result).P)); - Console.WriteLine(); - Console.WriteLine(" Ecliptic Heliocentric Coordinates - Dynamical Frame J2000"); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions X (au)", ((DynamicalXYZ)Result).X)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Y (au)", ((DynamicalXYZ)Result).Y)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Z (au)", ((DynamicalXYZ)Result).Z)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities X'(au/d)", ((DynamicalXYZ)Result).dX)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Y'(au/d)", ((DynamicalXYZ)Result).dY)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Z'(au/d)", ((DynamicalXYZ)Result).dZ)); - Console.WriteLine(); - Console.WriteLine(" Equatorial Heliocentric Coordinates - ICRS Frame J2000"); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions X (au)", ((ICRSXYZ)Result).X)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Y (au)", ((ICRSXYZ)Result).Y)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Z (au)", ((ICRSXYZ)Result).Z)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities X'(au/d)", ((ICRSXYZ)Result).dX)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Y'(au/d)", ((ICRSXYZ)Result).dY)); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Z'(au/d)", ((ICRSXYZ)Result).dZ)); - Console.WriteLine(); - } + if (Result.CoordinatesType == CoordinatesType.Elliptic) + { + WriteColorLine("Coordinates Type: ", ConsoleColor.Green, "\tHeliocentric Elliptic"); + WriteColorLine("Reference Frame: ", ConsoleColor.Green, "\tDynamical Equinox and Ecliptic J2000"); + WriteColorLine("Body: ", ConsoleColor.Green, $"\t\t\t{Enum.GetName(Result.Body)}"); + WriteColorLine("At UTC: ", ConsoleColor.Green, $"\t\t{Result.Time.UTC.ToUniversalTime().ToString("o")}"); + WriteColorLine("At TDB: ", ConsoleColor.Green, $"\t\t{Result.Time.TDB.ToString("o")}"); + Console.WriteLine("---------------------------------------------------------------"); + Console.WriteLine(String.Format("{0,-33}{1,30}", "semi-major axis (au)", (Result as VSOPResult_ELL).a)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "mean longitude (rad)", (Result as VSOPResult_ELL).l)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "k = e*cos(pi) (rad)", (Result as VSOPResult_ELL).k)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "h = e*sin(pi) (rad)", (Result as VSOPResult_ELL).h)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "q = sin(i/2)*cos(omega) (rad)", (Result as VSOPResult_ELL).q)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "p = sin(i/2)*sin(omega) (rad)", (Result as VSOPResult_ELL).p)); + Console.WriteLine("---------------------------------------------------------------"); + Console.WriteLine("e: eccentricity"); + Console.WriteLine("pi: perihelion longitude"); + Console.WriteLine("i: inclination"); + Console.WriteLine("omega: ascending node longitude"); + Console.WriteLine("==============================================================="); + Console.WriteLine(); + } + else if (Result.CoordinatesType == CoordinatesType.Rectangular) + { + WriteColorLine("Coordinates Type: ", ConsoleColor.Green, "\t Heliocentric Rectangular"); + if (Result.InertialFrame == InertialFrame.Dynamical) + { + WriteColorLine("Reference Frame: ", ConsoleColor.Green, "\tDynamical Ecliptic J2000"); + } + else + { + WriteColorLine("Reference Frame: ", ConsoleColor.Green, "\tEquatorial ICRS J2000"); + } - public static void PerformanceTest(int cycle) - { - Console.WriteLine(); - Console.WriteLine("=====================Start Performance Test====================="); - Stopwatch sw = new Stopwatch(); - sw.Start(); - int completedCycle = 0; - for (int i = 0; i < cycle; i++) + WriteColorLine("Body: ", ConsoleColor.Green, $"\t\t\t{Enum.GetName(Result.Body)}"); + WriteColorLine("At UTC: ", ConsoleColor.Green, $"\t\t{Result.Time.UTC.ToUniversalTime().ToString("o")}"); + WriteColorLine("At TDB: ", ConsoleColor.Green, $"\t\t{Result.Time.TDB.ToString("o")}"); + + Console.WriteLine("---------------------------------------------------------------"); + Console.WriteLine(String.Format("{0,-33}{1,30}", "position x (au)", (Result as VSOPResult_XYZ).x)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "position y (au)", (Result as VSOPResult_XYZ).y)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "position z (au)", (Result as VSOPResult_XYZ).z)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "velocity x (au/day)", (Result as VSOPResult_XYZ).dx)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "velocity y (au/day)", (Result as VSOPResult_XYZ).dy)); + Console.WriteLine(String.Format("{0,-33}{1,30}", "velocity z (au/day)", (Result as VSOPResult_XYZ).dz)); + Console.WriteLine("---------------------------------------------------------------"); + Console.WriteLine("==============================================================="); + Console.WriteLine(); + } + else if (Result.CoordinatesType == CoordinatesType.Spherical) { - DynamicalELL[] ResultAll = vsop.CalcAllPlanet(vTime); - completedCycle++; - if (completedCycle % 100 == 0 && completedCycleEasy to convert time by calling ```VSOPTime.UTC```, ```VSOPTime.TAI```, ```VSOPTime.TDB``` +3. Veryhigh performance per solution +
![Performance Test](https://github.com/kingsznhone/VSOP2013.NET/blob/main/PerformanceTest.png) +4. Useful Utility class. Convert Elliptic coordinates to cartesian or spherical +5. Async Api. +6. precalculation on φ in terms, which gives 10%+ speed of calculation. +7. Use [MessagePack](https://github.com/neuecc/MessagePack-CSharp#lz4-compression"MessagePack for C#") for binary serialize. +
Initialization time becomes less than 10% of previous. +7. Brotli compression on source data. ~300Mb -> ~50MB. + +
+ ## Introduction The VSOP2013 files contain the series of the elliptic elements for the 8 planets Mercury, Venus, Earth-Moon barycenter, Mars, Jupiter, Saturn, Uranus, and Neptune and for the dwarf planet Pluto of the solution VSOP2013. The planetary solution VSOP2013 is fitted to the numerical integration INPOP10a built at IMCCE, Paris Observatory over the time interval +1890...+2000. @@ -25,169 +47,735 @@ The precision is of a few 0.1″ for the telluric planets (1.6″ for Mars) over G. FRANCOU & J.-L. SIMON (MAY 2013) Ref: Simon J.-L., Francou G., Fienga A., Manche H., A&A 557, A49 (2013) + +The planetary solution VSOP2013 is fitted to the numerical integration INPOP10a built at IMCCE, Paris Observatory (https://ftp.imcce.fr/pub/ephem/planets/vsop2013/). -List of the data files: -VSOP2013p1.dat : Mercury +## How to use -VSOP2013p2.dat : Venus +* NuGet Package Manager + ``` + PM> NuGet\Install-Package VSOP2013.NET -Version 1.1.2 + ``` -VSOP2013p3.dat : Earth-Moon Barycenter +``` +using VSOP2013; -VSOP2013p4.dat : Mars +Calculator vsop = new Calculator(); -VSOP2013p5.dat : Jupiter +//Create VSOPTime using UTC . +DateTime Tinput = DateTime.Now; +VSOPTime vTime = new VSOPTime(Tinput.ToUniversalTime()); -VSOP2013p6.dat : Saturn +//Calculate EMB's present position +VSOPResult_ELL ell = vsop.GetPlanetPosition(VSOPBody.EMB, vTime); -VSOP2013p7.dat : Uranus +//Cast to diffirent coordinate system. +VSOPResult_XYZ xyz=(VSOPResult_XYZ)ell; +VSOPResult_LBR lbr=(VSOPResult_LBR)ell; -VSOP2013p8.dat : Neptune +//Print result +Console.WriteLine($"Coordinates Type: {Enum.GetName(ell.CoordinatesType)}"); +Console.WriteLine($"InertialFrame: {Enum.GetName(ell.InertialFrame)}"); +Console.WriteLine($"Time Frame Reference: {Enum.GetName(ell.TimeFrameReference)}"); +Console.WriteLine($"Body: {Enum.GetName(ell.Body)}"); +Console.WriteLine($"Time: {ell.Time.UTC.ToString("o")}"); +Console.WriteLine("---------------------------------------------------------------"); +Console.WriteLine(String.Format("{0,-33}{1,30}", "semi-major axis (au)", ell.a)); +Console.WriteLine(String.Format("{0,-33}{1,30}", "mean longitude (rad)", ell.l)); +Console.WriteLine(String.Format("{0,-33}{1,30}", "k = e*cos(pi) (rad)", ell.k)); +Console.WriteLine(String.Format("{0,-33}{1,30}", "h = e*sin(pi) (rad)", ell.h)); +Console.WriteLine(String.Format("{0,-33}{1,30}", "q = sin(i/2)*cos(omega) (rad)", ell.q)); +Console.WriteLine(String.Format("{0,-33}{1,30}", "p = sin(i/2)*sin(omega) (rad)", ell.p)); +Console.WriteLine("==============================================================="); +``` -VSOP2013p9.dat : Pluto - -The planetary solution VSOP2013 is fitted to the numerical integration INPOP10a built at IMCCE, Paris Observatory (https://ftp.imcce.fr/pub/ephem/planets/vsop2013/). +## Change Log +### 2023.7.05 v1.1.2 -## FILES DESCRIPTION - - Each VSOP2013 file corresponds to a planet and contains trigonometric series, functions of time (Periodic series and Poisson series), that represent the 6 elliptic elements of the planet: +Many improvements. -Variable 1 : a = semi-major axis (ua) +Some of them are from VSOP87.NET -Variable 2 : λ = mean longitude (radian) +### 2023.03.26 V1.1 -Variable 3 : k = e cos ω +Initialize performance upgrade. -Variable 4 : h = e sin ω +Add result classes. -Variable 5 : q = sin(i/2) cos Ω +Use MessagePack to compress original data. -Variable 6 : p = sin(i/2) sin Ω - -### with: +### 2022.06.02 -e : eccentricity +New features. -ω : perihelion longitude +Performance Optimization. -i : inclination +Upgrade to .NET 6 -Ω : ascending node longitude - -### VSOP2013 series are characterized by 3 parameters: +### 2020.11.14 -- the planet index 1-9 from Mercury to Pluto, +Upgrade to .NET 5 + +### 2020.07.06 -- the variable index 1-6 for a, λ, k, h, q, p, +Initial PR. -- the time power α. +
-## Feature +## Enviroment Require -1. Use ResultBase class to manage calculate results. Easy to cast diffirent frame. -2. Use VSOPTime class to manage time. Easy to convert time by calling VSOPTime.UTC VSOPTime.TAI VSOPTime.TDB -3. Extremely optimized. About 6ms per full planet result on 5900HX. +.NET6 Runtime -![Performance Test](https://github.com/kingsznhone/VSOP2013.NET/blob/main/Performance.png) +.NET7 Runtime -Use highly optimized binary serialize lib to load data. +.NET8 Runtime (Planned) + +Windows 10 64bit + +## Reference + + [MessagePack](https://github.com/neuecc/MessagePack-CSharp#lz4-compression"MessagePack for C#") -With LZ4 compress. Core DLL is less than 100m now. +
-## API +# API -### public double CalcIV(VSOPBody body,int iv, VSOPTime time) +## Class Calculator + +Provide methods to calculate planet position. + +### ```public double GetVariable(VSOPBody body,int iv, VSOPTime time)``` Calculate a specific variable of a planet. +
+ #### Parameters ```body``` VSOPBody Planet enum. +
+ ```iv``` Variable index 0-5 : a l k h q p +
+ ```time``` VSOPTime Exclusive time class in VSOP2013. +
+ #### Return ```double``` variable result. -### public VSOPResult CalcPlanet(VSOPBody body, VSOPTime time) +
+ +### ```public Task GetVariableAsync(VSOPBody body,int iv, VSOPTime time)``` + +Calculate a specific variable of a planet. + +
+ +#### Parameters + +```body``` VSOPBody + +Planet enum. + +
+ +```iv``` Variable index + +0-5 : a l k h q p + +
+ +```time``` VSOPTime + +Exclusive time class in VSOP2013. + +
+ +#### Return + +```Task``` + +variable result. + +
+ +### ```public VSOPResult_ELL GetPlanetPosition(VSOPBody body, VSOPTime time)``` Calculate all variable of a planet. +
+ #### Parameters ```body``` VSOPBody Planet enum. +
+ ```time``` VSOPTime Exclusive time class in VSOP2013. +
+ #### Return -```VSOPResult``` Full Result with 6 variable in 3 type of Coordinates. +```VSOPResult_ELL``` -- Elliptic Elements Dynamical Frame J2000 +Full Result with 6 variable of Elliptic Coordinates. -- Ecliptic Heliocentric Coordinates Dynamical Frame J2000 - -- Equatorial Heliocentric Coordinates ICRS Frame J2000 +Can be explicit cast to ```VSOPResult_XYZ``` and ```VSOPResult_LBR``` -### public VSOPResult[] CalcAllPlanet(VSOPTime time) +
+ +### ```public Task GetPlanetPositionAsync(VSOPBody body, VSOPTime time)``` Calculate all variable of a planet. +
+ #### Parameters +```body``` VSOPBody + +Planet enum. + +
+ ```time``` VSOPTime Exclusive time class in VSOP2013. +
+ #### Return -```VSOPResult[]``` Full result with 6 variable in 3 type of Coordinates of all 9 planet. +```Task``` -## Change Log +Full Result with 6 variable of Elliptic Coordinates. -### 2023.03.26 V1.1 +Can be explicit cast to ```VSOPResult_XYZ``` and ```VSOPResult_LBR``` -Initialize performance upgrade. +
-Add result classes. +## Static Class Utility -Use MessagePack to compress original data. +### Overview -### 2022.06.02 +This Class Provide some useful function. -New features. +
-Performance Optimization. +### ```static double[,] MultiplyMatrix(double[,] A, double[,] B)``` -Upgrade to .NET 6 +A handy matrix multiply function -### 2020.11.14 +#### Parameters -Upgrade to .NET 5 +```A``` double[,] -### 2020.07.06 +Matrix A -Initial PR. +
+ +```B``` double[,] + +Matrix B + +
+ +#### Return + +```double[,]``` + +Matrix C = AB. + +
+ +### ```static double[] XYZtoLBR(double[] xyz)``` + +#### Parameters + +```xyz``` double[] + +Array of cartesian coordinate elements + +
+ +#### Return + +```double[]``` + +Array of spherical coordinate elements + +
+ +### ```static double[] LBRtoXYZ(double[] lbr)``` + +#### Parameters + +```lbr``` double[] + +Array of spherical coordinate elements + +
+ +#### Return + +```double[]``` + +Array of cartesian coordinate elements + +
+ +### ```static double[] ELLtoXYZ(double[] ell)``` + +This is a magic function I directly copy from VSOP2013. + +It's way beyond my math level. + +So I can't find how to inverse XYZ elements to ELL elements. + +Need help. + +#### Parameters + +```ell``` double[] + +Array of elliptic coordinate elements + +
+ +#### Return + +```double[]``` + +Array of cartesian coordinate elements + +
+ +### ```static double[] ELLtoLBR(double[] ell)``` + +#### Parameters + +```ell``` double[] + +Array of elliptic coordinate elements + +
+ +#### Return + +```double[]``` + +Array of spherical coordinate elements + +
+ +### ```static double[] DynamicaltoICRS(double[] xyz)``` + +#### Parameters + +```xyz``` double[] + +Array of cartesian coordinate elements that inertial frame of dynamical equinox and ecliptic. + +
+ +#### Return + +```double[]``` + +Array of cartesian coordinate elements that inertial frame of ICRS equinox and ecliptic. + +
+ +### ```static double[] ICRStoDynamical(double[] xyz)``` + +#### Parameters + +```xyz``` double[] + +Array of cartesian coordinate elements that inertial frame of ICRS equinox and ecliptic. + +
+ +#### Return + +```double[]``` + +Array of cartesian coordinate elements that inertial frame of dynamical equinox and ecliptic. + +
+ +## Class VSOPResult_XYZ : VSOPResult + +### Properties + +```VSOPBody Body { get; }``` + +Planet of this result. + +
+ +```CoordinatesType CoordinatesType { get; }``` + +Coordinates type of this result. + +
+ +```InertialFrame InertialFrame { get; }``` + +Inertial frame of this result. + +```TimeFrameReference TimeFrameReference { get; }``` + +Time frame reference of this result. + +
+ +```VSOPTime Time { get; }``` + +Input time of this result. + +
+ +```double[] Variables_ELL { get;}``` + +Raw data of this result in elliptic coordinate. + +
+ +```double[] Variables_XYZ { get;}``` + +Raw data of this result in cartesian coordinate. + +
+ +```double x {get;}``` + +Position x (au) + +```double y {get;}``` + +Position y (au) + +```double z {get;}``` + +Position z (au) + +```double dx {get;}``` + +Velocity x (au/day) + +```double dy {get;}``` + +Velocity y (au/day) + +```double dz {get;}``` + +Velocity z (au/day) + +
+ +### Methods + +```void ToICRSFrame()``` + +Convert this result to ICRS Inertial Frame. + +
+ +```void ToDynamicalFrame()``` + +Convert this result to Dynamical Inertial Frame. + +
+ +## Class VSOPResult_ELL : VSOPResult + +### Properties + +```VSOPBody Body { get; }``` + +Planet of this result. + +
+ +```CoordinatesType CoordinatesType { get; }``` + +Coordinates type of this result. + +
+ +```InertialFrame InertialFrame { get; }``` + +Inertial frame of this result. + +```TimeFrameReference TimeFrameReference { get; }``` + +Time frame reference of this result. + +
+ +```VSOPTime Time { get; }``` + +Input time of this result. + +
+ +```double[] Variables_ELL { get;}``` + +Raw data of this result in elliptic coordinate. + +
+ +```double[] Variables_LBR { get;}``` + +Raw data of this result in spherical coordinate. + +
+ + +```double a {get;}``` + +Semi-major axis (au) + +```double l {get;}``` + +Mean longitude (rd) + +```double k {get;}``` + +e*cos(pi) (rd) + +```double h {get;}``` + +e*sin(pi) (rd) + +```double q {get;}``` +sin(i/2)*cos(omega) (rd) + +```double p {get;}``` + +sin(i/2)*sin(omega) (rd) + +
+ +## Class VSOPResult_LBR : VSOPResult + +### Properties + +```VSOPVersion Version { get; }``` + +VSOP version of this result. + +
+ +```VSOPBody Body { get; }``` + +Planet of this result. + +
+ +```CoordinatesReference CoordinatesRefrence { get; }``` + +Coordinates reference of this result. + +
+ +```CoordinatesType CoordinatesType { get; }``` + +Coordinates type of this result. + +
+ +```TimeFrameReference TimeFrameReference { get; }``` + +Time frame reference of this result. + +
+ +```VSOPTime Time { get; }``` + +Input time of this result. + +
+ +```double[] Variables { get;}``` + +Raw data of this result. + +
+ +```double l {get;}``` + +longitude (rd) + +```double b {get;}``` + +latitude (rd) + +```double r {get;}``` + +radius (au) + +```double dl {get;}``` + +longitude velocity (rd/day) + +```double db {get;}``` + +latitude velocity (rd/day) + +```double dr {get;}``` + +radius velocity (au/day) + +
+ +### Methods + +```void ToICRSFrame()``` + +Convert this result to ICRS Inertial Frame. + +
+ +```void ToDynamicalFrame()``` + +Convert this result to Dynamical Inertial Frame. + +
+ +## Class VSOPTime +#### summary + +This class provide time convert and management for VSOP87. + +
+ +#### Constructor + +```VSOPTime(DateTime UTC)``` + +Use UTC Time to initialize VSOPTime. + +
+ +#### Properties + +```DateTime UTC``` + +UTC Time frame. + +
+ +```DateTime TAI``` + +International Atomic Time + +
+ +```DateTime TT``` + +Terrestrial Time (aka. TDT) + +
+ +```DateTime TDB``` + +Barycentric Dynamical Time + +VSOP87 use this time frame. + +```double J2000``` + +Get J2000 from TDB. + +
+ +### Methods + +```DateTime ChangeFrame(DateTime dt, TimeFrame TargetFrame)``` + +#### Parameters + +```dt``` DateTime + +A Datetime of any frame. + +
+ +```SourceFrame``` TimeFrame + +Time frame of ```dt``` + +
+ +```TargetFrame``` TimeFrame + +Target time frame. + +
+ +#### Return + +```DateTime``` + +Datetime of target time Frame. + +
+ +```static double ToJ2000(DateTime dt)``` +#### Parameters + +```dt``` DateTime + +Datetime to convert + +
+ +#### Return + +```double``` + +Julian date. + +
+ +```static DateTime FromJ2000(double JD)``` + +#### Parameters + +```double``` + +Julian date to analyze. + +
+ +#### Return + +```dt``` DateTime + +Datetime Class. + +
-## Enviroment Require -.NET6 Runtime -## Reference - - [MessagePack](https://github.com/neuecc/MessagePack-CSharp#lz4-compression"MessagePack for C#") diff --git a/VSOP2013.sln b/VSOP2013.sln index 1be164f..aeb3341 100644 --- a/VSOP2013.sln +++ b/VSOP2013.sln @@ -3,7 +3,7 @@ Microsoft Visual Studio Solution File, Format Version 12.00 # Visual Studio Version 17 VisualStudioVersion = 17.2.32526.322 MinimumVisualStudioVersion = 10.0.40219.1 -Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "VSOP2013", "VSOP2013\VSOP2013.csproj", "{A36BC2BE-D01A-497C-B17A-3B15DA455D9F}" +Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "VSOP2013.NET", "VSOP2013\VSOP2013.NET.csproj", "{A36BC2BE-D01A-497C-B17A-3B15DA455D9F}" EndProject Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Demo", "Demo\Demo.csproj", "{9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}" EndProject @@ -12,7 +12,7 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution .editorconfig = .editorconfig EndProjectSection EndProject -Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "DataConverter", "DataConverter\DataConverter.csproj", "{2A5C5EFE-C49F-436A-9FC1-176AE11E0334}" +Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "DataConverter", "DataConverter\DataConverter.csproj", "{2A5C5EFE-C49F-436A-9FC1-176AE11E0334}" EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution diff --git a/VSOP2013/Calculator.cs b/VSOP2013/Calculator.cs index fa3940f..a870744 100644 --- a/VSOP2013/Calculator.cs +++ b/VSOP2013/Calculator.cs @@ -1,12 +1,9 @@ -using MessagePack; -using System.Diagnostics; +using System.IO.Compression; using System.Reflection; -using System.Runtime.Serialization.Formatters.Binary; -using VSOP2013.VSOPResult; +using MessagePack; namespace VSOP2013 { - public class Calculator { public List VSOP2013DATA; @@ -14,7 +11,7 @@ public class Calculator /// /// //Planetary frequency in longitude /// - readonly double[] freqpla = + private readonly double[] freqpla = { 0.2608790314068555e5, 0.1021328554743445e5, @@ -27,8 +24,7 @@ public class Calculator 0.2533566020437000e2 }; - const double dpi = 2 * Math.PI; - const double a1000 = 365250.0d; + private const double a1000 = 365250.0d; public Calculator() { @@ -36,43 +32,37 @@ public Calculator() var lz4Options = MessagePackSerializerOptions.Standard.WithCompression(MessagePackCompression.Lz4Block); var assembly = Assembly.GetExecutingAssembly(); var names = assembly.GetManifestResourceNames(); - VSOP2013DATA = new List(); + VSOP2013DATA = new List(9); ParallelLoopResult result = Parallel.For(0, 9, ip => { - string datafilename = $"VSOP2013.Resources.VSOP2013_{((VSOPBody)(ip)).ToString()}.BIN"; - using (Stream s = assembly.GetManifestResourceStream(datafilename)) - { - var data = MessagePackSerializer.Deserialize(s, lz4Options); - VSOP2013DATA.Add(data); - } + string datafilename = $"VSOP2013.NET.Resources.VSOP2013_{(VSOPBody)(ip)}.BIN"; + using Stream s = assembly.GetManifestResourceStream(datafilename); + using BrotliStream bs = new(s, CompressionMode.Decompress); + var data = MessagePackSerializer.Deserialize(bs); + VSOP2013DATA.Add(data); }); + VSOP2013DATA = VSOP2013DATA.OrderBy(x => x.body).ToList(); + GC.Collect(); GC.WaitForPendingFinalizers(); } - public DynamicalELL[] CalcAllPlanet(VSOPTime time) - { - DynamicalELL[] vsopresult = new DynamicalELL[9]; - - ParallelLoopResult result = Parallel.For(0, 9, ip => - { - vsopresult[ip] = CalcPlanet((VSOPBody)ip, time); - }); - return vsopresult; - } - - public DynamicalELL CalcPlanet(VSOPBody body, VSOPTime time) + public VSOPResult_ELL GetPlanetPosition(VSOPBody body, VSOPTime time) { - double[] ELL = new double[6]; ParallelLoopResult result = Parallel.For(0, 6, iv => { - ELL[iv] = CalcIV(body, iv, time); + ELL[iv] = GetVariable(body, iv, time); }); - DynamicalELL Coordinate = new DynamicalELL(body, time, ELL); + VSOPResult_ELL Coordinate = new(body, time, ELL); return Coordinate; } + public async Task GetPlanetPositionAsync(VSOPBody body, VSOPTime time) + { + return await Task.Run(() => GetPlanetPosition(body, time)); + } + /// /// Calculate a specific variable /// @@ -80,24 +70,29 @@ public DynamicalELL CalcPlanet(VSOPBody body, VSOPTime time) /// 0-5 : a l k h q p /// /// - public double CalcIV(VSOPBody body, int iv, VSOPTime time) + public double GetVariable(VSOPBody body, int iv, VSOPTime time) + { + return Calculate(VSOP2013DATA[(int)body].variables[iv], time.J2000); + } + + public async Task GetVariableAsync(VSOPBody body, int iv, VSOPTime time) { - return CalcIV(VSOP2013DATA[(int)body].variables[iv], VSOPTime.ToJulianDate2000(time.TDB)); + return await Task.Run(() => GetVariable(body, iv, time)); } /// - /// + /// /// - /// + /// /// /// Elliptic Elements - private double CalcIV(VariableTable vTable, double JD2000) + private double Calculate(VariableTable Table, double JD2000) { //Thousand of Julian Years double tj = JD2000 / a1000; - //Iteration on Time - double[] t = new double[21]; + //Iteration on Time + Span t = stackalloc double[21]; t[0] = 1.0d; t[1] = tj; for (int i = 2; i < 21; i++) @@ -106,27 +101,25 @@ private double CalcIV(VariableTable vTable, double JD2000) } double result = 0d; - double arg; - double sarg; - double carg; + double u, su, cu; double xl; - double tit; - for (int it = 0; it < vTable.PowerTables.Length; it++) + Term[] terms; + for (int it = 0; it < Table.PowerTables.Length; it++) { - tit = t[it]; - if (vTable.PowerTables[it].Terms == null) continue; - for (int n = 0; n < vTable.PowerTables[it].Terms.Length; n++) + if (Table.PowerTables[it].Terms == null) continue; + terms = Table.PowerTables[it].Terms; + for (int n = 0; n < terms.Length; n++) { - arg = vTable.PowerTables[it].Terms[n].aa + vTable.PowerTables[it].Terms[n].bb * tj; - (sarg, carg) = Math.SinCos(arg); - result += tit * (vTable.PowerTables[it].Terms[n].ss * sarg + vTable.PowerTables[it].Terms[n].cc * carg); + u = terms[n].aa + terms[n].bb * tj; + (su, cu) = Math.SinCos(u); + result += t[it] * (terms[n].ss * su + terms[n].cc * cu); } } - if (vTable.iv == 1) + if (Table.iv == 1) { - xl = result + freqpla[(int)vTable.Body] * tj; - xl = xl % dpi; - if (xl < 0) xl = xl + dpi; + xl = result + freqpla[(int)Table.Body] * tj; + xl %= Math.Tau; + if (xl < 0) xl += Math.Tau; result = xl; } return result; diff --git a/VSOP2013/DataStruct.cs b/VSOP2013/DataStruct.cs index d348f22..0894b51 100644 --- a/VSOP2013/DataStruct.cs +++ b/VSOP2013/DataStruct.cs @@ -1,10 +1,5 @@ -using MessagePack; -using System; -using System.Collections.Generic; -using System.Linq; -using System.Runtime.InteropServices; -using System.Text; -using System.Threading.Tasks; +using System.Runtime.InteropServices; +using MessagePack; namespace VSOP2013 { @@ -36,7 +31,7 @@ public enum VSOPFileName public struct PlanetEphemeris { - //Elliptic Elements + //Elliptic Elements //a (au), lambda (radian), k, h, q, p //Dynamical Frame J2000' public double[] DynamicalELL; @@ -52,17 +47,15 @@ public struct PlanetEphemeris public double[] ICRSXYZ; } - - [MessagePackObject] [Serializable] public struct PlanetTable { [Key(0)] public VSOPBody body; + [Key(1)] public VariableTable[] variables; - } [MessagePackObject] @@ -71,8 +64,10 @@ public struct VariableTable { [Key(0)] public VSOPBody Body; + [Key(1)] public int iv; + [Key(2)] public PowerTable[] PowerTables; } @@ -83,12 +78,13 @@ public struct PowerTable { [Key(0)] public VSOPBody Body; + [Key(1)] public int iv; + [Key(2)] public int it; - [Key(3)] - public Header header; + [Key(4)] public Term[] Terms; } @@ -99,10 +95,13 @@ public struct Header { [Key(0)] public int ip; + [Key(1)] public int iv; + [Key(2)] public int it; + [Key(3)] public int nt; } @@ -113,19 +112,23 @@ public struct Header public struct Term { [Key(0)] - [FieldOffset(0)] + [FieldOffset(0)] public long rank; + [Key(1)] [FieldOffset(8)] public double ss; + [Key(2)] [FieldOffset(16)] public double cc; + [Key(3)] [FieldOffset(24)] public double aa; + [Key(4)] [FieldOffset(32)] public double bb; } -} +} \ No newline at end of file diff --git a/VSOP2013/Resources/VSOP2013_EMB.BIN b/VSOP2013/Resources/VSOP2013_EMB.BIN index c2fd13c..d494151 100644 Binary files a/VSOP2013/Resources/VSOP2013_EMB.BIN and b/VSOP2013/Resources/VSOP2013_EMB.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_JUPITER.BIN b/VSOP2013/Resources/VSOP2013_JUPITER.BIN index e8354dd..26606a4 100644 Binary files a/VSOP2013/Resources/VSOP2013_JUPITER.BIN and b/VSOP2013/Resources/VSOP2013_JUPITER.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_MARS.BIN b/VSOP2013/Resources/VSOP2013_MARS.BIN index 59deb43..469c765 100644 Binary files a/VSOP2013/Resources/VSOP2013_MARS.BIN and b/VSOP2013/Resources/VSOP2013_MARS.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_MERCURY.BIN b/VSOP2013/Resources/VSOP2013_MERCURY.BIN index 6df3bf7..0dd66e3 100644 Binary files a/VSOP2013/Resources/VSOP2013_MERCURY.BIN and b/VSOP2013/Resources/VSOP2013_MERCURY.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN b/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN index 946a941..11d94f3 100644 Binary files a/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN and b/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_PLUTO.BIN b/VSOP2013/Resources/VSOP2013_PLUTO.BIN index 3beb4b2..ccf5769 100644 Binary files a/VSOP2013/Resources/VSOP2013_PLUTO.BIN and b/VSOP2013/Resources/VSOP2013_PLUTO.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_SATURN.BIN b/VSOP2013/Resources/VSOP2013_SATURN.BIN index f122054..f77bc7e 100644 Binary files a/VSOP2013/Resources/VSOP2013_SATURN.BIN and b/VSOP2013/Resources/VSOP2013_SATURN.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_URANUS.BIN b/VSOP2013/Resources/VSOP2013_URANUS.BIN index 43c89e0..7908313 100644 Binary files a/VSOP2013/Resources/VSOP2013_URANUS.BIN and b/VSOP2013/Resources/VSOP2013_URANUS.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_VENUS.BIN b/VSOP2013/Resources/VSOP2013_VENUS.BIN index f2f5b2b..3e58a71 100644 Binary files a/VSOP2013/Resources/VSOP2013_VENUS.BIN and b/VSOP2013/Resources/VSOP2013_VENUS.BIN differ diff --git a/VSOP2013/Utility.cs b/VSOP2013/Utility.cs new file mode 100644 index 0000000..9128d41 --- /dev/null +++ b/VSOP2013/Utility.cs @@ -0,0 +1,335 @@ +using System.Numerics; + +namespace VSOP2013 +{ + public static class Utility + { + public static readonly double[] gmp = {4.9125474514508118699e-11d, + 7.2434524861627027000e-10d, + 8.9970116036316091182e-10d, + 9.5495351057792580598e-11d, + 2.8253458420837780000e-07d, + 8.4597151856806587398e-08d, + 1.2920249167819693900e-08d, + 1.5243589007842762800e-08d, + 2.1886997654259696800e-12d}; + + public static readonly double gmsol = 2.9591220836841438269e-04d; + + public static double[,] MultiplyMatrix(double[,] A, double[,] B) + { + int rA = A.GetLength(0); + int cA = A.GetLength(1); + int rB = B.GetLength(0); + int cB = B.GetLength(1); + + if (cA != rB) + { + Console.WriteLine("Matrixes can't be multiplied!!"); + return null; + } + else + { + double temp; + double[,] kHasil = new double[rA, cB]; + + for (int i = 0; i < rA; i++) + { + for (int j = 0; j < cB; j++) + { + temp = 0; + for (int k = 0; k < cA; k++) + { + temp += A[i, k] * B[k, j]; + } + kHasil[i, j] = temp; + } + } + + return kHasil; + } + } + + /// + /// Convert cardinal coordinate to spherical coordinate + /// + /// x,y,z,dx,dy,dz + /// l,b,r,dl,db,dr + public static double[] XYZtoLBR(double[] xyz) + { + Span lbr = stackalloc double[6]; + + double x = xyz[0]; + double y = xyz[1]; + double z = xyz[2]; + double dx = xyz[3]; + double dy = xyz[4]; + double dz = xyz[5]; + double l, b, r, dl, db, dr; + + //r = sqrt(x^2+y^2+z^2) + r = Math.Sqrt(x * x + y * y + z * z); + + //L = sgn(y) * arccos(x/sqrt(x^2+y^2)) + if (y >= 0) + { + l = Math.Acos(x / Math.Sqrt(x * x + y * y)); + } + else + { + l = -Math.Acos(x / Math.Sqrt(x * x + y * y)); + } + + //θ = arccos(z/r) + //b = pi/2 - θ. + //Sin(θ) = Cos(b), Cos(θ) = Sin(b) + b = Math.Asin(z / r); + + // Inverse Jacobian matrix From Caridnal to Sperical + //https://en.wikipedia.org/wiki/Spherical_coordinate_system#Integration_and_differentiation_in_spherical_coordinates + double[,] J_1 = { + {x/r , y/r , z/r }, + {x*z/(r*r*Math.Sqrt(x*x+y*y)) , y*z/(r*r*Math.Sqrt(x*x+y*y)),-(x*x+y*y)/(r*r*Math.Sqrt(x*x+y*y)) }, + {-y/(x*x+y*y) , x/(x*x+y*y) , 0 }}; + double[,] Velocity = { { dx }, { dy }, { dz } }; + var C = MultiplyMatrix(J_1, Velocity); + dl = C[2, 0]; + db = C[1, 0]; + dr = C[0, 0]; + + //modulo r into [0,Tau) + l -= Math.Floor(l / Math.Tau) * Math.Tau; + if (l < 0) + l += Math.Tau; + + lbr[0] = l; + lbr[1] = b; + lbr[2] = r; + lbr[3] = dl; + lbr[4] = -db; + lbr[5] = dr; + return lbr.ToArray(); + } + + /// + /// convert spherical coordinate to cardinal coordinate + /// + /// l,b,r,dl,db,dr + /// x,y,z,dx,dy,dz + public static double[] LBRtoXYZ(double[] lbr) + { + Span xyz = stackalloc double[6]; + double x, y, z, dx, dy, dz; + double l, b, r, dl, db, dr; + + l = lbr[0]; + b = lbr[1]; + r = lbr[2]; + dl = lbr[3]; + db = lbr[4]; + dr = lbr[5]; + + x = r * Math.Cos(b) * Math.Cos(l); + y = r * Math.Cos(b) * Math.Sin(l); + z = r * Math.Sin(b); + + // Jacobian matrix From Sperical to Caridnal + //https://en.wikipedia.org/wiki/Spherical_coordinate_system#Integration_and_differentiation_in_spherical_coordinates + double[,] J = { + { Math.Cos(b) * Math.Cos(l), r * Math.Sin(b) * Math.Cos(l), -r * Math.Cos(b) * Math.Sin(l) }, + { Math.Cos(b) * Math.Sin(l), r * Math.Sin(b) * Math.Sin(l), r * Math.Cos(b) * Math.Cos(l) }, + { Math.Sin(b), -r * Math.Cos(b), 0 }}; + double[,] Velocity = { { dr }, { db }, { dl } }; + var C = MultiplyMatrix(J, Velocity); + dx = C[0, 0]; + dy = C[1, 0]; + dz = C[2, 0]; + + xyz[0] = x; + xyz[1] = y; + xyz[2] = z; + xyz[3] = dx; + xyz[4] = dy; + xyz[5] = -dz; + return xyz.ToArray(); + } + + /// + /// Convert Elliptic coordinate to cardinal coordinate. + /// This is kind of magic that I will never undersdand. + /// Directly translate from VSOP2013.f. + /// + /// planet + /// Elliptic Elements: a,l,k,h,q,p + /// Cardinal Heliocentric Coordinates + public static double[] ELLtoXYZ(VSOPBody body, double[] ell) + { + Span xyz = stackalloc double[6]; + double a, l, k, h, q, p; + a = ell[0]; + l = ell[1]; + k = ell[2]; + h = ell[3]; + q = ell[4]; + p = ell[5]; + Complex z, z1, z2, z3, zto, zteta; + double rgm, xfi, xki; + double u, ex, ex2, ex3, gl, gm, e, dl, rsa; + double xm, xr, xms, xmc, xn; + + //Initialization + rgm = Math.Sqrt(gmp[(int)body] + gmsol); + + //Computation + xfi = Math.Sqrt(1.0d - (k * k) - (h * h)); + xki = Math.Sqrt(1.0d - (q * q) - (p * p)); + u = 1.0d / (1.0d + xfi); + z = new Complex(k, h); + ex = z.Magnitude; + ex2 = ex * ex; + ex3 = ex2 * ex; + z1 = Complex.Conjugate(z); + gl = ell[1] % (Math.Tau); + gm = gl - Math.Atan2(h, k); + e = gl + (ex - 0.125d * ex3) * Math.Sin(gm) + + 0.5d * ex2 * Math.Sin(2.0d * gm) + + 0.375d * ex3 * Math.Sin(3.0d * gm); + + z2 = new Complex(0d, e); + zteta = Complex.Exp(z2); + while (true) + { + z3 = z1 * zteta; + dl = gl - e + z3.Imaginary; + rsa = 1.0d - z3.Real; + e += dl / rsa; + if (Math.Abs(dl) < Math.Pow(10, -15)) break; + } + + z1 = u * z * z3.Imaginary; + z2 = new Complex(z1.Imaginary, -z1.Real); + zto = (-z + zteta + z2) / rsa; + xm = p * zto.Real - q * zto.Imaginary; + xr = a * rsa; + + xyz[0] = xr * (zto.Real - 2.0d * p * xm); + xyz[1] = xr * (zto.Imaginary + 2.0d * q * xm); + xyz[2] = -2.0d * xr * xki * xm; + xms = a * (h + zto.Imaginary) / xfi; + xmc = a * (k + zto.Real) / xfi; + + xn = rgm / Math.Pow(a, 1.5d); + xyz[3] = xn * (((2.0d * p * p) - 1.0d) * xms + (2.0d * p * q * xmc)); + xyz[4] = xn * ((1.0d - (2.0d * q * q)) * xmc - (2.0d * p * q * xms)); + xyz[5] = xn * 2.0d * xki * (p * xms + q * xmc); + + return xyz.ToArray(); + } + + /// + /// Convert Elliptic coordinate to spherical coordinate. + /// + /// planet + /// Spherical Heliocentric Coordinates + public static double[] ELLtoLBR(VSOPBody body, double[] ell) + { + double[] xyz = ELLtoXYZ(body, ell); + return XYZtoLBR(xyz); + } + + /// + /// Another magic function + /// + /// Ecliptic Heliocentric Coordinates - Dynamical Frame J2000 + /// Equatorial Heliocentric Coordinates - ICRS Frame J2000 + public static double[] DynamicaltoICRS(double[] dynamical) + { + Span icrs = stackalloc double[6]; + + double dgrad = Math.PI / 180.0d; + double sdrad = Math.PI / 180.0d / 3600.0d; + double eps = (23.0d + 26.0d / 60.0d + 21.411360d / 3600.0d) * dgrad; + double phi = -0.051880d * sdrad; + + double Seps, Ceps, Sphi, Cphi; + Seps = Math.Sin(eps); + Ceps = Math.Cos(eps); + Sphi = Math.Sin(phi); + Cphi = Math.Cos(phi); + + //Rotation Matrix + double[,] R = new double[,] { {Cphi, -Sphi*Ceps, Sphi*Seps }, + {Sphi, Cphi*Ceps, -Cphi*Seps }, + {0, Seps, Ceps }}; + + // Vector for cardinal coordnate element + double[,] A = new double[,] { {dynamical[0]}, + {dynamical[1]}, + {dynamical[2]}}; + + double[,] C = MultiplyMatrix(R, A); + + dynamical[0] = C[0, 0]; + dynamical[1] = C[1, 0]; + dynamical[2] = C[2, 0]; + + A = new double[,] { {dynamical[3]}, + {dynamical[4]}, + {dynamical[5]}}; + + C = MultiplyMatrix(R, A); + dynamical[3] = C[0, 0]; + dynamical[4] = C[1, 0]; + dynamical[5] = C[2, 0]; + + return dynamical.ToArray(); + } + + /// + /// Another magic function + /// + /// Ecliptic Heliocentric Coordinates - Dynamical Frame J2000 + /// Equatorial Heliocentric Coordinates - ICRS Frame J2000 + public static double[] ICRStoDynamical(double[] icrs) + { + double[] dynamical = new double[6]; + + double dgrad = Math.PI / 180.0d; + double sdrad = Math.PI / 180.0d / 3600.0d; + double eps = (23.0d + 26.0d / 60.0d + 21.411360d / 3600.0d) * dgrad; + double phi = -0.051880d * sdrad; + + double Seps, Ceps, Sphi, Cphi; + Seps = Math.Sin(eps); + Ceps = Math.Cos(eps); + Sphi = Math.Sin(phi); + Cphi = Math.Cos(phi); + + //Reverse Matrix + double[,] R_1 = new double[,] {{ Cphi, Sphi, 0 }, + {-Sphi*Ceps, Cphi*Ceps, Seps }, + { Sphi*Seps, -Cphi*Seps, Ceps }}; + // Vector for cardinal coordnate element + double[,] A = new double[,] { {icrs[0]}, + {icrs[1]}, + {icrs[2]}}; + + double[,] C = MultiplyMatrix(R_1, A); + + dynamical[0] = C[0, 0]; + dynamical[1] = C[1, 0]; + dynamical[2] = C[2, 0]; + + A = new double[,] { {icrs[3]}, + {icrs[4]}, + {icrs[5]}}; + + C = MultiplyMatrix(R_1, A); + dynamical[3] = C[0, 0]; + dynamical[4] = C[1, 0]; + dynamical[5] = C[2, 0]; + + return dynamical.ToArray(); + } + } +} \ No newline at end of file diff --git a/VSOP2013/VSOP2013.csproj b/VSOP2013/VSOP2013.NET.csproj similarity index 72% rename from VSOP2013/VSOP2013.csproj rename to VSOP2013/VSOP2013.NET.csproj index 7208847..d975dff 100644 --- a/VSOP2013/VSOP2013.csproj +++ b/VSOP2013/VSOP2013.NET.csproj @@ -1,18 +1,20 @@ - + - net6.0 + net6.0;net7.0 x64 enable x64 - 1.1.0 + 1.1.2 VSOP2013.NET - KingsZNHONE + KingsZNHONE VSOP was developed and is maintained (updated with the latest data) by the scientists at the Bureau des Longitudes in Paris. - -VSOP2013,another version of VSOP, computed the positions of the planets directly at any moment, as well as their orbital elements with improved accuracy. +VSOP2013, computed the positions of the planets directly at any moment, as well as their orbital elements with improved accuracy. + https://github.com/kingsznhone/VSOP2013.NET.git + git https://github.com/kingsznhone/VSOP2013.NET - + True + VSOP2013.NET diff --git a/VSOP2013/VSOPResult/DynamicalELL.cs b/VSOP2013/VSOPResult/DynamicalELL.cs deleted file mode 100644 index 1cf5929..0000000 --- a/VSOP2013/VSOPResult/DynamicalELL.cs +++ /dev/null @@ -1,55 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Threading.Tasks; - -namespace VSOP2013.VSOPResult -{ - /// - /// Elliptic Elements - /// a (au), lambda (radian), k, h, q, p - /// Dynamical Frame J2000' - /// - public sealed class DynamicalELL : ResultBase - { - /// - /// a = semi-major axis (au) - /// - public double A { get => data_ell[0]; } - /// - /// l = mean longitude (rd) - /// - public double L { get => data_ell[1]; } - /// - /// k = e*cos(pi) (rd) - /// - public double K { get => data_ell[2]; } - /// - /// h = e*sin(pi) (rd) - /// - public double H { get => data_ell[3]; } - /// - /// q = sin(i/2)*cos(omega) (rd) - /// - public double Q { get => data_ell[4]; } - /// - /// p = sin(i/2)*sin(omega) (rd) - /// - public double P { get => data_ell[5]; } - - public DynamicalELL(VSOPBody body, VSOPTime time, double[] variables_ell) : base(body, time, variables_ell) { } - - - public static explicit operator DynamicalXYZ(DynamicalELL ELL) - { - return new DynamicalXYZ(ELL.Body, ELL.Time, ELL.data_ell); - } - - public static explicit operator ICRSXYZ(DynamicalELL ELL) - { - return new ICRSXYZ(ELL.Body, ELL.Time, ELL.data_ell); - } - - } -} diff --git a/VSOP2013/VSOPResult/DynamicalXYZ.cs b/VSOP2013/VSOPResult/DynamicalXYZ.cs deleted file mode 100644 index 9acedf9..0000000 --- a/VSOP2013/VSOPResult/DynamicalXYZ.cs +++ /dev/null @@ -1,60 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Numerics; -using System.Text; -using System.Threading.Tasks; - -namespace VSOP2013.VSOPResult -{ - /// - /// Ecliptic Heliocentric Coordinates - /// X,Y,Z (au) X'',Y'',Z'' (au/d) - /// Dynamical Frame J2000' - /// - public sealed class DynamicalXYZ : ResultBase - { - /// - /// Positions X (au) - /// - public double X { get => data_xyz[0]; } - /// - /// Positions Y (au) - /// - public double Y { get => data_xyz[1]; } - /// - /// Positions Z (au) - /// - public double Z { get => data_xyz[2]; } - /// - /// Velocities X'(au/d) - /// - public double dX { get => data_xyz[3]; } - /// - /// Velocities Y'(au/d) - /// - public double dY { get => data_xyz[4]; } - /// - /// Velocities Z'(au/d) - /// - public double dZ { get => data_xyz[5]; } - - private double[] data_xyz { get; set; } - - public DynamicalXYZ(VSOPBody body, VSOPTime time, double[] variables_ell) : base(body, time, variables_ell) - { - data_xyz = ELLtoXYZ(body,data_ell); - } - - public static explicit operator DynamicalELL(DynamicalXYZ xyz) - { - return new DynamicalELL(xyz.Body, xyz.Time, xyz.data_ell); - } - - public static explicit operator ICRSXYZ(DynamicalXYZ xyz) - { - return new ICRSXYZ(xyz.Body, xyz.Time, xyz.data_ell); - } - - } -} diff --git a/VSOP2013/VSOPResult/ICRSXYZ.cs b/VSOP2013/VSOPResult/ICRSXYZ.cs deleted file mode 100644 index e37656a..0000000 --- a/VSOP2013/VSOPResult/ICRSXYZ.cs +++ /dev/null @@ -1,55 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Threading.Tasks; - -namespace VSOP2013.VSOPResult -{ - //Equatorial Heliocentric Coordinates: - //X,Y,Z (au) X'',Y'',Z'' (au/d) - //ICRS Frame J2000 - public sealed class ICRSXYZ : ResultBase - { - /// - /// Positions X (au) - /// - public double X { get => data_ICRS[0]; } - /// - /// Positions Y (au) - /// - public double Y { get => data_ICRS[1]; } - /// - /// Positions Z (au) - /// - public double Z { get => data_ICRS[2]; } - /// - /// Velocities X'(au/d) - /// - public double dX { get => data_ICRS[3]; } - /// - /// Velocities Y'(au/d) - /// - public double dY { get => data_ICRS[4]; } - /// - /// Velocities Z'(au/d) - /// - public double dZ { get => data_ICRS[5]; } - - private double[] data_ICRS; - public ICRSXYZ(VSOPBody body, VSOPTime time, double[] variables_ell) : base(body, time, variables_ell) - { - data_ICRS = DynamicaltoICRS(ELLtoXYZ(body,data_ell)); - } - - public static explicit operator DynamicalELL(ICRSXYZ xyz) - { - return new DynamicalELL(xyz.Body, xyz.Time, xyz.data_ell); - } - - public static explicit operator DynamicalXYZ(ICRSXYZ xyz) - { - return new DynamicalXYZ(xyz.Body, xyz.Time, xyz.data_ell); - } - } -} diff --git a/VSOP2013/VSOPResult/ResultBase.cs b/VSOP2013/VSOPResult/ResultBase.cs deleted file mode 100644 index 428f7f9..0000000 --- a/VSOP2013/VSOPResult/ResultBase.cs +++ /dev/null @@ -1,152 +0,0 @@ -using System; -using System.Collections.Generic; -using System.ComponentModel; -using System.Linq; -using System.Numerics; -using System.Text; -using System.Threading.Tasks; - -namespace VSOP2013.VSOPResult -{ - public abstract class ResultBase - { - protected readonly double[] gmp = {4.9125474514508118699e-11d, - 7.2434524861627027000e-10d, - 8.9970116036316091182e-10d, - 9.5495351057792580598e-11d, - 2.8253458420837780000e-07d, - 8.4597151856806587398e-08d, - 1.2920249167819693900e-08d, - 1.5243589007842762800e-08d, - 2.1886997654259696800e-12d}; - - protected readonly double gmsol = 2.9591220836841438269e-04d; - - public VSOPBody Body { get; private set; } - public VSOPTime Time { get; private set; } - - //Elliptic Elements - //a (au), lambda (radian), k, h, q, p - //Dynamical Frame J2000' - protected double[] data_ell; - - public ResultBase(VSOPBody body, VSOPTime time, double[] variables) - { - Body = body; - Time = time; - data_ell = variables; - } - - /// - /// This is kind of magic that I can't undersdand - /// translate from FORTRAN code - /// - /// planet - /// Elliptic Elements - /// Ecliptic Heliocentric Coordinates - protected double[] ELLtoXYZ(VSOPBody body, double[] ELL) - { - - double[] xyz; - Complex z, z1, z2, z3, zto, zteta; - double rgm, xfi, xki; - double u, ex, ex2, ex3, gl, gm, e, dl, rsa; - double xm, xr, xms, xmc, xn; - - //Initialization - rgm = Math.Sqrt(gmp[(int)body] + gmsol); - xyz = new double[6]; - - //Computation - xfi = Math.Sqrt(1.0d - (ELL[2] * ELL[2]) - (ELL[3] * ELL[3])); - xki = Math.Sqrt(1.0d - (ELL[4] * ELL[4]) - (ELL[5] * ELL[5])); - u = 1.0d / (1.0d + xfi); - z = new Complex(ELL[2], ELL[3]); - ex = z.Magnitude; - ex2 = ex * ex; - ex3 = ex2 * ex; - z1 = Complex.Conjugate(z); - gl = ELL[1] % (2 * Math.PI); - gm = gl - Math.Atan2(ELL[3], ELL[2]); - e = gl + (ex - 0.125d * ex3) * Math.Sin(gm) - + 0.5d * ex2 * Math.Sin(2.0d * gm) - + 0.375d * ex3 * Math.Sin(3.0d * gm); - - z2 = new Complex(0d, e); - zteta = Complex.Exp(z2); - while (true) - { - z3 = z1 * zteta; - dl = gl - e + z3.Imaginary; - rsa = 1.0d - z3.Real; - e = e + dl / rsa; - if (Math.Abs(dl) < Math.Pow(10, -15)) break; - } - - z1 = u * z * z3.Imaginary; - z2 = new Complex(z1.Imaginary, -z1.Real); - zto = (-z + zteta + z2) / rsa; - xm = ELL[5] * zto.Real - ELL[4] * zto.Imaginary; - xr = ELL[0] * rsa; - - xyz[0] = xr * (zto.Real - 2.0d * ELL[5] * xm); - xyz[1] = xr * (zto.Imaginary + 2.0d * ELL[4] * xm); - xyz[2] = -2.0d * xr * xki * xm; - xms = ELL[0] * (ELL[3] + zto.Imaginary) / xfi; - xmc = ELL[0] * (ELL[2] + zto.Real) / xfi; - - xn = rgm / Math.Pow(ELL[0], 1.5d); - xyz[3] = xn * ((2.0d * ELL[5] * ELL[5] - 1.0d) * xms + 2.0d * ELL[5] * ELL[4] * xmc); - xyz[4] = xn * ((1.0d - 2.0d * ELL[4] * ELL[4]) * xmc - 2.0d * ELL[5] * ELL[4] * xms); - xyz[5] = 2.0d * xn * xki * (ELL[5] * xms + ELL[4] * xmc); - - return xyz; - } - - - /// - /// Another magic function - /// - /// Ecliptic Heliocentric Coordinates - Dynamical Frame J2000 - /// Equatorial Heliocentric Coordinates - ICRS Frame J2000 - protected double[] DynamicaltoICRS(double[] xyz) - { - double[] icrs = new double[6]; - //Rotation Matrix - double[,] Rot = new double[3, 3]; - double pi = Math.PI; - double dgrad = pi / 180.0d; - double sdrad = pi / 180.0d / 3600.0d; - double eps = (23.0d + 26.0d / 60.0d + 21.411360d / 3600.0d) * dgrad; - double phi = -0.051880d * sdrad; - - double seps, ceps, sphi, cphi; - (seps, ceps) = Math.SinCos(eps); - (sphi, cphi) = Math.SinCos(phi); - - Rot[0, 0] = cphi; - Rot[0, 1] = -sphi * ceps; - Rot[0, 2] = sphi * seps; - Rot[1, 0] = sphi; - Rot[1, 1] = cphi * ceps; - Rot[1, 2] = -cphi * seps; - Rot[2, 0] = 0.0d; - Rot[2, 1] = seps; - Rot[2, 2] = ceps; - - //Computation - for (int i = 0; i < 3; i++) - { - icrs[i] = 0.0d; - icrs[i + 3] = 0.0d; - for (int j = 0; j < 3; j++) - { - icrs[i] = icrs[i] + Rot[i, j] * xyz[j]; - icrs[i + 3] = icrs[i + 3] + Rot[i, j] * xyz[j + 3]; - } - } - return icrs; - } - - } -} diff --git a/VSOP2013/VSOPResult/VSOPResult.cs b/VSOP2013/VSOPResult/VSOPResult.cs new file mode 100644 index 0000000..0246d4f --- /dev/null +++ b/VSOP2013/VSOPResult/VSOPResult.cs @@ -0,0 +1,42 @@ +namespace VSOP2013 +{ + public enum CoordinatesType + { + Elliptic = 0, + Rectangular = 1, + Spherical = 2 + } + + public enum InertialFrame + { + /// + /// Ecliptic Heliocentric Dynamical Frame + /// + Dynamical = 0, + + /// + /// Equatorial Heliocentric ICRS Frame + /// + ICRS = 1 + } + + public enum TimeFrameReference + { + EclipticJ2000 = 0, + } + + public abstract class VSOPResult + { + public abstract VSOPBody Body { get; } + public abstract CoordinatesType CoordinatesType { get; } + public abstract InertialFrame InertialFrame { get; } + public abstract TimeFrameReference TimeFrameReference { get; } + + public abstract VSOPTime Time { get; } + + //Elliptic Elements + //a (au), lambda (radian), k, h, q, p + //Dynamical Frame J2000' + public abstract double[] Variables_ELL { get; set; } + } +} \ No newline at end of file diff --git a/VSOP2013/VSOPResult/VSOPResult_ELL.cs b/VSOP2013/VSOPResult/VSOPResult_ELL.cs new file mode 100644 index 0000000..ffddc40 --- /dev/null +++ b/VSOP2013/VSOPResult/VSOPResult_ELL.cs @@ -0,0 +1,69 @@ +namespace VSOP2013 +{ + /// + /// Elliptic Elements + /// a (au), lambda (radian), k, h, q, p + /// Dynamical Frame J2000' + /// + public sealed class VSOPResult_ELL : VSOPResult + { + public override VSOPBody Body { get; } + + public override CoordinatesType CoordinatesType => CoordinatesType.Elliptic; + + public override InertialFrame InertialFrame => InertialFrame.Dynamical; + + public override TimeFrameReference TimeFrameReference => TimeFrameReference.EclipticJ2000; + + public override VSOPTime Time { get; } + + public override double[] Variables_ELL { get; set; } + + public VSOPResult_ELL(VSOPBody body, VSOPTime time, double[] ell) + { + Body = body; + Time = time; + Variables_ELL = ell; + } + + /// + /// a = semi-major axis (au) + /// + public double a { get => Variables_ELL[0]; } + + /// + /// l = mean longitude (rd) + /// + public double l { get => Variables_ELL[1]; } + + /// + /// k = e*cos(pi) (rd) + /// + public double k { get => Variables_ELL[2]; } + + /// + /// h = e*sin(pi) (rd) + /// + public double h { get => Variables_ELL[3]; } + + /// + /// q = sin(i/2)*cos(omega) (rd) + /// + public double q { get => Variables_ELL[4]; } + + /// + /// p = sin(i/2)*sin(omega) (rd) + /// + public double p { get => Variables_ELL[5]; } + + public static explicit operator VSOPResult_XYZ(VSOPResult_ELL ELL) + { + return new VSOPResult_XYZ(ELL.Body, ELL.Time, ELL.Variables_ELL); + } + + public static explicit operator VSOPResult_LBR(VSOPResult_ELL ELL) + { + return new VSOPResult_LBR(ELL.Body, ELL.Time, ELL.Variables_ELL); + } + } +} \ No newline at end of file diff --git a/VSOP2013/VSOPResult/VSOPResult_LBR.cs b/VSOP2013/VSOPResult/VSOPResult_LBR.cs new file mode 100644 index 0000000..beed640 --- /dev/null +++ b/VSOP2013/VSOPResult/VSOPResult_LBR.cs @@ -0,0 +1,104 @@ +namespace VSOP2013 +{ + public class VSOPResult_LBR : VSOPResult + { + public override VSOPBody Body { get; } + + public override CoordinatesType CoordinatesType => CoordinatesType.Spherical; + + private InertialFrame _inertialFrame; + public override InertialFrame InertialFrame + { get { return _inertialFrame; } } + + public override TimeFrameReference TimeFrameReference => TimeFrameReference.EclipticJ2000; + + public override VSOPTime Time { get; } + + public override double[] Variables_ELL { get; set; } + + public double[] Variables_LBR { get; set; } + + /// + /// + /// + /// + /// + /// Elliptic elements + public VSOPResult_LBR(VSOPBody body, VSOPTime time, double[] ell) + { + Body = body; + Time = time; + _inertialFrame = InertialFrame.Dynamical; + Variables_ELL = ell; + Variables_LBR = Utility.ELLtoLBR(body, ell); + } + + /// + /// a = semi-major axis (au) + /// + public double l { get => Variables_LBR[0]; } + + /// + /// l = mean longitude (rd) + /// + public double b { get => Variables_LBR[1]; } + + /// + /// k = e*cos(pi) (rd) + /// + public double r { get => Variables_LBR[2]; } + + /// + /// h = e*sin(pi) (rd) + /// + public double dl { get => Variables_LBR[3]; } + + /// + /// q = sin(i/2)*cos(omega) (rd) + /// + public double db { get => Variables_LBR[4]; } + + /// + /// p = sin(i/2)*sin(omega) (rd) + /// + public double dr { get => Variables_LBR[5]; } + + public void ToICRSFrame() + { + if (InertialFrame == InertialFrame.Dynamical) + { + _inertialFrame = InertialFrame.ICRS; + Variables_LBR = Utility.DynamicaltoICRS(Utility.ELLtoXYZ(this.Body, Variables_ELL)); + } + } + + public void ToDynamicalFrame() + { + if (InertialFrame == InertialFrame.ICRS) + { + _inertialFrame = InertialFrame.Dynamical; + Variables_LBR = Utility.ELLtoXYZ(this.Body, Variables_ELL); + } + } + + public static explicit operator VSOPResult_XYZ(VSOPResult_LBR LBR) + { + double[] v = Utility.ELLtoLBR(LBR.Body, LBR.Variables_ELL); + var result = new VSOPResult_XYZ(LBR.Body, LBR.Time, v); + if (LBR.InertialFrame == InertialFrame.Dynamical) + { + return result; + } + else + { + result.ToICRSFrame(); + return result; + } + } + + public static explicit operator VSOPResult_ELL(VSOPResult_LBR LBR) + { + return new VSOPResult_ELL(LBR.Body, LBR.Time, LBR.Variables_ELL); + } + } +} \ No newline at end of file diff --git a/VSOP2013/VSOPResult/VSOPResult_XYZ.cs b/VSOP2013/VSOPResult/VSOPResult_XYZ.cs new file mode 100644 index 0000000..f94b2c9 --- /dev/null +++ b/VSOP2013/VSOPResult/VSOPResult_XYZ.cs @@ -0,0 +1,98 @@ +namespace VSOP2013 +{ + public class VSOPResult_XYZ : VSOPResult + { + public override VSOPBody Body { get; } + + public override CoordinatesType CoordinatesType => CoordinatesType.Rectangular; + + private InertialFrame _inertialFrame; + public override InertialFrame InertialFrame + { get { return _inertialFrame; } } + + public override TimeFrameReference TimeFrameReference => TimeFrameReference.EclipticJ2000; + + public override VSOPTime Time { get; } + + public override double[] Variables_ELL { get; set; } + + public double[] Variables_XYZ { get; set; } + + public VSOPResult_XYZ(VSOPBody body, VSOPTime time, double[] ell) + { + Body = body; + Time = time; + _inertialFrame = InertialFrame.Dynamical; + Variables_ELL = ell; + Variables_XYZ = Utility.ELLtoXYZ(body, ell); + } + + /// + /// a = semi-major axis (au) + /// + public double x { get => Variables_XYZ[0]; } + + /// + /// l = mean longitude (rd) + /// + public double y { get => Variables_XYZ[1]; } + + /// + /// k = e*cos(pi) (rd) + /// + public double z { get => Variables_XYZ[2]; } + + /// + /// h = e*sin(pi) (rd) + /// + public double dx { get => Variables_XYZ[3]; } + + /// + /// q = sin(i/2)*cos(omega) (rd) + /// + public double dy { get => Variables_XYZ[4]; } + + /// + /// p = sin(i/2)*sin(omega) (rd) + /// + public double dz { get => Variables_XYZ[5]; } + + public void ToICRSFrame() + { + if (InertialFrame == InertialFrame.Dynamical) + { + _inertialFrame = InertialFrame.ICRS; + Variables_XYZ = Utility.DynamicaltoICRS(Variables_XYZ); + } + } + + public void ToDynamicalFrame() + { + if (InertialFrame == InertialFrame.ICRS) + { + _inertialFrame = InertialFrame.Dynamical; + Variables_XYZ = Utility.ICRStoDynamical(Variables_XYZ); + } + } + + public static explicit operator VSOPResult_LBR(VSOPResult_XYZ XYZ) + { + double[] v = Utility.ELLtoLBR(XYZ.Body, XYZ.Variables_ELL); + var result = new VSOPResult_LBR(XYZ.Body, XYZ.Time, v); + if (XYZ.InertialFrame == InertialFrame.Dynamical) + { + return result; + } + else + { + result.ToICRSFrame(); + return result; + } + } + + public static explicit operator VSOPResult_ELL(VSOPResult_XYZ XYZ) + { + return new VSOPResult_ELL(XYZ.Body, XYZ.Time, XYZ.Variables_ELL); + } + } +} \ No newline at end of file diff --git a/VSOP2013/VSOPTime.cs b/VSOP2013/VSOPTime.cs index ab6293c..7c27ea8 100644 --- a/VSOP2013/VSOPTime.cs +++ b/VSOP2013/VSOPTime.cs @@ -1,7 +1,4 @@ -using System; -using System.Collections.Generic; - -namespace VSOP2013 +namespace VSOP2013 { public enum TimeFrame { @@ -18,25 +15,30 @@ public class VSOPTime /// /// UTC:Coordinated Universal Time /// - public DateTime UTC { get { return _dt; } } + public DateTime UTC => _dt; /// /// TAI:International Atomic Time /// - public DateTime TAI { get { return ChangeFrame(_dt, TimeFrame.TAI); } } + public DateTime TAI => ChangeFrame(_dt, TimeFrame.UTC, TimeFrame.TAI); /// /// TT :Terrestrial Time (aka. TDT) /// - public DateTime TT { get { return ChangeFrame(_dt, TimeFrame.TT); } } + public DateTime TT => ChangeFrame(_dt, TimeFrame.UTC, TimeFrame.TT); /// /// TDB:Barycentric Dynamical Time /// - public DateTime TDB { get { return ChangeFrame(_dt, TimeFrame.TDB); } } + public DateTime TDB => ChangeFrame(_dt, TimeFrame.UTC, TimeFrame.TDB); + + /// + /// Get JD from TDB + /// + public double J2000 => VSOPTime.ToJ2000(TDB); - private List> UpGradeFuncs; - private List> DownGradeFuncs; + private readonly List> UpGradeFuncs; + private readonly List> DownGradeFuncs; public VSOPTime(DateTime UTC) { @@ -51,107 +53,98 @@ public VSOPTime(DateTime UTC) DownGradeFuncs = new List>( new Func[] { TAItoUTC, TTtoTAI, TDBtoTT }); - } /// - /// + /// /// /// DateTime in UTC Frame /// /// - public DateTime ChangeFrame(DateTime dt, TimeFrame TargetFrame) + public DateTime ChangeFrame(DateTime dt, TimeFrame SourceFrame, TimeFrame TargetFrame) { - var buffer = TimeFrame.UTC; - while (buffer != TargetFrame) + while (SourceFrame != TargetFrame) { - if (TargetFrame > buffer) + if (TargetFrame > SourceFrame) { - dt = UpGradeFuncs[(int)buffer](dt); - buffer += 1; + dt = UpGradeFuncs[(int)SourceFrame](dt); + SourceFrame += 1; } - else if (TargetFrame < buffer) + else if (TargetFrame < SourceFrame) { - dt = DownGradeFuncs[(int)buffer - 1](dt); - buffer -= 1; + dt = DownGradeFuncs[(int)SourceFrame - 1](dt); + SourceFrame -= 1; } } return dt; } - #region UTC To TDB - public static DateTime UTCtoTAI(DateTime UTC) + private static DateTime UTCtoTAI(DateTime UTC) { return UTC.AddSeconds(37); } - public static DateTime TAItoTT(DateTime TAI) + private static DateTime TAItoTT(DateTime TAI) { return TAI.AddSeconds(32.184); } - public static DateTime TTtoTDB(DateTime TT) + private static DateTime TTtoTDB(DateTime TT) { //Error btw TT&TDB is so small that can be ignored. return TT; } - #endregion + + #endregion UTC To TDB #region TDB to UTC - public static DateTime TDBtoTT(DateTime TDB) + + private static DateTime TDBtoTT(DateTime TDB) { return TDB; } - public static DateTime TTtoTAI(DateTime TT) + private static DateTime TTtoTAI(DateTime TT) { return TT.AddSeconds(-32.184); } - public static DateTime TAItoUTC(DateTime TAI) + private static DateTime TAItoUTC(DateTime TAI) { return TAI.AddSeconds(-37); } - #endregion + #endregion TDB to UTC #region JulianDate Convert /// - /// + /// /// /// Barycentric Dynamical Time - /// JulianDate in double - public static double ToJulianDate(DateTime dt) - { - return dt.ToOADate() + 2415018.5d; - } - - /// - /// - /// - /// Barycentric Dynamical Time - /// JulianDate 2000 - public static double ToJulianDate2000(DateTime dt) + /// JulianDate 2000 in double + public static double ToJ2000(DateTime dt) { //Input TDB double j2000 = 2451545.0d; - //OADate + JD - j2000(JD) = TDB from J2000 + // TDB from J2000 = OADate + JD - j2000(JD) return dt.ToOADate() + 2415018.5d - j2000; } /// - /// + /// /// - /// + /// /// Datetime in TDB Frame - public static DateTime FromJulianDate(double JD) + public static DateTime FromJ2000(double J2000) { - return DateTime.FromOADate(JD - 2415018.5); + double j2000 = 2451545.0d; + + return DateTime.FromOADate(J2000 + j2000 - 2415018.5).ToUniversalTime(); } - #endregion + #endregion JulianDate Convert } -} +} \ No newline at end of file