diff --git a/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/IgrfMagneticFieldBenchmarks.cs b/FMA_Optimizations_Summary.md similarity index 100% rename from benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/IgrfMagneticFieldBenchmarks.cs rename to FMA_Optimizations_Summary.md diff --git a/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/MagneticFieldLineTraceBenchmarks.cs b/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/MagneticFieldLineTraceBenchmarks.cs deleted file mode 100644 index e69de29..0000000 diff --git a/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/results/v2/Issues_42_50_51_Linux.md b/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/results/v2/Issues_42_50_51_Linux.md index e576854..568e515 100644 --- a/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/results/v2/Issues_42_50_51_Linux.md +++ b/benchmarks/AuroraScienceHub.Geopack.Benchmarks/Geopack/results/v2/Issues_42_50_51_Linux.md @@ -9,62 +9,62 @@ AMD Ryzen 5 5500U with Radeon Graphics 0.40GHz, 1 CPU, 12 logical and 6 physical ``` -| Method | Job | Runtime | Mean | Error | StdDev | Median | Ratio | RatioSD | Gen0 | Allocated | Alloc Ratio | -|-------------------|-----------|-----------|----------------:|--------------:|--------------:|----------------:|--------:|--------:|--------:|----------:|------------:| -| ToSphericalVector | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| ToCartesianVector | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| ToSpherical | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0001 ns | 0.0001 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| ToCartesian | .NET 10.0 | .NET 10.0 | 0.0019 ns | 0.0025 ns | 0.0023 ns | 0.0012 ns | 0.000 | 0.00 | - | - | 0.00 | -| GeiToGeo | .NET 10.0 | .NET 10.0 | 0.0017 ns | 0.0014 ns | 0.0012 ns | 0.0015 ns | 0.000 | 0.00 | - | - | 0.00 | -| GeoToGei | .NET 10.0 | .NET 10.0 | 0.0009 ns | 0.0006 ns | 0.0005 ns | 0.0009 ns | 0.000 | 0.00 | - | - | 0.00 | -| ToGeocentric | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| ToGeodetic | .NET 10.0 | .NET 10.0 | 332.4512 ns | 0.3533 ns | 0.3132 ns | 332.3268 ns | 0.290 | 0.00 | - | - | 0.00 | -| GeoToGsw | .NET 10.0 | .NET 10.0 | 0.0007 ns | 0.0008 ns | 0.0007 ns | 0.0006 ns | 0.000 | 0.00 | - | - | 0.00 | -| GswToGeo | .NET 10.0 | .NET 10.0 | 0.0007 ns | 0.0006 ns | 0.0006 ns | 0.0008 ns | 0.000 | 0.00 | - | - | 0.00 | -| GeoToMag | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| MagToGeo | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| GswToGse | .NET 10.0 | .NET 10.0 | 0.0007 ns | 0.0013 ns | 0.0011 ns | 0.0001 ns | 0.000 | 0.00 | - | - | 0.00 | -| GseToGsw | .NET 10.0 | .NET 10.0 | 0.0002 ns | 0.0005 ns | 0.0005 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| Method | Job | Runtime | Mean | Error | StdDev | Median | Ratio | RatioSD | Gen0 | Allocated | Alloc Ratio | +|------------------ |---------- |---------- |----------------:|--------------:|--------------:|----------------:|--------:|--------:|--------:|----------:|------------:| +| ToSphericalVector | .NET 10.0 | .NET 10.0 | 0.0012 ns | 0.0014 ns | 0.0013 ns | 0.0006 ns | 0.000 | 0.00 | - | - | 0.00 | +| ToCartesianVector | .NET 10.0 | .NET 10.0 | 0.0017 ns | 0.0017 ns | 0.0015 ns | 0.0011 ns | 0.000 | 0.00 | - | - | 0.00 | +| ToSpherical | .NET 10.0 | .NET 10.0 | 0.0017 ns | 0.0011 ns | 0.0011 ns | 0.0018 ns | 0.000 | 0.00 | - | - | 0.00 | +| ToCartesian | .NET 10.0 | .NET 10.0 | 0.0001 ns | 0.0002 ns | 0.0002 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| GeiToGeo | .NET 10.0 | .NET 10.0 | 0.0003 ns | 0.0005 ns | 0.0004 ns | 0.0002 ns | 0.000 | 0.00 | - | - | 0.00 | +| GeoToGei | .NET 10.0 | .NET 10.0 | 0.0001 ns | 0.0002 ns | 0.0001 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| ToGeocentric | .NET 10.0 | .NET 10.0 | 0.0001 ns | 0.0002 ns | 0.0001 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| ToGeodetic | .NET 10.0 | .NET 10.0 | 329.3568 ns | 0.4036 ns | 0.3775 ns | 329.5001 ns | 0.290 | 0.00 | - | - | 0.00 | +| GeoToGsw | .NET 10.0 | .NET 10.0 | 0.0000 ns | 0.0001 ns | 0.0001 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| GswToGeo | .NET 10.0 | .NET 10.0 | 0.0001 ns | 0.0005 ns | 0.0004 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| GeoToMag | .NET 10.0 | .NET 10.0 | 0.0010 ns | 0.0019 ns | 0.0017 ns | 0.0002 ns | 0.000 | 0.00 | - | - | 0.00 | +| MagToGeo | .NET 10.0 | .NET 10.0 | 0.0009 ns | 0.0018 ns | 0.0016 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | +| GswToGse | .NET 10.0 | .NET 10.0 | 0.0005 ns | 0.0004 ns | 0.0004 ns | 0.0006 ns | 0.000 | 0.00 | - | - | 0.00 | +| GseToGsw | .NET 10.0 | .NET 10.0 | 0.0024 ns | 0.0016 ns | 0.0014 ns | 0.0019 ns | 0.000 | 0.00 | - | - | 0.00 | | MagToSm | .NET 10.0 | .NET 10.0 | 0.0007 ns | 0.0005 ns | 0.0005 ns | 0.0006 ns | 0.000 | 0.00 | - | - | 0.00 | -| SmToMag | .NET 10.0 | .NET 10.0 | 0.0009 ns | 0.0018 ns | 0.0016 ns | 0.0000 ns | 0.000 | 0.00 | - | - | 0.00 | -| SmToGsw | .NET 10.0 | .NET 10.0 | 0.0013 ns | 0.0014 ns | 0.0013 ns | 0.0008 ns | 0.000 | 0.00 | - | - | 0.00 | -| GswToSm | .NET 10.0 | .NET 10.0 | 0.0031 ns | 0.0006 ns | 0.0006 ns | 0.0033 ns | 0.000 | 0.00 | - | - | 0.00 | -| Recalc | .NET 10.0 | .NET 10.0 | 1,146.5225 ns | 3.5475 ns | 3.3183 ns | 1,146.0720 ns | 1.000 | 0.00 | 1.3828 | 2896 B | 1.00 | -| ShuMgnp | .NET 10.0 | .NET 10.0 | 248.2604 ns | 0.3242 ns | 0.2531 ns | 248.2565 ns | 0.217 | 0.00 | - | - | 0.00 | -| T96Mgnp | .NET 10.0 | .NET 10.0 | 0.0707 ns | 0.0327 ns | 0.0674 ns | 0.0701 ns | 0.000 | 0.00 | - | - | 0.00 | -| IgrfGeo | .NET 10.0 | .NET 10.0 | 455.8333 ns | 1.9287 ns | 1.7097 ns | 455.6252 ns | 0.398 | 0.00 | 0.1297 | 272 B | 0.09 | -| IgrfGsw | .NET 10.0 | .NET 10.0 | 515.2586 ns | 1.7177 ns | 1.5227 ns | 515.1673 ns | 0.449 | 0.00 | 0.1373 | 288 B | 0.10 | -| Dip | .NET 10.0 | .NET 10.0 | 0.2075 ns | 0.0042 ns | 0.0037 ns | 0.2064 ns | 0.000 | 0.00 | - | - | 0.00 | -| Sun | .NET 10.0 | .NET 10.0 | 153.8404 ns | 0.2023 ns | 0.1793 ns | 153.8521 ns | 0.134 | 0.00 | - | - | 0.00 | -| T89 | .NET 10.0 | .NET 10.0 | 124.3691 ns | 0.1685 ns | 0.1494 ns | 124.3627 ns | 0.108 | 0.00 | - | - | 0.00 | -| Trace_NS | .NET 10.0 | .NET 10.0 | 397,184.3338 ns | 691.4255 ns | 612.9307 ns | 397,275.4460 ns | 346.428 | 1.10 | 97.6563 | 204472 B | 70.60 | -| Trace_SN | .NET 10.0 | .NET 10.0 | 376,351.1130 ns | 1,677.9388 ns | 1,487.4491 ns | 376,693.0295 ns | 328.257 | 1.55 | 87.4023 | 183200 B | 63.26 | +| SmToMag | .NET 10.0 | .NET 10.0 | 0.0007 ns | 0.0006 ns | 0.0005 ns | 0.0006 ns | 0.000 | 0.00 | - | - | 0.00 | +| SmToGsw | .NET 10.0 | .NET 10.0 | 0.0005 ns | 0.0006 ns | 0.0005 ns | 0.0004 ns | 0.000 | 0.00 | - | - | 0.00 | +| GswToSm | .NET 10.0 | .NET 10.0 | 0.0007 ns | 0.0008 ns | 0.0007 ns | 0.0006 ns | 0.000 | 0.00 | - | - | 0.00 | +| Recalc | .NET 10.0 | .NET 10.0 | 1,135.7724 ns | 3.9808 ns | 3.5289 ns | 1,134.5793 ns | 1.000 | 0.00 | 1.3828 | 2896 B | 1.00 | +| ShuMgnp | .NET 10.0 | .NET 10.0 | 248.2120 ns | 0.2740 ns | 0.2288 ns | 248.2190 ns | 0.219 | 0.00 | - | - | 0.00 | +| T96Mgnp | .NET 10.0 | .NET 10.0 | 0.0004 ns | 0.0009 ns | 0.0007 ns | 0.0001 ns | 0.000 | 0.00 | - | - | 0.00 | +| IgrfGeo | .NET 10.0 | .NET 10.0 | 454.4769 ns | 1.6005 ns | 1.4188 ns | 453.9866 ns | 0.400 | 0.00 | 0.1297 | 272 B | 0.09 | +| IgrfGsw | .NET 10.0 | .NET 10.0 | 518.6890 ns | 8.4863 ns | 7.9381 ns | 516.7420 ns | 0.457 | 0.01 | 0.1373 | 288 B | 0.10 | +| Dip | .NET 10.0 | .NET 10.0 | 0.2110 ns | 0.0040 ns | 0.0035 ns | 0.2098 ns | 0.000 | 0.00 | - | - | 0.00 | +| Sun | .NET 10.0 | .NET 10.0 | 153.8762 ns | 0.4146 ns | 0.3675 ns | 153.9401 ns | 0.135 | 0.00 | - | - | 0.00 | +| T89 | .NET 10.0 | .NET 10.0 | 124.7912 ns | 0.5287 ns | 0.4945 ns | 124.5573 ns | 0.110 | 0.00 | - | - | 0.00 | +| Trace_NS | .NET 10.0 | .NET 10.0 | 394,777.6689 ns | 1,984.5582 ns | 1,549.4125 ns | 394,729.2681 ns | 347.588 | 1.67 | 97.6563 | 204472 B | 70.60 | +| Trace_SN | .NET 10.0 | .NET 10.0 | 376,446.2077 ns | 6,044.0743 ns | 5,653.6308 ns | 373,395.9683 ns | 331.448 | 4.92 | 87.4023 | 183200 B | 63.26 | | | | | | | | | | | | | | -| ToSphericalVector | .NET 8.0 | .NET 8.0 | 58.6147 ns | 0.1306 ns | 0.1090 ns | 58.6216 ns | 0.047 | 0.00 | - | - | 0.00 | -| ToCartesianVector | .NET 8.0 | .NET 8.0 | 35.8947 ns | 0.1011 ns | 0.0946 ns | 35.8679 ns | 0.029 | 0.00 | - | - | 0.00 | -| ToSpherical | .NET 8.0 | .NET 8.0 | 35.4515 ns | 0.0713 ns | 0.0632 ns | 35.4479 ns | 0.029 | 0.00 | - | - | 0.00 | -| ToCartesian | .NET 8.0 | .NET 8.0 | 0.0151 ns | 0.0003 ns | 0.0003 ns | 0.0150 ns | 0.000 | 0.00 | - | - | 0.00 | -| GeiToGeo | .NET 8.0 | .NET 8.0 | 5.5052 ns | 0.0129 ns | 0.0120 ns | 5.5073 ns | 0.004 | 0.00 | - | - | 0.00 | -| GeoToGei | .NET 8.0 | .NET 8.0 | 5.7522 ns | 0.0145 ns | 0.0135 ns | 5.7537 ns | 0.005 | 0.00 | - | - | 0.00 | -| ToGeocentric | .NET 8.0 | .NET 8.0 | 159.2611 ns | 0.1921 ns | 0.1797 ns | 159.2355 ns | 0.128 | 0.00 | - | - | 0.00 | -| ToGeodetic | .NET 8.0 | .NET 8.0 | 336.6362 ns | 0.4566 ns | 0.4047 ns | 336.6917 ns | 0.271 | 0.00 | - | - | 0.00 | -| GeoToGsw | .NET 8.0 | .NET 8.0 | 6.4986 ns | 0.0073 ns | 0.0065 ns | 6.4988 ns | 0.005 | 0.00 | - | - | 0.00 | -| GswToGeo | .NET 8.0 | .NET 8.0 | 6.2557 ns | 0.0096 ns | 0.0080 ns | 6.2534 ns | 0.005 | 0.00 | - | - | 0.00 | -| GeoToMag | .NET 8.0 | .NET 8.0 | 6.0182 ns | 0.0146 ns | 0.0130 ns | 6.0166 ns | 0.005 | 0.00 | - | - | 0.00 | -| MagToGeo | .NET 8.0 | .NET 8.0 | 6.0025 ns | 0.0039 ns | 0.0033 ns | 6.0035 ns | 0.005 | 0.00 | - | - | 0.00 | -| GswToGse | .NET 8.0 | .NET 8.0 | 6.5121 ns | 0.0172 ns | 0.0161 ns | 6.5089 ns | 0.005 | 0.00 | - | - | 0.00 | -| GseToGsw | .NET 8.0 | .NET 8.0 | 6.2485 ns | 0.0180 ns | 0.0168 ns | 6.2463 ns | 0.005 | 0.00 | - | - | 0.00 | -| MagToSm | .NET 8.0 | .NET 8.0 | 5.5219 ns | 0.0106 ns | 0.0099 ns | 5.5218 ns | 0.004 | 0.00 | - | - | 0.00 | -| SmToMag | .NET 8.0 | .NET 8.0 | 5.7446 ns | 0.0055 ns | 0.0046 ns | 5.7455 ns | 0.005 | 0.00 | - | - | 0.00 | -| SmToGsw | .NET 8.0 | .NET 8.0 | 5.5089 ns | 0.0154 ns | 0.0144 ns | 5.5125 ns | 0.004 | 0.00 | - | - | 0.00 | -| GswToSm | .NET 8.0 | .NET 8.0 | 5.7521 ns | 0.0055 ns | 0.0049 ns | 5.7498 ns | 0.005 | 0.00 | - | - | 0.00 | -| Recalc | .NET 8.0 | .NET 8.0 | 1,241.8564 ns | 4.5656 ns | 4.2707 ns | 1,242.5220 ns | 1.000 | 0.00 | 1.3828 | 2896 B | 1.00 | -| ShuMgnp | .NET 8.0 | .NET 8.0 | 508.6458 ns | 0.6105 ns | 0.4766 ns | 508.7414 ns | 0.410 | 0.00 | - | - | 0.00 | -| T96Mgnp | .NET 8.0 | .NET 8.0 | 112.3719 ns | 0.1178 ns | 0.1045 ns | 112.3408 ns | 0.090 | 0.00 | - | - | 0.00 | -| IgrfGeo | .NET 8.0 | .NET 8.0 | 466.2773 ns | 2.1515 ns | 2.0125 ns | 465.5527 ns | 0.375 | 0.00 | 0.1297 | 272 B | 0.09 | -| IgrfGsw | .NET 8.0 | .NET 8.0 | 519.4500 ns | 3.3455 ns | 3.1294 ns | 519.8844 ns | 0.418 | 0.00 | 0.1373 | 288 B | 0.10 | -| Dip | .NET 8.0 | .NET 8.0 | 129.0686 ns | 0.1325 ns | 0.1174 ns | 129.0692 ns | 0.104 | 0.00 | - | - | 0.00 | -| Sun | .NET 8.0 | .NET 8.0 | 155.7287 ns | 0.3621 ns | 0.3387 ns | 155.5598 ns | 0.125 | 0.00 | - | - | 0.00 | -| T89 | .NET 8.0 | .NET 8.0 | 126.7337 ns | 0.2672 ns | 0.2499 ns | 126.8159 ns | 0.102 | 0.00 | - | - | 0.00 | -| Trace_NS | .NET 8.0 | .NET 8.0 | 436,818.7748 ns | 2,004.9580 ns | 1,777.3430 ns | 437,057.5444 ns | 351.750 | 1.81 | 97.6563 | 204608 B | 70.65 | -| Trace_SN | .NET 8.0 | .NET 8.0 | 410,329.1146 ns | 1,226.9733 ns | 1,024.5783 ns | 410,627.9526 ns | 330.420 | 1.36 | 87.4023 | 183336 B | 63.31 | +| ToSphericalVector | .NET 8.0 | .NET 8.0 | 58.3548 ns | 0.0434 ns | 0.0363 ns | 58.3503 ns | 0.048 | 0.00 | - | - | 0.00 | +| ToCartesianVector | .NET 8.0 | .NET 8.0 | 35.8617 ns | 0.1579 ns | 0.1318 ns | 35.8871 ns | 0.029 | 0.00 | - | - | 0.00 | +| ToSpherical | .NET 8.0 | .NET 8.0 | 35.0748 ns | 0.0731 ns | 0.0648 ns | 35.0476 ns | 0.029 | 0.00 | - | - | 0.00 | +| ToCartesian | .NET 8.0 | .NET 8.0 | 0.0182 ns | 0.0019 ns | 0.0016 ns | 0.0178 ns | 0.000 | 0.00 | - | - | 0.00 | +| GeiToGeo | .NET 8.0 | .NET 8.0 | 5.5140 ns | 0.0152 ns | 0.0142 ns | 5.5104 ns | 0.004 | 0.00 | - | - | 0.00 | +| GeoToGei | .NET 8.0 | .NET 8.0 | 5.7500 ns | 0.0110 ns | 0.0097 ns | 5.7525 ns | 0.005 | 0.00 | - | - | 0.00 | +| ToGeocentric | .NET 8.0 | .NET 8.0 | 158.6567 ns | 0.2616 ns | 0.2447 ns | 158.6431 ns | 0.129 | 0.00 | - | - | 0.00 | +| ToGeodetic | .NET 8.0 | .NET 8.0 | 332.1483 ns | 0.3184 ns | 0.2486 ns | 332.0785 ns | 0.271 | 0.00 | - | - | 0.00 | +| GeoToGsw | .NET 8.0 | .NET 8.0 | 6.5120 ns | 0.0118 ns | 0.0099 ns | 6.5145 ns | 0.005 | 0.00 | - | - | 0.00 | +| GswToGeo | .NET 8.0 | .NET 8.0 | 6.2514 ns | 0.0187 ns | 0.0146 ns | 6.2569 ns | 0.005 | 0.00 | - | - | 0.00 | +| GeoToMag | .NET 8.0 | .NET 8.0 | 6.0040 ns | 0.0188 ns | 0.0157 ns | 5.9964 ns | 0.005 | 0.00 | - | - | 0.00 | +| MagToGeo | .NET 8.0 | .NET 8.0 | 6.0036 ns | 0.0144 ns | 0.0120 ns | 5.9971 ns | 0.005 | 0.00 | - | - | 0.00 | +| GswToGse | .NET 8.0 | .NET 8.0 | 6.6122 ns | 0.0398 ns | 0.0372 ns | 6.6150 ns | 0.005 | 0.00 | - | - | 0.00 | +| GseToGsw | .NET 8.0 | .NET 8.0 | 6.2298 ns | 0.0186 ns | 0.0156 ns | 6.2233 ns | 0.005 | 0.00 | - | - | 0.00 | +| MagToSm | .NET 8.0 | .NET 8.0 | 5.4997 ns | 0.0058 ns | 0.0054 ns | 5.5004 ns | 0.004 | 0.00 | - | - | 0.00 | +| SmToMag | .NET 8.0 | .NET 8.0 | 5.7476 ns | 0.0035 ns | 0.0032 ns | 5.7480 ns | 0.005 | 0.00 | - | - | 0.00 | +| SmToGsw | .NET 8.0 | .NET 8.0 | 5.5090 ns | 0.0069 ns | 0.0061 ns | 5.5096 ns | 0.004 | 0.00 | - | - | 0.00 | +| GswToSm | .NET 8.0 | .NET 8.0 | 5.7518 ns | 0.0091 ns | 0.0080 ns | 5.7499 ns | 0.005 | 0.00 | - | - | 0.00 | +| Recalc | .NET 8.0 | .NET 8.0 | 1,225.5996 ns | 5.8998 ns | 5.2300 ns | 1,225.1386 ns | 1.000 | 0.01 | 1.3828 | 2896 B | 1.00 | +| ShuMgnp | .NET 8.0 | .NET 8.0 | 509.3406 ns | 0.5872 ns | 0.5493 ns | 509.1909 ns | 0.416 | 0.00 | - | - | 0.00 | +| T96Mgnp | .NET 8.0 | .NET 8.0 | 111.8199 ns | 0.1666 ns | 0.1301 ns | 111.8400 ns | 0.091 | 0.00 | - | - | 0.00 | +| IgrfGeo | .NET 8.0 | .NET 8.0 | 462.8087 ns | 1.8867 ns | 1.7649 ns | 462.5069 ns | 0.378 | 0.00 | 0.1297 | 272 B | 0.09 | +| IgrfGsw | .NET 8.0 | .NET 8.0 | 519.6007 ns | 1.3307 ns | 1.2447 ns | 519.7432 ns | 0.424 | 0.00 | 0.1373 | 288 B | 0.10 | +| Dip | .NET 8.0 | .NET 8.0 | 129.1234 ns | 0.2076 ns | 0.1942 ns | 129.0977 ns | 0.105 | 0.00 | - | - | 0.00 | +| Sun | .NET 8.0 | .NET 8.0 | 158.6323 ns | 0.1673 ns | 0.1483 ns | 158.6875 ns | 0.129 | 0.00 | - | - | 0.00 | +| T89 | .NET 8.0 | .NET 8.0 | 126.9628 ns | 0.2185 ns | 0.2043 ns | 126.9991 ns | 0.104 | 0.00 | - | - | 0.00 | +| Trace_NS | .NET 8.0 | .NET 8.0 | 435,764.6175 ns | 1,343.5597 ns | 1,121.9332 ns | 435,625.8623 ns | 355.558 | 1.71 | 97.6563 | 204608 B | 70.65 | +| Trace_SN | .NET 8.0 | .NET 8.0 | 410,702.3057 ns | 2,027.2870 ns | 1,896.3255 ns | 410,279.8330 ns | 335.109 | 2.04 | 87.4023 | 183336 B | 63.31 | diff --git a/docs/ReleaseNotes_GeopackV2.md b/docs/ReleaseNotes_GeopackV2.md index a0d901b..d33b2d8 100644 --- a/docs/ReleaseNotes_GeopackV2.md +++ b/docs/ReleaseNotes_GeopackV2.md @@ -30,3 +30,34 @@ This PR adds dependency injection support for the Geopack library by introducing - Added AddExternalFieldModels() extension method to register IT89 implementation - Added Microsoft.Extensions.DependencyInjection.Abstractions package dependency (version 9.0.10) +## Issue #46 (PR #70): + +This PR implements comprehensive code optimizations focused on improving computational performance across the Geopack library. +The optimizations target mathematical operations, memory allocation, and leverage modern .NET features including SIMD vectorization. + +Key changes: + +- Replaced Math.Pow(x, 2) with direct multiplication x * x for better performance; +- Adopted Math.SinCos() to compute sine and cosine simultaneously, reducing redundant calculations; +- Extracted magic numbers into well-documented constants in GeopackConstants.cs; +- Implemented SIMD vectorization using Vector for IGRF coefficient interpolation and extrapolation; +- Refactored Newton's method iteration into a separate method for better code organization; +- Optimized list initialization with capacity hints to reduce allocations. + +| File | Description | +|------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------------| +| `src/Geopack/Geopack.Trace.cs` | Optimized list capacity, dot product calculation | +| `src/Geopack/Geopack.T96Mgnp.cs` | Replaced Math.Pow with multiplication, used Math.SinCos, extracted pressure factor constant | +| `src/Geopack/Geopack.Sun.cs` | Extracted constants, adopted Math.SinCos | +| `src/Geopack/Geopack.ShuMgnp.cs` | Extracted model coefficients as constants, refactored Newton's method to private method, optimized trigonometric calculations | +| `src/Geopack/Geopack.Recalc.cs` | Implemented SIMD vectorization for coefficient processing, eliminated zero multiplications, cached repeated calculations | +| `src/Geopack/Geopack.IgrfGsw.cs` | Replaced Math.Pow with direct multiplication | +| `src/Geopack/Geopack.IgrfGeo.cs` | Adopted Math.SinCos for simultaneous trigonometric calculations | +| `src/Geopack/Geopack.Dip.cs` | Replaced Math.Pow with multiplication, cached repeated array accesses | +| `src/ExternalFieldModels/T89/T89.cs` | Adopted Math.SinCos | +| `src/Contracts/Spherical/SphericalVector.cs` | Minor formatting adjustment | +| `src/Contracts/Spherical/SphericalLocation.cs` | Cached trigonometric calculations to avoid recomputation | +| `src/Contracts/Engine/GeopackConstants.cs` | Added numerous well-documented constants for physics calculations and algorithm parameters | +| `src/Contracts/Coordinates/GeodeticCoordinates.cs` | Extracted WGS84 constants, replaced Math.Pow, adopted Math.SinCos | +| `src/Contracts/Coordinates/GeocentricCoordinates.cs` | Optimized iterative calculation with constant extraction and Math.SinCos | +| `src/Contracts/Cartesian/CartesianVector.cs` | Replaced Math.Pow with direct multiplication | diff --git a/src/Contracts/Cartesian/CartesianLocation.cs b/src/Contracts/Cartesian/CartesianLocation.cs index 386725b..d5f4084 100644 --- a/src/Contracts/Cartesian/CartesianLocation.cs +++ b/src/Contracts/Cartesian/CartesianLocation.cs @@ -55,8 +55,8 @@ public SphericalLocation ToSpherical() { double phi; double theta; - double sq = X * X + Y * Y; - double r = Math.Sqrt(sq + Z * Z); + double sq = Math.FusedMultiplyAdd(X, X, Y * Y); + double r = Math.Sqrt(Math.FusedMultiplyAdd(Z, Z, sq)); if (Math.Abs(sq) > double.Epsilon) { diff --git a/src/Contracts/Cartesian/CartesianVector.cs b/src/Contracts/Cartesian/CartesianVector.cs index 32bf58e..aa0dd10 100644 --- a/src/Contracts/Cartesian/CartesianVector.cs +++ b/src/Contracts/Cartesian/CartesianVector.cs @@ -37,10 +37,10 @@ private CartesianVector(double x, double y, double z, CoordinateSystem coordinat /// public SphericalVector ToSphericalVector(CartesianLocation location) { - double rho2 = Math.Pow(location.X, 2.0D) + Math.Pow(location.Y, 2.0D); + double rho2 = location.X * location.X + location.Y * location.Y; double rho = Math.Sqrt(rho2); - double r = Math.Sqrt(rho2 + Math.Pow(location.Z, 2.0D)); + double r = Math.Sqrt(rho2 + location.Z * location.Z); if (Math.Abs(r) <= double.Epsilon) { @@ -64,9 +64,9 @@ public SphericalVector ToSphericalVector(CartesianLocation location) double ct = location.Z / r; double st = rho / r; - double br = (location.X * X + location.Y * Y + location.Z * Z) / r; - double btheta = (X * cphi + Y * sphi) * ct - Z * st; - double bphi = Y * cphi - X * sphi; + double br = Math.FusedMultiplyAdd(location.X, X, Math.FusedMultiplyAdd(location.Y, Y, location.Z * Z)) / r; + double btheta = Math.FusedMultiplyAdd(Math.FusedMultiplyAdd(X, cphi, Y * sphi), ct, -Z * st); + double bphi = Math.FusedMultiplyAdd(Y, cphi, -X * sphi); return SphericalVector.New(br, btheta, bphi, CoordinateSystem); } diff --git a/src/Contracts/Coordinates/GeocentricCoordinates.cs b/src/Contracts/Coordinates/GeocentricCoordinates.cs index 61115d0..0b1c99f 100644 --- a/src/Contracts/Coordinates/GeocentricCoordinates.cs +++ b/src/Contracts/Coordinates/GeocentricCoordinates.cs @@ -1,3 +1,5 @@ +using AuroraScienceHub.Geopack.Contracts.Engine; + namespace AuroraScienceHub.Geopack.Contracts.Coordinates; /// @@ -37,32 +39,33 @@ public GeocentricCoordinates(double r, double theta) /// public GeodeticCoordinates ToGeodetic() { - const double r_eq = 6378.137D; - const double beta = 6.73949674228e-3; - const double tol = 1e-6; - int n = 0; - double phi = 1.570796327D - Theta; + double phi = GeopackConstants.HalfPi - Theta; double phi1 = phi; + double r2 = R * R; - double xmus, rs, cosfims, h, z, x, rr, dphi; + double xmus, h, dphi; do { - double sp = Math.Sin(phi1); - double arg = sp * (1.0D + beta) / Math.Sqrt(1.0D + beta * (2.0D + beta) * Math.Pow(sp, 2)); + (double sinPhi, double cosPhi) = Math.SinCos(phi1); + double sp2 = sinPhi * sinPhi; + double arg = sinPhi * GeopackConstants.WGS84Ex / Math.Sqrt(1.0D + GeopackConstants.WGS84FirstEx * sp2); + double rs = GeopackConstants.REq / Math.Sqrt(1.0D + GeopackConstants.WGS84Beta * sp2); + double rs2 = rs * rs; xmus = Math.Asin(arg); - rs = r_eq / Math.Sqrt(1.0D + beta * Math.Pow(Math.Sin(phi1), 2)); - cosfims = Math.Cos(phi1 - xmus); - h = Math.Sqrt(rs * cosfims * rs * cosfims + R * R - rs * rs) - rs * cosfims; - z = rs * Math.Sin(phi1) + h * Math.Sin(xmus); - x = rs * Math.Cos(phi1) + h * Math.Cos(xmus); - rr = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(z, 2)); + (double sinXmus, double cosXmus) = Math.SinCos(xmus); + + double cosfims = Math.Cos(phi1 - xmus); + h = Math.Sqrt(rs2 * cosfims * cosfims + r2 - rs2) - rs * cosfims; + double z = rs * sinPhi + h * sinXmus; + double x = rs * cosPhi + h * cosXmus; + double rr = Math.Sqrt(x * x + z * z); dphi = Math.Asin(z / rr) - phi; phi1 -= dphi; n++; } - while (Math.Abs(dphi) > tol && n < 100); + while (Math.Abs(dphi) > GeopackConstants.CoordinateConvergenceTolerance && n < 100); return new GeodeticCoordinates(xmus, h); } diff --git a/src/Contracts/Coordinates/GeodeticCoordinates.cs b/src/Contracts/Coordinates/GeodeticCoordinates.cs index b5b4c63..162355d 100644 --- a/src/Contracts/Coordinates/GeodeticCoordinates.cs +++ b/src/Contracts/Coordinates/GeodeticCoordinates.cs @@ -1,3 +1,5 @@ +using AuroraScienceHub.Geopack.Contracts.Engine; + namespace AuroraScienceHub.Geopack.Contracts.Coordinates; /// @@ -34,18 +36,16 @@ public GeodeticCoordinates(double geodLatitude, double altitude) /// public GeocentricCoordinates ToGeocentric() { - const double r_eq = 6378.137D; - const double beta = 6.73949674228e-3; - - double cosxmu = Math.Cos(Latitude); double sinxmu = Math.Sin(Latitude); - double den = Math.Sqrt(Math.Pow(cosxmu, 2) + Math.Pow(sinxmu / (1.0D + beta), 2)); + double cosxmu = Math.Cos(Latitude); + double sinxmuBeta = sinxmu / GeopackConstants.WGS84Ex; + double den = Math.Sqrt(cosxmu * cosxmu + sinxmuBeta * sinxmuBeta); double coslam = cosxmu / den; - double sinlam = sinxmu / (den * (1.0D + beta)); - double rs = r_eq / Math.Sqrt(1.0D + beta * Math.Pow(sinlam, 2)); + double sinlam = sinxmu / (den * GeopackConstants.WGS84Ex); + double rs = GeopackConstants.REq / Math.Sqrt(1.0D + GeopackConstants.WGS84Beta * sinlam * sinlam); double x = rs * coslam + Altitude * cosxmu; double z = rs * sinlam + Altitude * sinxmu; - double r = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(z, 2)); + double r = Math.Sqrt(x * x + z * z); double theta = Math.Acos(z / r); return new GeocentricCoordinates(r, theta); diff --git a/src/Contracts/Engine/GeopackConstants.cs b/src/Contracts/Engine/GeopackConstants.cs index a3ef54b..c7cf168 100644 --- a/src/Contracts/Engine/GeopackConstants.cs +++ b/src/Contracts/Engine/GeopackConstants.cs @@ -10,6 +10,11 @@ public static class GeopackConstants /// public const double Pi = 3.141592654D; + /// + /// PI/2 value from original Geopack (lower precision for compatibility) + /// + public const double HalfPi = 1.570796327D; + /// /// 2π value from original Geopack /// @@ -19,4 +24,106 @@ public static class GeopackConstants /// Radians to degrees conversion factor (180/π) /// public const double Rad = 57.295779513D; + + /// + /// Equatorial Earth's radius in kilometers + /// + public const double REq = 6378.137D; + + /// + /// WGS84 ellipsoid flattening coefficient (Beta = (a-b)/a where a is equatorial radius, b is polar radius) + /// + public const double WGS84Beta = 6.73949674228e-3; + + /// + /// Convergence tolerance for iterative coordinate transformations + /// + public const double CoordinateConvergenceTolerance = 1e-6; + + /// + /// WGS84 ellipsoid parameter: 1 + Beta + /// + public const double WGS84Ex = 1.0 + WGS84Beta; + + /// + /// WGS84 ellipsoid parameter: Beta * (2 + Beta) + /// + public const double WGS84FirstEx = WGS84Beta * (2.0 + WGS84Beta); + + /// + /// Earth's axial tilt at J2000 epoch (obliquity of the ecliptic) in degrees + /// + public const double EarthObliquityJ2000 = 23.45229; + + /// + /// Rate of change of Earth's obliquity per Julian century in degrees + /// + public const double EarthObliquityRatePerCentury = 0.0130125; + + /// + /// Number of days in a Julian century + /// + public const double DaysPerJulianCentury = 36525.0; + + /// + /// Solar wind proton mass density conversion factor (converts n*v^2 to pressure in nPa) + /// Formula: 1.6726e-27 kg (proton mass) * 1e15 (km^3 to m^3) * 1e9 (Pa to nPa) = 1.6726e-6 + /// But empirical value used in models is 1.94e-6 + /// + public const double SolarWindDynamicPressureFactor = 1.94e-6; + + /// + /// Seconds per hour + /// + public const double SecondsPerHour = 3600.0; + + /// + /// Seconds per day + /// + public const double SecondsPerDay = 86400.0; + + /// + /// Degrees in a full circle + /// + public const double DegreesPerCircle = 360.0; + + /// + /// Degrees in a half circle (semicircle) + /// + public const double DegreesPerSemicircle = 180.0; + + /// + /// Average number of days in a year (accounting for leap years) + /// + public const double DaysPerYear = 365.25; + + /// + /// Reciprocal of IGRF interpolation interval (1/5 = 0.2) for optimization + /// + public const double IgrfInterpolationIntervalReciprocal = 0.2; + + /// + /// Total number of IGRF coefficients + /// + public const int IgrfCoefficientCount = 105; + + /// + /// Number of IGRF delta coefficients used in extrapolation + /// + public const int IgrfDeltaCoefficientCount = 45; + + /// + /// IGRF extrapolation base year + /// + public const int IgrfExtrapolationBaseYear = 2025; + + /// + /// Maximum iterations for Newton's method convergence + /// + public const int NewtonMaxIterations = 1000; + + /// + /// Convergence tolerance for Newton's iterative methods + /// + public const double NewtonConvergenceTolerance = 1e-4; } diff --git a/src/Contracts/Spherical/SphericalLocation.cs b/src/Contracts/Spherical/SphericalLocation.cs index c967b9c..367a155 100644 --- a/src/Contracts/Spherical/SphericalLocation.cs +++ b/src/Contracts/Spherical/SphericalLocation.cs @@ -35,10 +35,14 @@ public static SphericalLocation New(double r, double theta, double phi, Coordina /// public CartesianLocation ToCartesian() { - double sq = R * Math.Sin(Theta); - double x = sq * Math.Cos(Phi); - double y = sq * Math.Sin(Phi); - double z = R * Math.Cos(Theta); + double sinTheta = Math.Sin(Theta); + double cosTheta = Math.Cos(Theta); + double sinPhi = Math.Sin(Phi); + double cosPhi = Math.Cos(Phi); + double sq = R * sinTheta; + double x = sq * cosPhi; + double y = sq * sinPhi; + double z = R * cosTheta; return CartesianLocation.New(x, y, z, CoordinateSystem); } diff --git a/src/Contracts/Spherical/SphericalVector.cs b/src/Contracts/Spherical/SphericalVector.cs index f35408d..ac5f1ad 100644 --- a/src/Contracts/Spherical/SphericalVector.cs +++ b/src/Contracts/Spherical/SphericalVector.cs @@ -42,10 +42,11 @@ public CartesianVector ToCartesianVector(SphericalLocation location) double c = Math.Cos(location.Theta); double sf = Math.Sin(location.Phi); double cf = Math.Cos(location.Phi); - double be = R * s + Theta * c; - double bx = be * cf - Phi * sf; - double by = be * sf + Phi * cf; - double bz = R * c - Theta * s; + + double be = Math.FusedMultiplyAdd(R, s, Theta * c); + double bx = Math.FusedMultiplyAdd(be, cf, -Phi * sf); + double by = Math.FusedMultiplyAdd(be, sf, Phi * cf); + double bz = Math.FusedMultiplyAdd(R, c, -Theta * s); return CartesianVector.New(bx, by, bz, CoordinateSystem); } diff --git a/src/ExternalFieldModels/T89/T89.cs b/src/ExternalFieldModels/T89/T89.cs index 876874b..348712a 100644 --- a/src/ExternalFieldModels/T89/T89.cs +++ b/src/ExternalFieldModels/T89/T89.cs @@ -101,8 +101,7 @@ public CartesianVector Calculate(int IOPT, double[] PARMOD, doubl YND = 2.0D * YN; } - double SPS = Math.Sin(PS); - double CPS = Math.Cos(PS); + (double SPS, double CPS) = Math.SinCos(PS); double X2 = location.X * location.X; double Y2 = location.Y * location.Y; @@ -249,8 +248,8 @@ public CartesianVector Calculate(int IOPT, double[] PARMOD, doubl double FXMN = -location.X * FXYM; double FYPL = location.Y * FXYP; double FYMN = -location.Y * FXYM; - double FZPL = WCSP + XYWC * SZRP; - double FZMN = WCSM + XYWC * SZRM; + double FZPL = Math.FusedMultiplyAdd(WCSP, 1.0, XYWC * SZRP); + double FZMN = Math.FusedMultiplyAdd(WCSM, 1.0, XYWC * SZRM); double DER13 = FXPL + FXMN; double DER14 = (FXPL - FXMN) * SPS; double DER23 = FYPL + FYMN; diff --git a/src/Geopack/Geopack.Dip.cs b/src/Geopack/Geopack.Dip.cs index 639a758..3bb4e72 100644 --- a/src/Geopack/Geopack.Dip.cs +++ b/src/Geopack/Geopack.Dip.cs @@ -14,12 +14,15 @@ public CartesianVector Dip(ComputationContext context, CartesianL throw new InvalidOperationException("Location must be in GSW coordinate system."); } - double dipmom = Math.Sqrt(Math.Pow(context.G[1], 2) + Math.Pow(context.G[2], 2) + Math.Pow(context.H[2], 2)); + double g1 = context.G[1]; + double g2 = context.G[2]; + double h2 = context.H[2]; + double dipmom = Math.Sqrt(g1 * g1 + g2 * g2 + h2 * h2); - double p = Math.Pow(location.X, 2); - double u = Math.Pow(location.Z, 2); + double p = location.X * location.X; + double u = location.Z * location.Z; double v = 3.0D * location.Z * location.X; - double t = Math.Pow(location.Y, 2); + double t = location.Y * location.Y; if (p + t + u is 0D) { diff --git a/src/Geopack/Geopack.GeiGeo.cs b/src/Geopack/Geopack.GeiGeo.cs index 09355e0..efaa652 100644 --- a/src/Geopack/Geopack.GeiGeo.cs +++ b/src/Geopack/Geopack.GeiGeo.cs @@ -18,16 +18,16 @@ private static T GeiGeoInternal(ComputationContext context, T components, Ope { OperationType.Direct => components.CoordinateSystem is CoordinateSystem.GEI ? T.New( - components.X * context.CGST + components.Y * context.SGST, - components.Y * context.CGST - components.X * context.SGST, + Math.FusedMultiplyAdd(components.X, context.CGST, components.Y * context.SGST), + Math.FusedMultiplyAdd(components.Y, context.CGST, -components.X * context.SGST), components.Z, CoordinateSystem.GEO) : throw new InvalidOperationException("Input coordinates must be in GEI system."), OperationType.Reversed => components.CoordinateSystem is CoordinateSystem.GEO ? T.New( - components.X * context.CGST - components.Y * context.SGST, - components.Y * context.CGST + components.X * context.SGST, + Math.FusedMultiplyAdd(components.X, context.CGST, -components.Y * context.SGST), + Math.FusedMultiplyAdd(components.Y, context.CGST, components.X * context.SGST), components.Z, CoordinateSystem.GEI) : throw new InvalidOperationException("Input coordinates must be in GEO system."), diff --git a/src/Geopack/Geopack.GeoGsw.cs b/src/Geopack/Geopack.GeoGsw.cs index 1d578e3..4c3f23f 100644 --- a/src/Geopack/Geopack.GeoGsw.cs +++ b/src/Geopack/Geopack.GeoGsw.cs @@ -18,17 +18,17 @@ private static T GeoGswInternal(ComputationContext context, T components, Ope { OperationType.Direct => components.CoordinateSystem is CoordinateSystem.GEO ? T.New( - context.A11 * components.X + context.A12 * components.Y + context.A13 * components.Z, - context.A21 * components.X + context.A22 * components.Y + context.A23 * components.Z, - context.A31 * components.X + context.A32 * components.Y + context.A33 * components.Z, + Math.FusedMultiplyAdd(context.A11, components.X, Math.FusedMultiplyAdd(context.A12, components.Y, context.A13 * components.Z)), + Math.FusedMultiplyAdd(context.A21, components.X, Math.FusedMultiplyAdd(context.A22, components.Y, context.A23 * components.Z)), + Math.FusedMultiplyAdd(context.A31, components.X, Math.FusedMultiplyAdd(context.A32, components.Y, context.A33 * components.Z)), CoordinateSystem.GSW) : throw new InvalidOperationException("Input coordinates must be in GEO system."), OperationType.Reversed => components.CoordinateSystem is CoordinateSystem.GSW ? T.New( - context.A11 * components.X + context.A21 * components.Y + context.A31 * components.Z, - context.A12 * components.X + context.A22 * components.Y + context.A32 * components.Z, - context.A13 * components.X + context.A23 * components.Y + context.A33 * components.Z, + Math.FusedMultiplyAdd(context.A11, components.X, Math.FusedMultiplyAdd(context.A21, components.Y, context.A31 * components.Z)), + Math.FusedMultiplyAdd(context.A12, components.X, Math.FusedMultiplyAdd(context.A22, components.Y, context.A32 * components.Z)), + Math.FusedMultiplyAdd(context.A13, components.X, Math.FusedMultiplyAdd(context.A23, components.Y, context.A33 * components.Z)), CoordinateSystem.GEO) : throw new InvalidOperationException("Input coordinates must be in GSW system."), _ => throw new NotSupportedException($"Specify correct OperationType: {operation}. Available types are Direct and Reversed.") diff --git a/src/Geopack/Geopack.GeoMag.cs b/src/Geopack/Geopack.GeoMag.cs index 02ed88c..76c8b1c 100644 --- a/src/Geopack/Geopack.GeoMag.cs +++ b/src/Geopack/Geopack.GeoMag.cs @@ -18,17 +18,17 @@ private static T GeoMagInternal(ComputationContext context, T components, Ope { OperationType.Direct => components.CoordinateSystem is CoordinateSystem.GEO ? T.New( - components.X * context.CTCL + components.Y * context.CTSL - components.Z * context.ST0, - components.Y * context.CL0 - components.X * context.SL0, - components.X * context.STCL + components.Y * context.STSL + components.Z * context.CT0, + Math.FusedMultiplyAdd(components.X, context.CTCL, Math.FusedMultiplyAdd(components.Y, context.CTSL, -components.Z * context.ST0)), + Math.FusedMultiplyAdd(components.Y, context.CL0, -components.X * context.SL0), + Math.FusedMultiplyAdd(components.X, context.STCL, Math.FusedMultiplyAdd(components.Y, context.STSL, components.Z * context.CT0)), CoordinateSystem.MAG) : throw new InvalidOperationException("Input coordinates must be in GEO system."), OperationType.Reversed => components.CoordinateSystem is CoordinateSystem.MAG ? T.New( - components.X * context.CTCL - components.Y * context.SL0 + components.Z * context.STCL, - components.X * context.CTSL + components.Y * context.CL0 + components.Z * context.STSL, - components.Z * context.CT0 - components.X * context.ST0, + Math.FusedMultiplyAdd(components.X, context.CTCL, Math.FusedMultiplyAdd(-components.Y, context.SL0, components.Z * context.STCL)), + Math.FusedMultiplyAdd(components.X, context.CTSL, Math.FusedMultiplyAdd(components.Y, context.CL0, components.Z * context.STSL)), + Math.FusedMultiplyAdd(components.Z, context.CT0, -components.X * context.ST0), CoordinateSystem.GEO) : throw new InvalidOperationException("Input coordinates must be in MAG system."), _ => throw new NotSupportedException($"Specify correct OperationType: {operation}. Available types are Direct and Reversed.") diff --git a/src/Geopack/Geopack.GswGse.cs b/src/Geopack/Geopack.GswGse.cs index 7a505f5..4ca4b5e 100644 --- a/src/Geopack/Geopack.GswGse.cs +++ b/src/Geopack/Geopack.GswGse.cs @@ -18,17 +18,17 @@ private static T GswGseInternal(ComputationContext context, T components, Ope { OperationType.Direct => components.CoordinateSystem is CoordinateSystem.GSW ? T.New( - context.E11 * components.X + context.E12 * components.Y + context.E13 * components.Z, - context.E21 * components.X + context.E22 * components.Y + context.E23 * components.Z, - context.E31 * components.X + context.E32 * components.Y + context.E33 * components.Z, + Math.FusedMultiplyAdd(context.E11, components.X, Math.FusedMultiplyAdd(context.E12, components.Y, context.E13 * components.Z)), + Math.FusedMultiplyAdd(context.E21, components.X, Math.FusedMultiplyAdd(context.E22, components.Y, context.E23 * components.Z)), + Math.FusedMultiplyAdd(context.E31, components.X, Math.FusedMultiplyAdd(context.E32, components.Y, context.E33 * components.Z)), CoordinateSystem.GSE) : throw new InvalidOperationException("Input coordinates must be in GSW system."), OperationType.Reversed => components.CoordinateSystem is CoordinateSystem.GSE ? T.New( - context.E11 * components.X + context.E21 * components.Y + context.E31 * components.Z, - context.E12 * components.X + context.E22 * components.Y + context.E32 * components.Z, - context.E13 * components.X + context.E23 * components.Y + context.E33 * components.Z, + Math.FusedMultiplyAdd(context.E11, components.X, Math.FusedMultiplyAdd(context.E21, components.Y, context.E31 * components.Z)), + Math.FusedMultiplyAdd(context.E12, components.X, Math.FusedMultiplyAdd(context.E22, components.Y, context.E32 * components.Z)), + Math.FusedMultiplyAdd(context.E13, components.X, Math.FusedMultiplyAdd(context.E23, components.Y, context.E33 * components.Z)), CoordinateSystem.GSW) : throw new InvalidOperationException("Input coordinates must be in GSE system."), _ => throw new NotSupportedException($"Specify correct OperationType: {operation}. Available types are Direct and Reversed.") diff --git a/src/Geopack/Geopack.IgrfGeo.cs b/src/Geopack/Geopack.IgrfGeo.cs index 9c4f567..aa38196 100644 --- a/src/Geopack/Geopack.IgrfGeo.cs +++ b/src/Geopack/Geopack.IgrfGeo.cs @@ -19,10 +19,8 @@ public SphericalVector IgrfGeo(ComputationContext context, Spheri throw new InvalidOperationException("Location must be in GEO coordinate system."); } - double c = Math.Cos(location.Theta); - double s = Math.Sin(location.Theta); - double cf = Math.Cos(location.Phi); - double sf = Math.Sin(location.Phi); + (double s, double c) = Math.SinCos(location.Theta); + (double sf, double cf) = Math.SinCos(location.Phi); double pp = 1.0D / location.R; double p = pp; diff --git a/src/Geopack/Geopack.IgrfGsw.cs b/src/Geopack/Geopack.IgrfGsw.cs index 3662602..7c8fe3e 100644 --- a/src/Geopack/Geopack.IgrfGsw.cs +++ b/src/Geopack/Geopack.IgrfGsw.cs @@ -16,8 +16,8 @@ public CartesianVector IgrfGsw(ComputationContext context, Cartes CartesianLocation geoLocation = GswToGeo(context, location); - double rho2 = Math.Pow(geoLocation.X, 2) + Math.Pow(geoLocation.Y, 2); - double r = Math.Sqrt(rho2 + Math.Pow(geoLocation.Z, 2)); + double rho2 = geoLocation.X * geoLocation.X + geoLocation.Y * geoLocation.Y; + double r = Math.Sqrt(rho2 + geoLocation.Z * geoLocation.Z); if (Math.Abs(r) <= double.Epsilon) { @@ -142,10 +142,10 @@ public CartesianVector IgrfGsw(ComputationContext context, Cartes bf = bbf / s; } - double he = br * s + bt * c; - double bx = he * cf - bf * sf; - double by = he * sf + bf * cf; - double bz = br * c - bt * s; + double he = Math.FusedMultiplyAdd(br, s, bt * c); + double bx = Math.FusedMultiplyAdd(he, cf, -bf * sf); + double by = Math.FusedMultiplyAdd(he, sf, bf * cf); + double bz = Math.FusedMultiplyAdd(br, c, -bt * s); return GeoToGsw(context, CartesianVector.New(bx, by, bz, CoordinateSystem.GEO)); } diff --git a/src/Geopack/Geopack.MagSm.cs b/src/Geopack/Geopack.MagSm.cs index 9239ca6..b389b1f 100644 --- a/src/Geopack/Geopack.MagSm.cs +++ b/src/Geopack/Geopack.MagSm.cs @@ -18,16 +18,16 @@ private static T MagSmInternal(ComputationContext context, T components, Oper { OperationType.Direct => components.CoordinateSystem is CoordinateSystem.MAG ? T.New( - components.X * context.CFI - components.Y * context.SFI, - components.X * context.SFI + components.Y * context.CFI, + Math.FusedMultiplyAdd(components.X, context.CFI, -components.Y * context.SFI), + Math.FusedMultiplyAdd(components.X, context.SFI, components.Y * context.CFI), components.Z, CoordinateSystem.SM) : throw new InvalidOperationException("Input coordinates must be in MAG system."), OperationType.Reversed => components.CoordinateSystem is CoordinateSystem.SM ? T.New( - components.X * context.CFI + components.Y * context.SFI, - components.Y * context.CFI - components.X * context.SFI, + Math.FusedMultiplyAdd(components.X, context.CFI, components.Y * context.SFI), + Math.FusedMultiplyAdd(components.Y, context.CFI, -components.X * context.SFI), components.Z, CoordinateSystem.MAG) : throw new InvalidOperationException("Input coordinates must be in SM system."), diff --git a/src/Geopack/Geopack.Recalc.cs b/src/Geopack/Geopack.Recalc.cs index fc21638..d5a9174 100644 --- a/src/Geopack/Geopack.Recalc.cs +++ b/src/Geopack/Geopack.Recalc.cs @@ -1,3 +1,4 @@ +using System.Numerics; using AuroraScienceHub.Geopack.Contracts.Cartesian; using AuroraScienceHub.Geopack.Contracts.Coordinates; using AuroraScienceHub.Geopack.Contracts.Engine; @@ -13,35 +14,34 @@ public ComputationContext Recalc(DateTime dateTime, CartesianVector? s { swVelocity ??= CartesianVector.New(-400D, 0D, 0D, CoordinateSystem.GSE); - // if (swVelocity.Required().CoordinateSystem is not CoordinateSystem.GSE) if (swVelocity.Required().CoordinateSystem is not CoordinateSystem.GSE) { throw new InvalidOperationException("Solar wind velocity must be in GSE coordinate system."); } - int IY = dateTime.Year; - int IDAY = dateTime.DayOfYear; - int IHOUR = dateTime.Hour; - int MIN = dateTime.Minute; - int ISEC = dateTime.Second; + int year = dateTime.Year; + int doy = dateTime.DayOfYear; + int hour = dateTime.Hour; + int minutes = dateTime.Minute; + int seconds = dateTime.Second; - if (IY < 1965) + if (year < 1965) { - IY = 1965; - Console.WriteLine($"Warning: Year {dateTime.Year} is out of range. Using {IY} instead."); + year = 1965; + Console.WriteLine($"Warning: Year {dateTime.Year} is out of range. Using {year} instead."); } - if (IY > 2030) + if (year > 2030) { - IY = 2030; - Console.WriteLine($"Warning: Year {dateTime.Year} is out of range. Using {IY} instead."); + year = 2030; + Console.WriteLine($"Warning: Year {dateTime.Year} is out of range. Using {year} instead."); } double[] REC = new double[105]; for (int N = 1; N <= 14; N++) { int N2 = 2 * N - 1; - N2 *= (N2 - 2); + N2 *= N2 - 2; for (int M = 1; M <= N; M++) { int MN = N * (N - 1) / 2 + M; @@ -50,46 +50,46 @@ public ComputationContext Recalc(DateTime dateTime, CartesianVector? s } double[] G, H; - switch (IY) + switch (year) { case < 1970: - (G, H) = Interpolate(1965, IY, IDAY, IgrfCoefficients.G65, IgrfCoefficients.G70, IgrfCoefficients.H65, IgrfCoefficients.H70); + (G, H) = Interpolate(1965, year, doy, IgrfCoefficients.G65, IgrfCoefficients.G70, IgrfCoefficients.H65, IgrfCoefficients.H70); break; case < 1975: - (G, H) = Interpolate(1970, IY, IDAY, IgrfCoefficients.G70, IgrfCoefficients.G75, IgrfCoefficients.H70, IgrfCoefficients.H75); + (G, H) = Interpolate(1970, year, doy, IgrfCoefficients.G70, IgrfCoefficients.G75, IgrfCoefficients.H70, IgrfCoefficients.H75); break; case < 1980: - (G, H) = Interpolate(1975, IY, IDAY, IgrfCoefficients.G75, IgrfCoefficients.G80, IgrfCoefficients.H75, IgrfCoefficients.H80); + (G, H) = Interpolate(1975, year, doy, IgrfCoefficients.G75, IgrfCoefficients.G80, IgrfCoefficients.H75, IgrfCoefficients.H80); break; case < 1985: - (G, H) = Interpolate(1980, IY, IDAY, IgrfCoefficients.G80, IgrfCoefficients.G85, IgrfCoefficients.H80, IgrfCoefficients.H85); + (G, H) = Interpolate(1980, year, doy, IgrfCoefficients.G80, IgrfCoefficients.G85, IgrfCoefficients.H80, IgrfCoefficients.H85); break; case < 1990: - (G, H) = Interpolate(1985, IY, IDAY, IgrfCoefficients.G85, IgrfCoefficients.G90, IgrfCoefficients.H85, IgrfCoefficients.H90); + (G, H) = Interpolate(1985, year, doy, IgrfCoefficients.G85, IgrfCoefficients.G90, IgrfCoefficients.H85, IgrfCoefficients.H90); break; case < 1995: - (G, H) = Interpolate(1990, IY, IDAY, IgrfCoefficients.G90, IgrfCoefficients.G95, IgrfCoefficients.H90, IgrfCoefficients.H95); + (G, H) = Interpolate(1990, year, doy, IgrfCoefficients.G90, IgrfCoefficients.G95, IgrfCoefficients.H90, IgrfCoefficients.H95); break; case < 2000: - (G, H) = Interpolate(1995, IY, IDAY, IgrfCoefficients.G95, IgrfCoefficients.G00, IgrfCoefficients.H95, IgrfCoefficients.H00); + (G, H) = Interpolate(1995, year, doy, IgrfCoefficients.G95, IgrfCoefficients.G00, IgrfCoefficients.H95, IgrfCoefficients.H00); break; case < 2005: - (G, H) = Interpolate(2000, IY, IDAY, IgrfCoefficients.G00, IgrfCoefficients.G05, IgrfCoefficients.H00, IgrfCoefficients.H05); + (G, H) = Interpolate(2000, year, doy, IgrfCoefficients.G00, IgrfCoefficients.G05, IgrfCoefficients.H00, IgrfCoefficients.H05); break; case < 2010: - (G, H) = Interpolate(2005, IY, IDAY, IgrfCoefficients.G05, IgrfCoefficients.G10, IgrfCoefficients.H05, IgrfCoefficients.H10); + (G, H) = Interpolate(2005, year, doy, IgrfCoefficients.G05, IgrfCoefficients.G10, IgrfCoefficients.H05, IgrfCoefficients.H10); break; case < 2015: - (G, H) = Interpolate(2010, IY, IDAY, IgrfCoefficients.G10, IgrfCoefficients.G15, IgrfCoefficients.H10, IgrfCoefficients.H15); + (G, H) = Interpolate(2010, year, doy, IgrfCoefficients.G10, IgrfCoefficients.G15, IgrfCoefficients.H10, IgrfCoefficients.H15); break; case < 2020: - (G, H) = Interpolate(2015, IY, IDAY, IgrfCoefficients.G15, IgrfCoefficients.G20, IgrfCoefficients.H15, IgrfCoefficients.H20); + (G, H) = Interpolate(2015, year, doy, IgrfCoefficients.G15, IgrfCoefficients.G20, IgrfCoefficients.H15, IgrfCoefficients.H20); break; case < 2025: - (G, H) = Interpolate(2020, IY, IDAY, IgrfCoefficients.G20, IgrfCoefficients.G25, IgrfCoefficients.H20, IgrfCoefficients.H25); + (G, H) = Interpolate(2020, year, doy, IgrfCoefficients.G20, IgrfCoefficients.G25, IgrfCoefficients.H20, IgrfCoefficients.H25); break; case >= 2025: - (G, H) = Extrapolate(IY, IDAY); + (G, H) = Extrapolate(year, doy); break; } @@ -127,36 +127,39 @@ public ComputationContext Recalc(DateTime dateTime, CartesianVector? s double CTCL = CT0 * CL0; Sun sun = Sun(dateTime); - double S1 = Math.Cos(sun.Srasn) * Math.Cos(sun.Sdec); - double S2 = Math.Sin(sun.Srasn) * Math.Cos(sun.Sdec); - double S3 = Math.Sin(sun.Sdec); + (double sinSrasn, double cosSrasn) = Math.SinCos(sun.Srasn); + (double sinSdec, double cosSdec) = Math.SinCos(sun.Sdec); + double S1 = cosSrasn * cosSdec; + double S2 = sinSrasn * cosSdec; + double S3 = sinSdec; - double DJ = 365d * (IY - 1900) + (IY - 1901) / 4d + IDAY - 0.5d + (IHOUR * 3600 + MIN * 60 + ISEC) / 86400d; + double DJ = 365d * (year - 1900) + (year - 1901) / 4d + doy - 0.5d + (hour * 3600 + minutes * 60 + seconds) / 86400d; double T = DJ / 36525d; double OBLIQ = (23.45229d - 0.0130125d * T) / 57.2957795d; - double DZ1 = 0; - double DZ2 = -Math.Sin(OBLIQ); - double DZ3 = Math.Cos(OBLIQ); + (double sinOBLIQ, double cosOBLIQ) = Math.SinCos(OBLIQ); + double DZ2 = -sinOBLIQ; + double DZ3 = cosOBLIQ; double DY1 = DZ2 * S3 - DZ3 * S2; - double DY2 = DZ3 * S1 - DZ1 * S3; - double DY3 = DZ1 * S2 - DZ2 * S1; + double DY2 = DZ3 * S1; + double DY3 = -DZ2 * S1; - double V = Math.Sqrt( - Math.Pow(swVelocity.Required().X, 2) - + Math.Pow(swVelocity.Required().Y, 2) - + Math.Pow(swVelocity.Required().Z, 2)); + CartesianVector sw = swVelocity.Required(); + double vx = sw.X; + double vy = sw.Y; + double vz = sw.Z; + double V = Math.Sqrt(vx * vx + vy * vy + vz * vz); + double invV = 1.0d / V; - double DX1 = -swVelocity.Required().X / V; - double DX2 = -swVelocity.Required().Y / V; - double DX3 = -swVelocity.Required().Z / V; + double DX1 = -vx * invV; + double DX2 = -vy * invV; + double DX3 = -vz * invV; - double X1 = DX1 * S1 + DX2 * DY1 + DX3 * DZ1; + double X1 = DX1 * S1 + DX2 * DY1; double X2 = DX1 * S2 + DX2 * DY2 + DX3 * DZ2; double X3 = DX1 * S3 + DX2 * DY3 + DX3 * DZ3; - double CGST = Math.Cos(sun.Gst); - double SGST = Math.Sin(sun.Gst); + (double SGST, double CGST) = Math.SinCos(sun.Gst); double DIP1 = STCL * CGST - STSL * SGST; double DIP2 = STCL * SGST + STSL * CGST; double DIP3 = CT0; @@ -165,9 +168,10 @@ public ComputationContext Recalc(DateTime dateTime, CartesianVector? s double Y2 = DIP3 * X1 - DIP1 * X3; double Y3 = DIP1 * X2 - DIP2 * X1; double Y = Math.Sqrt(Y1 * Y1 + Y2 * Y2 + Y3 * Y3); - Y1 /= Y; - Y2 /= Y; - Y3 /= Y; + double invY = 1.0d / Y; + Y1 *= invY; + Y2 *= invY; + Y3 *= invY; double Z1 = X2 * Y3 - X3 * Y2; double Z2 = X3 * Y1 - X1 * Y3; @@ -179,31 +183,32 @@ public ComputationContext Recalc(DateTime dateTime, CartesianVector? s double E21 = DY1 * X1 + DY2 * X2 + DY3 * X3; double E22 = DY1 * Y1 + DY2 * Y2 + DY3 * Y3; double E23 = DY1 * Z1 + DY2 * Z2 + DY3 * Z3; - double E31 = DZ1 * X1 + DZ2 * X2 + DZ3 * X3; - double E32 = DZ1 * Y1 + DZ2 * Y2 + DZ3 * Y3; - double E33 = DZ1 * Z1 + DZ2 * Z2 + DZ3 * Z3; + double E31 = DZ2 * X2 + DZ3 * X3; + double E32 = DZ2 * Y2 + DZ3 * Y3; + double E33 = DZ2 * Z2 + DZ3 * Z3; - double SPS = DIP1 * X1 + DIP2 * X2 + DIP3 * X3; - double CPS = Math.Sqrt(1 - SPS * SPS); + double SPS = Math.FusedMultiplyAdd(DIP1, X1, Math.FusedMultiplyAdd(DIP2, X2, DIP3 * X3)); + double SPS2 = SPS * SPS; + double CPS = Math.Sqrt(1.0d - SPS2); double PSI = Math.Asin(SPS); - double A11 = X1 * CGST + X2 * SGST; - double A12 = -X1 * SGST + X2 * CGST; + double A11 = Math.FusedMultiplyAdd(X1, CGST, X2 * SGST); + double A12 = Math.FusedMultiplyAdd(-X1, SGST, X2 * CGST); double A13 = X3; - double A21 = Y1 * CGST + Y2 * SGST; - double A22 = -Y1 * SGST + Y2 * CGST; + double A21 = Math.FusedMultiplyAdd(Y1, CGST, Y2 * SGST); + double A22 = Math.FusedMultiplyAdd(-Y1, SGST, Y2 * CGST); double A23 = Y3; - double A31 = Z1 * CGST + Z2 * SGST; - double A32 = -Z1 * SGST + Z2 * CGST; + double A31 = Math.FusedMultiplyAdd(Z1, CGST, Z2 * SGST); + double A32 = Math.FusedMultiplyAdd(-Z1, SGST, Z2 * CGST); double A33 = Z3; - double EXMAGX = CT0 * (CL0 * CGST - SL0 * SGST); - double EXMAGY = CT0 * (CL0 * SGST + SL0 * CGST); + double EXMAGX = CT0 * (Math.FusedMultiplyAdd(CL0, CGST, -SL0 * SGST)); + double EXMAGY = CT0 * (Math.FusedMultiplyAdd(CL0, SGST, SL0 * CGST)); double EXMAGZ = -ST0; - double EYMAGX = -(SL0 * CGST + CL0 * SGST); - double EYMAGY = -(SL0 * SGST - CL0 * CGST); - double CFI = Y1 * EYMAGX + Y2 * EYMAGY; - double SFI = Y1 * EXMAGX + Y2 * EXMAGY + Y3 * EXMAGZ; + double EYMAGX = -(Math.FusedMultiplyAdd(SL0, CGST, CL0 * SGST)); + double EYMAGY = -(Math.FusedMultiplyAdd(SL0, SGST, -CL0 * CGST)); + double CFI = Math.FusedMultiplyAdd(Y1, EYMAGX, Y2 * EYMAGY); + double SFI = Math.FusedMultiplyAdd(Y1, EXMAGX, Math.FusedMultiplyAdd(Y2, EXMAGY, Y3 * EXMAGZ)); return new ComputationContext( ST0: ST0, CT0: CT0, SL0: SL0, CL0: CL0, @@ -216,19 +221,62 @@ public ComputationContext Recalc(DateTime dateTime, CartesianVector? s H: H, G: G, REC: REC); } - private static (double[], double[]) Interpolate( - int year1, int IY, int IDAY, - double[] G1, double[] G2, double[] H1, double[] H2) + private static (double[] G, double[] H) Interpolate( + int year1, int IY, int doy, + double[] G1, double[] G2, double[] H1, double[] H2) { - double[] G = new double[105]; - double[] H = new double[105]; + double[] G = new double[GeopackConstants.IgrfCoefficientCount]; + double[] H = new double[GeopackConstants.IgrfCoefficientCount]; - double F2 = (IY + (IDAY - 1) / 365.25d - year1) / 5d; + double F2 = (IY + (doy - 1) / GeopackConstants.DaysPerYear - year1) * GeopackConstants.IgrfInterpolationIntervalReciprocal; double F1 = 1.0d - F2; - for (int N = 0; N <= 104; N++) + + Vector vF1 = new(F1); + Vector vF2 = new(F2); + int vectorSize = Vector.Count; + + int i = 0; + + for (; i <= GeopackConstants.IgrfCoefficientCount - vectorSize; i += vectorSize) { - G[N] = G1[N] * F1 + G2[N] * F2; - H[N] = H1[N] * F1 + H2[N] * F2; + Vector vG1 = new(G1, i); + Vector vG2 = new(G2, i); + Vector vH1 = new(H1, i); + Vector vH2 = new(H2, i); + + (vG1 * vF1 + vG2 * vF2).CopyTo(G, i); + (vH1 * vF1 + vH2 * vF2).CopyTo(H, i); + } + + if (i >= GeopackConstants.IgrfCoefficientCount) + { + return (G, H); + } + + int remaining = GeopackConstants.IgrfCoefficientCount - i; + if (remaining >= vectorSize / 2) + { + Vector vG1 = new(G1, i); + Vector vG2 = new(G2, i); + Vector vH1 = new(H1, i); + Vector vH2 = new(H2, i); + + Vector vG = vG1 * vF1 + vG2 * vF2; + Vector vH = vH1 * vF1 + vH2 * vF2; + + for (int j = 0; j < remaining; j++) + { + G[i + j] = vG[j]; + H[i + j] = vH[j]; + } + } + else + { + for (; i < GeopackConstants.IgrfCoefficientCount; i++) + { + G[i] = G1[i] * F1 + G2[i] * F2; + H[i] = H1[i] * F1 + H2[i] * F2; + } } return (G, H); @@ -236,24 +284,51 @@ private static (double[], double[]) Interpolate( private static (double[], double[]) Extrapolate(int iy, int iday) { - double DT = iy + (iday - 1) / 365.25d - 2025; - double[] G = new double[105]; - double[] H = new double[105]; + double DT = iy + (iday - 1) / GeopackConstants.DaysPerYear - GeopackConstants.IgrfExtrapolationBaseYear; + double[] G = new double[GeopackConstants.IgrfCoefficientCount]; + double[] H = new double[GeopackConstants.IgrfCoefficientCount]; + + int vectorSize = Vector.Count; + Vector vDT = new(DT); + + int vectorizedLength = GeopackConstants.IgrfCoefficientCount; - for (int N = 0; N <= 104; N++) + for (int i = 0; i < vectorizedLength; i += vectorSize) { - G[N] = IgrfCoefficients.G25[N]; - H[N] = IgrfCoefficients.H25[N]; - if (N > 44) - { - continue; - } + Vector vG25 = new(IgrfCoefficients.G25, i); + Vector vH25 = new(IgrfCoefficients.H25, i); + vG25.CopyTo(G, i); + vH25.CopyTo(H, i); + } - G[N] += IgrfCoefficients.DG25[N] * DT; - H[N] += IgrfCoefficients.DH25[N] * DT; + for (int i = vectorizedLength; i < GeopackConstants.IgrfCoefficientCount; i++) + { + G[i] = IgrfCoefficients.G25[i]; + H[i] = IgrfCoefficients.H25[i]; + } + + int deltaVectorizedLength = (GeopackConstants.IgrfDeltaCoefficientCount / vectorSize) * vectorSize; + + for (int i = 0; i < deltaVectorizedLength; i += vectorSize) + { + Vector vG = new(G, i); + Vector vH = new(H, i); + Vector vDG25 = new(IgrfCoefficients.DG25, i); + Vector vDH25 = new(IgrfCoefficients.DH25, i); + + vG += vDG25 * vDT; + vG.CopyTo(G, i); + + vH += vDH25 * vDT; + vH.CopyTo(H, i); + } + + for (int i = deltaVectorizedLength; i < GeopackConstants.IgrfDeltaCoefficientCount; i++) + { + G[i] = Math.FusedMultiplyAdd(IgrfCoefficients.DG25[i], DT, G[i]); + H[i] = Math.FusedMultiplyAdd(IgrfCoefficients.DH25[i], DT, H[i]); } return (G, H); } - } diff --git a/src/Geopack/Geopack.ShuMgnp.cs b/src/Geopack/Geopack.ShuMgnp.cs index 67d3dfe..2d34b3a 100644 --- a/src/Geopack/Geopack.ShuMgnp.cs +++ b/src/Geopack/Geopack.ShuMgnp.cs @@ -1,11 +1,52 @@ using AuroraScienceHub.Geopack.Contracts.Cartesian; using AuroraScienceHub.Geopack.Contracts.Coordinates; +using AuroraScienceHub.Geopack.Contracts.Engine; using AuroraScienceHub.Geopack.Contracts.PhysicalObjects; namespace AuroraScienceHub.Geopack; internal sealed partial class Geopack { + /// + /// Shu magnetopause model coefficient 1 + /// + private const double ShuMgnpCoeff1 = 10.22D; + + /// + /// Shu magnetopause model coefficient 2 + /// + private const double ShuMgnpCoeff2 = 1.29D; + + /// + /// Shu magnetopause model coefficient 3 + /// + private const double ShuMgnpCoeff3 = 0.184D; + + /// + /// Shu magnetopause model coefficient 4 + /// + private const double ShuMgnpCoeff4 = 8.14D; + + /// + /// Shu magnetopause model pressure exponent (-1/6.6) + /// + private const double ShuMgnpPressureExponent = -0.15151515D; + + /// + /// Shu magnetopause model alpha coefficient 1 + /// + private const double ShuMgnpAlphaCoeff1 = 0.58D; + + /// + /// Shu magnetopause model alpha coefficient 2 + /// + private const double ShuMgnpAlphaCoeff2 = 0.007D; + + /// + /// Shu magnetopause model alpha coefficient 3 + /// + private const double ShuMgnpAlphaCoeff3 = 0.024D; + public Magnetopause ShuMgnp( double xnPd, double vel, @@ -31,7 +72,7 @@ public Magnetopause ShuMgnp( else { // Solar wind dynamic pressure in nPa - p = 1.94e-6 * xnPd * Math.Pow(vel, 2); + p = GeopackConstants.SolarWindDynamicPressureFactor * xnPd * vel * vel; } double phi; @@ -45,8 +86,8 @@ public Magnetopause ShuMgnp( } MagnetopausePosition id; - double r0 = (10.22D + 1.29D * Math.Tanh(0.184D * (bzImf + 8.14D))) * Math.Pow(p, -0.15151515D); - double alpha = (0.58D - 0.007D * bzImf) * (1.0D + 0.024D * Math.Log(p)); + double r0 = (ShuMgnpCoeff1 + ShuMgnpCoeff2 * Math.Tanh(ShuMgnpCoeff3 * (bzImf + ShuMgnpCoeff4))) * Math.Pow(p, ShuMgnpPressureExponent); + double alpha = (ShuMgnpAlphaCoeff1 - ShuMgnpAlphaCoeff2 * bzImf) * (1.0D + ShuMgnpAlphaCoeff3 * Math.Log(p)); double r = Math.Sqrt(location.X * location.X + location.Y * location.Y + location.Z * location.Z); double rm = r0 * Math.Pow(2.0D / (1.0D + location.X / r), alpha); @@ -71,6 +112,37 @@ public Magnetopause ShuMgnp( double ct = xmt96 / r; // Newton's iterative method to find nearest boundary point + (r, double sinT, double cosT) = FindMagnetopauseBoundaryPoint(r0, alpha, r, st, ct); + double xMgnp = r * cosT; + double rho = r * sinT; + (double sinPhi, double cosPhi) = Math.SinCos(phi); + double yMgnp = rho * sinPhi; + double zMgnp = rho * cosPhi; + + double dx = location.X - xMgnp; + double dy = location.Y - yMgnp; + double dz = location.Z - zMgnp; + double dist = Math.Sqrt(dx * dx + dy * dy + dz * dz); + + return new Magnetopause(CartesianLocation.New(xMgnp, yMgnp, zMgnp, CoordinateSystem.GSW), dist, id); + } + + /// + /// Finds the magnetopause boundary point using Newton's iterative method + /// + /// Base magnetopause standoff distance + /// Magnetopause shape parameter + /// Initial radial distance + /// Initial sine of theta + /// Initial cosine of theta + /// Tuple containing final r, sin(theta), and cos(theta) + private static (double r, double sinT, double cosT) FindMagnetopauseBoundaryPoint( + double r0, + double alpha, + double r, + double st, + double ct) + { int nit = 0; double t, rm_val, f, gradf_r, gradf_t, gradf, dr, dt, ds; @@ -89,8 +161,7 @@ public Magnetopause ShuMgnp( r += dr; t += dt; - st = Math.Sin(t); - ct = Math.Cos(t); + (st, ct) = Math.SinCos(t); ds = Math.Sqrt(dr * dr + (r * dt) * (r * dt)); nit++; @@ -103,16 +174,7 @@ public Magnetopause ShuMgnp( } while (ds > 1e-4); - double xMgnp = r * Math.Cos(t); - double rho = r * Math.Sin(t); - double yMgnp = rho * Math.Sin(phi); - double zMgnp = rho * Math.Cos(phi); - - double dist = Math.Sqrt( - Math.Pow(location.X - xMgnp, 2) + - Math.Pow(location.Y - yMgnp, 2) + - Math.Pow(location.Z - zMgnp, 2)); - - return new Magnetopause(CartesianLocation.New(xMgnp, yMgnp, zMgnp, CoordinateSystem.GSW), dist, id); + (double sinT, double cosT) = Math.SinCos(t); + return (r, sinT, cosT); } } diff --git a/src/Geopack/Geopack.SmGsw.cs b/src/Geopack/Geopack.SmGsw.cs index d16593e..61a84d4 100644 --- a/src/Geopack/Geopack.SmGsw.cs +++ b/src/Geopack/Geopack.SmGsw.cs @@ -18,17 +18,17 @@ private static T SmGswInternal(ComputationContext context, T components, Oper { OperationType.Direct => components.CoordinateSystem is CoordinateSystem.SM ? T.New( - components.X * context.CPS + components.Z * context.SPS, + Math.FusedMultiplyAdd(components.X, context.CPS, components.Z * context.SPS), components.Y, - components.Z * context.CPS - components.X * context.SPS, + Math.FusedMultiplyAdd(components.Z, context.CPS, -components.X * context.SPS), CoordinateSystem.GSW) : throw new InvalidOperationException("Input coordinates must be in SM system."), OperationType.Reversed => components.CoordinateSystem is CoordinateSystem.GSW ? T.New( - components.X * context.CPS - components.Z * context.SPS, + Math.FusedMultiplyAdd(components.X, context.CPS, -components.Z * context.SPS), components.Y, - components.X * context.SPS + components.Z * context.CPS, + Math.FusedMultiplyAdd(components.X, context.SPS, components.Z * context.CPS), CoordinateSystem.SM) : throw new InvalidOperationException("Input coordinates must be in GSW system."), _ => throw new NotSupportedException($"Specify correct OperationType: {operation}. Available types are Direct and Reversed.") diff --git a/src/Geopack/Geopack.Sun.cs b/src/Geopack/Geopack.Sun.cs index 8bc1f76..9a4ecb4 100644 --- a/src/Geopack/Geopack.Sun.cs +++ b/src/Geopack/Geopack.Sun.cs @@ -12,33 +12,34 @@ public Sun Sun(DateTime dateTime) return new Sun(dateTime); } - double fday = (dateTime.Hour * 3600 + dateTime.Minute * 60 + dateTime.Second) / 86400.0D; + double fday = (dateTime.Hour * GeopackConstants.SecondsPerHour + dateTime.Minute * 60 + dateTime.Second) / GeopackConstants.SecondsPerDay; //TODO Ask Tsyganenko if there 4 is really not 4.0D. Here we lose fraction double dj = 365 * (dateTime.Year - 1900) + (dateTime.Year - 1901) / 4 + dateTime.DayOfYear - 0.5D + fday; - double t = dj / 36525.0D; - double vl = (279.696678D + 0.9856473354D * dj) % 360.0D; - double gst = (279.690983D + 0.9856473354D * dj + 360.0D * fday + 180.0D) % 360.0D / GeopackConstants.Rad; - double g = (358.475845D + 0.985600267D * dj) % 360.0D / GeopackConstants.Rad; + double t = dj / GeopackConstants.DaysPerJulianCentury; + double vl = (279.696678D + 0.9856473354D * dj) % GeopackConstants.DegreesPerCircle; + double gst = (279.690983D + 0.9856473354D * dj + GeopackConstants.DegreesPerCircle * fday + GeopackConstants.DegreesPerSemicircle) % GeopackConstants.DegreesPerCircle / GeopackConstants.Rad; + double g = (358.475845D + 0.985600267D * dj) % GeopackConstants.DegreesPerCircle / GeopackConstants.Rad; double slong = (vl + (1.91946D - 0.004789D * t) * Math.Sin(g) + 0.020094D * Math.Sin(2.0D * g)) / GeopackConstants.Rad; - if (slong > 6.2831853D) + if (slong > GeopackConstants.TwoPi) { - slong -= 6.283185307D; + slong -= GeopackConstants.TwoPi; } if (slong < 0.0D) { - slong += 6.283185307D; + slong += GeopackConstants.TwoPi; } - double obliq = (23.45229D - 0.0130125D * t) / GeopackConstants.Rad; - double sob = Math.Sin(obliq); + double obliq = Math.FusedMultiplyAdd(-GeopackConstants.EarthObliquityRatePerCentury, t, GeopackConstants.EarthObliquityJ2000) / GeopackConstants.Rad; + (double sob, double cob) = Math.SinCos(obliq); double slp = slong - 9.924e-5D; - double sind = sob * Math.Sin(slp); - double cosd = Math.Sqrt(1.0D - sind * sind); + (double sinSlp, double cosSlp) = Math.SinCos(slp); + double sind = sob * sinSlp; + double cosd = Math.Sqrt(Math.FusedMultiplyAdd(-sind, sind, 1.0D)); double sc = sind / cosd; double sdec = Math.Atan2(sind, cosd); - double srasn = GeopackConstants.Pi - Math.Atan2(Math.Cos(obliq) / sob * sc, -Math.Cos(slp) / cosd); + double srasn = GeopackConstants.Pi - Math.Atan2(cob / sob * sc, -cosSlp / cosd); return new Sun(dateTime, gst, slong, srasn, sdec); } diff --git a/src/Geopack/Geopack.T96Mgnp.cs b/src/Geopack/Geopack.T96Mgnp.cs index fc9d80f..5a37693 100644 --- a/src/Geopack/Geopack.T96Mgnp.cs +++ b/src/Geopack/Geopack.T96Mgnp.cs @@ -1,5 +1,6 @@ using AuroraScienceHub.Geopack.Contracts.Cartesian; using AuroraScienceHub.Geopack.Contracts.Coordinates; +using AuroraScienceHub.Geopack.Contracts.Engine; using AuroraScienceHub.Geopack.Contracts.PhysicalObjects; namespace AuroraScienceHub.Geopack; @@ -22,7 +23,7 @@ public Magnetopause T96Mgnp( } else { - pd = 1.94e-6 * xnPd * vel * vel; + pd = GeopackConstants.SolarWindDynamicPressureFactor * xnPd * vel * vel; } if (pd is 0D) @@ -54,14 +55,15 @@ public Magnetopause T96Mgnp( phi = 0.0D; } - double rho = Math.Sqrt(Math.Pow(location.Y, 2) + Math.Pow(location.Z, 2)); + double rho = Math.Sqrt(location.Y * location.Y + location.Z * location.Z); if (location.X < xm) { double xMgnp = location.X; - double rhomGnp = a * Math.Sqrt(Math.Pow(s0, 2) - 1.0D); - double yMgnp = rhomGnp * Math.Sin(phi); - double zMgnp = rhomGnp * Math.Cos(phi); + double rhomGnp = a * Math.Sqrt(s0 * s0 - 1.0D); + (double sinPhi, double cosPhi) = Math.SinCos(phi); + double yMgnp = rhomGnp * sinPhi; + double zMgnp = rhomGnp * cosPhi; double dist = Math.Sqrt( (location.X - xMgnp) * (location.X - xMgnp) + (location.Y - yMgnp) * (location.Y - yMgnp) + @@ -91,8 +93,9 @@ public Magnetopause T96Mgnp( } double rhomGnpOut = a * Math.Sqrt(arg); - double yMgnpOut = rhomGnpOut * Math.Sin(phi); - double zMgnpOut = rhomGnpOut * Math.Cos(phi); + (double sinPhiOut, double cosPhiOut) = Math.SinCos(phi); + double yMgnpOut = rhomGnpOut * sinPhiOut; + double zMgnpOut = rhomGnpOut * cosPhiOut; double distOut = Math.Sqrt( (location.X - xMgnpOut) * (location.X - xMgnpOut) + diff --git a/src/Geopack/Geopack.Trace.cs b/src/Geopack/Geopack.Trace.cs index 646f1e0..f0768f6 100644 --- a/src/Geopack/Geopack.Trace.cs +++ b/src/Geopack/Geopack.Trace.cs @@ -23,7 +23,7 @@ public FieldLine Trace(ComputationContext context, throw new InvalidOperationException("Location must be in GSW system."); } - List points = new(); + List points = new(lMax + 1); double direction = (double)dir; int l = 0; @@ -39,7 +39,8 @@ public FieldLine Trace(ComputationContext context, FieldLineRhsVector initialRhs = Rhand(context, x, y, z, iopt, parmod, exName, inName, ds3); double ad = 0.01D; - if (x * initialRhs.R1 + y * initialRhs.R2 + z * initialRhs.R3 < 0.0D) + double dotProduct = x * initialRhs.R1 + y * initialRhs.R2 + z * initialRhs.R3; + if (dotProduct < 0.0D) { ad = -0.01D; } @@ -147,14 +148,16 @@ private FieldLineRhsVector Rhand(ComputationContext context, InternalFieldModel inName, double ds3) { - CartesianVector externalField = exName.Calculate(iopt, parmod, context.PSI, CartesianLocation.New(x, y, z, CoordinateSystem.GSW)); - CartesianVector internalField = inName(context, CartesianLocation.New(x, y, z, CoordinateSystem.GSW)); + CartesianLocation location = CartesianLocation.New(x, y, z, CoordinateSystem.GSW); + + CartesianVector externalField = exName.Calculate(iopt, parmod, context.PSI, location); + CartesianVector internalField = inName(context, location); double bx = externalField.X + internalField.X; double by = externalField.Y + internalField.Y; double bz = externalField.Z + internalField.Z; - double b = ds3 / Math.Sqrt(bx * bx + by * by + bz * bz); + double b = ds3 / Math.Sqrt(Math.FusedMultiplyAdd(bx, bx, Math.FusedMultiplyAdd(by, by, bz * bz))); double r1 = bx * b; double r2 = by * b; @@ -177,7 +180,10 @@ private FieldLineRhsVector Rhand(ComputationContext context, { ds3 = -currentDs / 3.0D; + // First Runge-Kutta stage FieldLineRhsVector r1 = Rhand(context, x, y, z, iopt, parmod, exName, inName, ds3); + + // Second stage FieldLineRhsVector r2 = Rhand(context, x + r1.R1, y + r1.R2, z + r1.R3, iopt, parmod, exName, inName, ds3); FieldLineRhsVector r3 = Rhand(context, x + 0.5D * (r1.R1 + r2.R1),