interpolating    mapping    searching

The Backtrack Orbit Search Algorithm

Worked Example: Pittsburgh

The picture below is a relatively ordinary orbit for a polar orbiter. The ascending equatorial crossing is at 77.75E. The satellite has an inclination of 98.78 degrees (declination of 8.78 degrees) and a period of 101 minutes. The swath width is 1400 km.

We want to find all the orbits for which this sensor sees Pittsburgh on the ascending pass. Given the location of Pittsburgh (40.50, -80.22) the algorithm looks at what the swath would look like if the groundtrack of the satellite went right through Pittsburgh on the ascending pass.

1) Find the orbit with nadir crossing the given point.

sin(theta) = sin(90 - Lat_inf) / sin(90 - Lat_p )
sin(theta) =  sin(8.78) / sin(49.5) = 0.15264 / 0.76041 = 0.2007
theta = 11.5779

Lon_n = asin(sin(theta) * sin(Lat_p) / Sin(Lat _inf ))  + Lon_p
Lon_n = asin(0.2007 * sin(40.5) / sin(81.22))  - 80.22
Lon_n = asin(0.2007 * 0.64945 / 0.98828)  - 80.22
Lon_n = asin(0.13189)  - 80.22
Lon_n = 7.57883  - 80.22 = -72.64117

Then correct for rotation

sin(length) / sin(90) = sin(lat_p) / sin(lat_inf )  so
length = asin(sin(lat_p) / sin(lat_inf ))
length = asin(sin(40.5) / sin(81.22))
length = asin(0.64945 / 0.98828)
length = asin(0.65715)
length = 41.0830 degrees.
 

Which means the true equatorial crossing can be found by adjusting for the rotation of the Earth during the time it takes the satellite to
travel 41.0830 degrees.

RealLon_n = Lon_n + [ period * length* ROTATION_RATE / circumference]
RealLon_n = -72.64117 + [ 101 * 41.0830 * (15/60) / 360 ]
RealLon_n = -72.64117 + [ 101 * 41.0830 * (15/60) / 360 ]
RealLon_n = -72.64117 + 2.8815
RealLon_n = -69.75965

And if that's all we wanted we'd be done.  But we need to the crossings for all the orbits that cover Pittsburgh, which means we need to take the swath width into account, which means we need to go back to the static sphere.

2) Convert to Cartesian 3-space and find the orbital plane.

The two points (lat_p, lon_p) and (lat_n , lon_n) lie on the Great Circle that is the orbit on the staic sphere.  Combined with the origin they define a plane.

x = r * cos(lon) * cos(lat);
y = r * sin(lon) * cos(lat);
z = r * sin(lat);

The radius is just a constant multiplier so we'll use a radius of 1.

x_p = cos(lon_p) * cos(lat_p) = cos(-80.22) * cos(40.5) = 0.16986 * 0.76041 = 0.12916
y_p = sin(lon_p) * cos(lat_p) = sin(-80.22) * cos(40.5) = -0.98547 * 0.76041 = -0.74936
z_p = sin(lat_p) = sin(40.5) = 0.64945

x_n = cos(lon_n) * cos(lat_n) = cos(-72.64117) * cos(0.0) = 0.29836
y_n = sin(lon_n) * cos(lat_n) = sin(-72.64117) * cos(0.0) = -0.95445
z_n = sin(lat_n) = sin(0.0) = 0

The plane of the Great Circle is defined by:
(y_pz_n - y_nz_p)x - (x_p z_n - x_nz_p)y + (x_py_n - x_ny_p)z = 0

Let

a = (y_pz_n - y_nz_p) = (-0.74936 * 0) - (-0.95445 * 0.64945) = 0.61987
b =  -(x_pz_n - x_nz_p ) = -((0.12916 * 0) - (0.29836 * 0.64945)) = 0.19377
c = (x_py_n - x_ny_p) = (0.12916 * -0.95445) - (0.29836 * -0.74936) = -0.123277 + 0.22358 = 0.10030

And the plane of the orbit on a static Earth is: ax + by + cz = 0

Sanity check:

(0.61987 * 0.12916) + (0.19377 * -0.74936) + (0.10030 * 0.64945) =
0.08006 - 0.14520 + 0.06514 = 0

(0.61987 * 0.29836 + (0.19377 * -0.95445) + (0.10030 * 0) =
0.18494 - 0.18494 + 0 = 0

3) Find points on the swath edges:

