ロケット打上げ軌道投入時の必要ΔVの簡易計算

deltaV.py

"""
軌道投入への必要ΔVの簡易計算スクリプト
"""

import sys
import os
import configparser
import numpy as np
from numpy import sqrt, cos, sin, arcsin, tan, deg2rad
import matplotlib.pyplot as plt


class Earth:
    def __init__(self, setting_file, reload=False):
        if reload:  # 再読込の際はself.settingの値をそのまま使う
            pass
        else:
            self.setting_file = setting_file
            self.setting = configparser.ConfigParser()
            self.setting.optionxform = str  # 大文字小文字を区別するおまじない
            self.setting.read(setting_file, encoding="utf8")
        setting = self.setting

        self.perigee = setting.getfloat("軌道", "近地点高度[km]")
        self.apogee = setting.getfloat("軌道", "遠地点高度[km]")
        self.inclination = setting.getfloat("軌道", "軌道傾斜角[deg]")

        self.longitude = setting.getfloat("打上げ", "射点緯度[deg]")
        self.vel_loss = setting.getfloat("打上げ", "損失速度[km/s]")

        self.gravity_const = 3.986e05  # 地球重力定数 [km3/sec2]
        self.radius = 6378  # 地球半径 [km]
        self.omega = 7.292e-05  # 地球時点角速度 [rad/s]

        self.vel_initial = self.radius * self.omega * cos(deg2rad(self.longitude))

    def calc_deltaV_2_Hohmann_transfer_orbit(self):
        """
        地表面から2ホーマン軌道で投入軌道に入れる、⊿Vを計算する。
        1回目のインパルス噴射で地表からペリジ高度まで上げ、2回目でアポジまで増速
        self.deltaV1 (float) : 1回目⊿V(km/s)
        self.deltaV2 (float) : 2回め⊿V(km/s)
        self.deltaV (float) : 必要⊿V（合計⊿V）(km/s)
        self.vel_require (float) : トータルで必要な⊿V(km/s)
        """
        r_earth = self.radius
        r_apogee = self.apogee + self.radius
        r_perigee = self.perigee + self.radius
        ak = (r_earth + r_perigee) / 2
        a2 = (r_apogee + r_perigee) / 2
        vel_0 = self.vel_initial
        vel_1 = sqrt(self.gravity_const * (2 / r_earth - 1 / ak))
        vel_2k = sqrt(self.gravity_const * (2 / r_perigee - 1 / ak))
        vel_2 = sqrt(self.gravity_const * (2 / r_perigee - 1 / a2))
        inc = deg2rad(self.inclination)  # 軌道傾斜角 [rad]
        lon = deg2rad(self.longitude)  # 打上射場緯度 [rad]
        vel_vec0 = np.array([vel_0, 0, 0])

        if self.longitude < self.inclination:
            alpha = arcsin(tan(lon) / tan(inc))
            vel_vec1 = np.array(
                [
                    sin(inc) * sin(lon) * sin(alpha) + cos(inc) * cos(lon),
                    -sin(inc) * sin(lon) * cos(alpha),
                    sin(inc) * cos(lon) * cos(alpha),
                ]
            )
            vel_vec1 = vel_vec1 * vel_1
            self.deltaV1 = np.linalg.norm(vel_vec1 - vel_vec0)
            self.deltaV2 = vel_2 - vel_2k
            self.deltaV = self.deltaV1 + self.deltaV2
            self.vel_require = self.deltaV + self.vel_loss
        else:
            phi = deg2rad(self.longitude - self.inclination)  # 緯度と軌道傾斜角の差分 [rad]
            vel_vec1 = np.array([1, 0, 0]) * vel_1
            self.deltaV1 = np.linalg.norm(vel_vec1 - vel_vec0)
            self.deltaV2 = vel_2 - vel_2k + 2 * vel_2 * sin(phi / 2)
            self.deltaV = self.deltaV1 + self.deltaV2
            self.vel_require = self.deltaV + self.vel_loss

    def display(self):
        print("ホーマン軌道 ⊿V1 = %.3f km/s" % (self.deltaV1))
        print("ホーマン軌道 ⊿V2 = %.3f km/s" % (self.deltaV2))
        print("ホーマン軌道 ⊿V = %.3f km/s" % (self.deltaV))
        print("損失速度 = %.3f km/s" % (self.vel_loss))
        print("初期速度 = %.3f km/s" % (self.vel_initial))
        print("必要⊿V = %.3f km/s" % (self.vel_require))

    def change_setting_value(self, section, key, value):
        """設定ファイルの中身を変更してファイルを保存しなおす
        Args:
            sectiojn (str) : 設定ファイルの[]で囲まれたセクション
            key (str) : 設定ファイルのkey
            value (str or float) : 書き換える値（中でstringに変換）
        """
        self.setting.set(section, key, str(value))
        self.__init__(self.setting_file, reload=True)


