Last active
May 13, 2018 09:17
-
-
Save komasaru/e233f884fe7335de8178082d0a80315a to your computer and use it in GitHub Desktop.
Python script to calc ICRS coordinates with JPL binary data DE430.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #! /usr/local/bin/python3.6 | |
| """ | |
| JPLEPH(JPL の DE430 バイナリデータ)読み込み、座標（位置・速度）を計算 | |
| date name version | |
| 2018.04.05 mk-mode.com 1.00 新規作成 | |
| Copyright(C) 2018 mk-mode.com All Rights Reserved. | |
| ---------------------------------------------------------------------- | |
| * 引数 | |
| [第１] 対象天体番号（必須） | |
| 1: 水星 (Mercury) | |
| 2: 金星 (Venus) | |
| 3: 地球 (Earth) | |
| 4: 火星 (Mars) | |
| 5: 木星 (Jupiter) | |
| 6: 土星 (Saturn) | |
| 7: 天王星 (Uranus) | |
| 8: 海王星 (Neptune) | |
| 9: 冥王星 (Pluto) | |
| 10: 月 (Moon) | |
| 11: 太陽 (Sun) | |
| 12: 太陽系重心 (Solar system Barycenter) | |
| 13: 地球 - 月の重心 (Earth-Moon Barycenter) | |
| 14: 地球の章動 (Earth Nutations) | |
| 15: 月の秤動 (Lunar mantle Librations) | |
| [第２] 基準天体番号（必須。 0, 1 - 13） | |
| （ 0 は、対象天体番号が 14, 15 のときのみ） | |
| [第３] ユリウス日（省略可。省略時は現在日時をUTCとみなしたユリウス日） | |
| * 注意事項 | |
| - 求める座標は「赤道直角座標(ICRS)」 | |
| - 天体番号は、係数データの番号（並び順）と若干異なるので注意する。 | |
| （特に、天体番号 3, 10, と 12 以降） | |
| 係数データの並び順： | |
| 水星(1)、金星(2)、地球 - 月の重心(3)、火星(4)、木星(5)、土星(6)、 | |
| 天王星(7)、海王星(8)、冥王星(9)、月（地心）(10)、太陽(11)、 | |
| 地球の章動(12)、月の秤動(13) | |
| - 時刻系は「太陽系力学時(TDB)」である。（≒地球時(TT)） | |
| - 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度（6要素）、 | |
| 天体番号が 14 の場合は、黄経における章動 Δψ, 黄道傾斜における章動 Δε の | |
| 角位置・角速度（4要素）、 | |
| 天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度（6要素）。 | |
| - 対象天体番号 = 基準天体番号 は、無意味なので処理しない。 | |
| - 天体番号が 12 の場合は、 x, y, z の位置・速度の値は全て 0.0 とする。 | |
| - その他、JPL 提供の FORTRAN プログラム "testeph.f" を参考にした。 | |
| """ | |
| import datetime | |
| import struct | |
| import sys | |
| import traceback | |
| import lib_read_de430 | |
| class JplCalcDe430: | |
| USAGE = ( | |
| "【使用方法】 ./jpl_calc_de430.rb 対象天体番号 基準天体番号 [ユリウス日]\n" | |
| "【天体番号】（対象: 1 - 15, 基準: 0 - 13）\n" | |
| " 1: 水星, 2: 金星, 3: 地球, 4: 火星, 5: 木星,\n" | |
| " 6: 土星, 7: 天王星, 8: 海王星, 9: 冥王星, 10: 月,\n" | |
| " 11: 太陽, 12: 太陽系重心, 13: 地球 - 月の重心,\n" | |
| " 14: 地球の章動, 15: 月の秤動,\n" | |
| " 0: 対象天体番号 14, 15 の場合に指定\n" | |
| " ※対象天体番号≠基準天体番号" | |
| ) | |
| ASTRS = [ | |
| "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", | |
| "Neptune", "Pluto", "Moon", "Sun", "Solar system Barycenter", | |
| "Earth-Moon barycenter", "Earth Nutations", "Lunar mantle Librations" | |
| ] | |
| KIND = 2 # 計算区分（0: 計算しない、1: 位置のみ計算、2: 位置・速度を計算） | |
| # （但し、現時点では 1 の場合の処理は実装していない） | |
| BARY = True # 基準フラグ（True: 太陽系重心が基準, False: 太陽が基準） | |
| KM = False # 単位フラグ（True: km, km/sec, False: AU, AU/day） | |
| def __init__(self): | |
| # コマンドライン引数取得 | |
| self.__get_args() | |
| # DE430 読み込み | |
| de = lib_read_de430.JplReadDe430(self.jd) | |
| de.read_de430() | |
| self.sss = de.sss | |
| self.au = de.au | |
| self.emrat = de.emrat | |
| self.ipts = de.ipts | |
| self.jds = de.jds | |
| self.coeffs = de.coeffs | |
| # 各種初期化 | |
| self.lists = [0 for _ in range(12)] # 計算対象フラグ一覧配列 | |
| self.pvs = [[0.0 for _ in range(6)] for _ in range(11)] | |
| # 位置・角度データ配列 | |
| self.pvs_2 = [None for _ in range(13)] # 位置・角度データ配列（対象 - 基準 算出用） | |
| self.rrds = [0.0 for _ in range(6)] # 算出データ（対象 - 基準）配列 | |
| def exec(self): | |
| """ Execution """ | |
| try: | |
| self.__check_jd() # 引数のユリウス日をチェック | |
| self.__get_list() # 計算対象フラグ一覧（係数データの並びに対応）取得 | |
| self.__calc() # 各種計算 | |
| self.__display() # 結果出力 | |
| except Exception as e: | |
| raise | |
| def __get_args(self): | |
| """ コマンドライン引数取得 | |
| * 第１引数の値をインスタンス変数 self.astrs[0] に設定する | |
| * 第２引数の値をインスタンス変数 self.astrs[1] に設定する | |
| * 第３引数の値をインスタンス変数 self.jd に設定する | |
| （但し、第３引数が存在しない場合は現在日時を設定する） | |
| """ | |
| try: | |
| # 対象、基準天体番号 | |
| if len(sys.argv) < 3: | |
| print(self.USAGE) | |
| sys.exit(0) | |
| target, center = int(sys.argv[1]), int(sys.argv[2]) | |
| if ((target < 1 or target > 15) or | |
| (center < 0 or center > 13) or | |
| target == center or | |
| (target == 14 and center != 0) or | |
| (target == 15 and center != 0) or | |
| (target != 14 and target != 15 and center == 0)): | |
| print(self.USAGE) | |
| sys.exit(0) | |
| self.astrs = [target, center] | |
| # ユリウス日 | |
| if len(sys.argv) > 3: | |
| self.jd = float(sys.argv[3]) | |
| else: | |
| self.jd = self.__gc_to_jd(datetime.datetime.now()) | |
| except Exception as e: | |
| raise | |
| def __gc_to_jd(self, t): | |
| """ 年月日(グレゴリオ暦)からユリウス日(JD)を計算する | |
| * フリーゲルの公式を使用する | |
| JD = int(365.25 * year) | |
| + int(year / 400) | |
| - int(year / 100) | |
| + int(30.59 (month - 2)) | |
| + day | |
| + 1721088.5 | |
| * 上記の int(x) は厳密には、 x を超えない最大の整数 | |
| :param object t: datetime | |
| :return float jd: ユリウス日 | |
| """ | |
| year, month, day = t.year, t.month, t.day | |
| hour, minute, second = t.hour, t.minute, t.second | |
| try: | |
| # 1月,2月は前年の13月,14月とする | |
| if month < 3: | |
| year -= 1 | |
| month += 12 | |
| # 日付(整数)部分計算 | |
| jd = int(365.25 * year) + year // 400 - year // 100 \ | |
| + int(30.59 * (month - 2)) + day + 1721088.5 | |
| # 時間(小数)部分計算 | |
| t = (second / 3600 + minute / 60 + hour) / 24 | |
| return jd + t | |
| except Exception as e: | |
| raise | |
| def __check_jd(self): | |
| """ 引数のユリウス日をチェック """ | |
| try: | |
| if self.jd < self.sss[0] or self.jd >= self.sss[1]: | |
| print("Please input JD s.t. {} <= JD < {}.".format(*self.sss)) | |
| sys.exit(0) | |
| except Exception as e: | |
| raise | |
| def __get_list(self): | |
| """ 計算対象フラグ一覧取得 | |
| * チェビシェフ多項式による計算が必要な天体の一覧をインスタンス変数 | |
| （リスト） self.lists に設定する | |
| * 配列の並び順（係数データの並び順から「太陽」を除外した13個） | |
| self.lists = [ | |
| 水星, 金星, 地球 - 月の重心, 火星, 木星, 土星, 天王星, 海王星, | |
| 冥王星, 月（地心）, 地球の章動, 月の秤動 | |
| ] | |
| """ | |
| try: | |
| if self.astrs[0] == 14: | |
| if self.ipts[11][1] > 0: | |
| self.lists[10] = self.KIND | |
| return | |
| if self.astrs[0] == 15: | |
| if self.ipts[12][1] > 0: | |
| self.lists[11] = self.KIND | |
| for k in self.astrs: | |
| if k <= 10: | |
| self.lists[k - 1] = self.KIND | |
| if k == 10: | |
| self.lists[2] = self.KIND | |
| if k == 3: | |
| self.lists[9] = self.KIND | |
| if k == 13: | |
| self.lists[2] = self.KIND | |
| except Exception as e: | |
| raise | |
| def __calc(self): | |
| """ Calclulation """ | |
| try: | |
| # 補間（11:太陽） | |
| self.pv_sun = self.__interpolate(11) | |
| # 補間（1:水星〜10:月） | |
| for i in range(10): | |
| if self.lists[i] == 0: | |
| continue | |
| self.pvs[i] = self.__interpolate(i + 1) | |
| if i > 8: | |
| continue | |
| if self.BARY: | |
| continue | |
| self.pvs[i] = [a - b for a, b in zip(self.pvs[i], self.pv_sun)] | |
| # 補間（14:地球の章動） | |
| if self.lists[10] > 0 and self.ipts[11][1] > 0: | |
| p_nut = self.__interpolate(14) | |
| # 補間（15:月の秤動） | |
| if self.lists[11] > 0 and self.ipts[12][1] > 0: | |
| self.pvs[10] = self.__interpolate(15) | |
| # 対象天体と基準天体の差 | |
| if self.astrs[0] == 14: | |
| if self.ipts[11][1] > 0: | |
| self.rrds = p_nut + [0.0, 0.0] | |
| elif self.astrs[0] == 15: | |
| if self.ipts[12][1] > 0: | |
| self.rrds = self.pvs[10] | |
| else: | |
| for i in range(10): | |
| self.pvs_2[i] = self.pvs[i] | |
| if 11 in self.astrs: | |
| self.pvs_2[10] = self.pv_sun | |
| if 12 in self.astrs: | |
| self.pvs_2[11] = [0.0 for _ in range(6)] | |
| if 13 in self.astrs: | |
| self.pvs_2[12] = self.pvs[2] | |
| if (self.astrs[0] * self.astrs[1] == 30 or \ | |
| self.astrs[0] + self.astrs[1] == 13): | |
| self.pvs_2[2] = [0.0 for _ in range(6)] | |
| else: | |
| if self.lists[2] != 0: | |
| self.pvs_2[2] = [ | |
| a - b / (1.0 + self.emrat) | |
| for a, b in zip(self.pvs[2], self.pvs[9]) | |
| ] | |
| if self.lists[9] != 0: | |
| self.pvs_2[9] = [ | |
| a + b | |
| for a, b in zip(self.pvs_2[2], self.pvs[9]) | |
| ] | |
| for i in range(6): | |
| self.rrds[i] = self.pvs_2[self.astrs[0] - 1][i] \ | |
| - self.pvs_2[self.astrs[1] - 1][i] | |
| except Exception as e: | |
| raise | |
| def __interpolate(self, astr): | |
| """ 補間 | |
| * 使用するチェビシェフ多項式の係数は、 | |
| * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度（6要素）、 | |
| 天体番号が 14 の場合は、 Δψ, Δε の角位置・角速度（4要素）、 | |
| 天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度（6要素）。 | |
| * 天体番号が 12 の場合は、 x, y, z の位置・速度の値は全て 0.0 とする。 | |
| :param int astr: 天体番号 | |
| :return list pvs: [ | |
| x 位置, y 位置, z 位置, | |
| x 速度, y 速度, z 速度 | |
| ] | |
| 但し、 | |
| 14（地球の章動）の場合は、 | |
| [ | |
| Δψ の角位置, Δε の角位置, | |
| Δψ の角速度, Δε の角速度 | |
| ] | |
| 15（月の秤動）の場合は、 | |
| [ | |
| φ の角位置, θ の角位置, ψ の角位置, | |
| φ の角速度, θ の角速度, ψ の角速度 | |
| ] | |
| """ | |
| pvs = [] | |
| try: | |
| tc, idx_sub = self.__norm_time(astr) | |
| n_item = 2 if astr == 14 else 3 # 要素数 | |
| i_ipt = astr - 3 if astr > 13 else astr - 1 | |
| i_coef = astr - 3 if astr > 13 else astr - 1 | |
| # 位置 | |
| ps = [1, tc] | |
| for _ in range(2, self.ipts[i_ipt][1]): | |
| ps.append(2 * tc * ps[-1] - ps[-2]) | |
| for i in range(n_item): | |
| val = 0 | |
| for j in range(self.ipts[i_ipt][1]): | |
| val += self.coeffs[i_coef][idx_sub][i][j] * ps[j] | |
| if not(self.KM) and astr < 14: | |
| val /= self.au | |
| pvs.append(val) | |
| # 速度 | |
| vs = [0, 1, 2 * 2 * tc] | |
| for i in range(3, self.ipts[i_ipt][1]): | |
| vs.append(2 * tc * vs[-1] + 2 * ps[i - 1] - vs[-2]) | |
| for i in range(n_item): | |
| val = 0 | |
| for j in range(self.ipts[i_ipt][1]): | |
| val += self.coeffs[i_coef][idx_sub][i][j] \ | |
| * vs[j] * 2 * self.ipts[i_ipt][2] / self.sss[2] | |
| if astr < 14: | |
| val /= 86400 if self.KM else self.au | |
| pvs.append(val) | |
| return pvs | |
| except Exception as e: | |
| raise | |
| def __norm_time(self, astr): | |
| """ チェビシェフ多項式用に時刻を正規化、サブ区間のインデックス算出 | |
| :param int astr: 天体番号 | |
| :return list: [チェビシェフ時間, サブ区間のインデックス] | |
| """ | |
| try: | |
| idx = astr - 2 if astr > 13 else astr | |
| jd_start = self.jds[0] | |
| tc = (self.jd - jd_start) / self.sss[2] | |
| temp = tc * self.ipts[idx - 1][2] | |
| idx = int(temp - int(tc)) # サブ区間のインデックス | |
| tc = (temp % 1 + int(tc)) * 2 - 1 # チェビシェフ時間 | |
| return [tc, idx] | |
| except Exception as e: | |
| raise | |
| def __display(self): | |
| """ Display """ | |
| try: | |
| s = ( | |
| " Target: {:2d} ({})\n" | |
| " Center: {:2d}" | |
| ).format( | |
| self.astrs[0], self.ASTRS[self.astrs[0] - 1], | |
| self.astrs[1] | |
| ) | |
| if self.astrs[1] != 0: | |
| s += " ({})".format(self.ASTRS[self.astrs[1] - 1]) | |
| s += ( | |
| "\n" | |
| " JD: {:16.8f} day\n" | |
| " 1AU: {:11.1f} km\n" | |
| ).format(self.jd, self.au) | |
| print(s) | |
| s = "" | |
| if self.astrs[0] == 14: | |
| s = ( | |
| " Position(Δψ) = {:32.20f} rad\n" | |
| " Position(Δε) = {:32.20f} rad\n" | |
| " Velocity(Δψ) = {:32.20f} rad/day\n" | |
| " Velocity(Δε) = {:32.20f} rad/day" | |
| ).format( | |
| self.rrds[0], self.rrds[1], | |
| self.rrds[2], self.rrds[3] | |
| ) | |
| elif self.astrs[0] == 15: | |
| s = ( | |
| " Position(φ) = {:32.20f} rad\n" | |
| " Position(θ) = {:32.20f} rad\n" | |
| " Position(ψ) = {:32.20f} rad\n" | |
| " Velocity(φ) = {:32.20f} rad/day\n" | |
| " Velocity(θ) = {:32.20f} rad/day\n" | |
| " Velocity(ψ) = {:32.20f} rad/day" | |
| ).format( | |
| self.rrds[0], self.rrds[1], self.rrds[2], | |
| self.rrds[3], self.rrds[4], self.rrds[5] | |
| ) | |
| else: | |
| unit_0, unit_1 = ["km", "sec"] if self.KM else ["AU", "day"] | |
| s = ( | |
| " Position(x) = {:32.20f} {}\n" | |
| " Position(y) = {:32.20f} {}\n" | |
| " Position(z) = {:32.20f} {}\n" | |
| " Velocity(x) = {:32.20f} {}/{}\n" | |
| " Velocity(y) = {:32.20f} {}/{}\n" | |
| " Velocity(z) = {:32.20f} {}/{}" | |
| ).format( | |
| self.rrds[0], unit_0, | |
| self.rrds[1], unit_0, | |
| self.rrds[2], unit_0, | |
| self.rrds[3], unit_0, unit_1, | |
| self.rrds[4], unit_0, unit_1, | |
| self.rrds[5], unit_0, unit_1 | |
| ) | |
| print(s) | |
| except Exception as e: | |
| raise | |
| if __name__ == '__main__': | |
| try: | |
| obj = JplCalcDe430() | |
| obj.exec() | |
| except Exception as e: | |
| traceback.print_exc() | |
| sys.exit(1) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| JPLEPH(JPL の DE430 バイナリデータ)読み込みライブラリ | |
| date name version | |
| 2018.04.05 mk-mode.com 1.00 新規作成 | |
| Copyright(C) 2018 mk-mode.com All Rights Reserved. | |
| ---------------------------------------------------------------------- | |
| 1: 水星 (Mercury) | |
| 2: 金星 (Venus) | |
| 3: 地球 - 月の重心 (Earth-Moon barycenter) | |
| 4: 火星 (Mars) | |
| 5: 木星 (Jupiter) | |
| 6: 土星 (Saturn) | |
| 7: 天王星 (Uranus) | |
| 8: 海王星 (Neptune) | |
| 9: 冥王星 (Pluto) | |
| 10: 月（地心） (Moon (geocentric)) | |
| 11: 太陽 (Sun) | |
| 12: 地球の章動（1980年IAUモデル） (Earth Nutations in longitude and obliquity(IAU 1980 model)) | |
| 13: 月の秤動 (Lunar mantle libration) | |
| """ | |
| import datetime | |
| import struct | |
| class JplReadDe430: | |
| FILE_BIN = "JPLEPH" | |
| KSIZE = 2036 | |
| RECL = 4 | |
| def __init__(self, jd): | |
| self.pos = 0 # レコード位置 | |
| self.jd = jd | |
| def read_de430(self): | |
| """ Execution """ | |
| try: | |
| with open(self.FILE_BIN, "rb") as f: | |
| self.__get_ttl(f) # TTL | |
| self.__get_cnam(f) # CNAM | |
| self.__get_ss(f) # SS | |
| self.__get_ncon(f) # NCON | |
| self.__get_au(f) # AU | |
| self.__get_emrat(f) # EMRAT | |
| self.__get_ipt(f) # IPT | |
| self.__get_numde(f) # NUMDE | |
| self.__get_ipt_13(f) # IPT | |
| self.__get_cval(f) # CVAL | |
| self.__get_coeff(f) # 係数取得 | |
| self.__get_jdepoc() # JDEPOC | |
| except Exception as e: | |
| raise | |
| def __get_ttl(self, f): | |
| """ TTL（タイトル）取得 | |
| - 84 byte * 3 | |
| - ASCII文字列(後続のnull文字やスペースを削除) | |
| :param file_object f | |
| """ | |
| len_rec = 84 | |
| self.ttl = "" | |
| try: | |
| for i in range(3): | |
| f.seek(self.pos + len_rec * i) | |
| a = struct.unpack(str(len_rec) + "s", f.read(len_rec))[0] | |
| if i != 0: | |
| self.ttl += "\n" | |
| self.ttl += a.decode("utf-8").rstrip() | |
| self.pos += len_rec * 3 | |
| except Exception as e: | |
| raise | |
| def __get_cnam(self, f): | |
| """ CNAM（定数名）取得 | |
| - 6 byte * 400 | |
| - ASCII文字列(後続のnull文字やスペースを削除) | |
| :param file_object f | |
| """ | |
| len_rec = 6 | |
| self.cnams = [] | |
| try: | |
| for i in range(400): | |
| f.seek(self.pos + len_rec * i) | |
| a = struct.unpack(str(len_rec) + "s", f.read(len_rec))[0] | |
| self.cnams.append(a.decode("utf-8").rstrip()) | |
| self.pos += len_rec * 400 | |
| except Exception as e: | |
| raise | |
| def __get_ss(self, f): | |
| """ SS（ユリウス日（開始、終了）、分割日数）取得 | |
| - 8 byte * 3 | |
| - 倍精度浮動小数点数(機種依存) | |
| :param file_object f | |
| """ | |
| len_rec = 8 | |
| self.sss = [] | |
| try: | |
| for i in range(3): | |
| f.seek(self.pos + len_rec * i) | |
| a = struct.unpack("d", f.read(len_rec))[0] | |
| self.sss.append(a) | |
| self.pos += len_rec * 3 | |
| except Exception as e: | |
| raise | |
| def __get_ncon(self, f): | |
| """ NCON（定数の数）取得 | |
| - 4 byte * 1 | |
| - unsigned int (符号なし整数, エンディアンと int のサイズに依存) | |
| :param file_object f | |
| """ | |
| len_rec = 4 | |
| try: | |
| f.seek(self.pos) | |
| a = struct.unpack("I", f.read(len_rec))[0] | |
| self.ncon = a | |
| self.pos += len_rec | |
| except Exception as e: | |
| raise | |
| def __get_au(self, f): | |
| """ AU（天文単位）取得 | |
| - 8 byte * 1 | |
| - 倍精度浮動小数点数(機種依存) | |
| :param file_object f | |
| """ | |
| len_rec = 8 | |
| try: | |
| f.seek(self.pos) | |
| a = struct.unpack("d", f.read(len_rec))[0] | |
| self.au = a | |
| self.pos += len_rec | |
| except Exception as e: | |
| raise | |
| def __get_emrat(self, f): | |
| """ EMRAT（地球と月の質量比）取得 | |
| - 8 byte * 1 | |
| - 倍精度浮動小数点数(機種依存) | |
| :param file_object f | |
| """ | |
| len_rec = 8 | |
| try: | |
| f.seek(self.pos) | |
| a = struct.unpack("d", f.read(len_rec))[0] | |
| self.emrat = a | |
| self.pos += len_rec | |
| except Exception as e: | |
| raise | |
| def __get_numde(self, f): | |
| """ NUMDE（DEバージョン番号）取得 | |
| - 4 byte * 1 | |
| - unsigned int (符号なし整数, エンディアンと int のサイズに依存) | |
| :param file_object f | |
| """ | |
| len_rec = 4 | |
| try: | |
| f.seek(self.pos) | |
| a = struct.unpack("I", f.read(len_rec))[0] | |
| self.numde = a | |
| self.pos += len_rec | |
| except Exception as e: | |
| raise | |
| def __get_ipt(self, f): | |
| """ IPT（オフセット、係数の数、サブ区間数）（水星〜月の章動）取得 | |
| - 4 byte * 12 * 3 | |
| - unsigned int (符号なし整数, エンディアンと int のサイズに依存) | |
| :param file_object f | |
| """ | |
| len_rec = 4 | |
| self.ipts = [] | |
| try: | |
| for i in range(12): | |
| l = [] | |
| for j in range(3): | |
| f.seek(self.pos + len_rec * j) | |
| a = struct.unpack("I", f.read(len_rec))[0] | |
| l.append(a) | |
| self.ipts.append(l) | |
| self.pos += len_rec * 3 | |
| except Exception as e: | |
| raise | |
| def __get_ipt_13(self, f): | |
| """ IPT_13（オフセット、係数の数、サブ区間数）（月の秤動）取得 | |
| - 4 byte * 1 * 3 | |
| - unsigned int (符号なし整数, エンディアンと int のサイズに依存) | |
| :param file_object f | |
| """ | |
| len_rec = 4 | |
| try: | |
| l = [] | |
| for i in range(3): | |
| f.seek(self.pos + len_rec * i) | |
| a = struct.unpack("I", f.read(len_rec))[0] | |
| l.append(a) | |
| self.ipts.append(l) | |
| self.pos += len_rec * 3 | |
| except Exception as e: | |
| raise | |
| def __get_cval(self, f): | |
| """ CVAL（定数値）取得 | |
| - 8 byte * @ncon | |
| - 倍精度浮動小数点数(機種依存) | |
| :param file_object f | |
| """ | |
| self.pos = self.KSIZE * self.RECL | |
| len_rec = 8 | |
| self.cvals = [] | |
| try: | |
| for i in range(self.ncon): | |
| f.seek(self.pos + len_rec * i) | |
| a = struct.unpack("d", f.read(len_rec))[0] | |
| self.cvals.append(a) | |
| except Exception as e: | |
| raise | |
| def __get_coeff(self, f): | |
| """ COEFF 取得 | |
| * レコード位置計算 | |
| * 対象区間のユリウス日（開始、終了）をインスタンス変数(list) | |
| self.jds に設定する | |
| * 全ての係数の値をインスタンス変数（配列） self.coeffs に設定する | |
| （x, y, z やサブ区間毎に分割して格納） | |
| * 8 byte * ? | |
| * 倍精度浮動小数点数(機種依存) | |
| * 最初の2要素は当該データの開始・終了ユリウス日 | |
| """ | |
| idx = (self.jd - self.sss[0]) // self.sss[2] | |
| pos = int(self.KSIZE * self.RECL * (2 + idx)) | |
| recl = 8 | |
| items = [] | |
| self.coeffs = [] | |
| try: | |
| for i in range(self.KSIZE // 2): | |
| f.seek(pos + recl * i) | |
| a = struct.unpack("d", f.read(recl))[0] | |
| items.append(a) | |
| self.jds = [items.pop(0), items.pop(0)] | |
| for i, ipt in enumerate(self.ipts): | |
| n = 2 if i == 11 else 3 | |
| l_1 = [] | |
| for _ in range(ipt[2]): | |
| l_0 = [] | |
| for _ in range(n): | |
| l_0.append(items[:ipt[1]]) | |
| items = items[ipt[1]:] | |
| l_1.append(l_0) | |
| self.coeffs.append(l_1) | |
| except Exception as e: | |
| raise | |
| def __get_jdepoc(self): | |
| """ JDEPOC 取得 | |
| * self.cvals（定数値配列）の中から「元期」をインスタンス変数 | |
| self.jdepoc に設定する | |
| """ | |
| try: | |
| self.jdepoc = self.cvals[4] | |
| except Exception as e: | |
| raise |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment