#!/usr/bin/python3
# tle2ssc
# Convert a TLE into an .ssc file for Celestia
#
# Usage:
# tle2ssc [filename]
# Example:
# tle2ssc foo-tle.txt

# Check:
# https://pypi.org/project/sgp4/
# https://rhodesmill.org/skyfield/earth-satellites.html
# https://github.com/ivanstan/tle-api
# https://github.com/TruSat/trusat-orbit

import os
import skyfield
import math
from datetime import datetime
from skyfield.api import load, wgs84
from skyfield.api import EarthSatellite
from sgp4.api import Satrec, WGS72

satnum=37846
satradius=0.005
satname = 'GSAT0101'
#satellites_url = 'http://celestrak.com/NORAD/elements/galileo.txt'
satellites_url = './extras/galileo-gnss/galileo.txt'
satellites = load.tle_file(satellites_url)
#print('Loaded', len(satellites), 'satellites')

ts = load.timescale()
#t = ts.now()
# 2022-05-20 02:17:30
t = ts.utc(2022, 5, 20, 2, 17, 30)

xpdotp = 1440.0 / (2.0 * math.pi)

by_number = {sat.model.satnum: sat for sat in satellites}
satellite = by_number[satnum]

# Two different ways to get Epoch...
#satepoch=(t.tdb)
satepoch=satellite.model.jdsatepoch

satellite_name=satname
satellite_radius=satradius
satellite_epoch=satepoch
# ALLLLL XXX
#    The unique satellite NORAD catalog number given in the TLE file.
satellite_number=satellite.model.satnum
#    Satellite classification, or else 'U' for “Unknown”
satellite_classification=satellite.model.classification
#    International designator
satellite_intldesg=satellite.model.intldesg
#    Full four-digit year of this element set’s epoch moment.
satellite_epochyr=satellite.model.epochyr
#    Fractional days into the year of the epoch moment.
satellite_epochdays=satellite.model.epochdays
#    Julian date of the epoch (computed from epochyr and epochdays).
satellite_jdsatepoch=satellite.model.jdsatepoch
#    First time derivative of the mean motion (ignored by SGP4).
satellite_ndot=satellite.model.ndot
#print('satellite_ndot', satellite_ndot)
#    Second time derivative of the mean motion (ignored by SGP4).
satellite_nddot=satellite.model.nddot
#print('satellite_nddot', satellite_nddot)
#    Ballistic drag coefficient B* in inverse earth radii.
satellite_bstar=satellite.model.bstar
#print('satellite_bstar', satellite_bstar)
#    Ephemeris type (ignored by SGP4 as determination now automatic)
satellite_ephtype=satellite.model.ephtype
#    Element number
satellite_elnum=satellite.model.elnum

#    Inclination in radians. Convert radians to degrees.
satellite_inclination=math.degrees(satellite.model.inclo)
satellite_obliquity=math.degrees(satellite.model.inclo)

#    Right ascension of ascending node in radians. Convert to degrees.
satellite_ascending_node=math.degrees(satellite.model.nodeo)
satellite_equator_ascending_node=math.degrees(satellite.model.nodeo)

#    Eccentricity.
satellite_eccentricity=satellite.model.ecco

#    Argument of perigee in radians.
satellite_arg_of_pericenter=math.degrees(satellite.model.argpo)

#    Mean anomaly in radians.
satellite_mean_anomoly=math.degrees(satellite.model.mo)

#    Mean motion in radians per minute.
satellite_no_kozai=satellite.model.no_kozai
satellite_period=(1 / (satellite.model.no_kozai * xpdotp))

#    Revolution number at epoch [Revs]
satellite_revnum=satellite.model.revnum


#    'i',             # 'a' = old AFSPC mode, 'i' = improved mode
#    5,               # satnum: Satellite number
#    18441.785,       # epoch: days since 1949 December 31 00:00 UT
#    2.8098e-05,      # bstar: drag coefficient (/earth radii)
#    6.969196665e-13, # ndot: ballistic coefficient (revs/day)
#    0.0,             # nddot: second derivative of mean motion (revs/day^3)
#    0.1859667,       # ecco: eccentricity
#    5.7904160274885, # argpo: argument of perigee (radians)
#    0.5980929187319, # inclo: inclination (radians)
#    0.3373093125574, # mo: mean anomaly (radians)
#    0.0472294454407, # no_kozai: mean motion (radians/minute)
#    6.0863854713832, # nodeo: right ascension of ascending node (radians)


# Create SSC
# From TLE:
# GSAT0101 (PRN E11)      
# 1 37846U 11060A   22140.09549104 -.00000093  00000+0  00000+0 0  9998
# 2 37846  56.9858  22.1062 0004117  28.0726 331.9949  1.70474933 65739

# Satellite name: GSAT0101
# Satellite number: 37846U
# Satellite number: 37846
# International designator: 11060A
# Inclination: 56.9858
# Epoch year and Julian day fraction: 22140.09549104
# Right ascension of ascending node: 22.1062
# Eccentricity: 0004117
# First derivative of mean motion or ballistic coefficient: -.00000093
# Argument of perigee: 28.0726
# Second derivative of mean motion: 00000+0
# Mean anomaly: 331.9949
# Drag term or radiation pressure coefficient: 00000+0
# Mean motion: 1.70474933
# Ephemeris type: 0
# Element number and checksum: 9998
# Revolution number at epoch and checksum: 65739


# Convert epoch from TLE to SSC
# Convert TLE to Julian Day
# C47 = Epoch Year = 2022
# =1721424.5-INT((C47-1)/100)+INT((C47-1)/400)+INT(365.25*(C47-1))+C49
print('"', satellite_name, '-', satellite_number, '" ','"Sol/Earth" {',sep="")
print('    Class "spacecraft"')
print('    # Mesh "foo.3ds XXX"')
print('    radius', satellite_radius)
print()
print('  EllipticalOrbit {')
print('    Epoch', satellite_epoch)
print('    Period', satellite_period)
print('    SemiMajorAxis  29600.181 XXX')
print('    Eccentricity', satellite_eccentricity)
print('    Inclination', satellite_inclination)
print('    AscendingNode', satellite_ascending_node)
print('    ArgOfPericenter', satellite_arg_of_pericenter)
print('    MeanAnomaly', satellite_mean_anomoly)
print('  }')
print('  Obliquity', satellite_obliquity)
print('  EquatorAscendingNode', satellite_equator_ascending_node)
print('  RotationOffset  312.7348 XXX')
print('  # Orientation  [  ]')
print('}')