We have the width in kilometers but we need to keep the units consistent.  For this example we'll use 111 km/degree but you may want to be more precise.

width = 1400 / 111 = 12.6126 degrees.

For the left edge of the swath we need to find the point (Lat_left , Lon_left) that is distance (width/2) from the point of interest (Lat_p , Lon_p) along the heading (180 + theta).  

lat_left =  asin((sin(lat_p)*cos(width/2)) + (cos(lat _p)*sin(width/2)*sin(180 + theta)) )
lat_left =  asin((sin(40.5)*cos(6.3063)) + (cos(40.5)*sin(6.3063)*sin(191.5779)) )
lat_left =  asin((0.64945 * 0.99395) + (0.7604 * 0.10984 * -0.2007) )
lat_left =  asin(0.6455 + (-0.01676) )
lat_left =  asin(0.62874) = 38.957

lon_left = lon_p - acos( (cos(width/2) - sin(lat _p) * sin(lat_left)) / (cos(lat_p)*cos(lat _left )) )
lon_left = -80.22 - acos( (cos(6.3063) - sin(40.5) * sin(38.957)) / (cos(40.5)*cos(38.957)) )
lon_left = -80.22 - acos( (0.99395 - 0.64945 * 0.62874) / (0.7604*0.77762) )
lon_left = -80.22 - acos( (0.99395- 0.40834) / (0.5913) )
lon_left = -80.22 - acos( (0.58561481) / (0.5913) )
lon_left = -80.22 - acos(0.990384)
lon_left = -80.22 - 7.952
lon_left = -88.1721
 

For the right edge of the swath we need to find the point (lat_right , lon_right) that is distance (width/2) from the point of interest (lat_p , lon_p) along the heading (theta).

lat_right =  asin( (sin(lat_p)*cos(width/2)) + (cos(lat_p)*sin(width/2)*sin(theta))
lat_right =  asin( (sin(40.5)*cos(6.3063)) + (cos(40.5)*sin(6.3063)*sin(11.5779))
lat_right =  asin( 0.64945 * 0.99395) + (0.7604 * 0.10984 * 0.2007)
lat_right =  asin(0.6455 + 0.01676)
lat_right =  asin(0.66226) = 41.472

lon_right = lon_p + acos( (cos(width/2) - sin(lat _p)*sin(lat_right)) / (cos(lat_p)*cos(lat _right )) )
lon_right = -80.22 + acos( (cos(6.3063) - sin(40.5)*sin(41.472)) / (cos(40.5)*cos(41.472)) )
lon_right = -80.22 + acos( (0.99395 - 0.64945 * 0.66226) / (0.7604 * 0.74927) )
lon_right = -80.22 + acos( (0.99395 - 0.4301) / (0.7604 * 0.74927) )
lon_right = -80.22 + acos( 0.56384524 / 0.5697)
lon_right = -80.22 + acos( 0.98964)
lon_right = -80.22 + 8.253 = -71.967

4) Find the planes of the swath edges:

Convert the edge points to Cartesian and then the planes of the swath edges are given by ax + by + cz = d.

x_left = cos(lon_left) * cos(lat_left) = 0.031887 * 0.77762 = 0.02479
y_left = sin(lon_left) * cos(lat_left) = -0.99949 * 0.77762 = -0.77722
z_left = sin(lat_left) = 0.62874

x_right = cos(lon_right) * cos(lat_right ) = 0.30956 * 0.74928 = 0.23195
y_right = sin(lon_right) * cos(lat_right ) = -0.95088 * 0.74928 = -0.71247
z_right = sin(lat_right) = 0.66225
 
Perform a sanity check.  The two points we found should be the same distance from the orbital plane.

ax_left +by_left + cz_left = d = -(ax _right +by_right + cz_right)
(0.61987 * 0.02479) + (0.19377 * -0.77722) + (0.10030 * 0.62874) = d = -((0.61987 * 0.23195) + (0.19377 * -0.71247) + (0.10030 * 0.66225))
(0.015367 - 0.15060 + 0.06306) = d = -(0.14378 - 0.138055 + 0.06642)
-0.07217 = d = -0.07214    

There's some cumulative rounding error, but it's pretty good.

5) Find the intersects with the given latitude:

Now we have those planes we can find the points where the planes intersect the latitude of our point of interest.  That is:
ax + by + cz = d    : swath edge
x² +y² + z² = r²     : sphere
z = z_p                    : lat = lat_p

