Skip to content

Commit 3edd31b

Browse files
oldrevoldrev
committed
merge
2 parents 69872bf + 41f1f19 commit 3edd31b

File tree

9 files changed

+283
-421
lines changed

9 files changed

+283
-421
lines changed

‎GeographicLib/Accumulator.cs

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -29,10 +29,6 @@ public struct Accumulator
2929
/// </summary>
3030
private double _s, _t;
3131

32-
public static Accumulator NaN => new Accumulator(double.NaN);
33-
34-
public static Accumulator Zero => new Accumulator(0D);
35-
3632
/// <summary>
3733
/// Construct from a double.
3834
/// </summary>

‎GeographicLib/Directory.Build.props

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,6 @@
55
<Title>This is an attempt to port the Java version of GeographicLib to C#.</Title>
66
<PackageTags>geographic;gis;geodesic;gnomonic;</PackageTags>
77
<VersionPrefix>1.49.3</VersionPrefix>
8-
<VersionSuffix>beta</VersionSuffix>
98
</PropertyGroup>
109

1110
<Import Project="../Shared.props" />

‎GeographicLib/GeoMath.cs

Lines changed: 0 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -187,22 +187,6 @@ public static double Polyval(int N, double[] p, int s, double x)
187187
return y;
188188
}
189189

