In [1]:
import astropy.units as u
from astropy.constants import c as speed_of_light
from astropy.coordinates import EarthLocation, solar_system_ephemeris
from astropy.coordinates import SkyCoord, AltAz, ICRS, get_moon, get_sun, get_body
from astropy.time import Time
from datetime import datetime
from astropy.table import Table
import numpy as np
import pandas as pd 

import astropy.constants
import astropy.units as u

from datetime import datetime

# SPICE
import spiceypy
In [2]:
spiceypy.furnsh("kernels/naif0012.tls")
spiceypy.furnsh("kernels/de440s.bsp")
spiceypy.furnsh("kernels/pck00011.tpc")
spiceypy.furnsh("kernels/earth_latest_high_prec.bpc")
In [3]:
[dim, radii] = spiceypy.bodvrd("EARTH", "RADII", 3)
flattening = (radii[0]-radii[2])/radii[0];

# Dwingeloo
h = 70.26/1000;
lon =  6.396346463227839;
lat =  52.81213723180477;

dwingeloo_obs = spiceypy.georec(lon*np.pi/180, lat*np.pi/180, h, radii[0], flattening)

# Stockert
h = 434/1000;
lon =  6.721943350231514;
lat =  50.56944039751571;

stockert_obs = spiceypy.georec(lon*np.pi/180, lat*np.pi/180, h, radii[0], flattening)
In [4]:
tx_obs = dwingeloo_obs
current_obs = stockert_obs
In [5]:
# object
body = "VENUS"
[dim, radii] = spiceypy.bodvrd(body, "RADII", 3)
body_radius = radii[0]*1000
body_radius
Out[5]:
6051800.0
In [6]:
tx_time = (["2025-03-22T12:01:00.000", 
            "2025-03-22T12:11:00.000", 
            "2025-03-22T12:21:00.000", 
            "2025-03-22T12:31:00.000",
            ])

l = len(tx_time)

time_astropy = Time(tx_time)
ttx = time_astropy
lighttime_up = np.zeros(l)
lighttime_down = np.zeros(l)
In [7]:
for t in range(int(len(tx_time))):
    et = spiceypy.str2et( Time(ttx[t].isot).utc.value )
    [azlsta,lt] = spiceypy.azlcpo("ELLIPSOID", body, et, "XCN+S", False, True, tx_obs, "EARTH", "ITRF93")
    lighttime_up[int(t)] = lt-body_radius/astropy.constants.c.value
In [8]:
trx = ttx + 2*lighttime_up*u.s
trx.precision = 9
ttx.precision = 9

for i in range(5):
    for t in range(int(len(tx_time))):
        et = spiceypy.str2et( trx[t].utc.value )
        [azlsta,lt] = spiceypy.azlcpo("ELLIPSOID", body, et, "CN+S", False, True, current_obs, "EARTH", "ITRF93")
        lighttime_down[t] = lt-body_radius/astropy.constants.c.value
        trx[t] = ttx[t] + lighttime_down[t]*u.s + lighttime_up[t]*u.s
In [9]:
print(lighttime_down+lighttime_up)

[279.9744859  279.97437085 279.97431935 279.97433093]
In [10]:
start_rx = time_astropy + (lighttime_down+lighttime_up)*u.s
start_rx.precision = 6
start_rx.isot
Out[10]:
array(['2025-03-22T12:05:39.974486', '2025-03-22T12:15:39.974371',
       '2025-03-22T12:25:39.974319', '2025-03-22T12:35:39.974331'],
      dtype='<U26')
In [11]:
# Doppler
In [12]:
time_astropy = Time("2025-03-22T11:59:59.000000")
et = spiceypy.str2et( Time(time_astropy.isot).utc.value )

duration = 50*60
f0 = 1299.5e6
v_down = np.zeros(int(duration))
v_up = np.zeros(int(duration))
d = np.zeros(int(duration))
lighttime_up = np.zeros(int(duration))
lighttime_down = np.zeros(int(duration))
In [13]:
trx = Time(time_astropy.isot) + np.arange(int(duration))*u.s

for t in range(int(duration)):
    timestamp = trx[t]
    et = spiceypy.str2et( timestamp.utc.value )
    [azlsta,lt] = spiceypy.azlcpo("ELLIPSOID", body, et, "CN+S", False, True, current_obs, "EARTH", "ITRF93")
    v_down[t] = 1000*azlsta[3]
    d[t] = 1000*azlsta[0]
    lighttime_down[t] = lt-body_radius/astropy.constants.c.value
In [14]:
ttx = trx - 2*lighttime_down*u.s

for i in range(5):
    for t in range(int(duration)):
        timestamp = ttx[t]
        et = spiceypy.str2et( timestamp.utc.value )
        [azlsta,lt] = spiceypy.azlcpo("ELLIPSOID", body, et, "XCN+S", False, True, dwingeloo_obs, "EARTH", "ITRF93")
        v_up[t] = 1000*azlsta[3]
        lighttime_up[t] = lt-body_radius/astropy.constants.c.value
    ttx = trx - lighttime_down*u.s - lighttime_up*u.s
In [15]:
def doppler(v,f):
    return (1 - v / astropy.constants.c.value) * f
In [16]:
dop_up = doppler(v_up,f0)-f0
dop_down = doppler(v_down,f0)-f0
In [17]:
freq_offsets = dop_up[1:]+dop_down[1:]
doppler_rates = np.diff(dop_up+dop_down)
In [18]:
# a = np.asarray([trx[1:], freq_offsets, doppler_rates ]).T
# df = pd.DataFrame(a)
# df.to_csv("stockert_venus_doppler.csv", index=False, header=("rx_time_utc","freq_offset_hz", "doppler_rate_hz_s"))
In [ ]: