Fix julian calculations.

CHANGELOG.md

@@ -1,3 +1,7 @@
+## [1.0.3]
+
+* Fix julian calculations.
+
 ## [1.0.2]
 
 * Topocentric, ECF and ECI calculations.

lib/latlng.dart

@@ -2,6 +2,7 @@
 library latlng;
 
 import 'dart:convert';
+import 'dart:core';
 import 'dart:math';
 
 import 'src/lerp.dart' as l;
@@ -16,6 +17,7 @@ part 'src/eci.dart';
 part 'src/topocentric.dart';
 part 'src/look_angle.dart';
 part 'src/julian.dart';
+part 'src/const.dart';
 
 part 'src/geojson/feature.dart';
 part 'src/geojson/feature_collection.dart';

lib/src/const.dart

@@ -0,0 +1,10 @@
+part of '../latlng.dart';
+
+const double _rad2deg = 180.0 / pi;
+
+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

lib/src/ecef.dart

@@ -13,4 +13,20 @@ class EarthCenteredEarthFixed {
 
   /// Z Coordinate.
   final double z;
+
+  EarthCenteredEarthFixed operator +(EarthCenteredEarthFixed other) {
+    final dx = x + other.x;
+    final dy = x + other.y;
+    final dz = x + other.z;
+
+    return EarthCenteredEarthFixed(dx, dy, dz);
+  }
+
+  EarthCenteredEarthFixed operator -(EarthCenteredEarthFixed other) {
+    final dx = x - other.x;
+    final dy = x - other.y;
+    final dz = x - other.z;
+
+    return EarthCenteredEarthFixed(dx, dy, dz);
+  }
 }

lib/src/julian.dart

@@ -1,22 +1,5 @@
 part of '../latlng.dart';
 
