Skip to content

Commit

Permalink
v1.1.6
Browse files Browse the repository at this point in the history 
SIMD Accel Utility. Bug Fix
  • Loading branch information
@kingsznhone
kingsznhone committed Sep 3, 2023
1 parent 6d01a02 commit 2dbf80c
Show file tree
Hide file tree
Showing 6 changed files with 101 additions and 8 deletions.
1 change: 1 addition & 0 deletions Demo/Program.cs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
using System;
using System.Globalization;
using BenchmarkDotNet.Configs;
using BenchmarkDotNet.Running;
using VSOP2013;

Expand Down
6 changes: 6 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
96 changes: 89 additions & 7 deletions VSOP2013.NET/Utility.cs
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
using System.Numerics;
using System.Runtime.Intrinsics;

namespace VSOP2013
{
Expand Down Expand Up @@ -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<double> v1 = Vector256.Create(x / r, y / r, z / r, 0);
Vector256<double> 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<double> v3 = Vector256.Create(-y / (x * x + y * y), x / (x * x + y * y), 0, 0);
Vector256<double> 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 = {
Expand Down Expand Up @@ -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<double> 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<double> 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<double> v3 = Vector256.Create(Math.Sin(b), -r * Math.Cos(b), 0, 0);
Vector256<double> 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 = {
Expand Down Expand Up @@ -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<double> v1 = Vector256.Create(Cphi, -Sphi * Ceps, Sphi * Seps, 0);
Vector256<double> v2 = Vector256.Create(Sphi, Cphi * Ceps, -Cphi * Seps, 0);
Vector256<double> v3 = Vector256.Create(0, Seps, Ceps, 0);
Vector256<double> vv = Vector256.Create(dynamical[0], dynamical[1], dynamical[2], 0);
Vector256<double> 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 },
Expand All @@ -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();
}

/// <summary>
Expand All @@ -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<double> v1 = Vector256.Create(Cphi, Sphi, 0, 0);
Vector256<double> v2 = Vector256.Create(-Sphi * Ceps, Cphi * Ceps, Seps, 0);
Vector256<double> v3 = Vector256.Create(Sphi * Seps, -Cphi * Seps, Ceps, 0);
Vector256<double> vv = Vector256.Create(icrs[0], icrs[1], icrs[2], 0);
Vector256<double> 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 },
Expand Down
2 changes: 1 addition & 1 deletion VSOP2013.NET/VSOP2013.NET.csproj
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
<PlatformTarget>x64</PlatformTarget>
<ImplicitUsings>enable</ImplicitUsings>
<Platforms>x64</Platforms>
<Version>1.1.5</Version>
<Version>1.1.6</Version>
<Title>VSOP2013.NET</Title>
<Authors>KingsZNHONE</Authors>
<Description>VSOP was developed and is maintained (updated with the latest data) by the scientists at the Bureau des Longitudes in Paris.
Expand Down
2 changes: 2 additions & 0 deletions VSOP2013.NET/VSOPResult/VSOPResult_LBR.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ public override ReferenceFrame ReferenceFrame
/// <param name="ell">Elliptic elements</param>
public VSOPResult_LBR(VSOPResult_ELL result)
{
Variables_ELL = result.Variables_ELL;
_coordinatesReference = result.CoordinatesReference;
_referenceFrame = result.ReferenceFrame;
Body = result.Body;
Expand All @@ -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;
Expand Down
2 changes: 2 additions & 0 deletions VSOP2013.NET/VSOPResult/VSOPResult_XYZ.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down

0 comments on commit 2dbf80c

Please sign in to comment.