diff --git a/DataConverter/DataReader.cs b/DataConverter/DataReader.cs index 478f2da..d37b3d1 100644 --- a/DataConverter/DataReader.cs +++ b/DataConverter/DataReader.cs @@ -69,7 +69,7 @@ public static List ReadData() } VSOP2013DATA.Add(planet); } - + ParallelLoopResult result = Parallel.For(0, 9, ip => { ReadPlanet(VSOP2013DATA[ip], ip); @@ -108,7 +108,7 @@ private static void ReadPlanet(PlanetTable Planet, int ip) Planet.body = (VSOPBody)H.ip; Planet.variables[H.iv].Body = (VSOPBody)H.ip; - Planet.variables[H.iv].iv = H.iv; + Planet.variables[H.iv].Variable = (VSOPVariable) H.iv; Term[] buffer = new Term[H.nt]; for (int i = 0; i < H.nt; i++) @@ -116,9 +116,8 @@ private static void ReadPlanet(PlanetTable Planet, int ip) line = sr.ReadLine(); ReadTerm(line, ref buffer[i]); } - - 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].Variable = (VSOPVariable)H.iv; + Planet.variables[H.iv].PowerTables[H.it].Power = H.it; Planet.variables[H.iv].PowerTables[H.it].Body = (VSOPBody)H.ip; Planet.variables[H.iv].PowerTables[H.it].Terms = buffer; } @@ -128,22 +127,24 @@ private static void ReadPlanet(PlanetTable Planet, int ip) private static Header ReadHeader(string line) { + ReadOnlySpan lineSpan = line.AsSpan(); Header H = new(); int lineptr = 9; - H.ip = Convert.ToInt32(line.Substring(lineptr, 3).Trim()) - 1; + H.ip = int.Parse(lineSpan[lineptr..(lineptr + 3)].Trim()) - 1; lineptr += 3; - H.iv = Convert.ToInt32(line.Substring(lineptr, 3).Trim()) - 1; + H.iv = int.Parse(lineSpan[lineptr..(lineptr + 3)].Trim()) - 1; lineptr += 3; - H.it = Convert.ToInt32(line.Substring(lineptr, 3).Trim()); + H.it = int.Parse(lineSpan[lineptr..(lineptr + 3)].Trim()); lineptr += 3; - H.nt = Convert.ToInt32(line.Substring(lineptr, 7).Trim()); + H.nt = int.Parse(lineSpan[lineptr..(lineptr + 7)].Trim()); return H; } private static Term ReadTerm(string line, ref Term T) { + ReadOnlySpan lineSpan = line.AsSpan(); int lineptr; - int[] Bufferiphi = new int[17]; + Span Bufferiphi = stackalloc int[17]; int index = 0; double ci; @@ -155,7 +156,8 @@ private static Term ReadTerm(string line, ref Term T) // for (int counter = 0; counter < 4; counter++) { - Bufferiphi[index] = Convert.ToInt32(line.Substring(lineptr, 3).Trim()); + var debug = lineSpan[lineptr..(lineptr + 3)].Trim(); + Bufferiphi[index] = int.Parse(lineSpan[lineptr..(lineptr + 3)].Trim()); index++; lineptr += 3; } @@ -164,7 +166,8 @@ private static Term ReadTerm(string line, ref Term T) // for (int counter = 0; counter < 5; counter++) { - Bufferiphi[index] = Convert.ToInt32(line.Substring(lineptr, 3).Trim()); + + Bufferiphi[index] = int.Parse(lineSpan[lineptr..(lineptr + 3)].Trim()); index++; lineptr += 3; } @@ -173,14 +176,14 @@ private static Term ReadTerm(string line, ref Term T) // for (int counter = 0; counter < 4; counter++) { - Bufferiphi[index] = Convert.ToInt32(line.Substring(lineptr, 4).Trim()); + Bufferiphi[index] = int.Parse(lineSpan[lineptr..(lineptr + 4)].Trim()); index++; lineptr += 4; } // lineptr++; // - Bufferiphi[index] = Convert.ToInt32(line.Substring(lineptr, 6).Trim()); + Bufferiphi[index] = int.Parse(lineSpan[lineptr..(lineptr + 6)].Trim()); index++; lineptr += 6; // @@ -188,24 +191,22 @@ private static Term ReadTerm(string line, ref Term T) // for (int counter = 0; counter < 3; counter++) { - Bufferiphi[index] = Convert.ToInt32(line.Substring(lineptr, 3).Trim()); + Bufferiphi[index] = int.Parse(lineSpan[lineptr..(lineptr + 3)].Trim()); index++; lineptr += 3; } - //T.iphi = Bufferiphi; - - ci = Convert.ToDouble(line.Substring(lineptr, 20).Trim()); + ci = double.Parse(lineSpan[lineptr..(lineptr + 20)].Trim()); lineptr += 20; lineptr++; - ci *= Math.Pow(10, Convert.ToDouble(line.Substring(lineptr, 3).Trim())); + ci *= Math.Pow(10, double.Parse(lineSpan[lineptr..(lineptr + 3)].Trim())); lineptr += 3; T.ss = ci; - ci = Convert.ToDouble(line.Substring(lineptr, 20).Trim()); + ci = double.Parse(lineSpan[lineptr..(lineptr + 20)].Trim()); lineptr += 20; lineptr++; - ci *= Math.Pow(10, Convert.ToDouble(line.Substring(lineptr, 3).Trim())); + ci *= Math.Pow(10, double.Parse(lineSpan[lineptr..(lineptr + 3)].Trim())); T.cc = ci; diff --git a/Demo/Demo.csproj b/Demo/Demo.csproj index ac355ff..435ad9d 100644 --- a/Demo/Demo.csproj +++ b/Demo/Demo.csproj @@ -2,7 +2,7 @@ Exe - net6.0;net7.0 + net6.0;net7.0;net8.0 x64 diff --git a/Demo/PerfTest.cs b/Demo/PerfTest.cs index f628137..ab577c9 100644 --- a/Demo/PerfTest.cs +++ b/Demo/PerfTest.cs @@ -1,4 +1,5 @@ using System; +using System.Threading.Tasks; using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Jobs; using VSOP2013; @@ -6,8 +7,8 @@ namespace Demo { [SimpleJob(RuntimeMoniker.Net60)] + [SimpleJob(RuntimeMoniker.Net70)] [SimpleJob(RuntimeMoniker.Net80)] - [MemoryDiagnoser] public class PerfTest { @@ -28,11 +29,18 @@ public void init() vTime = new VSOPTime(dt, TimeFrame.UTC); } - [Benchmark] + [Benchmark(Baseline = true)] public VSOPResult Compute() { var ell = vsop.GetPlanetPosition(VSOPBody.JUPITER, vTime); return ell; } + + [Benchmark] + public VSOPResult Compute_Native() + { + var ell = vsop.GetPlanetPosition_Native(VSOPBody.JUPITER, vTime); + return ell; + } } } \ No newline at end of file diff --git a/Demo/Program.cs b/Demo/Program.cs index 9a5d27c..34f2e53 100644 --- a/Demo/Program.cs +++ b/Demo/Program.cs @@ -37,24 +37,24 @@ private static void Main(string[] args) ell = vsop.GetPlanetPosition(VSOPBody.EMB, vTime); FormattedPrint(ell, vTime); + ell = vsop.GetPlanetPosition_Native(VSOPBody.EMB, vTime); + FormattedPrint(ell, vTime); + //xyz = (VSOPResult_XYZ)ell; + //FormattedPrint(xyz, vTime); + //xyz.ReferenceFrame = ReferenceFrame.ICRSJ2000; + //FormattedPrint(xyz, vTime); - xyz = (VSOPResult_XYZ)ell; - FormattedPrint(xyz, vTime); - xyz.ReferenceFrame = ReferenceFrame.ICRSJ2000; - FormattedPrint(xyz, vTime); - - lbr = (VSOPResult_LBR)ell; - FormattedPrint(lbr, vTime); - lbr.ReferenceFrame = ReferenceFrame.ICRSJ2000; - FormattedPrint(lbr, vTime); + //lbr = (VSOPResult_LBR)ell; + //FormattedPrint(lbr, vTime); + //lbr.ReferenceFrame = ReferenceFrame.ICRSJ2000; + //FormattedPrint(lbr, vTime); - //for(int i = 0; i < 1000; i++) - { - foreach (VSOPBody body in Enum.GetValues(typeof(VSOPBody))) - { - ell = vsop.GetPlanetPosition(body, vTime); - } - } + //{ + // foreach (VSOPBody body in Enum.GetValues(typeof(VSOPBody))) + // { + // ell = vsop.GetPlanetPosition(body, vTime); + // } + //} Console.WriteLine("Press Enter to Start Performance Test..."); Console.ReadLine(); @@ -63,6 +63,7 @@ private static void Main(string[] args) #else var summary = BenchmarkRunner.Run(); #endif + Console.ReadLine(); } public static void FormattedPrint(VSOPResult Result, VSOPTime vtime) @@ -106,8 +107,8 @@ public static void FormattedPrint(VSOPResult Result, VSOPTime vtime) WriteColorLine("Reference Frame: ", ConsoleColor.Green, $"\tICRS equinox and ecliptic J2000"); break; } - 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")}"); + WriteColorLine("At UTC: ", ConsoleColor.Green, $"\t\t{Result.Time.UTC.ToUniversalTime():o}"); + WriteColorLine("At TDB: ", ConsoleColor.Green, $"\t\t{Result.Time.TDB:o}"); if (Result.CoordinatesType == CoordinatesType.Elliptic) { @@ -158,8 +159,8 @@ private static void WriteColorLine(params object[] oo) foreach (var o in oo) if (o == null) Console.ResetColor(); - else if (o is ConsoleColor) - Console.ForegroundColor = (ConsoleColor)o; + else if (o is ConsoleColor color) + Console.ForegroundColor = color; else Console.Write(o.ToString()); Console.WriteLine(); diff --git a/Demo/Properties/launchSettings.json b/Demo/Properties/launchSettings.json new file mode 100644 index 0000000..e84c5b7 --- /dev/null +++ b/Demo/Properties/launchSettings.json @@ -0,0 +1,8 @@ +{ + "profiles": { + "Demo": { + "commandName": "Project", + "nativeDebugging": true + } + } +} \ No newline at end of file diff --git a/NativeAccelerator/NativeAccelerator.c b/NativeAccelerator/NativeAccelerator.c new file mode 100644 index 0000000..b1e0187 --- /dev/null +++ b/NativeAccelerator/NativeAccelerator.c @@ -0,0 +1,13 @@ +#include "NativeAccelerator.h" +#include "math.h" + +double StartIteration(struct Term* terms, int length, double tj, double tit) { + double result = 0; + for (int n = 0; n < length; n++) { + double u = terms[n].aa + terms[n].bb * tj; + double su = sin(u); + double cu = cos(u); + result += tit * (terms[n].ss * su + terms[n].cc * cu); + } + return result; +} diff --git a/NativeAccelerator/NativeAccelerator.h b/NativeAccelerator/NativeAccelerator.h new file mode 100644 index 0000000..b42b8ba --- /dev/null +++ b/NativeAccelerator/NativeAccelerator.h @@ -0,0 +1,26 @@ +#pragma once + + +#ifdef DLLEXPORT// DLLEXPORT +#define DLL_EXPORT __declspec(dllexport) +#else +#define DLL_EXPORT __declspec(dllimport) +#endif + +struct Term +{ + double ss; + double cc; + double aa; + double bb; +}; + +#ifdef __cplusplus +extern "C" { +#endif + +DLL_EXPORT double StartIteration(struct Term* terms,int length, double tj, double tit); + +#ifdef __cplusplus +} +#endif \ No newline at end of file diff --git a/NativeAccelerator/NativeAccelerator.vcxproj b/NativeAccelerator/NativeAccelerator.vcxproj new file mode 100644 index 0000000..2b9ce1d --- /dev/null +++ b/NativeAccelerator/NativeAccelerator.vcxproj @@ -0,0 +1,165 @@ + + + + + Debug + Win32 + + + Release + Win32 + + + Debug + x64 + + + Release + x64 + + + + 17.0 + Win32Proj + {1f19706f-2068-4359-a307-ff6620d57100} + NativeMethod + 10.0.19041.0 + NativeAccelerator + + + + Application + true + v143 + Unicode + false + + + Application + false + v143 + true + Unicode + false + + + DynamicLibrary + true + v143 + Unicode + false + + + DynamicLibrary + false + v143 + true + Unicode + false + + + + + + + + + + + + + + + + + + + + + + Level3 + true + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + + + Console + true + + + + + Level3 + true + true + true + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + + + Console + true + true + true + + + + + Level3 + true + _DEBUG;_CONSOLE;DLLEXPORT;%(PreprocessorDefinitions) + true + stdcpp20 + stdc17 + AdvancedVectorExtensions2 + Fast + + + Console + true + C:\Program Files (x86)\Intel\oneAPI\ipp\latest\lib;C:\Program Files (x86)\Intel\oneAPI\mkl\latest\lib;C:\Program Files (x86)\Intel\oneAPI\tbb\latest\lib;%(AdditionalLibraryDirectories) + + + copy /Y "$(TargetPath)" "$(SolutionDir)VSOP2013.NET\Resources" + + + + + Level3 + true + true + true + NDEBUG;_CONSOLE;DLLEXPORT;%(PreprocessorDefinitions) + + + stdcpp20 + stdc17 + AdvancedVectorExtensions2 + Fast + Speed + + All + true + true + + + Console + true + true + true + C:\Program Files (x86)\Intel\oneAPI\ipp\latest\lib;C:\Program Files (x86)\Intel\oneAPI\mkl\latest\lib;C:\Program Files (x86)\Intel\oneAPI\tbb\latest\lib;%(AdditionalLibraryDirectories) + + + copy /Y "$(TargetPath)" "$(SolutionDir)VSOP2013.NET\Resources" + + + + + + + + + + + + \ No newline at end of file diff --git a/NativeAccelerator/NativeAccelerator.vcxproj.filters b/NativeAccelerator/NativeAccelerator.vcxproj.filters new file mode 100644 index 0000000..8f68bbe --- /dev/null +++ b/NativeAccelerator/NativeAccelerator.vcxproj.filters @@ -0,0 +1,21 @@ + + + + + {d9c5ed58-0657-4da5-b32b-ed112b925ae8} + + + {51430353-5ca6-44cb-92ca-53dbc6673792} + + + + + Header Files + + + + + Source Files + + + \ No newline at end of file diff --git a/README.md b/README.md index 61dd6a8..19edf14 100644 --- a/README.md +++ b/README.md @@ -23,21 +23,27 @@ I promise it will be much faster than origin algorithm. This is the best VSOP2013 library ever. -![Demo](https://github.com/kingsznhone/VSOP2013.NET/blob/main/Demo.png) +![Demo](./README/Demo.png) + +## Performance + +
![Performance Test](./README/NativeAccelerate.png) + +Note: .NET 8 occurs performance regression due to RyuJIT bugs. [Detail Here](https://github.com/dotnet/runtime/issues/95954#issuecomment-1956661569) ## Features 1. Use VSOPResult class to manage calculate results. 2. Use VSOPTime class to manage time.
Easy 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. +3. Very high performance per solution +4. Useful Utility class. Convert Elliptic coordinates to cartesian and spherical +5. Async Api included. 6. precalculation on φ in terms, which gives 20%+ speed up 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. +
Initialization time becomes less than 10% of previous version. 8. Brotli compression on source data. ~300Mb -> ~50MB. +9. Optional Native Side library accelerate. (Only on Windows)
@@ -58,7 +64,7 @@ The planetary solution VSOP2013 is fitted to the numerical integration INPOP10a * NuGet Package Manager ``` - PM> NuGet\Install-Package VSOP2013.NET -Version 1.1.8 + PM> NuGet\Install-Package VSOP2013.NET ``` ``` @@ -100,10 +106,23 @@ Console.WriteLine("============================================================= ## Change Log + +### 2024.04.24 v1.2.0 + +Add .NET 8. which cause performance regression. + +Add Native Accelerate method ```GetPlanetPosition_Native``` to accelerate calculation. with 30%+ performance Improvment. (Experimental) + +Native CPP code is only for Windows x64 AVX2 enviroment. + +Using fast floating-point compile options on native libraries can result in a decrease in precision and is difficult to estimate. + + ### 2024.01.14 v1.1.8 Fix critical error in ELL to XYZ convertion. + ### 2023.12.13 v1.1.7 Bug fix. @@ -210,6 +229,38 @@ Exclusive time class in VSOP2013.
+### ```public double GetVariable_Native(VSOPBody body,int iv, VSOPTime time)``` + +Calculate a specific variable of a planet. Using Native Accelerate. + +
+ +#### 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 Task GetVariableAsync(VSOPBody body,int iv, VSOPTime time)``` Calculate a specific variable of a planet. @@ -244,6 +295,40 @@ variable result.
+### ```public Task GetVariableAsync_Native(VSOPBody body,int iv, VSOPTime time)``` + +Calculate a specific variable of a planet. Using Native Accelerate. + +
+ +#### 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. @@ -274,6 +359,36 @@ Can be explicit cast to ```VSOPResult_XYZ``` and ```VSOPResult_LBR```
+### ```public VSOPResult_ELL GetPlanetPosition_Native(VSOPBody body, VSOPTime time)``` + +Calculate all variable of a planet. Using Native Accelerate. + +
+ +#### Parameters + +```body``` VSOPBody + +Planet enum. + +
+ +```time``` VSOPTime + +Exclusive time class in VSOP2013. + +
+ +#### Return + +```VSOPResult_ELL``` + +Full Result with 6 variable of Elliptic Coordinates. + +Can be explicit cast to ```VSOPResult_XYZ``` and ```VSOPResult_LBR``` + +
+ ### ```public Task GetPlanetPositionAsync(VSOPBody body, VSOPTime time)``` Calculate all variable of a planet. @@ -304,6 +419,37 @@ Can be explicit cast to ```VSOPResult_XYZ``` and ```VSOPResult_LBR```
+### ```public Task GetPlanetPositionAsync_Native(VSOPBody body, VSOPTime time)``` + +Calculate all variable of a planet. Using Native Accelerate. + +
+ +#### Parameters + +```body``` VSOPBody + +Planet enum. + +
+ +```time``` VSOPTime + +Exclusive time class in VSOP2013. + +
+ +#### Return + +```Task``` + +Full Result with 6 variable of Elliptic Coordinates. + +Can be explicit cast to ```VSOPResult_XYZ``` and ```VSOPResult_LBR``` + +
+ + ## Static Class Utility This Class Provide some useful function. @@ -376,7 +522,7 @@ Array of cartesian coordinate elements ### ```static double[] ELLtoXYZ(double[] ell)``` -This is a magic function I directly copy from VSOP2013. +This is a magic function I directly copy from original VSOP2013. It's way beyond my math level. diff --git a/Demo.png b/README/Demo.png similarity index 100% rename from Demo.png rename to README/Demo.png diff --git a/README/NativeAccelerate.png b/README/NativeAccelerate.png new file mode 100644 index 0000000..bfc65f2 Binary files /dev/null and b/README/NativeAccelerate.png differ diff --git a/PerformanceTest.png b/README/PerformanceTest.png similarity index 100% rename from PerformanceTest.png rename to README/PerformanceTest.png diff --git a/VSOP2013.NET.sln b/VSOP2013.NET.sln index 406531a..c3347f2 100644 --- a/VSOP2013.NET.sln +++ b/VSOP2013.NET.sln @@ -14,23 +14,47 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution EndProject Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "DataConverter", "DataConverter\DataConverter.csproj", "{2A5C5EFE-C49F-436A-9FC1-176AE11E0334}" EndProject +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "NativeAccelerator", "NativeAccelerator\NativeAccelerator.vcxproj", "{1F19706F-2068-4359-A307-FF6620D57100}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|x64 = Debug|x64 + Debug|x86 = Debug|x86 Release|x64 = Release|x64 + Release|x86 = Release|x86 EndGlobalSection GlobalSection(ProjectConfigurationPlatforms) = postSolution {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}.Debug|x86.ActiveCfg = Debug|x64 + {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Debug|x86.Build.0 = Debug|x64 {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|x64.ActiveCfg = Release|x64 {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|x64.Build.0 = Release|x64 + {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|x86.ActiveCfg = Release|x64 + {A36BC2BE-D01A-497C-B17A-3B15DA455D9F}.Release|x86.Build.0 = Release|x64 {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}.Debug|x86.ActiveCfg = Debug|x64 + {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Debug|x86.Build.0 = Debug|x64 {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|x64.ActiveCfg = Release|x64 {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|x64.Build.0 = Release|x64 + {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|x86.ActiveCfg = Release|x64 + {9CC39AF5-6179-4903-8F3C-ED8ACEFE2399}.Release|x86.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}.Debug|x86.ActiveCfg = Debug|x64 + {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Debug|x86.Build.0 = Debug|x64 {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Release|x64.ActiveCfg = Release|x64 + {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Release|x86.ActiveCfg = Release|x64 + {2A5C5EFE-C49F-436A-9FC1-176AE11E0334}.Release|x86.Build.0 = Release|x64 + {1F19706F-2068-4359-A307-FF6620D57100}.Debug|x64.ActiveCfg = Debug|x64 + {1F19706F-2068-4359-A307-FF6620D57100}.Debug|x64.Build.0 = Debug|x64 + {1F19706F-2068-4359-A307-FF6620D57100}.Debug|x86.ActiveCfg = Debug|Win32 + {1F19706F-2068-4359-A307-FF6620D57100}.Debug|x86.Build.0 = Debug|Win32 + {1F19706F-2068-4359-A307-FF6620D57100}.Release|x64.ActiveCfg = Release|x64 + {1F19706F-2068-4359-A307-FF6620D57100}.Release|x64.Build.0 = Release|x64 + {1F19706F-2068-4359-A307-FF6620D57100}.Release|x86.ActiveCfg = Release|Win32 + {1F19706F-2068-4359-A307-FF6620D57100}.Release|x86.Build.0 = Release|Win32 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/VSOP2013.NET/Calculator.cs b/VSOP2013.NET/Calculator.cs index e575c02..78dfb52 100644 --- a/VSOP2013.NET/Calculator.cs +++ b/VSOP2013.NET/Calculator.cs @@ -1,12 +1,17 @@ using MessagePack; using System.IO.Compression; using System.Reflection; +using System.Runtime.InteropServices; namespace VSOP2013 { public class Calculator { - public List VSOP2013DATA; + + [DllImport("Resources/NativeAccelerator.dll")] + private static extern double StartIteration(Term[] terms,int length, double tj, double tit); + + private List VSOP2013DATA; /// /// //Planetary frequency in longitude @@ -40,6 +45,7 @@ public Calculator() var data = MessagePackSerializer.Deserialize(bs); VSOP2013DATA.Add(data); }); + VSOP2013DATA = VSOP2013DATA.OrderBy(x => x.body).ToList(); GC.Collect(); @@ -56,11 +62,26 @@ public VSOPResult_ELL GetPlanetPosition(VSOPBody body, VSOPTime time) VSOPResult_ELL Coordinate = new(body, time, ELL); return Coordinate; } + public VSOPResult_ELL GetPlanetPosition_Native(VSOPBody body, VSOPTime time) + { + double[] ELL = new double[6]; + ParallelLoopResult result = Parallel.For(0, 6, iv => + { + ELL[iv] = GetVariable_Native(body, iv, time); + }); + VSOPResult_ELL Coordinate = new(body, time, ELL); + return Coordinate; + } public async Task GetPlanetPositionAsync(VSOPBody body, VSOPTime time) { return await Task.Run(() => GetPlanetPosition(body, time)); } + public async Task GetPlanetPositionAsync_Native(VSOPBody body, VSOPTime time) + { + return await Task.Run(() => GetPlanetPosition_Native(body, time)); + } + /// /// Calculate a specific variable @@ -73,11 +94,19 @@ public double GetVariable(VSOPBody body, int iv, VSOPTime time) { return Calculate(VSOP2013DATA[(int)body].variables[iv], time.J2000); } + public double GetVariable_Native(VSOPBody body, int iv, VSOPTime time) + { + return Calculate_Native(VSOP2013DATA[(int)body].variables[iv], time.J2000); + } public async Task GetVariableAsync(VSOPBody body, int iv, VSOPTime time) { return await Task.Run(() => GetVariable(body, iv, time)); } + public async Task GetVariableAsync_Native(VSOPBody body, int iv, VSOPTime time) + { + return await Task.Run(() => GetVariable_Native(body, iv, time)); + } /// /// @@ -109,18 +138,54 @@ private double Calculate(VariableTable Table, double JD2000) for (int n = 0; n < terms.Length; n++) { u = terms[n].aa + terms[n].bb * tj; +#if NET8_0 + su = Math.Sin(u); + cu = Math.Cos(u); +#else (su, cu) = Math.SinCos(u); +#endif result += t[it] * (terms[n].ss * su + terms[n].cc * cu); } } - if (Table.iv == 1) + if (Table.Variable == VSOPVariable.A) + { + xl = result + freqpla[(int)Table.Body] * tj; + xl = (xl % Math.Tau + Math.Tau) % Math.Tau; + result = xl; + } + return result; + } + + private unsafe double Calculate_Native(VariableTable Table, double JD2000) + { + //Thousand of Julian Years + double tj = JD2000 / a1000; + + //Iteration on Time + Span t = stackalloc double[21]; + t[0] = 1.0d; + t[1] = tj; + for (int i = 2; i < 21; i++) + { + t[i] = t[1] * t[i - 1]; + } + + double result = 0d; + double xl; + for (int it = 0; it < Table.PowerTables.Length; it++) + { + if (Table.PowerTables[it].Terms == null) continue; + Term[] terms = Table.PowerTables[it].Terms; + result += StartIteration(Table.PowerTables[it].Terms, terms.Length, tj, t[it]); + } + if (Table.Variable == (VSOPVariable)1) { xl = result + freqpla[(int)Table.Body] * tj; - //modulu into [0,tau) xl = (xl % Math.Tau + Math.Tau) % Math.Tau; result = xl; } return result; } + } } \ No newline at end of file diff --git a/VSOP2013.NET/DataStruct.cs b/VSOP2013.NET/DataStruct.cs index 6c626e5..249c54b 100644 --- a/VSOP2013.NET/DataStruct.cs +++ b/VSOP2013.NET/DataStruct.cs @@ -16,35 +16,14 @@ public enum VSOPBody PLUTO = 8 } - public enum VSOPFileName + public enum VSOPVariable { - VSOP2013p1 = 0, - VSOP2013p2 = 1, - VSOP2013p3 = 2, - VSOP2013p4 = 3, - VSOP2013p5 = 4, - VSOP2013p6 = 5, - VSOP2013p7 = 6, - VSOP2013p8 = 7, - VSOP2013p9 = 8 - } - - public struct PlanetEphemeris - { - //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; + A = 0, + L = 1, + K = 2, + H = 3, + Q = 4, + P = 5, } [MessagePackObject] @@ -66,7 +45,7 @@ public struct VariableTable public VSOPBody Body; [Key(1)] - public int iv; + public VSOPVariable Variable; [Key(2)] public PowerTable[] PowerTables; @@ -80,13 +59,16 @@ public struct PowerTable public VSOPBody Body; [Key(1)] - public int iv; + public VSOPVariable Variable; [Key(2)] - public int it; + public int Power; + + [Key(3)] + public int TermsCount; [Key(4)] - public Term[] Terms; + public Term[] Terms; } [MessagePackObject] diff --git a/VSOP2013.NET/Resources/NativeAccelerator.dll b/VSOP2013.NET/Resources/NativeAccelerator.dll new file mode 100644 index 0000000..c655071 Binary files /dev/null and b/VSOP2013.NET/Resources/NativeAccelerator.dll differ diff --git a/VSOP2013.NET/Resources/VSOP2013_EMB.BIN b/VSOP2013.NET/Resources/VSOP2013_EMB.BIN index bbad015..6457ae0 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_EMB.BIN and b/VSOP2013.NET/Resources/VSOP2013_EMB.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_JUPITER.BIN b/VSOP2013.NET/Resources/VSOP2013_JUPITER.BIN index d33ee33..e0bc339 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_JUPITER.BIN and b/VSOP2013.NET/Resources/VSOP2013_JUPITER.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_MARS.BIN b/VSOP2013.NET/Resources/VSOP2013_MARS.BIN index d4a23fc..69ef036 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_MARS.BIN and b/VSOP2013.NET/Resources/VSOP2013_MARS.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_MERCURY.BIN b/VSOP2013.NET/Resources/VSOP2013_MERCURY.BIN index dc6a5f6..e0b811b 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_MERCURY.BIN and b/VSOP2013.NET/Resources/VSOP2013_MERCURY.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_NEPTUNE.BIN b/VSOP2013.NET/Resources/VSOP2013_NEPTUNE.BIN index 1cc8cb2..89ba300 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_NEPTUNE.BIN and b/VSOP2013.NET/Resources/VSOP2013_NEPTUNE.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_PLUTO.BIN b/VSOP2013.NET/Resources/VSOP2013_PLUTO.BIN index c48e6e9..19864b8 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_PLUTO.BIN and b/VSOP2013.NET/Resources/VSOP2013_PLUTO.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_SATURN.BIN b/VSOP2013.NET/Resources/VSOP2013_SATURN.BIN index 937650d..0ea9fd7 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_SATURN.BIN and b/VSOP2013.NET/Resources/VSOP2013_SATURN.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_URANUS.BIN b/VSOP2013.NET/Resources/VSOP2013_URANUS.BIN index 163eb55..72021ce 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_URANUS.BIN and b/VSOP2013.NET/Resources/VSOP2013_URANUS.BIN differ diff --git a/VSOP2013.NET/Resources/VSOP2013_VENUS.BIN b/VSOP2013.NET/Resources/VSOP2013_VENUS.BIN index f09934b..f3d8093 100644 Binary files a/VSOP2013.NET/Resources/VSOP2013_VENUS.BIN and b/VSOP2013.NET/Resources/VSOP2013_VENUS.BIN differ diff --git a/VSOP2013.NET/Utility.cs b/VSOP2013.NET/Utility.cs index 239d360..8aa48b7 100644 --- a/VSOP2013.NET/Utility.cs +++ b/VSOP2013.NET/Utility.cs @@ -1,19 +1,22 @@ -using System.Numerics; +using System.Numerics; using System.Runtime.Intrinsics; 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[] 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; @@ -92,16 +95,16 @@ public static double[] XYZtoLBR(double[] xyz) if (Vector256.IsHardwareAccelerated) { - Vector256 v1 = Vector256.Create(x / r, y / r, z / r, 0); + Vector256 v1 = Vector256.Create(x / r , y / r , z / r , 0); Vector256 v2 = Vector256.Create(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)), 0); - Vector256 v3 = Vector256.Create(-y / (x * x + y * y), x / (x * x + y * y), 0, 0); - Vector256 vv = Vector256.Create(dx, dy, dz, 0); + Vector256 v3 = Vector256.Create(-y / (x * x + y * y) , x / (x * x + y * y) , 0 , 0); + Vector256 vv = Vector256.Create(dx , dy , dz , 0); lbr[0] = l; lbr[1] = b; lbr[2] = r; lbr[3] = Vector256.Sum(vv * v1); - lbr[4] = Vector256.Sum(vv * v2); + lbr[4] =-Vector256.Sum(vv * v2); lbr[5] = Vector256.Sum(vv * v3); return lbr.ToArray(); } @@ -131,7 +134,7 @@ public static double[] XYZtoLBR(double[] xyz) lbr[1] = b; lbr[2] = r; lbr[3] = dl; - lbr[4] = -db; + lbr[4] =-db; lbr[5] = dr; return lbr.ToArray(); } @@ -165,16 +168,16 @@ public static double[] LBRtoXYZ(double[] lbr) if (Vector256.IsHardwareAccelerated) { Vector256 m1 = Vector256.Create(Math.Cos(b) * Math.Cos(l), r * Math.Sin(b) * Math.Cos(l), -r * Math.Cos(b) * Math.Sin(l), 0); - Vector256 m2 = Vector256.Create(Math.Cos(b) * Math.Sin(l), r * Math.Sin(b) * Math.Sin(l), r * Math.Cos(b) * Math.Cos(l), 0); - Vector256 m3 = Vector256.Create(Math.Sin(b), -r * Math.Cos(b), 0, 0); - Vector256 vv = Vector256.Create(dr, db, dl, 0); + Vector256 m2 = Vector256.Create(Math.Cos(b) * Math.Sin(l), r * Math.Sin(b) * Math.Sin(l), r * Math.Cos(b) * Math.Cos(l) , 0); + Vector256 m3 = Vector256.Create(Math.Sin(b) , -r * Math.Cos(b) , 0 , 0); + Vector256 vv = Vector256.Create(dr , db , dl , 0); xyz[0] = x; xyz[1] = y; xyz[2] = z; xyz[3] = Vector256.Sum(vv * m1); xyz[4] = Vector256.Sum(vv * m2); - xyz[5] = Vector256.Sum(vv * m3); + xyz[5] =-Vector256.Sum(vv * m3); return xyz.ToArray(); } @@ -198,13 +201,13 @@ public static double[] LBRtoXYZ(double[] lbr) xyz[2] = z; xyz[3] = dx; xyz[4] = dy; - xyz[5] = -dz; + xyz[5] =-dz; return xyz.ToArray(); } /// /// Convert Elliptic coordinate to cardinal coordinate. - /// This is kind of magic that I will never undersdand. + /// This is kind of magic that I will never understand. /// Directly translate from VSOP2013.f. /// /// planet @@ -235,14 +238,13 @@ public static double[] ELLtoXYZ(VSOPBody body, double[] ell) z = new Complex(k, h); ex = z.Magnitude; ex2 = ex * ex; - ex3 = ex2 * ex; + ex3 = ex * ex * 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); - while (true) { z2 = new Complex(0d, e); @@ -263,12 +265,13 @@ public static double[] ELLtoXYZ(VSOPBody body, double[] ell) 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[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(); @@ -306,9 +309,6 @@ public static double[] DynamicaltoICRS(double[] dynamical) Cphi = Math.Cos(phi); #if NET7_0 - - #region vector matrix mul - if (Vector256.IsHardwareAccelerated) { Vector256 m1 = Vector256.Create(Cphi, -Sphi * Ceps, Sphi * Seps, 0); @@ -326,9 +326,6 @@ public static double[] DynamicaltoICRS(double[] dynamical) icrs[5] = Vector256.Sum(vdv * m3); return icrs.ToArray(); } - - #endregion vector matrix mul - #endif //Rotation Matrix diff --git a/VSOP2013.NET/VSOP2013.NET.csproj b/VSOP2013.NET/VSOP2013.NET.csproj index 6753f14..b44af1d 100644 --- a/VSOP2013.NET/VSOP2013.NET.csproj +++ b/VSOP2013.NET/VSOP2013.NET.csproj @@ -1,11 +1,11 @@  - net6.0;net7.0 + net6.0;net7.0;net8.0 x64 enable x64 - 1.1.8 + 1.2.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. @@ -15,6 +15,7 @@ VSOP2013, computed the positions of the planets directly at any moment, as well https://github.com/kingsznhone/VSOP2013.NET True VSOP2013.NET + true @@ -37,4 +38,16 @@ VSOP2013, computed the positions of the planets directly at any moment, as well + + + Always + + + Always + + + Always + + +