orbits-lib/resources/orbits/orbit.rb
##### equator_crossing
#
# Find the points at which the orbit passing through coord
# crosses the equator on its ascending pass (if ascending
# is true) or its descending pass (if ascending is false).
# If corrected is true, we will correct for the rotation of
# the Earth, otherwise we won't.
equator_crossing
Find the points at which the orbit passing through coord
crosses the equator on its ascending pass (if ascending
is true) or its descending pass (if ascending is false).
If corrected is true, we will correct for the rotation of
the Earth, otherwise we won't.
equator_crossing Find the points at which the orbit passing through coord crosses the equator on its ascending pass ( if ascending is true ) or its descending pass ( if ascending is false ) . If corrected is true we will correct for the rotation of the Earth otherwise we won ' t .
def equator_crossing(coord, ascending, corrected=true)
       [This image](http://geospatialmethods.org/bosa/triangles.gif)
       is also helpful.  (`Lat_p`, `Lon_p`) is what we call `C`.
       `(0, Lon_p)` is what we call `B`.  `(0, Lon_n)` is
       what we call `A`.
       In order to figure out the crossing point, we need to
       set up a spherical triangle formed by:
       point `C` - the given coordinate
       point `B` - the point at (0, coord.theta)
       point `A` - the (unknown) point where the orbit crosses the equator
       To remain consistent with the Wolfram and similar
       references, we'll use A, B, and C notation.  "`A`" will
       denote the angle at point `A`.  "`a`" will denote the arc-length
       of the side across from point `A`.  Our capital variable names will
       use underscores to remain valid in Ruby, e.g. "_A_"
       We know two angles:
       Because side `A` is on a latitude line and side `c` is on
       the equator
      _B_ = PI / 2
       `A` is an angle between the orbit and the equator, which
       is the inclination if it's a retrograde orbit in the
       descending pass or a prograde orbit in the ascending
       pass.  Otherwise it's the complement of the inclination
      if retrograde? != ascending
        _A_ = inclination_rad
      else
        _A_ = PI - inclination_rad
      end
       We also know one side's arc length:
      a = coord.phi
       Coord could lie above or below the orbit.  We find the nearest
       inflection point in that case.
      a = max_orbit_phi if a > max_orbit_phi
      a = min_orbit_phi if a < min_orbit_phi
       We need to find:
       1. `c`, which we can add to coord's longitude to find
           the crossing longitude
       2. `b`, which is the arc length traveled since crossing the
           equator.  We can use this to calculate the amount
           the Earth has rotated since the satellite crossed
           the equator and correct accordingly.
       Since we know `A`, `a`, and `B`, we can use the law of sines
       to calculate `b`
       `sin(A) / sin(a) = sin(B) / sin(b) = sin(\C) / sin(c)`
      sin_b = sin(_B_) * sin(a) / sin(_A_)
      b = asin(sin_b)
       Finding `c` is trickier because we don't know `C`.  We use [Napier's
       first analogy](http://mathworld.wolfram.com/NapiersAnalogies.html):
       `sin((A - B)/2) / sin((A + B)/2) = tan((a - b)/2) / tan(c/2)`
       Rearranging:
       `tan(c/2) = (tan((a - b)/2) * sin((A + B)/2)) / sin((A - B)/2)`
       Since `sin((A - B)/2)` could be 0, we'll need to use atan2 and
       keep the numerator and denominator separate
      numerator = tan((a - b) / 2) * sin(( _A_ + _B_ ) / 2)
      denominator = sin((_A_ - _B_) / 2)
      c = 2 * atan2(numerator, denominator)
       We now know where the orbit crosses
      crossing_theta = coord.theta + c
       For the descending pass, we ended up calculating the point half
       way around the globe, so we need to correct
      crossing_theta -= PI unless ascending
       Correct for the rotation of the Earth, if necessary
      if corrected
         Recall, `b` is the angle of the orbit travelled since
         crossing the equator on the ascending pass (in radians).
         We can use this to calculate the angular distance covered
         by the orbit on the requested pass.
        if ascending
          distance_rad = b
        else
          distance_rad = PI - b
        end
         If the distance is below the starting latitude, the orbit will
         need to make an additional pass before reaching the target point,
         so we add `2pi` to the distance.
        distance_rad += 2 * PI while distance_rad < start_clat_rad
        distance_rad -= 2 * PI while distance_rad > start_clat_rad + 2 * PI
         `"time" = "distance" / "rate"`
         We know distance from above and rate is given by the orbit parameters,
         so we can calculate the time since crossing the equator.
        time_since_crossing_s = distance_rad / angular_velocity_rad_s
         Now we want to know how far the Earth has rotated since we crossed
         the equator.  We know the Earth's rotation rate and the elapsed time
        earth_rotation_rad = EARTH_ANGULAR_VELOCITY_RAD_S * time_since_crossing_s
         Since we're working with a unit sphere, we can make the
         correction directly
        crossing_theta += earth_rotation_rad
      end
       Return the coordinate
      Coordinate.phi_theta(0, crossing_theta)
    end
def equator_crossing ( coord , ascending , corrected = true ) _B_ = PI / 2 if retrograde? != ascending _A_ = inclination_rad else _A_ = PI - inclination_rad end a = coord . phi a = max_orbit_phi if a > max_orbit_phi a = min_orbit_phi if a < min_orbit_phi sin_b = sin ( _B_ ) * sin ( a ) / sin ( _A_ ) b = asin ( sin_b ) numerator = tan ( ( a - b ) / 2 ) * sin ( ( _A_ + _B_ ) / 2 ) denominator = sin ( ( _A_ - _B_ ) / 2 ) c = 2 * atan2 ( numerator , denominator ) crossing_theta = coord . theta + c crossing_theta -= PI unless ascending if corrected if ascending distance_rad = b else distance_rad = PI - b end distance_rad += 2 * PI while distance_rad < start_clat_rad distance_rad -= 2 * PI while distance_rad > start_clat_rad + 2 * PI time_since_crossing_s = distance_rad / angular_velocity_rad_s earth_rotation_rad = EARTH_ANGULAR_VELOCITY_RAD_S * time_since_crossing_s crossing_theta += earth_rotation_rad end Coordinate . phi_theta ( 0 , crossing_theta ) end
equator_crossing
Find the points at which the orbit passing through coord
crosses the equator on its ascending pass (if ascending
is true) or its descending pass (if ascending is false).
equator_crossing Find the points at which the orbit passing through coord crosses the equator on its ascending pass ( if ascending is true ) or its descending pass ( if ascending is false ) .
# [This page](http://mathworld.wolfram.com/SphericalTrigonometry.html) # has a good diagram and a summary of the math we'll need to use. # # [This image](http://geospatialmethods.org/bosa/triangles.gif) # is also helpful.  (`Lat_p`, `Lon_p`) is what we call `C`. # `(0, Lon_p)` is what we call `B`.  `(0, Lon_n)` is # what we call `A`. # # In order to figure out the crossing point, we need to # set up a spherical triangle formed by: # # point `C` - the given coordinate # point `B` - the point at (0, coord.theta) # point `A` - the (unknown) point where the orbit crosses the equator # # To remain consistent with the Wolfram and similar # references, we'll use A, B, and C notation.  "`A`" will # denote the angle at point `A`.  "`a`" will denote the arc-length # of the side across from point `A`.  Our capital variable names will # use underscores to remain valid in Ruby, e.g. "_A_" # # We know two angles: # Because side `A` is on a latitude line and side `c` is on # the equator # `A` is an angle between the orbit and the equator, which # is the inclination if it's a retrograde orbit in the # descending pass or a prograde orbit in the ascending # pass.  Otherwise it's the complement of the inclination # We also know one side's arc length: # Coord could lie above or below the orbit.  We find the nearest # inflection point in that case. # We need to find: # # 1. `c`, which we can add to coord's longitude to find #     the crossing longitude # 2. `b`, which is the arc length traveled since crossing the #     equator.  We can use this to calculate the amount #     the Earth has rotated since the satellite crossed #     the equator and correct accordingly. # Since we know `A`, `a`, and `B`, we can use the law of sines # to calculate `b` # # `sin(A) / sin(a) = sin(B) / sin(b) = sin(\C) / sin(c)` # Finding `c` is trickier because we don't know `C`.  We use [Napier's # first analogy](http://mathworld.wolfram.com/NapiersAnalogies.html): # # `sin((A - B)/2) / sin((A + B)/2) = tan((a - b)/2) / tan(c/2)` # # Rearranging: # # `tan(c/2) = (tan((a - b)/2) * sin((A + B)/2)) / sin((A - B)/2)` # # Since `sin((A - B)/2)` could be 0, we'll need to use atan2 and # keep the numerator and denominator separate # We now know where the orbit crosses # For the descending pass, we ended up calculating the point half # way around the globe, so we need to correct # Correct for the rotation of the Earth, if necessary # Recall, `b` is the angle of the orbit travelled since # crossing the equator on the ascending pass (in radians). # We can use this to calculate the angular distance covered # by the orbit on the requested pass. # If the distance is below the starting latitude, the orbit will # need to make an additional pass before reaching the target point, # so we add `2pi` to the distance. # `"time" = "distance" / "rate"` # # We know distance from above and rate is given by the orbit parameters, # so we can calculate the time since crossing the equator. # Now we want to know how far the Earth has rotated since we crossed # the equator.  We know the Earth's rotation rate and the elapsed time # Since we're working with a unit sphere, we can make the # correction directly # Return the coordinate
coord
ascending
corrected