Three equations, three unknowns.  We're keeping the latitude constant so we only care about the longitude which is atan2(y, x).  Substituting z = z_p and rearranging we get:

ax + by + cz_p -d = 0  and
x² -y² + z_p² - r ² = 0

Let

planar_const = cz_p -d = (0.10030 * 0.64945 ) - 0.07215 = 0.0651398 - 0.07215 = -0.007010
sphere_const = z_p² - r² = 0.64945 ² - 1 = 0.42179 - 1 = -0.57821
scale = (a² + b²) = (0.61987² + 0.19377 ²) = 0.42179

radical = SQRT((-a⁴ * sphere_const) - (a² * b ² * sphere_const) - (a² * planar_const²))
radical = SQRT((-0.61987⁴ * -0.57821) - (0.61987² * 0.19377² * -0.57821) - (0.61987² * 0.007010 ² ))
radical = SQRT((0.147639 * -0.57821) - (0.38424 * 0.0375468 * -0.57821) - (0.384239 * 0.00004914))
radical = SQRT(0.0853666 + 0.0083418 - 0.00001888)
radical = SQRT(0.093689) = 0.3060

firstXterm = (b² * planar_const)/scale = (0.19377² * -0.007010) / 0.42179 = -0.000624
secondXterm = ((b*radical)/scale) = (0.19377 * 0.30609) / 0.42179 = 0.14062

So

x = (-planar_const + firstXterm - secondXterm) / a = (0.007010 + -0.000624 - 0.14062) / 0.61987 = -0.22555
y = (-(b*planar_const) + rad) / scale = (-(0.19377 * -0.007010) + 0.30609) / 0.42179 = 0.7289
lon = atan2(y, x) = atan2(0.7289, -0.22555) = atan(-3.2316) + 180 = -72.806 +180 = 107.194

OR              

x = (-planar_const + firstXterm + secondXterm) / a = (0.007010 + -0.000624 + 0.14062) / 0.61987 = 0.237156
y = (-(b*planar_const) - rad) / scale = (-(0.19377 * -0.007010) - 0.30609) / 0.42179 = -0.72247
lon = atan2(y, x) = atan2(-0.72247, 0.237156) = atan(-3.0464) = -71.827

On the other side of the orbit let:

planar_const = cz_p +d = (0.10030 * 0.64945 ) + 0.0727 = 0.06514 + 0.0727 = 0.13784

radical = SQRT((-a⁴ * sphere_const) - (a² * b ² * sphere_const) - (a² * planar_const²))
radical = SQRT((-0.61987⁴ * -0.57821) - (0.61987² * 0.19377² * -0.57821) - (0.61987² * 0.13784 ² ))
radical = SQRT((0.147639 * -0.57821) - (0.38424 * 0.0375468 * -0.57821) - (0.384239 * 0.0189998))
radical = SQRT(0.0853666 + 0.0083418 - 0.0073005)
radical = SQRT(0.086408) = 0.29395

firstXterm = (b² * planar_const)/scale = (0.19377² * 0.13784) / 0.42179 = 0.01227
secondXterm = ((b*radical)/scale) = (0.19377 * 0.29395) / 0.42179 = 0.13504

So

x = (-planar_const + firstXterm - secondXterm) / a = (-0.13784 + 0.01227 - 0.13504) / 0.61987 = -0.42043
y = (-(b*planar_const) + radical) / scale = (-(0.19377 * 0.13784) + 0.29395) / 0.42179 = 0.633587
lon = atan2(y, x) = atan2(0.633587, -0.42043) = atan(-1.506998) + 180 = -56.43 + 180 =  123.567

OR              

x = (-planar_const + firstX + secondX) / a = (-0.13784 + 0.01227 + 0.13504) / 0.61987 = 0.015277396
y = (-(b*planar_const) - radical) / scale = (-(0.19377 * 0.13784) - 0.29395) / 0.42179 =  -0.76023
lon = atan2(y, x) = atan2(-0.) = atan(-49.762) = -88.849

The points we're looking for are the closest two.  The other two are where the edge planes intersect the latitude plane on the other side of the sphere.  The closest points are (40.5, -88.85) to the West and (40.5, -71.83) to the East. The sensor sees those points on the edge of the swath sometime prior to, or after, it sees the point of interest and a correction for the rotation of the Earth could be introduced here for increased accuracy.  But the error is small so we skip that step here.

6) Find crossings for the extreme orbits:

