apparent motion of a satellite¶

In [1]:

import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import astropy.units as u
import astropy.constants as c
import numpy as np

import matplotlib.pyplot as plt from matplotlib.patches import Circle import astropy.units as u import astropy.constants as c import numpy as np

assumption

• satellite oribting around the Earth at the altitude of $h$.
• $\theta$ is the orbital phase of the satellite.
• $\phi$ is the viewing angle of the satellite from the observatory.
• $\theta = 0$ when the satellite is at the zenith of the observatory.
• Define $\phi = 0$ when $\theta = 0$

In [2]:

v = lambda h: \
    np.sqrt(c.G * c.M_earth / (c.R_earth + h))
v(600 * u.km).to('km/s')

v = lambda h: \ np.sqrt(c.G * c.M_earth / (c.R_earth + h)) v(600 * u.km).to('km/s')

Out[2]:

$7.5578848 \; \mathrm{\frac{km}{s}}$

In [3]:

fig, ax = plt.subplots(figsize=(9, 5))


h = np.linspace(550, 35000, 201) * u.km
l1, = ax.semilogx(h, v(h).to('km/s'), label='Velocity')
ax.set_xlabel('Height (km)', fontsize=16)
ax.set_ylabel('Velocity (km/s)', fontsize=16)
ax.set_ylim([0, 8])

ax2 = ax.twinx()
l2, = ax2.plot(h, (2 * np.pi * (c.R_earth + h)/v(h)).to('hour'),
  color='C1', label='Period')
ax2.set_ylabel('Oribital Period (hour)', fontsize=16)
ax2.set_ylim([0, 28])

ax.legend(handles=[l1, l2], loc='lower right', fontsize=14, frameon=False)


fig.tight_layout()
# fig.savefig('../satellite_velocity.png', dpi=144)
plt.show()

fig, ax = plt.subplots(figsize=(9, 5)) h = np.linspace(550, 35000, 201) * u.km l1, = ax.semilogx(h, v(h).to('km/s'), label='Velocity') ax.set_xlabel('Height (km)', fontsize=16) ax.set_ylabel('Velocity (km/s)', fontsize=16) ax.set_ylim([0, 8]) ax2 = ax.twinx() l2, = ax2.plot(h, (2 * np.pi * (c.R_earth + h)/v(h)).to('hour'), color='C1', label='Period') ax2.set_ylabel('Oribital Period (hour)', fontsize=16) ax2.set_ylim([0, 28]) ax.legend(handles=[l1, l2], loc='lower right', fontsize=14, frameon=False) fig.tight_layout() # fig.savefig('../satellite_velocity.png', dpi=144) plt.show()

In [4]:

d = lambda h, theta: \
    np.sqrt(c.R_earth**2 + (c.R_earth + h)**2 \
    - 2 * c.R_earth * (c.R_earth + h) * np.cos(theta))
sinphi = lambda h, theta: \
    (c.R_earth + h)/d(h, theta) * np.sin(theta)
cosphi = lambda h, theta: \
    ((c.R_earth + h)**2 - c.R_earth**2 - d(h, theta)**2) \
    / 2 / c.R_earth / d(h, theta)

d = lambda h, theta: \ np.sqrt(c.R_earth**2 + (c.R_earth + h)**2 \ - 2 * c.R_earth * (c.R_earth + h) * np.cos(theta)) sinphi = lambda h, theta: \ (c.R_earth + h)/d(h, theta) * np.sin(theta) cosphi = lambda h, theta: \ ((c.R_earth + h)**2 - c.R_earth**2 - d(h, theta)**2) \ / 2 / c.R_earth / d(h, theta)

In [5]:

dphi = lambda h, theta: \
    v(h) / d(h, theta) / cosphi(h, theta) \
    * (np.cos(theta) - c.R_earth / d(h, theta) \
    * sinphi(h, theta) * np.sin(theta)) * u.rad

phi = lambda h, theta: \
    np.sign(theta) * np.arccos(cosphi(h, theta))

dphi = lambda h, theta: \ v(h) / d(h, theta) / cosphi(h, theta) \ * (np.cos(theta) - c.R_earth / d(h, theta) \ * sinphi(h, theta) * np.sin(theta)) * u.rad phi = lambda h, theta: \ np.sign(theta) * np.arccos(cosphi(h, theta))

In [6]:

theta = np.linspace(-60, 60, 21) * u.deg
earth = Circle([0, 0], radius=1, fill=False, color='green')

h = 1000 * u.km
x = d(h, theta) * sinphi(h, theta) / c.R_earth
y = d(h, theta) * cosphi(h, theta) / c.R_earth + 1

fig, ax = plt.subplots(figsize=(9, 4))

ax.scatter([0], [1], s=200, marker='x', color='C1')
ax.plot(x, y, marker='o', label=f'h = {h}')
ax.add_patch(earth)
ax.set_ylim([0.45, 1.35])
ax.set_aspect(1.0)
ax.set_xlabel(r'Earth Coordinate X ($R_{\oplus}$)', fontsize=16)
ax.set_ylabel(r'Earth Coordinate Y ($R_{\oplus}$)', fontsize=16)
ax.legend()
ax.grid()
fig.tight_layout()

# fig.savefig('../satellite_orbit.png', dpi=144)
plt.show()

theta = np.linspace(-60, 60, 21) * u.deg earth = Circle([0, 0], radius=1, fill=False, color='green') h = 1000 * u.km x = d(h, theta) * sinphi(h, theta) / c.R_earth y = d(h, theta) * cosphi(h, theta) / c.R_earth + 1 fig, ax = plt.subplots(figsize=(9, 4)) ax.scatter([0], [1], s=200, marker='x', color='C1') ax.plot(x, y, marker='o', label=f'h = {h}') ax.add_patch(earth) ax.set_ylim([0.45, 1.35]) ax.set_aspect(1.0) ax.set_xlabel(r'Earth Coordinate X ($R_{\oplus}$)', fontsize=16) ax.set_ylabel(r'Earth Coordinate Y ($R_{\oplus}$)', fontsize=16) ax.legend() ax.grid() fig.tight_layout() # fig.savefig('../satellite_orbit.png', dpi=144) plt.show()

In [12]:

theta = np.linspace(-30, 30, 1001) * u.deg
elev = lambda h, theta: 90 * u.deg - np.abs(phi(h, theta))

fig, ax = plt.subplots(figsize=(9, 5))

h = 550 * u.km
ax.plot(theta, dphi(h, theta).to('arcmin/s'), label=f'h = {h}')
h = 1000 * u.km
ax.plot(theta, dphi(h, theta).to('arcmin/s'), label=f'h = {h}')
h = 1500 * u.km
ax.plot(theta, dphi(h, theta).to('arcmin/s'), label=f'h = {h}')

ax.set_ylim([0, None])
ax.set_xlabel('Orbital Phase (deg)', fontsize=16)
ax.set_ylabel('Apparent Motion (arcmin/sec)', fontsize=16)
ax.legend()
ax.grid()
fig.tight_layout()

# fig.savefig('../velocity_orbital_phase.png', dpi=144)
plt.show()

theta = np.linspace(-30, 30, 1001) * u.deg elev = lambda h, theta: 90 * u.deg - np.abs(phi(h, theta)) fig, ax = plt.subplots(figsize=(9, 5)) h = 550 * u.km ax.plot(theta, dphi(h, theta).to('arcmin/s'), label=f'h = {h}') h = 1000 * u.km ax.plot(theta, dphi(h, theta).to('arcmin/s'), label=f'h = {h}') h = 1500 * u.km ax.plot(theta, dphi(h, theta).to('arcmin/s'), label=f'h = {h}') ax.set_ylim([0, None]) ax.set_xlabel('Orbital Phase (deg)', fontsize=16) ax.set_ylabel('Apparent Motion (arcmin/sec)', fontsize=16) ax.legend() ax.grid() fig.tight_layout() # fig.savefig('../velocity_orbital_phase.png', dpi=144) plt.show()

In [8]:

theta = np.linspace(0, 25, 1001) * u.deg
elev = lambda h, theta: 90 * u.deg - np.abs(phi(h, theta))

fig, ax = plt.subplots(figsize=(9, 5))

h = 550 * u.km
ax.plot(elev(h, theta), dphi(h, theta).to('arcmin/s'), label=f'h = {h}')
h = 1000 * u.km
ax.plot(elev(h, theta), dphi(h, theta).to('arcmin/s'), label=f'h = {h}')
h = 1500 * u.km
ax.plot(elev(h, theta), dphi(h, theta).to('arcmin/s'), label=f'h = {h}')

