#!/usr/bin/env python3

# For use in GNU Radio w/ gr-satellites.   See: https://destevez.net/2022/07/real-time-doppler-correction-with-gnu-radio/

# Based on:  https://github.com/daniestevez/jupyter_notebooks/blob/master/Vega-C_MEO/Doppler.ipynb

print("\n")
print("Working...")
print("\n")

import numpy as np
from skyfield.api import EarthSatellite, load, wgs84, utc
from scipy.constants import c
import datetime
import calendar, time;

# ------------------------------------------------------
# ENTER ALL CHANGES FOR EACH RECORDING BELOW THIS LINE
# ------------------------------------------------------

# IQ file = gqrx_20220825_031838_435525000_78125_picsat.raw   (for reference only)

# PICSAT     435.525             
# 1 43132U 18004X   22236.21879458  .00006336  00000+0  20172-3 0  9997
# 2 43132  97.3692 312.2838 0008813 115.2004 245.0152 15.32892351256787

# set recording start date/time values (UTC):
rec_year = 2022
rec_month = 8
rec_day = 25
rec_hour = 3
rec_minute = 18
rec_second = 38

# set recording center RF frequency (Hz) in format:  f_carrier = 2232235000  (true RF freq, NOT the same as SDR record freq if downconverter used)
f_carrier = 435525000

# enter CURRENT TLE info for satellite tracked in this recording (don't forget ' marks before & after each line)
line1 = '1 43132U 18004X   22236.21879458  .00006336  00000+0  20172-3 0  9997'
line2 = '2 43132  97.3692 312.2838 0008813 115.2004 245.0152 15.32892351256787'

# enter lat/lon/elevation(meters) of your receive location
lat = 37.7790
lon = -77.6097
el = 61

# ------------------------------------------------------
# ENTER ALL CHANGES FOR EACH RECORDING ABOVE THIS LINE
# ------------------------------------------------------

ts = load.timescale()

sat = EarthSatellite(line1, line2, 'sat', ts)

satellites = [sat]

f_carrier = np.array([f_carrier])

groundstation = wgs84.latlon(lat, lon, el)

t0_dt = datetime.datetime(rec_year, rec_month, rec_day, rec_hour, rec_minute, rec_second, tzinfo=utc)
t0_ts = ts.utc(t0_dt.year, t0_dt.month, t0_dt.day, t0_dt.hour, t0_dt.minute, t0_dt.second)

# calculate unix epoch time at start of recording
unix_epoch_record_start_time = calendar.timegm(t0_dt.utctimetuple())

doppler_T = 0.1
T = 2 * 3600
dts = np.arange(int(T / doppler_T)) * doppler_T
times = ts.from_datetimes([t0_dt + datetime.timedelta(seconds=dt) for dt in dts])

range_rates = np.array([
    (sat - groundstation).at(times).frame_latlon_and_rates(groundstation)[-1].m_per_s
    for sat in satellites])
dopplers = -range_rates/c*f_carrier

for freqs, sat in zip(dopplers, satellites):
    with open(f'/tmp/doppler.txt', 'w') as file:
        for t, f in zip(times, freqs):
            unix_t = (t.utc_datetime()
                      - datetime.datetime(1970, 1, 1, tzinfo=utc)).total_seconds()
            print(unix_t, f, file=file)

print("\n")
print("Created file '/tmp/doppler.txt' for use in the 'Doppler Correction' block from gr-satellites")
print("\n")
print("... for the 'Start time' field in your 'Doppler Correction' block, enter: " + str(unix_epoch_record_start_time))
print("\n")