The longitudinal distance from the point of interest (40.5, -80.22) to the swath edges is 8.63 degrees to the west and 8.39 degrees to the east.  So if the statellite were to pass over the point 8.63 degrees east of Pittsburgh (40.5, -71.59) the sensor would see Pittsburgh on the west edge of the swath.  And if the statellite were to pass over the point 8.39 degrees west of Pittsburgh (40.5, -88.61) the sensor would see Pittsburgh on the east edge of the swath.  We have only to find the equatorial crossings for those two orbits and we'll know the range of crossings for which the sensor is guaranteed to see Pittsburgh somewhere in the scan on the ascending pass.

To find those crossings we use the same method as above:

sin(theta) = sin(90 - Lat_inf) / sin(90 - Lat_p )
sin(theta) =  sin(8.78) / sin(49.5) = 0.15264 / 0.76041 = 0.2007
theta = 11.5779

Lon_ne = asin(sin(theta) * sin(Lat_e) / sin(Lat _inf ))  + Lon_e
Lon_ne = asin(0.2007 * sin(40.5) / sin(81.22))  - 71.59
Lon_ne = asin(0.2007 * 0.64945 / 0.98828)  - 71.59
Lon_ne = asin(0.13189)  - 71.59
Lon_ne = 7.57883  - 71.59 = -64.01

RealLon_ne = Lon_ne + 2.88 = -64.01 + 2.88 = -61.13

Lon_nw = asin(sin(theta) * sin(Lat_w) / sin(Lat _inf ))  + Lon_w
Lon_nw = asin(0.2007 * sin(40.5) / sin(81.22))  - 88.61
Lon_nw = asin(0.2007 * 0.64945 / 0.98828)  - 88.61
Lon_nw = asin(0.13189)  - 88.61
Lon_nw = 7.57883  - 88.61 = -81.03

RealLon_nw = Lon_nw + 2.88 = -81.03 + 2.88 = -78.14

And we're done. For any orbit with an ascending crossing in the range [-78.14, -61,1] the sensor will see Pittsburgh on the ascending pass.

Special Cases: 

The radius of the small circle that defines the swath edge is proportional to the radius of the sphere by the cos of the distance from the parallel great circle.  That is: radius_sc = radius * cos(distance)  In this example the swath width is 12.6126 degrees and we're using a radius of 1.  So radius_sc = 1 * cos(6.3063) = 0.99395.

The distance between two points in cartesian space is sqrt((x₁ -x₂)² + (y₁-y₂)² + (z₁-z₂)²) so the distance between left edge point and the west edge point is:
sqrt((x_left-x_west)² + (y_left -y_west)² + (z_left-z_west) ²) =
sqrt((0.02479 - 0.01528)² + (-0.77722 + 0.76023)² + (0.62874 - 0.64945)²) =
sqrt((0.00951)² + (-0.01699)² + (-0.02071)² ) =
sqrt(0.00009044 + 0.00028866 + 0.0004289) =
sqrt(0.000808) = 0.028425

Those two points, and the center of the small circle, create an icosceles triangle in the plane of the small circle.  We can use ordinary Euclidean trigonometry to find the vertex angle of that triangle and hence the length of the arc on the small circle.

cos(v) = (2r² - D²)/ 2r² = 1 - (D² /2r²) = 1 - (0.000808/1.97587) = 1 - 0.0004089 = 0.99959
v = acos(0.99959) = 1.63862 degrees

The correction would be the distance the earth rotated in the time is took the satellite to go 1.63862 degrees

101 * 1.63862 * (15/60) / 360 = 0.11493 degrees

Which would give us a west distance of 8.74 instead of 8.63

Similarly to the right/east the distance between the right edge point and the east edge point is:

sqrt((x_right-x_east)² + (y_right -y_east)² + (z_right-z_east) ²) =
sqrt((0.23195 - 0.237156)² + (-0.71247 + 0.72247)² + (0.66225- 0.64945)²) =
sqrt((0.005206)² + (0.01)² + (0.0128)²) =
sqrt(0.0000271 + 0.0001 + 0.00016384) =
sqrt(0.00029094) = 0.017057

cos(v) = (2r² - D²)/ 2r² = 1 - (D ²/2r²) = 1 - (0.00029094/1.97587) = 1 - 0.000147246 = 0.999853
v = acos(0.999853) = 0.9832 degrees

The correction would be the distance the earth rotated in the time is took the satellite to go 0.9832 degrees

101 * 0.9832 * (15/60) / 360 = 0.06896 degrees

Which would give us an east distance of  8.32 instead of 8.39