orbits-lib/resources/orbits/orbit.rb
def horizontal_swath_edges(coord, ascending)
      # Candidate solutions.  We may end up finding more points than we need.
      # Before returning, we'll take the two points closest to the given
      # coordinate
      candidates = []

      # We need to be able to assume the coordinate isn't on the equator, so
      # we'll return early if we find that it is.  We know the swath width
      # from the orbit parameters.  This calculation is slightly off, since
      # the swath width is measured perpendicular to the orbit, but the
      # Java code makes this assumption.  This case should almost never
      # happen.  When it does, the numbers are likely close enough.

      if coord.phi.abs <= EPSILON
        return [Coordinate.phi_theta(0, coord.theta - swath_width_rad / 2),
                Coordinate.phi_theta(0, coord.theta + swath_width_rad / 2)]
      end

      # We ignore the rotation of the Earth for the purposes of this method.
      # The values we return (the swath edges) are not influenced by the
      # Earth's rotation.  Intermediate values are, but they're not visible
      # outside of this method.
      #
      # Ignoring the rotation of the Earth, we can imagine three circles
      # traced by the orbit.  The first is the path on the ground directly
      # beneath the orbit.  This is a great circle whose center is the origin.
      # The other two are the paths traced by the swath edges, which are not
      # great circles and not centered at the origin.
      #
      # These three circles each lie in a plane, and the three planes are
      # parallel to each other.
      #
      # There's another circle we're interested in, which is the circle
      # traced by the latitude line at coord.  This circle lies in its own
      # plane which is not parallel to the other three.
      #
      # The four original circles are the intersections of their respective
      # planes and the sphere of the Earth.
      #
      # So, to find the swath's edge coordinates, we need to find the
      # intersection between three things:
      #
      #   1. `x^2 + y^2 + z^2 = 1`
      #
      #      The sphere of the Earth.  We've scaled everything so the Earth
      #      has radius 1, and it's centered on the origin.  Convenient.
      #
      #   2. `ax + by + cz = d`
      #
      #      The swath's plane.  `a`, `b`, `c` and `d` are constants which
      #      we'll calculate.  Since they're parallel, the swath edges share
      #      `a`, `b`, and `c` with the orbital plane, but each edge has a
      #      different `d` value.
      #
      #   3. `z = "coord.z"`
      #
      #      The plane slicing through latitude = coord.lat

      # We already know z

      z = coord.z

      # We'll find `a`, `b`, and `c` for the orbit's plane, which are the same
      # constants for the swath's plane.  To find them, we need 3 non-colinear
      # points in the orbit's plane.  The origin is one, but we need two more.

      # First, we'll get a point that the orbit can pass over which is close
      # to coord.  Usually we'll just use coord, but if coord's latitude is
      # above the orbit's inflection point, we need to find a point with a
      # lower latitude

      target_phi = coord.phi
      target_phi = max_orbit_phi if target_phi > max_orbit_phi
      target_phi = min_orbit_phi if target_phi < min_orbit_phi
      target = Coordinate.phi_theta(target_phi, coord.theta)

      # For the second point, we'll use the equatorial crossing of the orbit
      # through the target.  This is safe, because we bailed earlier if
      # the target was on the equator.  The "false" parameter tells the method
      # to ignore the Earth's rotation.

      cross = equator_crossing(target, ascending, false)

      # For a plane passing through the origin, `ax + by + cz = 0`.  With our
      # two other points, we can solve for `a`, `b`, and `c`:

      a = target.y * cross.z - target.z * cross.y
      b = target.z * cross.x - target.x * cross.z
      c = target.x * cross.y - target.y * cross.x

      # If we can find a point on the swath edge, we can calculate `d`.  We
      # can find a point on each edge by finding the orbit's inflection point
      # and adding/subtracting half the swath width from its latitude.

      inf = Coordinate.phi_theta(-inclination_rad, cross.theta - PI / 2)

      # Calculate solutions for the point p on each swath edge

      [Coordinate.phi_theta(inf.phi + swath_width_rad / 2, inf.theta),
       Coordinate.phi_theta(inf.phi - swath_width_rad / 2, inf.theta)].each do |p|

        # Calculate `d` from the plane equation

        d = a * p.x + b * p.y + c * p.z

        # The fact that `z` is constant is very handy.  It means we really only
        # care about two dimensions.  We're really trying to find the intersection
        # of a line (the edge's plane when `z = "coord.z"`) and a circle (the part
        # of the Earth intersecting `z = "coord.z"`).
        #
        # Note: that simplified math could make it easier to switch this
        # algorithm to use an ellipsoidal Earth (WGS84) in the future, since we
        # would just be intersecting a line with a slightly different circle.
        #
        # The swath plane's equation is:
        #
        # `ax + by + cz = d`
        #
        # `a`, `b`, `c`, `d`, and `z` are all known at this point, so we just
        # need to figure out `x` and `y`.  Solve the equation for `y`
        #
        # `y = -a/b x + (d - cz) / b`
        #
        # Which is the line we need to intersect with the latitude circle.

        # We can't just compute constants for the line, because `b` could be near 0.
        # Let's deal with that special case first.

        if b.abs < 0.000001

          # If `b = 0`, then the plane's equation is:
          #
          # `ax + 0y + cz = d`
          #
          # Solving for `x`:
          #
          # `x = (d - cz) / a`

          # Ok, another problem.  It's possible that both `a` and `b` are 0.  In that
          # case, we have an orbit that follows the equator.  If coord lies within the
          # swath, then every orbit will cover the point, otherwise no orbits will.

          if a.abs < 0.000001
            if coord.phi.abs < swath_width_rad / 2
              return [Coordinate.phi_theta(0, -PI), Coordinate.phi_theta(0, PI)]
            else
              return []
            end
          end

          # Now we can handle the case where `b` is 0 and `a` isn't

          x = (d - c * z) / a

          # We now know `x` and `z`, so we can use the sphere equation to solve for `y`.
          #
          # `y = +- sqrt(1 - z^2 - x^2)`
          #
          # In this case, `1 - z^2 - x^2` could be negative.  This would indicate
          # that the latitude does not intersect the swath edge.  Perhaps the point is
          # completely above one swath edge but below the other.  We won't worry
          # about duplicate solutions (`y = +- 0`), since the LongitudeCoverage
          # class will combine them correctly.

          dis = 1 - z**2 - x**2

          if dis >= 0
            y = sqrt(dis)
            candidates << Coordinate.xyz(x, y, z)
            candidates << Coordinate.xyz(x, -y, z)
          end

        # Now the case when b is not 0.
        else
          # Recall, the line's equation is:
          #
          # `y = -a/b x + (d - cz) / b`
          #
          # We'll call those constants `m` and `k`

          m = -a / b
          k = (d - c * z) / b

          # Which means
          #
          # `y = mx + k`
          #
          # (We won't use the standard `b` as the `y`-intercept name, since
          # it means something different in the plane equation).

          # We need to find the intersection of that line with the latitude
          # circle at `z = "coord.z"`.  Using the sphere equation with a
          # constant z, we get the latitude circle's equation:
          #
          #  `x^2 + y^2 = 1 - z^2`

          # It's centered at the origin.  `1 - z^2` is the radius squared.
          #
          # Taking the square root is safe because we calculated z from a
          # lat/lon coordinate.  Further, a negative radius doesn't make sense.

          r = sqrt(1 - z**2)

          # So,
          #
          # `x^2 + y^2 = r^2`
          #
          # Substituting `y` with the value from our line equation:
          #
          # `x^2 + (mx + k)^2 = r^2`
          #
          # Multiplying out, we get:
          #
          # `x^2 + m^2x^2 + 2mkx + k^2 = r^2`

          # Grouping constants, we get:
          #
          # `(m^2 + 1)x^2 + (2mk)x + (k^2 - r^2) = 0`
          #
          # We call those constants `a_q`, `b_q`, and `c_q`

          a_q = m**2 + 1
          b_q = 2 * m * k
          c_q = k**2 - r**2

          # So:
          #
          # `a_q x^2 + b_qx + c_q = 0`

          # We can solve this using the quadratic equation:
          #
          # `x = (-b +- sqrt(b^2 - 4ac)) / (2a)`

          discriminant = b_q**2 - 4 * a_q * c_q

          # The discriminant could be negative if the line and circle do
          # not intersect
          if discriminant >= 0
            # There are two possible `x` values
            [(-b_q + sqrt(discriminant)) / (2*a_q),
             (-b_q - sqrt(discriminant)) / (2*a_q)].each do |x2|
              # Now we need `y`
              y = m * x2 + k

              candidates << Coordinate.xyz(x2, y, z)
            end
          end
        end
      end

      # Return if we haven't found any valid points
      return [] if candidates.empty?

      # At this point, we have either 2 or 4 candidate points that are swath
      # edge intersections with the target latitude.  Now we need to find
      # the two that are closest to the coordinate the caller passed.  We
      # want them in west-to-east order to make things simpler for the caller.

      # Find the index of the first candidate to the east of the target
      candidates.sort_by!(&:theta)
      east_index = candidates.find_index { |c2| coord.theta < c2.theta }

      # Handle wrapping around the date line
      east_index = 0 if east_index.nil?

      # The west index is the one to the left of the east index, wrapping around
      # the date line if needed
      west_index = (east_index - 1) % candidates.size

      west_edge = candidates[west_index]
      east_edge = candidates[east_index]

      # There is a special case here.  If we only found two candidates, that
      # indicates that the point is above one swath edge, but below the other,
      # which happens at latitudes near the inflection point.
      #
      # This method returns swath edges for either the ascending or descending
      # pass, not both.  If both of our candidates are on the same edge, then
      # one is seen on the ascending pass and the other on the descending. We
      # solve this problem by moving one edge point to the inflection longitude.
      #
      # Note: the rest of this method will produce the same results regardless
      #       of the value of "ascending."  Factoring this out could slightly
      #       improve performance.

      xprod = target * cross
      north_hemisphere = coord.phi > 0
      if candidates.size == 2
        inf_theta = north_hemisphere ? inf.theta + PI : inf.theta
        inflection = Coordinate.phi_theta(coord.phi, inf_theta)
        if north_hemisphere && retrograde? == ascending ||
            !north_hemisphere && retrograde? != ascending
          west_edge = inflection
        else
          east_edge = inflection
        end
      end

      # Return the edges
      [west_edge, east_edge]
    end