-/// Calculates the Julian date.
-double _julian(int year, double doy) {
-  // Now calculate Julian date
-  // Ref: "Astronomical Formulae for Calculators", Jean Meeus, pages 23-25
-
-  year--;
-
-  // Centuries are not leap years unless they divide by 400
-  int A = year ~/ 100;
-  int B = 2 - A + (A ~/ 4);
-
-  double jan01 =
-      (365.25 * year).toInt() + (30.6001 * 14).toInt() + 1720994.5 + B;
-
-  return jan01 + doy;
-}
-
 double _dayOfYear(DateTime someDate) {
   final diff = someDate.difference(DateTime(someDate.year, 1, 1, 0, 0));
   final diffInDays = diff.inMilliseconds / 86400000.0;
@@ -26,32 +9,195 @@ double _dayOfYear(DateTime someDate) {
 
 /// A set of extensions on [DateTime] object.
 extension DateTimeExtensions on DateTime {
-  /// Calculates the Julian date.
-  double get julian {
-    final doy = _dayOfYear(this);
-    final julianDay = _julian(year, doy);
+  double get dayOfYear {
+    final diff = difference(DateTime(year, 1, 1, 0, 0));
+    final diffInDays = diff.inMilliseconds / 86400000.0;
 
-    return julianDay;
+    return diffInDays;
   }
 
+  /// Calculates the Julian date.
+  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 j = julian;
+    final result = julian.toGmst();
+
+    return result;
+  }
+}
+
+/// Encapsulates a Julian date.
+class Julian {
+  // Jan  1.5 2000 = Jan  1 2000 12h UTC
 
-    final d1 = j - 2451545.0;
-    final d2 = (j + 0.5) % 1.0;
-    final d3 = (d1 - d2) / 36525.0;
-    var d4 =
-        24110.54841 + d3 * (8640184.812866 + d3 * (0.093104 - d3 * 6.2E-06));
+  /// Create a Julian date object from a DateTime object. The time
+  /// contained in the DateTime object is assumed to be UTC.
+  ///
+  /// [utc] The UTC time to convert.
+  factory Julian.fromDateTime(DateTime utc) {
+    final v = _dayOfYear(utc) +
+        ((utc.hour +
+                ((utc.minute +
+                        ((utc.second + (utc.millisecond / 1000.0)) / 60.0)) /
+                    60.0)) /
+            24.0);
 
-    d4 = (d4 + 86636.555366976 * d2) % 86400.0;
+    return Julian._(v);
+  }
+
+  /// Initialize the Julian date object.
+  ///
+  /// The first day of the year, Jan 1, is day 1.0. Noon on Jan 1 is
+  /// represented by the day value of 1.5, etc.
+
+  const Julian._(this.value);
+
+  /// Create a Julian date object given a year and day-of-year.
+  ///
+  /// [year] The year, including the century (i.e., 2012).
+  /// [doy] Day of year (1 means January 1, etc.).
 
-    if (d4 < 0.0) {
-      d4 += 86400.0;
+  /// The fractional part of the day value is the fractional portion of
+  /// the day.
+  /// Examples:
+  ///    day = 1.0  Jan 1 00h
+  ///    day = 1.5  Jan 1 12h
+  ///    day = 2.0  Jan 2 00h
+  factory Julian.fromYearAndDoy(int year, double doy) {
+    // Arbitrary years used for error checking
+    if (year < 1900 || year > 2100) {
+      throw Exception('Year (1900, 2100)');
     }
 
-    final result = pi * 2.0 * (d4 / 86400.0);
+    // The last day of a leap year is day 366
+    if (doy < 1.0 || doy >= 367.0) {
+      throw Exception('Day (1, 367)');
+    }
 
-    return result;
+    // Now calculate Julian date
+    // Ref: "Astronomical Formulae for Calculators", Jean Meeus, pages 23-25
+
+    year--;
+
+    // Centuries are not leap years unless they divide by 400
+    int A = year ~/ 100;
+    int B = 2 - A + (A ~/ 4);
+
+    final jan01 =
+        (365.25 * year).floor() + (30.6001 * 14).floor() + 1720994.5 + B;
+
+    return Julian._(jan01 + doy);
+  }
+
+  /// Calculates the time difference between two Julian dates.
+  ///
+  /// Returns a [Duration] representing the time difference between the two dates.
+  Duration diffrence(Julian other) {
+    return toDateTime().difference(other.toDateTime());
+  }
+
+  /// The Julian date.
+  final double value;
+
+  /// Dec 31.5 1899 = Dec 31 1899 12h UTC
+  double fromJan0_12h_1900() => value - _j0H12Y1900;
+
+  /// Jan 1.0 1900 = Jan  1 1900 00h UTC
+  double fromJan1_00h_1900() => value - _j1H00Y1900;
+
+  /// Jan 1.5 1900 = Jan  1 1900 12h UTC
+  double fromJan1_12h_1900() => value - _j1H12Y1900;
+
+  /// Jan 1.5 2000 = Jan  1 2000 12h UTC
+  double fromJan1_12h_2000() => value - _j1H12Y2000;
+
+  /// Dec 31.5 1899 = Dec 31 1899 12h UTC
+  static const _j0H12Y1900 = 2415020.0;
+
+  /// Jan 1.0 1900 = Jan  1 1900 00h UTC
+  static const _j1H00Y1900 = 2415020.5;
+
+  /// Jan 1.5 1900 = Jan  1 1900 12h UTC
+  static const _j1H12Y1900 = 2415021.0;
+
+  /// Jan 1.5 2000 = Jan  1 2000 12h UTC
+  static const _j1H12Y2000 = 2451545.0;
+
+  /// Adds days.
+  Julian addDays(double day) => Julian._(value + day);
+
+  /// Adds hours.
+  Julian addHours(double hours) => Julian._(value + hours / _hoursPerDay);
+
+  /// Adds minutes.
+  Julian addMinutes(double min) => Julian._(value + min / _minutesPerDay);
+
+  /// Adds seconds.
+  Julian addSeconds(double sec) => Julian._(value + _secondsPerDay);
+
+  /// Calculate Greenwich Mean Sidereal Time for the Julian date.
+  ///
+  /// 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);
+  }
+
+  /// Calculate Local Mean Sidereal Time for this Julian date at the given
+  /// longitude.
+  ///
+  /// [longitude] The longitude, in radians, measured west from Greenwich.
+  ///
+  /// The angle, in radians, measuring eastward from the Vernal Equinox to
+  /// the given longitude.
+  double toLmst(double longitude) {
+    return (toGmst() + longitude) % _twopi;
+  }
+
+  /// Returns a UTC DateTime object that corresponds to this Julian date.
+  DateTime toDateTime() {
+    double d2 = value + 0.5;
+    int Z = d2.floor();
+    int alpha = ((Z - 1867216.25) / 36524.25).floor();
+    int A = Z + 1 + alpha - (alpha ~/ 4);
+    int B = A + 1524;
+    int C = ((B - 122.1) / 365.25).floor();
+    int D = (365.25 * C).floor();
+    int E = ((B - D) / 30.6001).floor();
+
+    // For reference: the fractional day of the month can be
+    // calculated as follows:
+    //
+    // double day = B - D - (int)(30.6001 * E) + F;
+
+    int month = (E <= 13) ? (E - 1) : (E - 13);
+    int year = (month >= 3) ? (C - 4716) : (C - 4715);
+
+    //final jdJan01 = Julian.fromYearAndDoy(year, 1.0);
+    //double doy = Value - jdJan01.Value; // zero-relative
+
+    final dtJan01 = DateTime(year, 1, 1, 0, 0, 0);
+
+    return dtJan01; //.add(Duration(days:  doy));
   }
 }

lib/src/planet.dart

@@ -30,6 +30,61 @@ abstract class Planet {
   /// Flattening of the planet.
   final double flattening;
 
+  /// Gets a polygon on the surface of the planet where a viewer can see.
+  List<LatLng> getGroundTrack(
+    LatLngAlt satellite, {
+    double precesion = 1,
+  }) {
+    final zone = <LatLng>[];
+
+    final latitude = satellite.latitude;
+    final longitude = satellite.longitude;
+    final altitude = satellite.altitude;
+
+    final cosLat = cos(latitude);
+    final sinLat = sin(latitude);
+
+    double num4 = acos(radius / altitude);
+    if (num4.isNaN) {
+      num4 = 0.0;
+    }
+    final cosNum4 = cos(num4);
+    final sinNum4 = sin(num4);
+    int i = 0;
+    do {
+      final angle = pi / 180.0 * i;
+      final lat = asin(sinLat * cosNum4 + cos(angle) * sinNum4 * cosLat);
+      final num9 = (cosNum4 - sinLat * sin(lat)) / (cosLat * cos(lat));
+      final lng = (((i != 0 || !(num4 > pi / 2.0 - latitude)) && 0 == 0)
+          ? (((i == 180 && num4 > pi / 2.0 + latitude))
+              ? (longitude + pi)
+              : ((num9.abs() > 1.0)
+                  ? longitude
+                  : ((i > 180)
+                      ? (longitude - acos(num9))
+                      : (longitude + acos(num9)))))
+          : (longitude + pi));
+
+      final z = LatLng(
+        lat * _rad2deg,
+        lng * _rad2deg,
+      );
+
+      zone.add(z);
+
+      i++;
+    } while (i <= 359);
+
+    zone.add(
+      LatLng(
+        zone[0].latitude,
+        zone[0].longitude,
+      ),
+    );
+
+    return zone;
+  }
+
   /// Calculates Topocentrics.
   Topocentric topocentric(
       LatLngAlt observer, EarthCenteredEarthFixed satellite) {

pubspec.yaml

@@ -1,6 +1,6 @@
 name: latlng
 description: GeoJSON, Geodesy and Geographical calculations for Dart. Provides LatLong and Mercator projection (EPSG4326).
-version: 1.0.2
+version: 1.0.3
 repository: https://github.com/xclud/dart_latlng
 homepage: https://pwa.ir