Skip to content

Commit

Permalink
Fix GMST calculations.
Browse files Browse the repository at this point in the history
  • Loading branch information
xclud committed Feb 9, 2024
1 parent 0457d76 commit e670d8f
Show file tree
Hide file tree
Showing 8 changed files with 243 additions and 32 deletions.
4 changes: 4 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
## [2.0.3]

* Fix GMST calculations.

## [2.0.2]

* Adds ground track (Footprint) calculations.
Expand Down
1 change: 1 addition & 0 deletions lib/latlng.dart
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ part 'src/look_angle.dart';
part 'src/julian.dart';
part 'src/angle.dart';
part 'src/const.dart';
part 'src/point2d.dart';

part 'src/geojson/feature.dart';
part 'src/geojson/feature_collection.dart';
Expand Down
2 changes: 1 addition & 1 deletion lib/src/const.dart
Original file line number Diff line number Diff line change
Expand Up @@ -8,4 +8,4 @@ const double _twopi = 2.0 * pi;
const double _hoursPerDay = 24.0; // Hours per day (solar)
const double _minutesPerDay = 1440.0; // Minutes per day (solar)
const double _secondsPerDay = 86400.0; // Seconds per day (solar)
const double _omegaE = 1.00273790934; // Earth rotation per sidereal day
//const double _omegaE = 1.00273790934; // Earth rotation per sidereal day
12 changes: 6 additions & 6 deletions lib/src/eci.dart
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@ class EarthCenteredInertial {
final double z;

/// Converts this ECI object to ECF.
EarthCenteredEarthFixed toEcf(double gmst) {
double x = this.x * cos(gmst) + this.y * sin(gmst);
double y = this.x * (-sin(gmst)) + this.y * cos(gmst);
EarthCenteredEarthFixed toEcf(Angle gmst) {
final x = this.x * cos(gmst.radians) + this.y * sin(gmst.radians);
final y = this.x * (-sin(gmst.radians)) + this.y * cos(gmst.radians);

return EarthCenteredEarthFixed(x, y, z);
}
Expand All @@ -28,15 +28,15 @@ class EarthCenteredInertial {
return toEcf(utc.toUtc().gsmt);
}

LatLngAlt toGeodetic(Planet planet, double gmst) {
LatLngAlt toGeodetic(Planet planet, Angle gmst) {
// http://www.celestrak.com/columns/v02n03/
final a = planet.radius;
final f = planet.flattening;
final e2 = (2 * f) - (f * f);

final R = sqrt((x * x) + (y * y));

var longitude = atan2(y, x) - gmst;
var longitude = atan2(y, x) - gmst.radians;
while (longitude < -pi) {
longitude += pi * 2;
}
Expand Down Expand Up @@ -64,7 +64,7 @@ class EarthCenteredInertial {
}

LatLngAlt toGeodeticByDateTime(Planet planet, DateTime utc) {
final gmst = Julian.fromDateTime(utc).toGmst();
final gmst = Julian.fromDateTime(utc).gmst;
return toGeodetic(planet, gmst);
}
}
38 changes: 14 additions & 24 deletions lib/src/julian.dart
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,8 @@ extension DateTimeExtensions on DateTime {
Julian get julian => Julian.fromDateTime(this);

/// Calculate Greenwich Mean Sidereal Time according to http://aa.usno.navy.mil/faq/docs/GAST.php
double get gsmt {
final result = julian.toGmst();
Angle get gsmt {
final result = julian.gmst;

return result;
}
Expand Down Expand Up @@ -141,26 +141,14 @@ class Julian {
/// Returns the angle, in radians, measuring eastward from the Vernal Equinox to
/// the prime meridian. This angle is also referred to as "ThetaG"
/// (Theta GMST).
double toGmst() {
// References:
// The 1992 Astronomical Almanac, page B6.
// Explanatory Supplement to the Astronomical Almanac, page 50.
// Orbital Coordinate Systems, Part III, Dr. T.S. Kelso,
// Satellite Times, Nov/Dec 1995

final ut = (value + 0.5) % 1.0;
final tu = (fromJan1_12h_2000() - ut) / 36525.0;

double gmst = 24110.54841 +
(tu * (8640184.812866 + (tu * (0.093104 - (tu * 6.2e-06)))));

gmst = (gmst + (_secondsPerDay * _omegaE * ut)) % _secondsPerDay;

if (gmst < 0.0) {
gmst += _secondsPerDay; // "wrap" negative modulo value
}

return _twopi * (gmst / _secondsPerDay);
Angle get gmst {
/* Calculate Greenwich Mean Sidereal Time according to
http://aa.usno.navy.mil/faq/docs/GAST.php */
final d = value - 2451545.0;
// Low precision equation is good enough for our purposes.
final rad = (18.697374558 + 24.06570982441908 * d) % 24;

return Angle.degree(rad);
}

/// Calculate Local Mean Sidereal Time for this Julian date at the given
Expand All @@ -170,8 +158,10 @@ class Julian {
///
/// The angle, in radians, measuring eastward from the Vernal Equinox to
/// the given longitude.
double toLmst(double longitude) {
return (toGmst() + longitude) % _twopi;
Angle lmst(Angle longitude) {
final rad = (gmst.radians + longitude.radians) % _twopi;

return Angle.radian(rad);
}

/// Returns a UTC DateTime object that corresponds to this Julian date.
Expand Down
209 changes: 209 additions & 0 deletions lib/src/point2d.dart
Original file line number Diff line number Diff line change
@@ -0,0 +1,209 @@
part of '../latlng.dart';

/// Represents a 2D point on screen-space in pixels.
class Point2D {
const Point2D(this.x, this.y);
final double x;
final double y;

/// Offsets a line by given [amount] in pixels.
List<Point2D> offset(List<Point2D> points, double amount) {
final offsetSegments = _offsetPointLine(points, amount);
return _joinLineSegments(offsetSegments, amount);
}
}

void _forEachPair<T>(List<T> list, Function(T a, T b) callback) {
if (list.isEmpty) {
return;
}
for (var i = 1, l = list.length; i < l; i++) {
callback(list[i - 1], list[i]);
}
}

/// Find the coefficients (a,b) of a line of equation y = a.x + b,
/// or the constant x for vertical lines
/// Return null if there's no equation possible
class _LineEquation {
const _LineEquation({this.x, this.a, this.b});

final double? x;
final double? a;
final double? b;
}

class _OffsetSegment {
const _OffsetSegment({
required this.offsetAngle,
required this.original,
required this.offset,
});
final double offsetAngle;
final List<Point2D> original;
final List<Point2D> offset;
}

_LineEquation? _lineEquation(Point2D pt1, Point2D pt2) {
if (pt1.x == pt2.x) {
return pt1.y == pt2.y ? null : _LineEquation(x: pt1.x);
}

var a = (pt2.y - pt1.y) / (pt2.x - pt1.x);
return _LineEquation(
a: a,
b: pt1.y - a * pt1.x,
);
}

/// Return the intersection point of two lines defined by two points each
/// Return null when there's no unique intersection
Point2D? _intersection(Point2D l1a, Point2D l1b, Point2D l2a, Point2D l2b) {
var line1 = _lineEquation(l1a, l1b);
var line2 = _lineEquation(l2a, l2b);

if (line1 == null || line2 == null) {
return null;
}

if (line1.x != null) {
return line2.x != null
? null
: Point2D(
line1.x!,
line2.a! * line1.x! + line2.b!,
);
}
if (line2.x != null) {
return Point2D(
line2.x!,
line1.a! * line2.x! + line1.b!,
);
}

if (line1.a == line2.a) {
return null;
}

var x = (line2.b! - line1.b!) / (line1.a! - line2.a!);
return Point2D(
x,
line1.a! * x + line1.b!,
);
}

Point2D _translatePoint(Point2D pt, double dist, double heading) {
return Point2D(
pt.x + dist * cos(heading),
pt.y + dist * sin(heading),
);
}

List<_OffsetSegment> _offsetPointLine(List<Point2D> points, double distance) {
var offsetSegments = <_OffsetSegment>[];

_forEachPair<Point2D>(points, (a, b) {
if (a.x == b.x && a.y == b.y) {
return;
}

// angles in (-PI, PI]
var segmentAngle = atan2(a.y - b.y, a.x - b.x);
var offsetAngle = segmentAngle - pi / 2;

offsetSegments.add(_OffsetSegment(offsetAngle: offsetAngle, original: [
a,
b
], offset: [
_translatePoint(a, distance, offsetAngle),
_translatePoint(b, distance, offsetAngle)
]));
});

return offsetSegments;
}

/// Join 2 line segments defined by 2 points each with a circular arc
List<Point2D> _joinSegments(
_OffsetSegment s1,
_OffsetSegment s2,
double offset,
) {
// TO DO: different join styles
return _circularArc(s1, s2, offset)
.where((x) {
return x != null;
})
.map((e) => e!)
.toList();
}

List<Point2D> _joinLineSegments(List<_OffsetSegment> segments, double offset) {
var joinedPoints = <Point2D>[];

if (segments.isNotEmpty) {
var first = segments.first;
var last = segments.last;

joinedPoints.add(first.offset[0]);
_forEachPair<_OffsetSegment>(segments, (s1, s2) {
joinedPoints = [...joinedPoints, ..._joinSegments(s1, s2, offset)];
});
joinedPoints.add(last.offset[1]);
}

return joinedPoints;
}

Point2D _segmentAsVector(List<Point2D> s) {
return Point2D(
s[1].x - s[0].x,
s[1].y - s[0].y,
);
}

double _getSignedAngle(List<Point2D> s1, List<Point2D> s2) {
final a = _segmentAsVector(s1);
final b = _segmentAsVector(s2);
return atan2(a.x * b.y - a.y * b.x, a.x * b.x + a.y * b.y);
}

/// Interpolates points between two offset segments in a circular form
List<Point2D?> _circularArc(
_OffsetSegment s1, _OffsetSegment s2, double distance) {
// if the segments are the same angle,
// there should be a single join point
if (s1.offsetAngle == s2.offsetAngle) {
return [s1.offset[1]];
}

final signedAngle = _getSignedAngle(s1.offset, s2.offset);
// for inner angles, just find the offset segments intersection
if ((signedAngle * distance > 0) &&
(signedAngle * _getSignedAngle(s1.offset, [s1.offset[0], s2.offset[1]]) >
0)) {
return [
_intersection(s1.offset[0], s1.offset[1], s2.offset[0], s2.offset[1])
];
}

// draws a circular arc with R = offset distance, C = original meeting point
var points = <Point2D>[];
var center = s1.original[1];
// ensure angles go in the anti-clockwise direction
var rightOffset = distance > 0;
var startAngle = rightOffset ? s2.offsetAngle : s1.offsetAngle;
var endAngle = rightOffset ? s1.offsetAngle : s2.offsetAngle;
// and that the end angle is bigger than the start angle
if (endAngle < startAngle) {
endAngle += pi * 2;
}
var step = pi / 8;
for (var alpha = startAngle; alpha < endAngle; alpha += step) {
points.add(_translatePoint(center, distance, alpha));
}
points.add(_translatePoint(center, distance, endAngle));

return rightOffset ? points.reversed.toList() : points;
}
2 changes: 1 addition & 1 deletion pubspec.yaml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name: latlng
description: GeoJSON, Geodesy and Geographical calculations for Dart. Provides LatLong and Mercator projection (EPSG4326).
version: 2.0.2
version: 2.0.3
repository: https://github.com/xclud/dart_latlng
homepage: https://pwa.ir

Expand Down
7 changes: 7 additions & 0 deletions test/latlng_test.dart
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,13 @@ import 'package:test/test.dart';
import 'package:latlng/latlng.dart';

void main() {
test('Julian', () {
final date = DateTime.utc(2024, 2, 9, 22, 4, 11);
final julian = date.julian;

expect(julian.gmst.degrees, 8.943690756103024);
});

test('Lat, Long test', () {
final location = LatLng(Angle.degree(10), Angle.degree(20));
expect(location.latitude.degrees, 10);
Expand Down

0 comments on commit e670d8f

Please sign in to comment.