衛星などの上空にある物体の地上からの可視範囲のグラフ化

visible_range.py

# -*- coding: utf-8 -*-
# 衛星高度と可視範囲のグラフ化
#
# cf. 半揚 稔雄(2014) 「ミッション解析と軌道設計の基礎」
# Copyright (c) 2016 Takahiro Inagawa
# Released under the MIT license
import sys
reload(sys)
import platform
sys.setdefaultencoding('utf-8')

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.font_manager
from matplotlib.font_manager import FontProperties

if 'Windows' == platform.system():
    font_path = r'C:\Windows\Fonts\ipaexg.ttf'

if 'Darwin' == platform.system(): # for Mac
    font_path = '/Library/Fonts/Osaka.ttf'

font_prop = FontProperties(fname=font_path)
mpl.rcParams['font.family'] = font_prop.get_name()
plt.close('all')
plt.style.use('ggplot')
mpl.rcParams['axes.grid'] = True

# plt.ion()

def radius_visible(altitude, invalid_angle_deg = 3):
    # 衛星やロケットなどの対象の高度と観測者側の無効角度を入力して可視範囲の半径を計算
    # return 可視半径 [m]
    re = 6378137 # 地球半径 m
    epsilon = np.deg2rad(invalid_angle_deg)
    phi = np.arccos(re/(re+altitude) * np.cos(epsilon)) - epsilon
    return phi * re

if __name__ == '__main__':
    alt = np.linspace(0,1500000)
    for r in [0, 3, 10, 15, 30, 45, 60]:
        radius = radius_visible(alt, r)
        plt.plot(alt/1000, radius/1000, label=u"仰角 = %d deg" % (r))
    plt.title(u"衛星高度と衛星可視範囲円半径")
    plt.xlabel(u"高度 km")
    plt.ylabel(u"可視範囲の円半径 km")
    plt.legend(loc="best")
    plt.ylim([0,4000])
    plt.xlim([min(alt/1000), max(alt/1000)])
    plt.vlines(x=[400, 1150], ymin = 0, ymax = 5000, colors="b", linestyles="dashed")