Add next pass time prediction

Author: RX14

README.md

@@ -50,6 +50,16 @@ satellite_position.to_lat_long_alt(time) # => Predict::LatLongAlt(@latitude=0.89
 # Predict look angles from an observer (200m altitude)
 observer = Predict::LatLongAlt.from_degrees(52.9, -2.24, 200.0)
 observer.look_at(satellite_position, time) # => Predict::LookAngle(@azimuth=0.50118495648349193, @elevation=-0.55812612380797122, @range=7501.0178628601843)
+
+# Predict next pass time
+start_time, end_time = satellite.next_pass(at: observer, after: time)
+start_time # => 2016-12-05 16:23:55
+end_time # => 2016-12-05 16:32:29
+
+# And the one after that...
+start_time, end_time = satellite.next_pass(at: observer, after: end_time)
+start_time # => 2016-12-05 17:58:46
+end_time # => 2016-12-05 18:09:15
 ```
 
 ## Development

spec/predict/satellite_spec.cr

@@ -197,4 +197,56 @@ describe Predict::Satellite do
     look_angle.elevation.should be_close(0.0574066, 0.5e-7)
     look_angle.range.should be_close(40812.259, 0.5e-3)
   end
+
+  describe "#next_pass" do
+    it "finds the next satellite pass" do
+      tle = Predict::TLE.parse_three_line <<-TLE
+        AO-51 [+]
+        1 28375U 04025K   09105.66391970  .00000003  00000-0  13761-4 0  3643
+        2 28375 098.0551 118.9086 0084159 315.8041 043.6444 14.40638450251959
+        TLE
+      time = Time.new(2009, 1, 5, 0, 0, 0)
+
+      satellite = Predict::Satellite.new(tle)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: time)
+      pass_start.should be_close(Time.new(2009, 1, 5, 4, 28, 10), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 4, 32, 15), 5.seconds)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: pass_end)
+      pass_start.should be_close(Time.new(2009, 1, 5, 6, 4, 0), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 6, 18, 0), 5.seconds)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: pass_end)
+      pass_start.should be_close(Time.new(2009, 1, 5, 7, 42, 45), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 7, 57, 50), 5.seconds)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: pass_end)
+      pass_start.should be_close(Time.new(2009, 1, 5, 9, 22, 5), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 9, 34, 20), 5.seconds)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: pass_end)
+      pass_start.should be_close(Time.new(2009, 1, 5, 11, 2, 5), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 11, 7, 35), 5.seconds)
+    end
+
+    it "finds the current pass with find_occuring_pass" do
+      tle = Predict::TLE.parse_three_line <<-TLE
+        AO-51 [+]
+        1 28375U 04025K   09105.66391970  .00000003  00000-0  13761-4 0  3643
+        2 28375 098.0551 118.9086 0084159 315.8041 043.6444 14.40638450251959
+        TLE
+      time = Time.new(2009, 1, 5, 4, 30, 0)
+
+      satellite = Predict::Satellite.new(tle)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: time)
+      pass_start.should be_close(Time.new(2009, 1, 5, 6, 4, 0), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 6, 18, 0), 5.seconds)
+
+      pass_start, pass_end = satellite.next_pass(at: GROUND_STATION, after: time, find_occuring_pass: true)
+      pass_start.should be_close(Time.new(2009, 1, 5, 4, 28, 10), 5.seconds)
+      pass_end.should be_close(Time.new(2009, 1, 5, 4, 32, 15), 5.seconds)
+    end
+  end
 end

src/predict/satellite.cr

@@ -1,12 +1,45 @@
+private def converge_root(a, b, *, accuracy = 1e-6, max_iterations = 40)
+  xa = yield a
+  xb = yield b
+
+  raise "Root is not bracketed" unless (xa * xb) <= 0
+
+  if xa < 0
+    # Function is increasing a -> b
+    dx = b - a
+    root = a
+  else
+    # Function is decreasing a -> b
+    dx = a - b
+    root = b
+  end
+
+  # In this algorithm, the search section is described by the lower bound of
+  # *root*, and the upper bound of *root + dx*. Halving dx is equivalent to
+  # taking the left section, halving dx and setting root to be the midsection is
+  # equivalent to taking the right section.
+
+  max_iterations.times do |i|
+    dx *= 0.5                        # halve dx
+    xmid = root + dx                 # calculate x value of midsection
+    fmid = yield xmid                # find function value at midsection
+    root = xmid if fmid <= 0         # midsection is below root, take right section
+    return root if dx.abs < accuracy # the root is always within (root + dx), so if dx < accuracy we can return
+  end
+
+  raise "Root finding failed: too many iterations"
+end
+
 module Predict
   # Struct containing an initialised SGP4 model of a satellite.
-  struct Satellite
+  class Satellite
+    getter tle : TLE
     getter satrec = SGP4::Elset.new
     getter constants : GravityConstants
     @epoch : Time
 
     # Initialises the `Satellite` with *tle*, using *grav_type*
-    def initialize(tle : TLE, grav_type = GravityConstants::WGS72)
+    def initialize(@tle : TLE, grav_type = GravityConstants::WGS72)
       @constants = GravityConstants[grav_type]
       @epoch = tle.epoch
       initialize_satrec(tle)
@@ -64,6 +97,89 @@ module Predict
       {position, velocity}
     end
 
+    private def look_angles(location, time)
+      satellite_position, _ = predict(time)
+      location.look_at(satellite_position, time)
+    end
+
+    # Finds the next pass at *location* occuring after the given time. When
+    # *find_occuring_pass* is true, a satellite pass currently in progress at
+    # the supplied time may be returned. The returned value is a tuple of pass
+    # start time followed by pass end time.
+    def next_pass(*, at : LatLongAlt, after : Time, find_occuring_pass = false) : {Time, Time}
+      time = after
+      location = at
+
+      raise "Satellite will never be visible" unless pass_possible? location
+
+      orbit_time = (MINS_PER_DAY / tle.mean_motion).minutes
+
+      # Wind back 1/4 orbit if we want to find currently occurring passes
+      time -= orbit_time / 4 if find_occuring_pass
+
+      look_angles = look_angles(location, time)
+
+      # Ensure satellite is not above the horizon
+      if look_angles.elevation > 0
+        # Move forward in 60 second intervals until the satellite goes below the horizon
+        while look_angles.elevation > 0
+          time += 60.seconds
+          look_angles = look_angles(location, time)
+        end
+
+        # Move forward 3/4 orbit
+        time += orbit_time * 0.75
+        look_angles = look_angles(location, time)
+      end
+
+      # Find the time it comes over the horizon
+      while look_angles.elevation < 0
+        time += 60.seconds
+        look_angles = look_angles(location, time)
+      end
+
+      # Find pass start
+      start_time = converge_root(time, time - 60.seconds, accuracy: 1.millisecond) do |time|
+        look_angles(location, time).elevation
+      end
+
+      # Find the time when it goes below
+      while look_angles.elevation > 0
+        time += 30.seconds
+        look_angles = look_angles(location, time)
+      end
+
+      # p look_angles.elevation
+      # p look_angles(location, time - 30.seconds).elevation
+
+      # Find los time
+      end_time = converge_root(time, time - 30.seconds, accuracy: 1.millisecond) do |time|
+        look_angles(location, time).elevation
+      end
+
+      {start_time, end_time}
+    end
+
+    # Returns true if it's possible for the satellite to be visible from
+    # *location*.
+    def pass_possible?(location : LatLongAlt)
+      return false if tle.mean_motion < 1e-8
+
+      incl = tle.inclination
+      incl = 180 - incl if incl >= 90.0
+      incl *= DEG2RAD
+
+      # Cube root of Kepler's third law constant for the earth and satellite
+      # with negligible mass, expressed in (km^3/s^2).
+      k = 331.25
+
+      semi_major_axis = k * (MINS_PER_DAY / tle.mean_motion)**TWO_THIRDS
+      apogee = semi_major_axis * (1.0 + tle.eccentricity) # From center of earth
+      ground_track_angle = Math.acos(@constants.radiusearthkm / apogee)
+
+      ground_track_angle + incl > location.latitude.abs
+    end
+
     private def check_error
       case @satrec.error
       when 1