Skip to content

Commit 2af78f2

Browse files
oldrevoldrev
committed
optimization code, update to 1.49.2
1 parent 43cbe25 commit 2af78f2

File tree

5 files changed

+74
-64
lines changed

5 files changed

+74
-64
lines changed

‎GeographicLib/Directory.Build.props

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
<PackageId>Sandwych.GeographicLib</PackageId>
55
<Title>This is an attempt to port the Java version of GeographicLib to C#.</Title>
66
<PackageTags>geographic;gis;geodesic;gnomonic;</PackageTags>
7-
<VersionPrefix>1.49.1</VersionPrefix>
7+
<VersionPrefix>1.49.2</VersionPrefix>
88
</PropertyGroup>
99

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

‎GeographicLib/Geodesic.cs

Lines changed: 52 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -210,21 +210,21 @@ public sealed class Geodesic
210210
/**
211211
* The order of the expansions used by Geodesic.
212212
**********************************************************************/
213-
internal static readonly int GEOGRAPHICLIB_GEODESIC_ORDER = 6;
214-
215-
internal static readonly int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
216-
internal static readonly int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
217-
internal static readonly int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
218-
internal static readonly int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
219-
internal static readonly int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
220-
internal static readonly int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
221-
internal static readonly int nA3x_ = nA3_;
222-
internal static readonly int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
223-
internal static readonly int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
224-
internal static readonly int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
225-
internal static readonly int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
226-
private static readonly int maxit1_ = 20;
227-
private static readonly int maxit2_ = maxit1_ + GeoMath.Digits + 10;
213+
internal const int GEOGRAPHICLIB_GEODESIC_ORDER = 6;
214+
215+
internal const int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
216+
internal const int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
217+
internal const int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
218+
internal const int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
219+
internal const int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
220+
internal const int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
221+
internal const int nA3x_ = nA3_;
222+
internal const int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
223+
internal const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
224+
internal const int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
225+
internal const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
226+
private const int maxit1_ = 20;
227+
private const int maxit2_ = maxit1_ + GeoMath.Digits + 10;
228228