##### horizontal_swath_edges
#
# Given coord, a coordinate within the orbital swath, returns an array
# of 2 points with the same latitude as coord which lie on the swath edges.
# The first point lies on the edge to the west of coord, the second to the
# east.
#
# If the coordinate cannot fall within the swath, returns an empty array
#
# Parameters:
# * **coord** - A coordinate (Orbit::Coordinate) within the swath
# * **ascending** - A boolean value.  If true, we return edges for the
#   orbit's ascending pass.  If false, we return edges for the descending
#   pass.
#
# Returns an array containing 0 or 2 Orbit::Coordinate instances
horizontal_swath_edges
Given coord, a coordinate within the orbital swath, returns an array
of 2 points with the same latitude as coord which lie on the swath edges.
The first point lies on the edge to the west of coord, the second to the
east.

If the coordinate cannot fall within the swath, returns an empty array

coord** - A coordinate (Orbit::Coordinate) within the swath
ascending** - A boolean value.  If true, we return edges for the
orbit's ascending pass.  If false, we return edges for the descending
pass.

Returns an array containing 0 or 2 Orbit::Coordinate instances
horizontal_swath_edges Given coord a coordinate within the orbital swath returns an array of 2 points with the same latitude as coord which lie on the swath edges . The first point lies on the edge to the west of coord the second to the east . If the coordinate cannot fall within the swath returns an empty array coord ** - A coordinate ( Orbit :: Coordinate ) within the swath ascending ** - A boolean value . If true we return edges for the orbit ' s ascending pass . If false we return edges for the descending pass . Returns an array containing 0 or 2 Orbit :: Coordinate instances
def horizontal_swath_edges(coord, ascending)
      candidates = []
      if coord.phi.abs <= EPSILON
        return [Coordinate.phi_theta(0, coord.theta - swath_width_rad / 2),
                Coordinate.phi_theta(0, coord.theta + swath_width_rad / 2)]
      end
       Ignoring the rotation of the Earth, we can imagine three circles
       traced by the orbit.  The first is the path on the ground directly
       beneath the orbit.  This is a great circle whose center is the origin.
       The other two are the paths traced by the swath edges, which are not
       great circles and not centered at the origin.
       These three circles each lie in a plane, and the three planes are
       parallel to each other.
       There's another circle we're interested in, which is the circle
       traced by the latitude line at coord.  This circle lies in its own
       plane which is not parallel to the other three.
       The four original circles are the intersections of their respective
       planes and the sphere of the Earth.
       So, to find the swath's edge coordinates, we need to find the
       intersection between three things:
         1. `x^2 + y^2 + z^2 = 1`
            The sphere of the Earth.  We've scaled everything so the Earth
            has radius 1, and it's centered on the origin.  Convenient.
         2. `ax + by + cz = d`
            The swath's plane.  `a`, `b`, `c` and `d` are constants which
            we'll calculate.  Since they're parallel, the swath edges share
            `a`, `b`, and `c` with the orbital plane, but each edge has a
            different `d` value.
         3. `z = "coord.z"`
            The plane slicing through latitude = coord.lat
       We already know z
      z = coord.z
       We'll find `a`, `b`, and `c` for the orbit's plane, which are the same
       constants for the swath's plane.  To find them, we need 3 non-colinear
       points in the orbit's plane.  The origin is one, but we need two more.
       First, we'll get a point that the orbit can pass over which is close
       to coord.  Usually we'll just use coord, but if coord's latitude is
       above the orbit's inflection point, we need to find a point with a
       lower latitude
      target_phi = coord.phi
      target_phi = max_orbit_phi if target_phi > max_orbit_phi
      target_phi = min_orbit_phi if target_phi < min_orbit_phi
      target = Coordinate.phi_theta(target_phi, coord.theta)
       For the second point, we'll use the equatorial crossing of the orbit
       through the target.  This is safe, because we bailed earlier if
       the target was on the equator.  The "false" parameter tells the method
       to ignore the Earth's rotation.
      cross = equator_crossing(target, ascending, false)
       For a plane passing through the origin, `ax + by + cz = 0`.  With our
       two other points, we can solve for `a`, `b`, and `c`:
      a = target.y * cross.z - target.z * cross.y
      b = target.z * cross.x - target.x * cross.z
      c = target.x * cross.y - target.y * cross.x
       If we can find a point on the swath edge, we can calculate `d`.  We
       can find a point on each edge by finding the orbit's inflection point
       and adding/subtracting half the swath width from its latitude.
      inf = Coordinate.phi_theta(-inclination_rad, cross.theta - PI / 2)
       Calculate solutions for the point p on each swath edge
      [Coordinate.phi_theta(inf.phi + swath_width_rad / 2, inf.theta),
       Coordinate.phi_theta(inf.phi - swath_width_rad / 2, inf.theta)].each do |p|
         Calculate `d` from the plane equation
        d = a * p.x + b * p.y + c * p.z
         The fact that `z` is constant is very handy.  It means we really only
         care about two dimensions.  We're really trying to find the intersection
         of a line (the edge's plane when `z = "coord.z"`) and a circle (the part
         of the Earth intersecting `z = "coord.z"`).
         Note: that simplified math could make it easier to switch this
         algorithm to use an ellipsoidal Earth (WGS84) in the future, since we
         would just be intersecting a line with a slightly different circle.
         The swath plane's equation is:
         `ax + by + cz = d`
         `a`, `b`, `c`, `d`, and `z` are all known at this point, so we just
         need to figure out `x` and `y`.  Solve the equation for `y`
         `y = -a/b x + (d - cz) / b`
         Which is the line we need to intersect with the latitude circle.
         We can't just compute constants for the line, because `b` could be near 0.
         Let's deal with that special case first.
        if b.abs < 0.000001
           If `b = 0`, then the plane's equation is:
           `ax + 0y + cz = d`
           Solving for `x`:
           `x = (d - cz) / a`
           Ok, another problem.  It's possible that both `a` and `b` are 0.  In that
           case, we have an orbit that follows the equator.  If coord lies within the
           swath, then every orbit will cover the point, otherwise no orbits will.
          if a.abs < 0.000001
            if coord.phi.abs < swath_width_rad / 2
              return [Coordinate.phi_theta(0, -PI), Coordinate.phi_theta(0, PI)]
            else
              return []
            end
          end
           Now we can handle the case where `b` is 0 and `a` isn't
          x = (d - c * z) / a
           We now know `x` and `z`, so we can use the sphere equation to solve for `y`.
           `y = +- sqrt(1 - z^2 - x^2)`
           In this case, `1 - z^2 - x^2` could be negative.  This would indicate
           that the latitude does not intersect the swath edge.  Perhaps the point is
           completely above one swath edge but below the other.  We won't worry
           about duplicate solutions (`y = +- 0`), since the LongitudeCoverage
           class will combine them correctly.
          dis = 1 - z**2 - x**2
          if dis >= 0
            y = sqrt(dis)
            candidates << Coordinate.xyz(x, y, z)
            candidates << Coordinate.xyz(x, -y, z)
          end
         Now the case when b is not 0.
        else
           Recall, the line's equation is:
           `y = -a/b x + (d - cz) / b`
           We'll call those constants `m` and `k`
          m = -a / b
          k = (d - c * z) / b
           Which means
           `y = mx + k`
           (We won't use the standard `b` as the `y`-intercept name, since
           it means something different in the plane equation).
           We need to find the intersection of that line with the latitude
           circle at `z = "coord.z"`.  Using the sphere equation with a
           constant z, we get the latitude circle's equation:
            `x^2 + y^2 = 1 - z^2`
           It's centered at the origin.  `1 - z^2` is the radius squared.
           Taking the square root is safe because we calculated z from a
           lat/lon coordinate.  Further, a negative radius doesn't make sense.
          r = sqrt(1 - z**2)
           So,
           `x^2 + y^2 = r^2`
           Substituting `y` with the value from our line equation:
           `x^2 + (mx + k)^2 = r^2`
           Multiplying out, we get:
           `x^2 + m^2x^2 + 2mkx + k^2 = r^2`
           Grouping constants, we get:
           `(m^2 + 1)x^2 + (2mk)x + (k^2 - r^2) = 0`
           We call those constants `a_q`, `b_q`, and `c_q`
          a_q = m**2 + 1
          b_q = 2 * m * k
          c_q = k**2 - r**2
           So:
           `a_q x^2 + b_qx + c_q = 0`
           We can solve this using the quadratic equation:
           `x = (-b +- sqrt(b^2 - 4ac)) / (2a)`
          discriminant = b_q**2 - 4 * a_q * c_q
           The discriminant could be negative if the line and circle do
           not intersect
          if discriminant >= 0
             There are two possible `x` values
            [(-b_q + sqrt(discriminant)) / (2*a_q),
             (-b_q - sqrt(discriminant)) / (2*a_q)].each do |x2|
               Now we need `y`
              y = m * x2 + k
              candidates << Coordinate.xyz(x2, y, z)
            end
          end
        end
      end
       Return if we haven't found any valid points
      return [] if candidates.empty?
       At this point, we have either 2 or 4 candidate points that are swath
       edge intersections with the target latitude.  Now we need to find
       the two that are closest to the coordinate the caller passed.  We
       want them in west-to-east order to make things simpler for the caller.
       Find the index of the first candidate to the east of the target
      candidates.sort_by!(&:theta)
      east_index = candidates.find_index { |c2| coord.theta < c2.theta }
       Handle wrapping around the date line
      east_index = 0 if east_index.nil?
       The west index is the one to the left of the east index, wrapping around
       the date line if needed
      west_index = (east_index - 1) % candidates.size
      west_edge = candidates[west_index]
      east_edge = candidates[east_index]
       There is a special case here.  If we only found two candidates, that
       indicates that the point is above one swath edge, but below the other,
       which happens at latitudes near the inflection point.
       This method returns swath edges for either the ascending or descending
       pass, not both.  If both of our candidates are on the same edge, then
       one is seen on the ascending pass and the other on the descending. We
       solve this problem by moving one edge point to the inflection longitude.
       Note: the rest of this method will produce the same results regardless
             of the value of "ascending."  Factoring this out could slightly
             improve performance.
      xprod = target * cross
      north_hemisphere = coord.phi > 0
      if candidates.size == 2
        inf_theta = north_hemisphere ? inf.theta + PI : inf.theta
        inflection = Coordinate.phi_theta(coord.phi, inf_theta)
        if north_hemisphere && retrograde? == ascending ||
            !north_hemisphere && retrograde? != ascending
          west_edge = inflection
        else
          east_edge = inflection
        end
      end
       Return the edges
      [west_edge, east_edge]
    end
