diff --git a/DataConverter/DataConverter.csproj b/DataConverter/DataConverter.csproj new file mode 100644 index 0000000..2f3dd48 --- /dev/null +++ b/DataConverter/DataConverter.csproj @@ -0,0 +1,37 @@ + + + + Exe + net6.0 + enable + enable + x64 + x64 + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/VSOP2013/DataReader.cs b/DataConverter/DataReader.cs similarity index 92% rename from VSOP2013/DataReader.cs rename to DataConverter/DataReader.cs index 44c947b..58fcf7a 100644 --- a/VSOP2013/DataReader.cs +++ b/DataConverter/DataReader.cs @@ -1,15 +1,13 @@ -using System.Diagnostics; -using System.Reflection; +using System.Reflection; -namespace VSOP2013 +namespace VSOP2013.DataConverter { public static class DataReader { - /// /// Mean Longitude J2000 (radian) /// - static readonly double[] ci0 = + private static readonly double[] ci0 = { 0.4402608631669000e1d, 0.3176134461576000e1d, @@ -33,7 +31,7 @@ public static class DataReader /// /// //Mean Motions in longitude (radian/cy) /// - static readonly double[] ci1 = + private static readonly double[] ci1 = { 0.2608790314068555e5d, 0.1021328554743445e5d, @@ -56,8 +54,6 @@ public static class DataReader public static List ReadData() { - Stopwatch sw = new Stopwatch(); - sw.Start(); List VSOP2013DATA = new List(); for (int ip = 0; ip < 9; ip++) { @@ -85,13 +81,8 @@ public static List ReadData() } } - sw.Stop(); GC.Collect(); GC.WaitForPendingFinalizers(); - double ticks = sw.ElapsedTicks; - double Freq = Stopwatch.Frequency; - double milliseconds = (ticks / Freq) * 1000; - Console.WriteLine($"Data Analyze OK...Elapsed milliseconds: {milliseconds} ms"); return VSOP2013DATA; } @@ -105,8 +96,9 @@ private static void ReadPlanet(PlanetTable Planet, int ip) //var debug = Assembly.GetExecutingAssembly().GetManifestResourceNames(); var assembly = Assembly.GetExecutingAssembly(); - string datafilename = $"VSOP2013.Resources.VSOP2013p{ip + 1}.dat"; + string datafilename = $"DataConverter.Resources.VSOP2013p{ip + 1}.dat"; Stream s = assembly.GetManifestResourceStream(datafilename); + sr = new StreamReader(s); while ((line = sr.ReadLine()) != null) { @@ -123,7 +115,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; @@ -203,7 +194,6 @@ private static Term ReadTerm(string line, ref Term T) //T.iphi = Bufferiphi; - ci = Convert.ToDouble(line.Substring(lineptr, 20).Trim()); lineptr += 20; lineptr++; @@ -232,7 +222,6 @@ private static Term ReadTerm(string line, ref Term T) private static PowerTable[] DataPruning(PowerTable[] tables) { - for (int i = 0; i < tables.Length; i++) { if (tables[i].Terms is null) diff --git a/DataConverter/Program.cs b/DataConverter/Program.cs new file mode 100644 index 0000000..d7116ac --- /dev/null +++ b/DataConverter/Program.cs @@ -0,0 +1,87 @@ +using MessagePack; +using System.Diagnostics; +using System.Runtime.Serialization.Formatters.Binary; + +namespace VSOP2013.DataConverter +{ + internal class Program + { + private static void Main(string[] args) + { + var lz4Options = MessagePackSerializerOptions.Standard.WithCompression(MessagePackCompression.Lz4Block); + + #region Read Original Data + + Stopwatch sw = new Stopwatch(); + sw.Start(); + + List VSOP2013DATA = DataReader.ReadData(); + + 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(); + + #endregion Read Original Data + + #region Dump Data + + string OutputDirPath = Directory.GetCurrentDirectory() + @"\Data"; + if (!Directory.Exists(OutputDirPath)) + { + Directory.CreateDirectory(OutputDirPath); + } + DirectoryInfo OutputDir = new DirectoryInfo(Directory.GetCurrentDirectory() + @"\Data"); + + sw.Restart(); + + ParallelLoopResult result = Parallel.For(0, 9, ip => + { + string filename = Path.Combine(OutputDir.FullName, string.Format("VSOP2013_{0}.BIN", VSOP2013DATA[ip].body.ToString())); + if (File.Exists(filename)) + { + File.Delete(filename); + } + using (FileStream fs = new FileStream(filename, FileMode.OpenOrCreate)) + { + BinaryFormatter bf = new BinaryFormatter(); + MessagePackSerializer.Serialize(fs, VSOP2013DATA[ip], lz4Options); + } + }); + + 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(); + + 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)); + } + }); + + sw.Stop(); + ticks = sw.ElapsedTicks; + milliseconds = (ticks / Freq) * 1000; + Console.WriteLine($"Data Reload Test OK. Elapsed: {milliseconds}ms"); + Console.ReadLine(); + + #endregion Test + } + } +} \ No newline at end of file diff --git a/VSOP2013/Resources/README.pdf b/DataConverter/Resources/README.pdf similarity index 100% rename from VSOP2013/Resources/README.pdf rename to DataConverter/Resources/README.pdf diff --git a/VSOP2013/Resources/VSOP2013-secular.dat b/DataConverter/Resources/VSOP2013-secular.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013-secular.dat rename to DataConverter/Resources/VSOP2013-secular.dat diff --git a/VSOP2013/Resources/VSOP2013.ctl b/DataConverter/Resources/VSOP2013.ctl similarity index 91% rename from VSOP2013/Resources/VSOP2013.ctl rename to DataConverter/Resources/VSOP2013.ctl index 0dc18cb..130216a 100644 --- a/VSOP2013/Resources/VSOP2013.ctl +++ b/DataConverter/Resources/VSOP2013.ctl @@ -1,11 +1,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 0.3870979635 6.2605163414 0.0452614144 0.2005681842 0.0405434321 0.0457750937 0.3493879042 -0.1615770401 -0.0453430160 0.0063187162 0.0268317850 0.0016062487 0.3493878714 -0.1302077267 -0.1058730361 0.0063187222 0.0239787530 0.0121467712 @@ -52,11 +52,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 0.7233268460 3.0850544129 -0.0045575903 0.0051297811 0.0066724242 0.0288666759 -0.7178452078 0.0333612116 0.0419440586 -0.0010551049 -0.0202953755 -0.0002102008 -0.7178452043 0.0139241146 0.0517532468 -0.0010551096 -0.0185370311 -0.0082658890 @@ -103,11 +103,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 1.0000096358 4.8188777642 -0.0036242360 0.0163481300 0.0001248730 -0.0000110805 0.1117529336 -1.0104931629 -0.0002498901 0.0168189379 0.0018288000 0.0000008295 0.1117527004 -0.9270100498 -0.4021802015 0.0168189384 0.0016775571 0.0007282156 @@ -154,11 +154,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 1.5236841626 4.7846953863 0.0850047012 -0.0386037157 0.0104503533 0.0124027403 -0.1474458094 -1.4589014057 -0.0268451993 0.0144622332 -0.0002104338 -0.0003632842 -0.1474461434 -1.3278375334 -0.6049474049 0.0144622332 -0.0000485668 -0.0004170125 @@ -205,11 +205,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 5.2027787025 5.4273729856 0.0474307493 0.0119309207 -0.0020314597 0.0112167899 2.9837893541 -4.1334108921 -0.0501531373 0.0060327773 0.0047801335 -0.0001547883 2.9837884053 -3.7723816270 -1.6901903627 0.0060327784 0.0044472568 0.0017594117 @@ -256,11 +256,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 9.5223420509 2.6372558976 -0.0013261890 0.0534078734 -0.0088168973 0.0198122559 -8.5151107303 3.6958970183 0.2724266220 -0.0025150365 -0.0051284514 0.0001902254 -8.5151099046 3.2825565780 1.7200893604 -0.0025150377 -0.0047809292 -0.0018654516 @@ -307,11 +307,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 19.2176223319 3.5628257887 -0.0473181207 0.0060839286 0.0018734886 0.0065084828 -16.4159077405 -8.4156336905 0.1821646518 0.0017679079 -0.0036818680 -0.0000368112 -16.4159097008 -7.7936503250 -3.1804126500 0.0017679071 -0.0033634059 -0.0014983360 @@ -358,11 +358,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 30.0385847963 1.1425611229 0.0069759741 0.0045558233 -0.0102865506 0.0115258164 12.1323737542 27.2358156923 -0.8402970634 -0.0028808365 0.0013000400 0.0000396763 12.1323801234 25.3226220555 10.0628233249 -0.0028808362 0.0011769819 0.0005535283 @@ -409,11 +409,11 @@ 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 - 3/ Equatorial Heliocentric Coordinates: X,Y,Z (au) X',Y',Z' (au/d) - ICRS Frame J2000 + 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 +3/ Equatorial Heliocentric Coordinates: X, Y, Z(au) X',Y',Z' (au/d) - ICRS Frame J2000 - Julian Date JD 2411545.0 (TDB) +Julian Date JD 2411545.0 (TDB) 39.4227219159 1.3910260961 -0.1777508423 -0.1751661709 -0.0515906603 0.1400924565 17.6212944088 43.6567983387 -9.7715435955 -0.0020295942 0.0006425559 0.0005199070 17.6213054610 43.9412232438 8.4004533071 -0.0020295941 0.0003827270 0.0007325993 diff --git a/VSOP2013/Resources/VSOP2013.f b/DataConverter/Resources/VSOP2013.f similarity index 100% rename from VSOP2013/Resources/VSOP2013.f rename to DataConverter/Resources/VSOP2013.f diff --git a/VSOP2013/Resources/VSOP2013p1.dat b/DataConverter/Resources/VSOP2013p1.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p1.dat rename to DataConverter/Resources/VSOP2013p1.dat diff --git a/VSOP2013/Resources/VSOP2013p2.dat b/DataConverter/Resources/VSOP2013p2.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p2.dat rename to DataConverter/Resources/VSOP2013p2.dat diff --git a/VSOP2013/Resources/VSOP2013p3.dat b/DataConverter/Resources/VSOP2013p3.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p3.dat rename to DataConverter/Resources/VSOP2013p3.dat diff --git a/VSOP2013/Resources/VSOP2013p4.dat b/DataConverter/Resources/VSOP2013p4.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p4.dat rename to DataConverter/Resources/VSOP2013p4.dat diff --git a/VSOP2013/Resources/VSOP2013p5.dat b/DataConverter/Resources/VSOP2013p5.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p5.dat rename to DataConverter/Resources/VSOP2013p5.dat diff --git a/VSOP2013/Resources/VSOP2013p6.dat b/DataConverter/Resources/VSOP2013p6.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p6.dat rename to DataConverter/Resources/VSOP2013p6.dat diff --git a/VSOP2013/Resources/VSOP2013p7.dat b/DataConverter/Resources/VSOP2013p7.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p7.dat rename to DataConverter/Resources/VSOP2013p7.dat diff --git a/VSOP2013/Resources/VSOP2013p8.dat b/DataConverter/Resources/VSOP2013p8.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p8.dat rename to DataConverter/Resources/VSOP2013p8.dat diff --git a/VSOP2013/Resources/VSOP2013p9.dat b/DataConverter/Resources/VSOP2013p9.dat similarity index 100% rename from VSOP2013/Resources/VSOP2013p9.dat rename to DataConverter/Resources/VSOP2013p9.dat diff --git a/DataConverter/copy.bat b/DataConverter/copy.bat new file mode 100644 index 0000000..3d778f8 --- /dev/null +++ b/DataConverter/copy.bat @@ -0,0 +1 @@ +copy bin\Debug\net6.0\Data ..\VSOP2013\Resources \ No newline at end of file diff --git a/Demo/Demo.csproj b/Demo/Demo.csproj index 300f8fe..1ad34c1 100644 --- a/Demo/Demo.csproj +++ b/Demo/Demo.csproj @@ -3,7 +3,7 @@ Exe net6.0 - AnyCPU;x64 + x64 diff --git a/Demo/Program.cs b/Demo/Program.cs index 786709a..502d51c 100644 --- a/Demo/Program.cs +++ b/Demo/Program.cs @@ -4,6 +4,8 @@ using System.IO; using System.Threading.Tasks; using VSOP2013; +using VSOP2013.VSOPResult; + namespace Demo { class Program @@ -45,9 +47,9 @@ static void Main(string[] args) Console.WriteLine("Press Enter To Start Substitution..."); Console.ReadLine(); - VSOPResult[] ResultAll = vsop.CalcAllPlanet(vTime); + DynamicalELL[] ResultAll = vsop.CalcAllPlanet(vTime); - foreach(VSOPResult result in ResultAll) + foreach(DynamicalELL result in ResultAll) { FormattedPrint(result, vTime); } @@ -60,7 +62,7 @@ static void Main(string[] args) Console.ReadLine(); } - public static void FormattedPrint(VSOPResult Result, VSOPTime vtime) + public static void FormattedPrint(DynamicalELL Result, VSOPTime vtime) { Console.WriteLine("==============================================================="); Console.WriteLine("PLANETARY EPHEMERIS VSOP2013"); @@ -69,28 +71,28 @@ public static void FormattedPrint(VSOPResult Result, VSOPTime vtime) 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.DynamicalELL[0])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "mean longitude (rd)", Result.DynamicalELL[1])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "k = e*cos(pi) (rd)", Result.DynamicalELL[2])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "h = e*sin(pi) (rd)", Result.DynamicalELL[3])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "q = sin(i/2)*cos(omega) (rd)", Result.DynamicalELL[4])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "p = sin(i/2)*sin(omega) (rd)", Result.DynamicalELL[5])); + 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)", Result.DynamicalXYZ[0])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Y (au)", Result.DynamicalXYZ[1])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Z (au)", Result.DynamicalXYZ[2])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities X'(au/d)", Result.DynamicalXYZ[3])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Y'(au/d)", Result.DynamicalXYZ[4])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Z'(au/d)", Result.DynamicalXYZ[5])); + 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)", Result.ICRSXYZ[0])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Y (au)", Result.ICRSXYZ[1])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Positions Z (au)", Result.ICRSXYZ[2])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities X'(au/d)", Result.ICRSXYZ[3])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Y'(au/d)", Result.ICRSXYZ[4])); - Console.WriteLine(String.Format("{0,-30} : {1,-30}", "Velocities Z'(au/d)", Result.ICRSXYZ[5])); + 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(); } @@ -102,44 +104,17 @@ public static void PerformanceTest(int cycle) Stopwatch sw = new Stopwatch(); sw.Start(); int completedCycle = 0; - var result = Parallel.For(0, cycle, (i) => - { - { - VSOPResult[] ResultAll = vsop.CalcAllPlanet(vTime); - completedCycle++; - if (completedCycle % 1000 == 0) - { - Console.WriteLine($"Cycle: {completedCycle,-10} {sw.Elapsed.TotalMilliseconds,10} ms"); - } - } - }); - - sw.Stop(); - Console.WriteLine($"Cycle: {cycle,-10} {sw.Elapsed.TotalMilliseconds,10} ms"); - } - - - public static void PerformanceTestSync(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++) + for (int i = 0; i < cycle; i++) { - double Result = vsop.CalcIV(VSOPBody.EMB,1,vTime); + DynamicalELL[] ResultAll = vsop.CalcAllPlanet(vTime); completedCycle++; - if (completedCycle % 1000 == 0) + if (completedCycle % 100 == 0 && completedCycle4s. It's hard to choose which is better. - -Finally I decide to use original data file. +With LZ4 compress. Core DLL is less than 100m now. ## API @@ -171,7 +159,15 @@ Exclusive time class in VSOP2013. ## Change Log -## 2022.06.02 +### 2023.03.26 V1.1 + +Initialize performance upgrade. + +Add result classes. + +Use MessagePack to compress original data. + +### 2022.06.02 New features. @@ -179,14 +175,19 @@ Performance Optimization. Upgrade to .NET 6 -## 2020.11.14 +### 2020.11.14 Upgrade to .NET 5 -## 2020.07.06 +### 2020.07.06 Initial PR. ## 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 f1ee5ed..1be164f 100644 --- a/VSOP2013.sln +++ b/VSOP2013.sln @@ -12,30 +12,25 @@ 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}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution - Debug|Any CPU = Debug|Any CPU Debug|x64 = Debug|x64 - Release|Any CPU = Release|Any CPU Release|x64 = Release|x64 EndGlobalSection GlobalSection(ProjectConfigurationPlatforms) = postSolution - {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Debug|Any CPU.ActiveCfg = Debug|Any CPU - {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Debug|Any CPU.Build.0 = Debug|Any CPU {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Debug|x64.ActiveCfg = Debug|x64 {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Debug|x64.Build.0 = Debug|x64 - {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|Any CPU.ActiveCfg = Release|x64 - {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|Any CPU.Build.0 = Release|x64 {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|x64.ActiveCfg = Release|x64 {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|x64.Build.0 = Release|x64 - {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Debug|Any CPU.ActiveCfg = Debug|Any CPU - {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Debug|Any CPU.Build.0 = Debug|Any CPU {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Debug|x64.ActiveCfg = Debug|x64 {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Debug|x64.Build.0 = Debug|x64 - {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|Any CPU.ActiveCfg = Release|x64 - {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|Any CPU.Build.0 = Release|x64 {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|x64.ActiveCfg = Release|x64 {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|x64.Build.0 = Release|x64 + {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Debug|x64.ActiveCfg = Debug|x64 + {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Debug|x64.Build.0 = Debug|x64 + {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Release|x64.ActiveCfg = Release|x64 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/VSOP2013/Calculator.cs b/VSOP2013/Calculator.cs index 91b6b64..fa3940f 100644 --- a/VSOP2013/Calculator.cs +++ b/VSOP2013/Calculator.cs @@ -1,15 +1,15 @@ -using System.Diagnostics; +using MessagePack; +using System.Diagnostics; +using System.Reflection; +using System.Runtime.Serialization.Formatters.Binary; +using VSOP2013.VSOPResult; namespace VSOP2013 { public class Calculator { - public readonly List VSOP2013DATA; - - public TimeSpan TimeUsed; - - Stopwatch sw; + public List VSOP2013DATA; /// /// //Planetary frequency in longitude @@ -33,24 +33,26 @@ public class Calculator public Calculator() { //Import Planet Data - sw = new Stopwatch(); - TimeUsed = new TimeSpan(0); - - sw.Start(); - VSOP2013DATA = DataReader.ReadData(); - sw.Stop(); + var lz4Options = MessagePackSerializerOptions.Standard.WithCompression(MessagePackCompression.Lz4Block); + var assembly = Assembly.GetExecutingAssembly(); + var names = assembly.GetManifestResourceNames(); + VSOP2013DATA = new List(); + 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); + } + }); GC.Collect(); GC.WaitForPendingFinalizers(); - double ticks = sw.ElapsedTicks; - double Freq = Stopwatch.Frequency; - double milliseconds = (ticks / Freq) * 1000; - Console.WriteLine($"Data Load OK...Elapsed milliseconds: {milliseconds} ms"); - } - public VSOPResult[] CalcAllPlanet(VSOPTime time) + public DynamicalELL[] CalcAllPlanet(VSOPTime time) { - VSOPResult[] vsopresult = new VSOPResult[9]; + DynamicalELL[] vsopresult = new DynamicalELL[9]; ParallelLoopResult result = Parallel.For(0, 9, ip => { @@ -59,7 +61,7 @@ public VSOPResult[] CalcAllPlanet(VSOPTime time) return vsopresult; } - public VSOPResult CalcPlanet(VSOPBody body, VSOPTime time) + public DynamicalELL CalcPlanet(VSOPBody body, VSOPTime time) { double[] ELL = new double[6]; @@ -67,7 +69,7 @@ public VSOPResult CalcPlanet(VSOPBody body, VSOPTime time) { ELL[iv] = CalcIV(body, iv, time); }); - VSOPResult Coordinate = new VSOPResult(body, time, ELL); + DynamicalELL Coordinate = new DynamicalELL(body, time, ELL); return Coordinate; } @@ -91,8 +93,9 @@ public double CalcIV(VSOPBody body, int iv, VSOPTime time) /// Elliptic Elements private double CalcIV(VariableTable vTable, double JD2000) { - + //Thousand of Julian Years double tj = JD2000 / a1000; + //Iteration on Time double[] t = new double[21]; t[0] = 1.0d; diff --git a/VSOP2013/DataStruct.cs b/VSOP2013/DataStruct.cs index 2f941e7..d348f22 100644 --- a/VSOP2013/DataStruct.cs +++ b/VSOP2013/DataStruct.cs @@ -1,4 +1,5 @@ -using System; +using MessagePack; +using System; using System.Collections.Generic; using System.Linq; using System.Runtime.InteropServices; @@ -51,69 +52,79 @@ public struct PlanetEphemeris public double[] ICRSXYZ; } - + + [MessagePackObject] [Serializable] public struct PlanetTable { - public VariableTable[] variables; + [Key(0)] public VSOPBody body; + [Key(1)] + public VariableTable[] variables; + } + [MessagePackObject] [Serializable] public struct VariableTable { - - public PowerTable[] PowerTables; + [Key(0)] public VSOPBody Body; + [Key(1)] public int iv; + [Key(2)] + public PowerTable[] PowerTables; } + [MessagePackObject] [Serializable] public struct PowerTable { - public Header header; - public Term[] Terms; + [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; } + [MessagePackObject] [Serializable] public struct Header { - - /// - /// body - /// + [Key(0)] public int ip; - /// - /// variable - /// + [Key(1)] public int iv; - /// - /// power - /// + [Key(2)] public int it; - /// - /// sum of terms - /// + [Key(3)] public int nt; } + [MessagePackObject] [Serializable] [StructLayout(LayoutKind.Explicit)] public struct Term { + [Key(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; } diff --git a/VSOP2013/Resources/VSOP2013_EMB.BIN b/VSOP2013/Resources/VSOP2013_EMB.BIN new file mode 100644 index 0000000..c2fd13c Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_EMB.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_JUPITER.BIN b/VSOP2013/Resources/VSOP2013_JUPITER.BIN new file mode 100644 index 0000000..e8354dd Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_JUPITER.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_MARS.BIN b/VSOP2013/Resources/VSOP2013_MARS.BIN new file mode 100644 index 0000000..59deb43 Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_MARS.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_MERCURY.BIN b/VSOP2013/Resources/VSOP2013_MERCURY.BIN new file mode 100644 index 0000000..6df3bf7 Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_MERCURY.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN b/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN new file mode 100644 index 0000000..946a941 Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_NEPTUNE.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_PLUTO.BIN b/VSOP2013/Resources/VSOP2013_PLUTO.BIN new file mode 100644 index 0000000..3beb4b2 Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_PLUTO.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_SATURN.BIN b/VSOP2013/Resources/VSOP2013_SATURN.BIN new file mode 100644 index 0000000..f122054 Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_SATURN.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_URANUS.BIN b/VSOP2013/Resources/VSOP2013_URANUS.BIN new file mode 100644 index 0000000..43c89e0 Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_URANUS.BIN differ diff --git a/VSOP2013/Resources/VSOP2013_VENUS.BIN b/VSOP2013/Resources/VSOP2013_VENUS.BIN new file mode 100644 index 0000000..f2f5b2b Binary files /dev/null and b/VSOP2013/Resources/VSOP2013_VENUS.BIN differ diff --git a/VSOP2013/VSOP2013.csproj b/VSOP2013/VSOP2013.csproj index 4612642..7208847 100644 --- a/VSOP2013/VSOP2013.csproj +++ b/VSOP2013/VSOP2013.csproj @@ -4,44 +4,35 @@ net6.0 x64 enable - AnyCPU;x64 - + x64 + 1.1.0 + VSOP2013.NET + 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. + https://github.com/kingsznhone/VSOP2013.NET + + - - - - - - - - - + + + + + + + + + - + + + all + runtime; build; native; contentfiles; analyzers; buildtransitive + diff --git a/VSOP2013/VSOPResult.cs b/VSOP2013/VSOPResult.cs deleted file mode 100644 index 8119fd4..0000000 --- a/VSOP2013/VSOPResult.cs +++ /dev/null @@ -1,171 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Numerics; -using System.Text; -using System.Threading.Tasks; - -namespace VSOP2013 -{ - public class VSOPResult - { - //Masses system - 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}; - - double gmsol = 2.9591220836841438269e-04d; - - public VSOPBody Body { get; } - - public VSOPTime Time { get; } - - //Elliptic Elements - //a (au), lambda (radian), k, h, q, p - //Dynamical Frame J2000' - public double[] DynamicalELL; - - //Ecliptic Heliocentric Coordinates - //X,Y,Z (au) X'',Y'',Z'' (au/d) - //Dynamical Frame J2000' - public double[] DynamicalXYZ; - - //Equatorial Heliocentric Coordinates: - //X,Y,Z (au) X'',Y'',Z'' (au/d) - //ICRS Frame J2000 - public double[] ICRSXYZ; - - - public VSOPResult(VSOPBody body, VSOPTime time, double[] variables) - { - Body = body; - Time = time; - DynamicalELL = variables; - DynamicalXYZ = ELLtoXYZ((int)body, DynamicalELL); - ICRSXYZ = DynamicaltoICRS(DynamicalXYZ); - } - - /// - /// This is kind of magic that I can't undersdand - /// translate from FORTRAN code - /// - /// planet - /// Elliptic Elements - /// Ecliptic Heliocentric Coordinates - private double[] ELLtoXYZ(int ip, double[] ELL) - { - - double[] w; - Complex z, z1, z2, z3, zto, zteta; - double rgm, xa, xl, xk, xh, xq, xp, xfi, xki; - double u, ex, ex2, ex3, gl, gm, e, dl, rsa; - double xcw, xsw, xm, xr, xms, xmc, xn; - - //Initialization - rgm = Math.Sqrt(gmp[ip] + gmsol); - w = new double[6]; - xa = ELL[0]; - xl = ELL[1]; - xk = ELL[2]; - xh = ELL[3]; - xq = ELL[4]; - xp = ELL[5]; - //Computation - xfi = Math.Sqrt(1.0d - (xk * xk) - (xh * xh)); - xki = Math.Sqrt(1.0d - (xq * xq) - (xp * xp)); - u = 1.0d / (1.0d + xfi); - z = new Complex(xk, xh); - ex = z.Magnitude; - ex2 = ex * ex; - ex3 = ex2 * ex; - z1 = Complex.Conjugate(z); - gl = xl % (2 * Math.PI); - gm = gl - Math.Atan2(xh, xk); - 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); - - while (true) - { - z2 = new Complex(0d, e); - zteta = Complex.Exp(z2); - 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; - xcw = zto.Real; - xsw = zto.Imaginary; - xm = xp * xcw - xq * xsw; - xr = xa * rsa; - - w[0] = xr * (xcw - 2.0d * xp * xm); - w[1] = xr * (xsw + 2.0d * xq * xm); - w[2] = -2.0d * xr * xki * xm; - xms = xa * (xh + xsw) / xfi; - xmc = xa * (xk + xcw) / xfi; - - xn = rgm / Math.Pow(xa, 1.5d); - w[3] = xn * ((2.0d * xp * xp - 1.0d) * xms + 2.0d * xp * xq * xmc); - w[4] = xn * ((1.0d - 2.0d * xq * xq) * xmc - 2.0d * xp * xq * xms); - w[5] = 2.0d * xn * xki * (xp * xms + xq * xmc); - - return w; - } - - /// - /// Another magic function - /// - /// Ecliptic Heliocentric Coordinates - Dynamical Frame J2000 - /// Equatorial Heliocentric Coordinates - ICRS Frame J2000 - public double[] DynamicaltoICRS(double[] w) - { - double[] w2 = 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 ceps = Math.Cos(eps); - double seps = Math.Sin(eps); - double cphi = Math.Cos(phi); - double sphi = Math.Sin(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++) - { - w2[i] = 0.0d; - w2[i + 3] = 0.0d; - for (int j = 0; j < 3; j++) - { - w2[i] = w2[i] + rot[i, j] * w[j]; - w2[i + 3] = w2[i + 3] + rot[i, j] * w[j + 3]; - } - } - return w2; - } - - } -} diff --git a/VSOP2013/VSOPResult/DynamicalELL.cs b/VSOP2013/VSOPResult/DynamicalELL.cs new file mode 100644 index 0000000..1cf5929 --- /dev/null +++ b/VSOP2013/VSOPResult/DynamicalELL.cs @@ -0,0 +1,55 @@ +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 new file mode 100644 index 0000000..9acedf9 --- /dev/null +++ b/VSOP2013/VSOPResult/DynamicalXYZ.cs @@ -0,0 +1,60 @@ +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 new file mode 100644 index 0000000..e37656a --- /dev/null +++ b/VSOP2013/VSOPResult/ICRSXYZ.cs @@ -0,0 +1,55 @@ +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 new file mode 100644 index 0000000..428f7f9 --- /dev/null +++ b/VSOP2013/VSOPResult/ResultBase.cs @@ -0,0 +1,152 @@ +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; + } + + } +}