if __name__ == "__main__":
    if len(sys.argv) == 1:
        setting_file = "setting.ini"
    else:
        setting_file = sys.argv[1]
        assert os.path.exists(setting_file), "ファイルが存在しません"
    earth = Earth(setting_file)
    earth.calc_deltaV_2_Hohmann_transfer_orbit()
    earth.display()

    def plot_deltaV_at_lon(longitude_deg, alttitude_km=500):
        inclination_deg = np.arange(0, 101, 0.5)
        deltaV = []
        for i in inclination_deg:
            earth.change_setting_value("軌道", "軌道傾斜角[deg]", i)
            earth.change_setting_value("打上げ", "射点緯度[deg]", longitude_deg)
            earth.change_setting_value("軌道", "近地点高度[km]", alttitude_km)
            earth.change_setting_value("軌道", "遠地点高度[km]", alttitude_km)
            earth.calc_deltaV_2_Hohmann_transfer_orbit()
            deltaV.append(earth.vel_require)
        return inclination_deg, deltaV

    plt.close("all")

    plt.figure()
    alt = 36000
    lon_array = [5, 28, 30, 42, 45]
    for i in lon_array:
        inc, delV = plot_deltaV_at_lon(i, alt)
        plt.plot(inc, delV, label="打上緯度: %d deg" % (i))
        plt.plot(inc[0], delV[0], "ro")
    plt.vlines(x=0, ymin=9, ymax=15.7, colors="black")
    plt.xlabel("軌道傾斜角 [deg]")
    plt.ylabel("必要ΔV [km/s]")
    plt.title("静止軌道への必要ΔV計算 高度%dkm (損失速度=%.1fkm/s)" % (alt, earth.vel_loss))
    plt.grid()
    plt.ylim(ymin=9)
    plt.yticks(np.arange(9, 16.1, 0.25))
    plt.legend(loc="best")

    plt.figure()
    alt = 500
    lon_array = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90]
    for i in lon_array:
        inc, delV = plot_deltaV_at_lon(i, alt)
        plt.plot(inc, delV, label="緯度: %d deg" % (i))
    plt.xlabel("軌道傾斜角 [deg]")
    plt.ylabel("必要ΔV [km/s]")
    plt.title("必要ΔV計算, 円軌道高度 %d km (損失速度=%.1fkm/s)" % (alt, earth.vel_loss))
    plt.grid()
    plt.ylim([9, 11])
    plt.legend(loc="best", fontsize=7)

    plt.figure()
    alt_array = [200, 500, 1000]
    lon = 42
    for i in alt_array:
        inc, delV = plot_deltaV_at_lon(lon, i)
        plt.plot(inc, delV, label="高度: %d km" % (i))
    plt.xlabel("軌道傾斜角 [deg]")
    plt.ylabel("必要ΔV [km/s]")
    plt.title("必要ΔV計算, 打上緯度 %d deg (損失速度=%.1fkm/s)" % (lon, earth.vel_loss))
    plt.grid()
    plt.ylim([9, 11])
    plt.legend(loc="best")

    plt.show()

setting.ini

[軌道]
近地点高度[km] = 561
遠地点高度[km] = 561
軌道傾斜角[deg] = 98.1

[打上げ]
射点緯度[deg] = 41
損失速度[km/s] = 1.7