def horizontal_swath_edges ( coord , ascending ) candidates = [ ] if coord . phi . abs <= EPSILON return [ Coordinate . phi_theta ( 0 , coord . theta - swath_width_rad / 2 ) , Coordinate . phi_theta ( 0 , coord . theta + swath_width_rad / 2 ) ] end z = coord . z target_phi = coord . phi target_phi = max_orbit_phi if target_phi > max_orbit_phi target_phi = min_orbit_phi if target_phi < min_orbit_phi target = Coordinate . phi_theta ( target_phi , coord . theta ) cross = equator_crossing ( target , ascending , false ) a = target . y * cross . z - target . z * cross . y b = target . z * cross . x - target . x * cross . z c = target . x * cross . y - target . y * cross . x inf = Coordinate . phi_theta ( - inclination_rad , cross . theta - PI / 2 ) [ Coordinate . phi_theta ( inf . phi + swath_width_rad / 2 , inf . theta ) , Coordinate . phi_theta ( inf . phi - swath_width_rad / 2 , inf . theta ) ] . each do | p | d = a * p . x + b * p . y + c * p . z if b . abs < 0.000001 if a . abs < 0.000001 if coord . phi . abs < swath_width_rad / 2 return [ Coordinate . phi_theta ( 0 , - PI ) , Coordinate . phi_theta ( 0 , PI ) ] else return [ ] end end x = ( d - c * z ) / a dis = 1 - z ** 2 - x ** 2 if dis >= 0 y = sqrt ( dis ) candidates << Coordinate . xyz ( x , y , z ) candidates << Coordinate . xyz ( x , - y , z ) end else m = - a / b k = ( d - c * z ) / b r = sqrt ( 1 - z ** 2 ) a_q = m ** 2 + 1 b_q = 2 * m * k c_q = k ** 2 - r ** 2 discriminant = b_q ** 2 - 4 * a_q * c_q if discriminant >= 0 [ ( - b_q + sqrt ( discriminant ) ) / ( 2 * a_q ) , ( - b_q - sqrt ( discriminant ) ) / ( 2 * a_q ) ] . each do | x2 | y = m * x2 + k candidates << Coordinate . xyz ( x2 , y , z ) end end end end return [ ] if candidates . empty? candidates . sort_by! ( & :theta ) east_index = candidates . find_index { | c2 | coord . theta < c2 . theta } east_index = 0 if east_index . nil? west_index = ( east_index - 1 ) % candidates . size west_edge = candidates [ west_index ] east_edge = candidates [ east_index ] xprod = target * cross north_hemisphere = coord . phi > 0 if candidates . size == 2 inf_theta = north_hemisphere ? inf . theta + PI : inf . theta inflection = Coordinate . phi_theta ( coord . phi , inf_theta ) if north_hemisphere && retrograde? == ascending || ! north_hemisphere && retrograde? != ascending west_edge = inflection else east_edge = inflection end end [ west_edge , east_edge ] end
horizontal_swath_edges
Given coord, a coordinate within the orbital swath, returns an array
of 2 points with the same latitude as coord which lie on the swath edges.
horizontal_swath_edges Given coord a coordinate within the orbital swath returns an array of 2 points with the same latitude as coord which lie on the swath edges .
# Candidate solutions.  We may end up finding more points than we need. # Before returning, we'll take the two points closest to the given # coordinate # We need to be able to assume the coordinate isn't on the equator, so # we'll return early if we find that it is.  We know the swath width # from the orbit parameters.  This calculation is slightly off, since # the swath width is measured perpendicular to the orbit, but the # Java code makes this assumption.  This case should almost never # happen.  When it does, the numbers are likely close enough. # We ignore the rotation of the Earth for the purposes of this method. # The values we return (the swath edges) are not influenced by the # Earth's rotation.  Intermediate values are, but they're not visible # outside of this method. # # Ignoring the rotation of the Earth, we can imagine three circles # traced by the orbit.  The first is the path on the ground directly # beneath the orbit.  This is a great circle whose center is the origin. # The other two are the paths traced by the swath edges, which are not # great circles and not centered at the origin. # # These three circles each lie in a plane, and the three planes are # parallel to each other. # # There's another circle we're interested in, which is the circle # traced by the latitude line at coord.  This circle lies in its own # plane which is not parallel to the other three. # # The four original circles are the intersections of their respective # planes and the sphere of the Earth. # # So, to find the swath's edge coordinates, we need to find the # intersection between three things: # #   1. `x^2 + y^2 + z^2 = 1` # #      The sphere of the Earth.  We've scaled everything so the Earth #      has radius 1, and it's centered on the origin.  Convenient. # #   2. `ax + by + cz = d` # #      The swath's plane.  `a`, `b`, `c` and `d` are constants which #      we'll calculate.  Since they're parallel, the swath edges share #      `a`, `b`, and `c` with the orbital plane, but each edge has a #      different `d` value. # #   3. `z = "coord.z"` # #      The plane slicing through latitude = coord.lat # We already know z # We'll find `a`, `b`, and `c` for the orbit's plane, which are the same # constants for the swath's plane.  To find them, we need 3 non-colinear # points in the orbit's plane.  The origin is one, but we need two more. # First, we'll get a point that the orbit can pass over which is close # to coord.  Usually we'll just use coord, but if coord's latitude is # above the orbit's inflection point, we need to find a point with a # lower latitude # For the second point, we'll use the equatorial crossing of the orbit # through the target.  This is safe, because we bailed earlier if # the target was on the equator.  The "false" parameter tells the method # to ignore the Earth's rotation. # For a plane passing through the origin, `ax + by + cz = 0`.  With our # two other points, we can solve for `a`, `b`, and `c`: # If we can find a point on the swath edge, we can calculate `d`.  We # can find a point on each edge by finding the orbit's inflection point # and adding/subtracting half the swath width from its latitude. # Calculate solutions for the point p on each swath edge # Calculate `d` from the plane equation # The fact that `z` is constant is very handy.  It means we really only # care about two dimensions.  We're really trying to find the intersection # of a line (the edge's plane when `z = "coord.z"`) and a circle (the part # of the Earth intersecting `z = "coord.z"`). # # Note: that simplified math could make it easier to switch this # algorithm to use an ellipsoidal Earth (WGS84) in the future, since we # would just be intersecting a line with a slightly different circle. # # The swath plane's equation is: # # `ax + by + cz = d` # # `a`, `b`, `c`, `d`, and `z` are all known at this point, so we just # need to figure out `x` and `y`.  Solve the equation for `y` # # `y = -a/b x + (d - cz) / b` # # Which is the line we need to intersect with the latitude circle. # We can't just compute constants for the line, because `b` could be near 0. # Let's deal with that special case first. # If `b = 0`, then the plane's equation is: # # `ax + 0y + cz = d` # # Solving for `x`: # # `x = (d - cz) / a` # Ok, another problem.  It's possible that both `a` and `b` are 0.  In that # case, we have an orbit that follows the equator.  If coord lies within the # swath, then every orbit will cover the point, otherwise no orbits will. # Now we can handle the case where `b` is 0 and `a` isn't # We now know `x` and `z`, so we can use the sphere equation to solve for `y`. # # `y = +- sqrt(1 - z^2 - x^2)` # # In this case, `1 - z^2 - x^2` could be negative.  This would indicate # that the latitude does not intersect the swath edge.  Perhaps the point is # completely above one swath edge but below the other.  We won't worry # about duplicate solutions (`y = +- 0`), since the LongitudeCoverage # class will combine them correctly. # Now the case when b is not 0. # Recall, the line's equation is: # # `y = -a/b x + (d - cz) / b` # # We'll call those constants `m` and `k` # Which means # # `y = mx + k` # # (We won't use the standard `b` as the `y`-intercept name, since # it means something different in the plane equation). # We need to find the intersection of that line with the latitude # circle at `z = "coord.z"`.  Using the sphere equation with a # constant z, we get the latitude circle's equation: # #  `x^2 + y^2 = 1 - z^2` # It's centered at the origin.  `1 - z^2` is the radius squared. # # Taking the square root is safe because we calculated z from a # lat/lon coordinate.  Further, a negative radius doesn't make sense. # So, # # `x^2 + y^2 = r^2` # # Substituting `y` with the value from our line equation: # # `x^2 + (mx + k)^2 = r^2` # # Multiplying out, we get: # # `x^2 + m^2x^2 + 2mkx + k^2 = r^2` # Grouping constants, we get: # # `(m^2 + 1)x^2 + (2mk)x + (k^2 - r^2) = 0` # # We call those constants `a_q`, `b_q`, and `c_q` # So: # # `a_q x^2 + b_qx + c_q = 0` # We can solve this using the quadratic equation: # # `x = (-b +- sqrt(b^2 - 4ac)) / (2a)` # The discriminant could be negative if the line and circle do # not intersect # There are two possible `x` values # Now we need `y` # Return if we haven't found any valid points # At this point, we have either 2 or 4 candidate points that are swath # edge intersections with the target latitude.  Now we need to find # the two that are closest to the coordinate the caller passed.  We # want them in west-to-east order to make things simpler for the caller. # Find the index of the first candidate to the east of the target # Handle wrapping around the date line # The west index is the one to the left of the east index, wrapping around # the date line if needed # There is a special case here.  If we only found two candidates, that # indicates that the point is above one swath edge, but below the other, # which happens at latitudes near the inflection point. # # This method returns swath edges for either the ascending or descending # pass, not both.  If both of our candidates are on the same edge, then # one is seen on the ascending pass and the other on the descending. We # solve this problem by moving one edge point to the inflection longitude. # # Note: the rest of this method will produce the same results regardless #       of the value of "ascending."  Factoring this out could slightly #       improve performance. # Return the edges
coord
ascending