ax.set_ylim([0, None])
ax.set_xlabel('Elevation (deg)', fontsize=16)
ax.set_ylabel('Apparent Motion (arcmin/sec)', fontsize=16)
ax.legend()
ax.grid()
fig.tight_layout()

# fig.savefig('../velocity_elevation.png', dpi=144)
plt.show()

theta = np.linspace(0, 25, 1001) * u.deg elev = lambda h, theta: 90 * u.deg - np.abs(phi(h, theta)) fig, ax = plt.subplots(figsize=(9, 5)) h = 550 * u.km ax.plot(elev(h, theta), dphi(h, theta).to('arcmin/s'), label=f'h = {h}') h = 1000 * u.km ax.plot(elev(h, theta), dphi(h, theta).to('arcmin/s'), label=f'h = {h}') h = 1500 * u.km ax.plot(elev(h, theta), dphi(h, theta).to('arcmin/s'), label=f'h = {h}') ax.set_ylim([0, None]) ax.set_xlabel('Elevation (deg)', fontsize=16) ax.set_ylabel('Apparent Motion (arcmin/sec)', fontsize=16) ax.legend() ax.grid() fig.tight_layout() # fig.savefig('../velocity_elevation.png', dpi=144) plt.show()

In [9]:

theta = np.linspace(0, 25, 1001) * u.deg
elev = lambda h, theta: 90 * u.deg - np.abs(phi(h, theta))
pixel = 1.18 * u.arcsec

def conv(x):
  """Vectorized 1/x, treating x==0 manually"""
  x = np.array(x, float)
  near_zero = np.isclose(x, 0)
  x[near_zero] = 0
  x[~near_zero] = 1000 / x[~near_zero]
  return x

fig, ax = plt.subplots(figsize=(9, 5))

h = 550 * u.km
ax.plot(elev(h, theta), pixel / dphi(h, theta).to('arcmin/ms'), label=f'h = {h}')
h = 1000 * u.km
ax.plot(elev(h, theta), pixel / dphi(h, theta).to('arcmin/ms'), label=f'h = {h}')
h = 1500 * u.km
ax.plot(elev(h, theta), pixel / dphi(h, theta).to('arcmin/ms'), label=f'h = {h}')

ax.set_ylim([0, 350])
ax.set_xlabel('Elevation (deg)', fontsize=16)
ax.set_ylabel('Pixel Crossing Time (msec)', fontsize=16)
ax.legend()
ax.grid()
ax2 = ax.secondary_yaxis('right', functions=(conv, conv))
ax2.set_yticks([2, 3, 5, 10, 20, 50])
ax2.set_ylabel('Frame Rate (fps)', fontsize=16)
fig.tight_layout()

# fig.savefig('../pixel_time.png', dpi=144)
plt.show()

theta = np.linspace(0, 25, 1001) * u.deg elev = lambda h, theta: 90 * u.deg - np.abs(phi(h, theta)) pixel = 1.18 * u.arcsec def conv(x): """Vectorized 1/x, treating x==0 manually""" x = np.array(x, float) near_zero = np.isclose(x, 0) x[near_zero] = 0 x[~near_zero] = 1000 / x[~near_zero] return x fig, ax = plt.subplots(figsize=(9, 5)) h = 550 * u.km ax.plot(elev(h, theta), pixel / dphi(h, theta).to('arcmin/ms'), label=f'h = {h}') h = 1000 * u.km ax.plot(elev(h, theta), pixel / dphi(h, theta).to('arcmin/ms'), label=f'h = {h}') h = 1500 * u.km ax.plot(elev(h, theta), pixel / dphi(h, theta).to('arcmin/ms'), label=f'h = {h}') ax.set_ylim([0, 350]) ax.set_xlabel('Elevation (deg)', fontsize=16) ax.set_ylabel('Pixel Crossing Time (msec)', fontsize=16) ax.legend() ax.grid() ax2 = ax.secondary_yaxis('right', functions=(conv, conv)) ax2.set_yticks([2, 3, 5, 10, 20, 50]) ax2.set_ylabel('Frame Rate (fps)', fontsize=16) fig.tight_layout() # fig.savefig('../pixel_time.png', dpi=144) plt.show()