diff --git a/Demo/Program.cs b/Demo/Program.cs index 833127b..5144c02 100644 --- a/Demo/Program.cs +++ b/Demo/Program.cs @@ -1,5 +1,6 @@ using System; using System.Globalization; +using BenchmarkDotNet.Configs; using BenchmarkDotNet.Running; using VSOP2013; diff --git a/README.md b/README.md index f5a638b..92eda9b 100644 --- a/README.md +++ b/README.md @@ -96,6 +96,12 @@ Console.WriteLine("============================================================= ## Change Log +### 2023.9.3 v1.1.6 + +SIMD Accel Utility. + +Bug Fix + ### 2023.7.7 v1.1.5 Bug Fix. diff --git a/VSOP2013.NET/Utility.cs b/VSOP2013.NET/Utility.cs index 21d315f..f5d642e 100644 --- a/VSOP2013.NET/Utility.cs +++ b/VSOP2013.NET/Utility.cs @@ -1,4 +1,5 @@ using System.Numerics; +using System.Runtime.Intrinsics; namespace VSOP2013 { @@ -85,6 +86,24 @@ public static double[] XYZtoLBR(double[] xyz) //Sin(θ) = Cos(b), Cos(θ) = Sin(b) b = Math.Asin(z / r); + +#if NET7_0 + #region vector matrix mul + if (Vector256.IsHardwareAccelerated) + { + 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); + + lbr[3] = Vector256.Sum(vv * v1); + lbr[4] = Vector256.Sum(vv * v2); + lbr[5] = Vector256.Sum(vv * v3); + return lbr.ToArray(); + } + #endregion +#endif + // Inverse Jacobian matrix From Caridnal to Sperical //https://en.wikipedia.org/wiki/Spherical_coordinate_system#Integration_and_differentiation_in_spherical_coordinates double[,] J_1 = { @@ -133,6 +152,23 @@ public static double[] LBRtoXYZ(double[] lbr) y = r * Math.Cos(b) * Math.Sin(l); z = r * Math.Sin(b); + +#if NET7_0 + #region vector matrix mul + if (Vector256.IsHardwareAccelerated) + { + Vector256 v1 = Vector256.Create(Math.Cos(b) * Math.Cos(l), r * Math.Sin(b) * Math.Cos(l), -r * Math.Cos(b) * Math.Sin(l), 0); + Vector256 v2 = Vector256.Create(Math.Cos(b) * Math.Sin(l), r * Math.Sin(b) * Math.Sin(l), r * Math.Cos(b) * Math.Cos(l), 0); + Vector256 v3 = Vector256.Create(Math.Sin(b), -r * Math.Cos(b), 0, 0); + Vector256 vv = Vector256.Create(dr, db, dl, 0); + + xyz[3] = Vector256.Sum(vv * v1); + xyz[4] = Vector256.Sum(vv * v2); + xyz[5] = Vector256.Sum(vv * v3); + return xyz.ToArray(); + } + #endregion +#endif // Jacobian matrix From Sperical to Caridnal //https://en.wikipedia.org/wiki/Spherical_coordinate_system#Integration_and_differentiation_in_spherical_coordinates double[,] J = { @@ -257,6 +293,29 @@ public static double[] DynamicaltoICRS(double[] dynamical) Sphi = Math.Sin(phi); Cphi = Math.Cos(phi); + +#if NET7_0 + #region vector matrix mul + if (Vector256.IsHardwareAccelerated) + { + Vector256 v1 = Vector256.Create(Cphi, -Sphi * Ceps, Sphi * Seps, 0); + Vector256 v2 = Vector256.Create(Sphi, Cphi * Ceps, -Cphi * Seps, 0); + Vector256 v3 = Vector256.Create(0, Seps, Ceps, 0); + Vector256 vv = Vector256.Create(dynamical[0], dynamical[1], dynamical[2], 0); + Vector256 vv2 = Vector256.Create(dynamical[3], dynamical[4], dynamical[5], 0); + + icrs[0] = Vector256.Sum(vv * v1); + icrs[1] = Vector256.Sum(vv * v2); + icrs[2] = Vector256.Sum(vv * v3); + + icrs[3] = Vector256.Sum(vv2 * v1); + icrs[4] = Vector256.Sum(vv2 * v2); + icrs[5] = Vector256.Sum(vv2 * v3); + return icrs.ToArray(); + } + #endregion +#endif + //Rotation Matrix double[,] R = new double[,] { {Cphi, -Sphi*Ceps, Sphi*Seps }, {Sphi, Cphi*Ceps, -Cphi*Seps }, @@ -269,20 +328,20 @@ public static double[] DynamicaltoICRS(double[] dynamical) double[,] C = MultiplyMatrix(R, A); - dynamical[0] = C[0, 0]; - dynamical[1] = C[1, 0]; - dynamical[2] = C[2, 0]; + icrs[0] = C[0, 0]; + icrs[1] = C[1, 0]; + icrs[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]; + icrs[3] = C[0, 0]; + icrs[4] = C[1, 0]; + icrs[5] = C[2, 0]; - return dynamical.ToArray(); + return icrs.ToArray(); } /// @@ -305,6 +364,29 @@ public static double[] ICRStoDynamical(double[] icrs) Sphi = Math.Sin(phi); Cphi = Math.Cos(phi); + +#if NET7_0 + #region vector matrix mul + if (Vector256.IsHardwareAccelerated) + { + Vector256 v1 = Vector256.Create(Cphi, Sphi, 0, 0); + Vector256 v2 = Vector256.Create(-Sphi * Ceps, Cphi * Ceps, Seps, 0); + Vector256 v3 = Vector256.Create(Sphi * Seps, -Cphi * Seps, Ceps, 0); + Vector256 vv = Vector256.Create(icrs[0], icrs[1], icrs[2], 0); + Vector256 vv2 = Vector256.Create(icrs[3], icrs[4], icrs[5], 0); + + dynamical[0] = Vector256.Sum(vv * v1); + dynamical[1] = Vector256.Sum(vv * v2); + dynamical[2] = Vector256.Sum(vv * v3); + + dynamical[3] = Vector256.Sum(vv2 * v1); + dynamical[4] = Vector256.Sum(vv2 * v2); + dynamical[5] = Vector256.Sum(vv2 * v3); + return dynamical.ToArray(); + } + #endregion +#endif + //Reverse Matrix double[,] R_1 = new double[,] {{ Cphi, Sphi, 0 }, {-Sphi*Ceps, Cphi*Ceps, Seps }, diff --git a/VSOP2013.NET/VSOP2013.NET.csproj b/VSOP2013.NET/VSOP2013.NET.csproj index 667df63..5f40d7c 100644 --- a/VSOP2013.NET/VSOP2013.NET.csproj +++ b/VSOP2013.NET/VSOP2013.NET.csproj @@ -5,7 +5,7 @@ x64 enable x64 - 1.1.5 + 1.1.6 VSOP2013.NET KingsZNHONE VSOP was developed and is maintained (updated with the latest data) by the scientists at the Bureau des Longitudes in Paris. diff --git a/VSOP2013.NET/VSOPResult/VSOPResult_LBR.cs b/VSOP2013.NET/VSOPResult/VSOPResult_LBR.cs index d06ca68..826a119 100644 --- a/VSOP2013.NET/VSOPResult/VSOPResult_LBR.cs +++ b/VSOP2013.NET/VSOPResult/VSOPResult_LBR.cs @@ -44,6 +44,7 @@ public override ReferenceFrame ReferenceFrame /// Elliptic elements public VSOPResult_LBR(VSOPResult_ELL result) { + Variables_ELL = result.Variables_ELL; _coordinatesReference = result.CoordinatesReference; _referenceFrame = result.ReferenceFrame; Body = result.Body; @@ -53,6 +54,7 @@ public VSOPResult_LBR(VSOPResult_ELL result) public VSOPResult_LBR(VSOPResult_XYZ result) { + Variables_ELL = result.Variables_ELL; _coordinatesReference = result.CoordinatesReference; _referenceFrame = result.ReferenceFrame; Body = result.Body; diff --git a/VSOP2013.NET/VSOPResult/VSOPResult_XYZ.cs b/VSOP2013.NET/VSOPResult/VSOPResult_XYZ.cs index 0308499..afda0a9 100644 --- a/VSOP2013.NET/VSOPResult/VSOPResult_XYZ.cs +++ b/VSOP2013.NET/VSOPResult/VSOPResult_XYZ.cs @@ -38,6 +38,7 @@ public override ReferenceFrame ReferenceFrame public VSOPResult_XYZ(VSOPResult_LBR result) { + Variables_ELL = result.Variables_ELL; _coordinatesReference = result.CoordinatesReference; _referenceFrame = result.ReferenceFrame; Body = result.Body; @@ -47,6 +48,7 @@ public VSOPResult_XYZ(VSOPResult_LBR result) public VSOPResult_XYZ(VSOPResult_ELL result) { + Variables_ELL = result.Variables_ELL; _coordinatesReference = result.CoordinatesReference; _referenceFrame = result.ReferenceFrame; Body = result.Body;