190-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
191-
public static double Polyval(int N, in ReadOnlySpan<double> p, int s, double x)
192-
{
193-
double y = N < 0 ? 0 : p[s++];
194-
while (--N >= 0) y = y * x + p[s++];
195-
return y;
196-
}
197-
198-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
199-
public unsafe static double Polyval(int N, double* p, int s, double x)
200-
{
201-
double y = N < 0 ? 0 : p[s++];
202-
while (--N >= 0) y = y * x + p[s++];
203-
return y;
204-
}
205-
206190
[MethodImpl(MethodImplOptions.AggressiveInlining)]
207191
public static double AngRound(double x)
208192
{
@@ -283,7 +267,6 @@ public static Pair AngDiff(double x, double y)
283267
* The results obey exactly the elementary properties of the trigonometric
284268
* functions, e.g., sin 9&deg; = cos 81&deg; = &minus; sin 123456789&deg;.
285269
**********************************************************************/
286-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
287270
public static Pair Sincosd(double x)
288271
{
289272
// In order to minimize round-off errors, this function exactly reduces
@@ -320,7 +303,6 @@ public static Pair Sincosd(double x)
320303
* &minus;1) = &minus;180&deg;, for &epsilon; positive and tiny;
321304
* atan2d(&plusmn;0, 1) = &plusmn;0&deg;.
322305
**********************************************************************/
323-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
324306
public static double Atan2d(double y, double x)
325307
{
326308
// In order to minimize round-off errors, this function rearranges the

‎GeographicLib/Geodesic.cs

Lines changed: 200 additions & 234 deletions
Large diffs are not rendered by default.

‎GeographicLib/GeodesicData.cs

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ namespace GeographicLib
2020
* Geodesic#Inverse(double, double, double, double) Geodesic.Inverse} and it
2121
* always includes the field <i>a12</i>.
2222
**********************************************************************/
23-
public struct GeodesicData
23+
public sealed class GeodesicData
2424
{
2525
/// <summary>
2626
/// latitude of point 1 (degrees).

‎GeographicLib/GeodesicLine.cs

Lines changed: 66 additions & 137 deletions
Original file line numberDiff line numberDiff line change
@@ -94,7 +94,7 @@ namespace GeographicLib
9494
* }
9595
* }}</pre>
9696
**********************************************************************/
97-
public unsafe struct GeodesicLine
97+
public sealed class GeodesicLine
9898
{
9999

100100
private const int nC1_ = Geodesic.nC1_;
@@ -103,23 +103,15 @@ public unsafe struct GeodesicLine
103103
private const int nC3_ = Geodesic.nC3_;
104104
private const int nC4_ = Geodesic.nC4_;
105105

106-
private const int C1aSize = nC1_ + 1;
107-
private const int C1paSize = nC1p_ + 1;
108-
private const int C2aSize = nC2_ + 1;
109-
110-
private readonly double _lat1, _lon1, _azi1;
111-
private readonly double _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
106+
private double _lat1, _lon1, _azi1;
107+
private double _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
112108
_salp1, _calp1, _ssig1, _csig1, _dn1, _stau1, _ctau1, _somg1, _comg1,
113109
_A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41;
114110
private double _a13, _s13;
115111
// index zero elements of _C1a, _C1pa, _C2a, _C3a are unused
116-
117-
private fixed double _C1a[C1aSize];
118-
private fixed double _C1pa[C1paSize];
119-
private fixed double _C2a[C2aSize];
120-
private fixed double _C3a[nC3_];
121-
private fixed double _C4a[nC4_]; // all the elements of _C4a are used
122-
private readonly int _caps;
112+
private double[] _C1a, _C1pa, _C2a, _C3a,
113+
_C4a; // all the elements of _C4a are used
114+
private int _caps;
123115

124116
/**
125117
* Constructor for a geodesic line staring at latitude <i>lat1</i>, longitude
@@ -137,7 +129,7 @@ public unsafe struct GeodesicLine
137129
* fixed, writing <i>lat1</i> = &plusmn;(90&deg; &minus; &epsilon;), and
138130
* taking the limit &epsilon; &rarr; 0+.
139131
**********************************************************************/
140-
public GeodesicLine(in Geodesic g,
132+
public GeodesicLine(Geodesic g,
141133
double lat1, double lon1, double azi1) :
142134
this(g, lat1, lon1, azi1, GeodesicMask.ALL)
143135
{
@@ -189,31 +181,25 @@ public GeodesicLine(in Geodesic g,
189181
* <i>caps</i> |= {@link GeodesicMask#ALL} for all of the above.
190182
* </ul>
191183
**********************************************************************/
192-
public GeodesicLine(in Geodesic g,
184+
public GeodesicLine(Geodesic g,
193185
double lat1, double lon1, double azi1,
194-
int caps) : this(g, lat1, lon1, azi1, double.NaN, double.NaN, caps, true)
186+
int caps)
195187
{
188+
azi1 = GeoMath.AngNormalize(azi1);
189+
double salp1, calp1;
190+
// Guard against underflow in salp0
191+
{
192+
Pair p = GeoMath.Sincosd(GeoMath.AngRound(azi1));
193+
salp1 = p.First; calp1 = p.Second;
194+
}
195+
LineInit(g, lat1, lon1, azi1, salp1, calp1, caps);
196196
}
197197

198-
private GeodesicLine(in Geodesic g,
198+
private void LineInit(Geodesic g,
199199
double lat1, double lon1,
200200
double azi1, double salp1, double calp1,
201-
int caps, bool angNorm = false)
201+
int caps)
202202
{
203-
_a13 = double.NaN;
204-
_s13 = double.NaN;
205-
206-
if (angNorm)
207-
{
208-
azi1 = GeoMath.AngNormalize(azi1);
209-
// Guard against underflow in salp0
210-
{
211-
Pair p = GeoMath.Sincosd(GeoMath.AngRound(azi1));
212-
salp1 = p.First;
213-
calp1 = p.Second;
214-
}
215-
}
216-
217203
_a = g._a;
218204
_f = g._f;
219205
_b = g._b;
@@ -266,103 +252,67 @@ private GeodesicLine(in Geodesic g,
266252
if ((_caps & GeodesicMask.CAP_C1) != 0)
267253
{
268254
_A1m1 = Geodesic.A1m1f(eps);
269-
//_C1a = new double[nC1_ + 1];
270-
fixed (double* pC1a = _C1a)
271-
{
272-
var spanC1a = new Span<double>(pC1a, C1aSize);
273-
Geodesic.C1f(eps, spanC1a);
274-
_B11 = Geodesic.SinCosSeries(true, _ssig1, _csig1, spanC1a);
275-
}
255+
_C1a = new double[nC1_ + 1];
256+
Geodesic.C1f(eps, _C1a);
257+
_B11 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C1a);
276258
double s = Math.Sin(_B11), c = Math.Cos(_B11);
277259
// tau1 = sig1 + B11
278260
_stau1 = _ssig1 * c + _csig1 * s;
279261
_ctau1 = _csig1 * c - _ssig1 * s;
280262
// Not necessary because C1pa reverts C1a
281263
// _B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa, nC1p_);
282264
}
283-
else
284-
{
285-
_A1m1 = double.NaN;
286-
_B11 = double.NaN;
287-
_stau1 = double.NaN;
288-
_ctau1 = double.NaN;
289-
}
290265

291266
if ((_caps & GeodesicMask.CAP_C1p) != 0)
292267
{
293-
fixed (double* pC1pa = _C1pa)
294-
{
295-
Geodesic.C1pf(eps, new Span<double>(pC1pa, C1paSize));
296-
}
268+
_C1pa = new double[nC1p_ + 1];
269+
Geodesic.C1pf(eps, _C1pa);
297270
}
298271

299272
if ((_caps & GeodesicMask.CAP_C2) != 0)
300273
{
301-
fixed (double* pC2a = _C2a)
302-
{
303-
var spanC2a = new Span<double>(pC2a, C2aSize);
304-
_A2m1 = Geodesic.A2m1f(eps);
305-
Geodesic.C2f(eps, spanC2a);
306-
_B21 = Geodesic.SinCosSeries(true, _ssig1, _csig1, spanC2a);
307-
}
308-
}
309-
else
310-
{
311-
_A2m1 = double.NaN;
312-
_B21 = double.NaN;
274+
_C2a = new double[nC2_ + 1];
275+
_A2m1 = Geodesic.A2m1f(eps);
276+
Geodesic.C2f(eps, _C2a);
277+
_B21 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C2a);
313278
}
314279

315280
if ((_caps & GeodesicMask.CAP_C3) != 0)
316281
{
317-
fixed (double* pC3a = _C3a)
318-
{
319-
var spanC3a = new Span<double>(pC3a, nC3_);
320-
g.C3f(eps, spanC3a);
321-
_A3c = -_f * _salp0 * g.A3f(eps);
322-
_B31 = Geodesic.SinCosSeries(true, _ssig1, _csig1, spanC3a);
323-
}
324-
}
325-
else
326-
{
327-
_A3c = double.NaN;
328-
_B31 = double.NaN;
282+
_C3a = new double[nC3_];
283+
g.C3f(eps, _C3a);
284+
_A3c = -_f * _salp0 * g.A3f(eps);
285+
_B31 = Geodesic.SinCosSeries(true, _ssig1, _csig1, _C3a);
329286
}
330287

331288
if ((_caps & GeodesicMask.CAP_C4) != 0)
332289
{
333-
fixed (double* pC4a = _C4a)
334-
{
335-
var spanC4a = new Span<double>(pC4a, nC4_);
336-
g.C4f(eps, spanC4a);
337-
// Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
338-
_A4 = GeoMath.Sq(_a) * _calp0 * _salp0 * g._e2;
339-
_B41 = Geodesic.SinCosSeries(false, _ssig1, _csig1, spanC4a);
340-
}
341-
}
342-
else
343-
{
344-
_A4 = double.NaN;
345-
_B41 = double.NaN;
290+
_C4a = new double[nC4_];
291+
g.C4f(eps, _C4a);
292+
// Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
293+
_A4 = GeoMath.Sq(_a) * _calp0 * _salp0 * g._e2;
294+
_B41 = Geodesic.SinCosSeries(false, _ssig1, _csig1, _C4a);
346295
}
347296
}
348297

349-
internal GeodesicLine(in Geodesic g,
298+
internal GeodesicLine(Geodesic g,
350299
double lat1, double lon1,
351300
double azi1, double salp1, double calp1,
352-
int caps, bool arcmode, double s13_a13) : this(g, lat1, lon1, azi1, salp1, calp1, caps)
301+
int caps, bool arcmode, double s13_a13)
353302
{
303+
LineInit(g, lat1, lon1, azi1, salp1, calp1, caps);
354304
GenSetDistance(arcmode, s13_a13);
355305
}
356306

357-
/*
307+
/**
358308
* A default constructor. If GeodesicLine.Position is called on the
359309
* resulting object, it returns immediately (without doing any calculations).
360310
* The object can be set with a call to {@link Geodesic.Line}. Use {@link
361311
* Init()} to test whether object is still in this uninitialized state.
362312
* (This constructor was useful in C++, e.g., to allow vectors of
363313
* GeodesicLine objects. It may not be needed in Java, so make it private.)
364314
**********************************************************************/
365-
315+
private GeodesicLine() { _caps = 0; }
366316

367317
/**
368318
* Compute the position of point 2 which is a distance <i>s12</i> (meters)
@@ -532,10 +482,10 @@ public GeodesicData Position(bool arcmode, double s12_a12, int outmask)
532482
s = Math.Sin(tau12),
533483
c = Math.Cos(tau12);
534484
// tau2 = tau1 + tau12
535-
fixed (double* pC1pa = _C1pa)
536-
{
537-
B12 = -Geodesic.SinCosSeries(true, _stau1 * c + _ctau1 * s, _ctau1 * c - _stau1 * s, new Span<double>(pC1pa, C1paSize));
538-
}
485+
B12 = -Geodesic.SinCosSeries(true,
486+
_stau1 * c + _ctau1 * s,
487+
_ctau1 * c - _stau1 * s,
488+
_C1pa);
539489
sig12 = tau12 - (B12 - _B11);
540490
ssig12 = Math.Sin(sig12); csig12 = Math.Cos(sig12);
541491
if (Math.Abs(_f) > 0.01)
@@ -565,10 +515,7 @@ public GeodesicData Position(bool arcmode, double s12_a12, int outmask)
565515
double
566516
ssig2_ = _ssig1 * csig12 + _csig1 * ssig12,
567517
csig2_ = _csig1 * csig12 - _ssig1 * ssig12;
568-
fixed (double* pC1a = _C1a)
569-
{
570-
B12 = Geodesic.SinCosSeries(true, ssig2_, csig2_, new Span<double>(pC1a, C1aSize));
571-
}
518+
B12 = Geodesic.SinCosSeries(true, ssig2_, csig2_, _C1a);
572519
double serr = (1 + _A1m1) * (sig12 + (B12 - _B11)) - s12_a12 / _b;
573520
sig12 = sig12 - serr / Math.Sqrt(1 + _k2 * GeoMath.Sq(ssig2_));
574521
ssig12 = Math.Sin(sig12); csig12 = Math.Cos(sig12);
@@ -586,12 +533,7 @@ public GeodesicData Position(bool arcmode, double s12_a12, int outmask)
586533
GeodesicMask.GEODESICSCALE)) != 0)
587534
{
588535
if (arcmode || Math.Abs(_f) > 0.01)
589-
{
590-
fixed (double* pC1a = _C1a)
591-
{
592-
B12 = Geodesic.SinCosSeries(true, ssig2, csig2, new Span<double>(pC1a, C1aSize));
593-
}
594-
}
536+
B12 = Geodesic.SinCosSeries(true, ssig2, csig2, _C1a);
595537
AB1 = (1 + _A1m1) * (B12 - _B11);
596538
}
597539
// sin(bet2) = cos(alp0) * sin(sig2)
@@ -619,14 +561,9 @@ public GeodesicData Position(bool arcmode, double s12_a12, int outmask)
619561
+ (Math.Atan2(E * somg2, comg2) - Math.Atan2(E * _somg1, _comg1)))
620562
: Math.Atan2(somg2 * _comg1 - comg2 * _somg1,
621563
comg2 * _comg1 + somg2 * _somg1);
622-
623-
double lam12 = double.NaN;
624-
fixed (double* pC3a = _C3a)
625-
{
626-
var spanC3a = new Span<double>(pC3a, nC3_);
627-
lam12 = omg12 + _A3c *
628-
(sig12 + (Geodesic.SinCosSeries(true, ssig2, csig2, spanC3a) - _B31));
629-
}
564+
double lam12 = omg12 + _A3c *
565+
(sig12 + (Geodesic.SinCosSeries(true, ssig2, csig2, _C3a)
566+
- _B31));
630567
double lon12 = GeoMath.ToDegrees(lam12);
631568
r.lon2 = ((outmask & GeodesicMask.LONG_UNROLL) != 0) ? _lon1 + lon12 :
632569
GeoMath.AngNormalize(r.lon1 + GeoMath.AngNormalize(lon12));
@@ -641,35 +578,27 @@ public GeodesicData Position(bool arcmode, double s12_a12, int outmask)
641578
if ((outmask &
642579
(GeodesicMask.REDUCEDLENGTH | GeodesicMask.GEODESICSCALE)) != 0)
643580
{
644-
fixed (double* pC2a = _C2a)
581+
double
582+
B22 = Geodesic.SinCosSeries(true, ssig2, csig2, _C2a),
583+
AB2 = (1 + _A2m1) * (B22 - _B21),
584+
J12 = (_A1m1 - _A2m1) * sig12 + (AB1 - AB2);
585+
if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
586+
// Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
587+
// accurate cancellation in the case of coincident points.
588+
r.m12 = _b * ((dn2 * (_csig1 * ssig2) - _dn1 * (_ssig1 * csig2))
589+
- _csig1 * csig2 * J12);
590+
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
645591
{
646-
var spanC2a = new Span<double>(pC2a, C2aSize);
647-
double
648-
B22 = Geodesic.SinCosSeries(true, ssig2, csig2, spanC2a),
649-
AB2 = (1 + _A2m1) * (B22 - _B21),
650-
J12 = (_A1m1 - _A2m1) * sig12 + (AB1 - AB2);
651-
if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
652-
// Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
653-
// accurate cancellation in the case of coincident points.
654-
r.m12 = _b * ((dn2 * (_csig1 * ssig2) - _dn1 * (_ssig1 * csig2))
655-
- _csig1 * csig2 * J12);
656-
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
657-
{
658-
double t = _k2 * (ssig2 - _ssig1) * (ssig2 + _ssig1) / (_dn1 + dn2);
659-
r.M12 = csig12 + (t * ssig2 - csig2 * J12) * _ssig1 / _dn1;
660-
r.M21 = csig12 - (t * _ssig1 - _csig1 * J12) * ssig2 / dn2;
661-
}
592+
double t = _k2 * (ssig2 - _ssig1) * (ssig2 + _ssig1) / (_dn1 + dn2);
593+
r.M12 = csig12 + (t * ssig2 - csig2 * J12) * _ssig1 / _dn1;
594+
r.M21 = csig12 - (t * _ssig1 - _csig1 * J12) * ssig2 / dn2;
662595
}
663596
}
664597

665598
if ((outmask & GeodesicMask.AREA) != 0)
666599
{
667-
var B42 = double.NaN;
668-
fixed (double* pC4a = _C4a)
669-
{
670-
var spanC4a = new Span<double>(pC4a, nC4_);
671-
B42 = Geodesic.SinCosSeries(false, ssig2, csig2, spanC4a);
672-
}
600+
double
601+
B42 = Geodesic.SinCosSeries(false, ssig2, csig2, _C4a);
673602
double salp12, calp12;
674603
if (_calp0 == 0 || _salp0 == 0)
675604
{

0 commit comments

Comments
 (0)