229229
// Underflow guard. We require
230230
// tiny_ * epsilon() > 0
@@ -638,22 +638,18 @@ public GeodesicData Inverse(double lat1, double lon1,
638638

639639
internal struct InverseData
640640
{
641-
internal GeodesicData g;
642-
internal double salp1, calp1, salp2, calp2;
641+
private GeodesicData _g;
643642

644-
[MethodImpl(MethodImplOptions.AggressiveInlining)]
645-
private void Init()
646-
{
647-
g = new GeodesicData();
648-
salp1 = calp1 = salp2 = calp2 = Double.NaN;
649-
}
643+
internal GeodesicData Geodesic => _g;
644+
internal double salp1, calp1, salp2, calp2;
650645

651646
internal static InverseData NaN
652647
{
653648
get
654649
{
655650
var self = new InverseData();
656-
self.Init();
651+
self._g = GeodesicData.NaN;
652+
self.salp1 = self.calp1 = self.salp2 = self.calp2 = Double.NaN;
657653
return self;
658654
}
659655
}
@@ -663,11 +659,11 @@ private InverseData InverseInt(double lat1, double lon1,
663659
double lat2, double lon2, int outmask)
664660
{
665661
InverseData result = InverseData.NaN;
666-
GeodesicData r = result.g;
667662
// Compute longitude difference (AngDiff does this carefully). Result is
668663
// in [-180, 180] but -180 is only for west-going geodesics. 180 is for
669664
// east-going and meridional geodesics.
670-
r.lat1 = lat1 = GeoMath.LatFix(lat1); r.lat2 = lat2 = GeoMath.LatFix(lat2);
665+
result.Geodesic.lat1 = lat1 = GeoMath.LatFix(lat1);
666+
result.Geodesic.lat2 = lat2 = GeoMath.LatFix(lat2);
671667
// If really close to the equator, treat as on equator.
672668
lat1 = GeoMath.AngRound(lat1);
673669
lat2 = GeoMath.AngRound(lat2);
@@ -678,11 +674,13 @@ private InverseData InverseInt(double lat1, double lon1,
678674
}
679675
if ((outmask & GeodesicMask.LONG_UNROLL) != 0)
680676
{
681-
r.lon1 = lon1; r.lon2 = (lon1 + lon12) + lon12s;
677+
result.Geodesic.lon1 = lon1;
678+
result.Geodesic.lon2 = (lon1 + lon12) + lon12s;
682679
}
683680
else
684681
{
685-
r.lon1 = GeoMath.AngNormalize(lon1); r.lon2 = GeoMath.AngNormalize(lon2);
682+
result.Geodesic.lon1 = GeoMath.AngNormalize(lon1);
683+
result.Geodesic.lon2 = GeoMath.AngNormalize(lon2);
686684
}
687685
// Make longitude difference positive.
688686
int lonsign = lon12 >= 0 ? 1 : -1;
@@ -798,7 +796,8 @@ private InverseData InverseInt(double lat1, double lon1,
798796
s12x = v.s12b; m12x = v.m12b;
799797
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
800798
{
801-
r.M12 = v.M12; r.M21 = v.M21;
799+
result.Geodesic.M12 = v.M12;
800+
result.Geodesic.M21 = v.M21;
802801
}
803802
}
804803
// Add the check for sig12 since zero length geodesics might yield m12 <
@@ -835,7 +834,7 @@ private InverseData InverseInt(double lat1, double lon1,
835834
sig12 = omg12 = lam12 / _f1;
836835
m12x = _b * Math.Sin(sig12);
837836
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
838-
r.M12 = r.M21 = Math.Cos(sig12);
837+
result.Geodesic.M12 = result.Geodesic.M21 = Math.Cos(sig12);
839838
a12 = lon12 / _f1;
840839

841840
}
@@ -864,7 +863,7 @@ private InverseData InverseInt(double lat1, double lon1,
864863
s12x = sig12 * _b * dnm;
865864
m12x = GeoMath.Sq(dnm) * _b * Math.Sin(sig12 / dnm);
866865
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
867-
r.M12 = r.M21 = Math.Cos(sig12 / dnm);
866+
result.Geodesic.M12 = result.Geodesic.M21 = Math.Cos(sig12 / dnm);
868867
a12 = GeoMath.ToDegrees(sig12);
869868
omg12 = lam12 / (_f1 * dnm);
870869
}
@@ -966,7 +965,8 @@ private InverseData InverseInt(double lat1, double lon1,
966965
s12x = v.s12b; m12x = v.m12b;
967966
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
968967
{
969-
r.M12 = v.M12; r.M21 = v.M21;
968+
result.Geodesic.M12 = v.M12;
969+
result.Geodesic.M21 = v.M21;
970970
}
971971
}
972972
m12x *= _b;
@@ -983,10 +983,10 @@ private InverseData InverseInt(double lat1, double lon1,
983983
}
984984

985985
if ((outmask & GeodesicMask.DISTANCE) != 0)
986-
r.s12 = 0 + s12x; // Convert -0 to 0
986+
result.Geodesic.s12 = 0 + s12x; // Convert -0 to 0
987987

988988
if ((outmask & GeodesicMask.REDUCEDLENGTH) != 0)
989-
r.m12 = 0 + m12x; // Convert -0 to 0
989+
result.Geodesic.m12 = 0 + m12x; // Convert -0 to 0
990990

991991
if ((outmask & GeodesicMask.AREA) != 0)
992992
{
@@ -1015,14 +1015,13 @@ private InverseData InverseInt(double lat1, double lon1,
10151015
}
10161016
double[] C4a = new double[nC4_];
10171017
C4f(eps, C4a);
1018-
double
1019-
B41 = SinCosSeries(false, ssig1, csig1, C4a),
1018+
double B41 = SinCosSeries(false, ssig1, csig1, C4a),
10201019
B42 = SinCosSeries(false, ssig2, csig2, C4a);
1021-
r.S12 = A4 * (B42 - B41);
1020+
result.Geodesic.S12 = A4 * (B42 - B41);
10221021
}
10231022
else
10241023
// Avoid problems with indeterminate sig1, sig2 on equator
1025-
r.S12 = 0;
1024+
result.Geodesic.S12 = 0;
10261025

10271026
if (!meridian && somg12 > 1)
10281027
{
@@ -1058,10 +1057,10 @@ private InverseData InverseInt(double lat1, double lon1,
10581057
}
10591058
alp12 = Math.Atan2(salp12, calp12);
10601059
}
1061-
r.S12 += _c2 * alp12;
1062-
r.S12 *= swapp * lonsign * latsign;
1060+
result.Geodesic.S12 += _c2 * alp12;
1061+
result.Geodesic.S12 *= swapp * lonsign * latsign;
10631062
// Convert -0 to 0
1064-
r.S12 += 0;
1063+
result.Geodesic.S12 += 0;
10651064
}
10661065

10671066
// Convert calp, salp to azimuth accounting for lonsign, swapp, latsign.
@@ -1070,16 +1069,22 @@ private InverseData InverseInt(double lat1, double lon1,
10701069
{ double t = salp1; salp1 = salp2; salp2 = t; }
10711070
{ double t = calp1; calp1 = calp2; calp2 = t; }
10721071
if ((outmask & GeodesicMask.GEODESICSCALE) != 0)
1073-
{ double t = r.M12; r.M12 = r.M21; r.M21 = t; }
1072+
{
1073+
double t = result.Geodesic.M12;
1074+
result.Geodesic.M12 = result.Geodesic.M21;
1075+
result.Geodesic.M21 = t;
1076+
}
10741077
}
10751078

10761079
salp1 *= swapp * lonsign; calp1 *= swapp * latsign;
10771080
salp2 *= swapp * lonsign; calp2 *= swapp * latsign;
10781081

10791082
// Returned value in [0, 180]
1080-
r.a12 = a12;
1081-
result.salp1 = salp1; result.calp1 = calp1;
1082-
result.salp2 = salp2; result.calp2 = calp2;
1083+
result.Geodesic.a12 = a12;
1084+
result.salp1 = salp1;
1085+
result.calp1 = calp1;
1086+
result.salp2 = salp2;
1087+
result.calp2 = calp2;
10831088
return result;
10841089
}
10851090

@@ -1130,7 +1135,7 @@ public GeodesicData Inverse(double lat1, double lon1,
11301135
{
11311136
outmask &= GeodesicMask.OUT_MASK;
11321137
InverseData result = InverseInt(lat1, lon1, lat2, lon2, outmask);
1133-
GeodesicData r = result.g;
1138+
GeodesicData r = result.Geodesic;
11341139
if ((outmask & GeodesicMask.AZIMUTH) != 0)
11351140
{
11361141
r.azi1 = GeoMath.Atan2d(result.salp1, result.calp1);
@@ -1186,7 +1191,7 @@ public GeodesicLine InverseLine(double lat1, double lon1,
11861191
{
11871192
InverseData result = InverseInt(lat1, lon1, lat2, lon2, 0);
11881193
double salp1 = result.salp1, calp1 = result.calp1,
1189-
azi1 = GeoMath.Atan2d(salp1, calp1), a12 = result.g.a12;
1194+
azi1 = GeoMath.Atan2d(salp1, calp1), a12 = result.Geodesic.a12;
11901195
// Ensure that a12 can be converted to a distance
11911196
if ((caps & (GeodesicMask.OUT_MASK & GeodesicMask.DISTANCE_IN)) != 0)
11921197
caps |= GeodesicMask.DISTANCE;

‎GeographicLib/GeodesicData.cs

Lines changed: 9 additions & 4 deletions
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 sealed class GeodesicData
23+
public class GeodesicData
2424
{
2525
/// <summary>
2626
/// latitude of point 1 (degrees).
@@ -85,10 +85,15 @@ public sealed class GeodesicData
8585
/// <summary>
8686
/// Initialize all the fields to Double.NaN.
8787
/// </summary>
88-
public GeodesicData()
88+
public static GeodesicData NaN
8989
{
90-
lat1 = lon1 = azi1 = lat2 = lon2 = azi2 =
91-
s12 = a12 = m12 = M12 = M21 = S12 = Double.NaN;
90+
get
91+
{
92+
var g = new GeodesicData();
93+
g.lat1 = g.lon1 = g.azi1 = g.lat2 = g.lon2 = g.azi2 = g.s12
94+
= g.a12 = g.m12 = g.M12 = g.M21 = g.S12 = Double.NaN;
95+
return g;
96+
}
9297
}
9398
}
9499
}

‎GeographicLib/GeodesicLine.cs

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -97,11 +97,11 @@ namespace GeographicLib
9797
public sealed class GeodesicLine
9898
{
9999

100-
private static readonly int nC1_ = Geodesic.nC1_;
101-
private static readonly int nC1p_ = Geodesic.nC1p_;
102-
private static readonly int nC2_ = Geodesic.nC2_;
103-
private static readonly int nC3_ = Geodesic.nC3_;
104-
private static readonly int nC4_ = Geodesic.nC4_;
100+
private const int nC1_ = Geodesic.nC1_;
101+
private const int nC1p_ = Geodesic.nC1p_;
102+
private const int nC2_ = Geodesic.nC2_;
103+
private const int nC3_ = Geodesic.nC3_;
104+
private const int nC4_ = Geodesic.nC4_;
105105

106106
private double _lat1, _lon1, _azi1;
107107
private double _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
@@ -447,17 +447,17 @@ public GeodesicData ArcPosition(double a12, int outmask)
447447
* Requesting a value which the GeodesicLine object is not capable of
448448
* computing is not an error; Double.NaN is returned instead.
449449
**********************************************************************/
450-
public GeodesicData Position(bool arcmode, double s12_a12,
451-
int outmask)
450+
public GeodesicData Position(bool arcmode, double s12_a12, int outmask)
452451
{
453452
outmask &= _caps & GeodesicMask.OUT_MASK;
454-
GeodesicData r = new GeodesicData();
453+
GeodesicData r = GeodesicData.NaN;
455454
if (!(Init && (arcmode || (_caps & (GeodesicMask.OUT_MASK & GeodesicMask.DISTANCE_IN)) != 0)))
456455
{
457456
// Uninitialized or impossible distance calculation requested
458457
return r;
459458
}
460-
r.lat1 = _lat1; r.azi1 = _azi1;
459+
r.lat1 = _lat1;
460+
r.azi1 = _azi1;
461461
r.lon1 = ((outmask & GeodesicMask.LONG_UNROLL) != 0) ? _lon1 :
462462
GeoMath.AngNormalize(_lon1);
463463

‎GeographicLib/Gnomonic.cs

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -241,11 +241,11 @@ public GnomonicData Reverse(double lat0, double lon0, double x, double y)
241241
| GeodesicMask.GEODESICSCALE);
242242

243243
int count = numit_, trip = 0;
244-
GeodesicData pos = null;
244+
GeodesicData pos = GeodesicData.NaN;
245245

246246
while (count-- > 0)
247247
{
248-
pos = line.Position(s,
248+
pos = line.Position(s,
249249
GeodesicMask.LONGITUDE | GeodesicMask.LATITUDE
250250
| GeodesicMask.AZIMUTH | GeodesicMask.DISTANCE_IN
251251
| GeodesicMask.REDUCEDLENGTH
@@ -256,7 +256,7 @@ public GnomonicData Reverse(double lat0, double lon0, double x, double y)
256256
break;
257257
}
258258

259-
double ds = little ? ((pos.m12 / pos.M12) - rho) * pos.M12 * pos.M12
259+
double ds = little ? ((pos.m12 / pos.M12) - rho) * pos.M12 * pos.M12
260260
: (rho - (pos.M12 / pos.m12)) * pos.m12 * pos.m12;
261261
s -= ds;
262262

0 commit comments

Comments
 (0)