Note
Click here to download the full example code
Optimizing with interpolationΒΆ
Out:
Space object 1: <Time object: scale='utc' format='mjd' value=53005.0>:
a : 7.2000e+06 x : -1.3922e+06
e : 2.0000e-02 y : 3.4519e+06
i : 7.5000e+01 z : 6.0816e+06
omega: 0.0000e+00 vx: -1.3986e+03
Omega: 8.6000e+01 vy: -6.4897e+03
anom : 6.0000e+01 vz: 3.5175e+03
Parameters: C_D=2.3, m=1.0, C_R=1.0, d=0.1, A=1.0
2020-10-15 13:29:52.074 INFO ; Scanner:init
Temporal points obj 0: 1440
--------------------------------------- Performance analysis --------------------------------------
Name | Executions | Mean time | Total time
--------------------------------------------------------+--------------+---------------+--------------
get_state | 1 | 2.03223e+00 s | 1.27 %
find_passes | 1 | 3.71742e-03 s | 0.00 %
Obs.Param.:calculate_observation:get_state | 36 | 1.11971e+00 s | 25.23 %
Obs.Param.:calculate_observation:enus,range,range_rate | 36 | 5.28826e-04 s | 0.01 %
Scanner:generator:point_radar:_point_station[tx] | 78036 | 1.10661e-04 s | 5.40 %
Scanner:generator:point_radar:_point_station[rx] | 234108 | 1.98539e-04 s | 29.09 %
Scanner:generator:point_radar | 78036 | 8.53267e-04 s | 41.67 %
Obs.Param.:calculate_observation:snr-step:gain | 76274 | 6.22534e-04 s | 29.72 %
Obs.Param.:calculate_observation:snr-step:snr | 76274 | 2.17526e-05 s | 1.04 %
Obs.Param.:calculate_observation:snr-step | 78036 | 6.45173e-04 s | 31.51 %
Obs.Param.:calculate_observation:generator | 36 | 3.25961e+00 s | 73.44 %
Obs.Param.:calculate_observation | 36 | 4.38008e+00 s | 98.69 %
observe_passes | 1 | 9.89323e+01 s | 61.92 %
observe_passes_interpolator | 1 | 5.88130e+01 s | 36.81 %
total | 1 | 1.59781e+02 s | 100.00 %
---------------------------------------------------------------------------------------------------
/home/danielk/IRF/IRF_GITLAB/SORTS/examples/scheduler_interpolation_optimizing.py:152: RuntimeWarning: divide by zero encountered in log10
SNRdB = 10*np.log10(dat['snr'])
/home/danielk/IRF/IRF_GITLAB/SORTS/examples/scheduler_interpolation_optimizing.py:152: RuntimeWarning: divide by zero encountered in log10
SNRdB = 10*np.log10(dat['snr'])
/home/danielk/IRF/IRF_GITLAB/SORTS/examples/scheduler_interpolation_optimizing.py:152: RuntimeWarning: divide by zero encountered in log10
SNRdB = 10*np.log10(dat['snr'])
/home/danielk/IRF/IRF_GITLAB/SORTS/examples/scheduler_interpolation_optimizing.py:152: RuntimeWarning: divide by zero encountered in log10
SNRdB = 10*np.log10(dat['snr'])
/home/danielk/IRF/IRF_GITLAB/SORTS/examples/scheduler_interpolation_optimizing.py:152: RuntimeWarning: divide by zero encountered in log10
SNRdB = 10*np.log10(dat['snr'])
/home/danielk/IRF/IRF_GITLAB/SORTS/examples/scheduler_interpolation_optimizing.py:152: RuntimeWarning: divide by zero encountered in log10
SNRdB = 10*np.log10(dat['snr'])
import numpy as np
import matplotlib.pyplot as plt
import sorts
eiscat3d = sorts.radars.eiscat3d_interp
from sorts.scheduler import StaticList, ObservedParameters
from sorts.controller import Scanner
from sorts import SpaceObject
from sorts.profiling import Profiler
from sorts.radar.scans import Fence
from sorts.interpolation import Legendre8
from sorts.propagator import Orekit
orekit_data = '/home/danielk/IRF/IRF_GITLAB/orekit_build/orekit-data-master.zip'
Prop_cls = Orekit
Prop_opts = dict(
orekit_data = orekit_data,
settings = dict(
in_frame='GCRS',
out_frame='ITRS',
),
)
end_t = 12*3600.0
scan = Fence(azimuth=90, num=40, dwell=0.1, min_elevation=30)
p = Profiler()
logger = sorts.profiling.get_logger('scanning')
objs = [
SpaceObject(
Prop_cls,
propagator_options = Prop_opts,
a = 7200e3,
e = 0.02,
i = 75,
raan = 86,
aop = 0,
mu0 = 60,
epoch = 53005.0,
parameters = dict(
d = 0.1,
A = 1.0,
),
),
]
for obj in objs: print(obj)
class ObservedScanning(StaticList, ObservedParameters):
pass
scanner_ctrl = Scanner(eiscat3d, scan, profiler=p, logger=logger)
scanner_ctrl.t = np.arange(0, end_t, scan.dwell())
p.start('total')
scheduler = ObservedScanning(
radar = eiscat3d,
controllers = [scanner_ctrl],
logger = logger,
profiler = p,
)
t = np.arange(0.0,end_t,30.0)
datas = []
passes = []
states = []
for ind in range(len(objs)):
print(f'Temporal points obj {ind}: {len(t)}')
p.start('get_state')
states += [objs[ind].get_state(t)]
p.stop('get_state')
interpolator = Legendre8(states[ind], t)
p.start('find_passes')
#rename cache_data to something more descriptive
passes += [eiscat3d.find_passes(t, states[ind], cache_data = False)]
p.stop('find_passes')
p.start('observe_passes')
data = scheduler.observe_passes(passes[ind], space_object = objs[ind], interpolator = None, snr_limit=False)
p.stop('observe_passes')
p.start('observe_passes_interpolator')
data = scheduler.observe_passes(passes[ind], space_object = objs[ind], interpolator = interpolator, snr_limit=False)
p.stop('observe_passes_interpolator')
datas.append(data)
p.stop('total')
print(p.fmt(normalize='total'))
fig = plt.figure(figsize=(15,15))
axes = [
[
fig.add_subplot(221, projection='3d'),
fig.add_subplot(222),
],
[
fig.add_subplot(223),
fig.add_subplot(224),
],
]
sorts.plotting.grid_earth(axes[0][0])
for tx in eiscat3d.tx:
axes[0][0].plot([tx.ecef[0]],[tx.ecef[1]],[tx.ecef[2]],'or')
for rx in eiscat3d.rx:
axes[0][0].plot([rx.ecef[0]],[rx.ecef[1]],[rx.ecef[2]],'og')
for radar, meta in scanner_ctrl(np.arange(0,scan.cycle(),scan.dwell())):
for tx in radar.tx:
point_tx = tx.pointing_ecef/np.linalg.norm(tx.pointing_ecef, axis=0)*scanner_ctrl.r.max() + tx.ecef
axes[0][0].plot([tx.ecef[0], point_tx[0]], [tx.ecef[1], point_tx[1]], [tx.ecef[2], point_tx[2]], 'r-', alpha=0.15)
for rx in radar.rx:
pecef = rx.pointing_ecef/np.linalg.norm(rx.pointing_ecef, axis=0)
for ri in range(pecef.shape[1]):
point_tx = tx.pointing_ecef/np.linalg.norm(tx.pointing_ecef, axis=0)*scanner_ctrl.r[ri] + tx.ecef
point = pecef[:,ri]*np.linalg.norm(rx.ecef - point_tx) + rx.ecef
axes[0][0].plot([rx.ecef[0], point[0]], [rx.ecef[1], point[1]], [rx.ecef[2], point[2]], 'g-', alpha=0.05)
for ind in range(len(objs)):
for pi in range(len(passes[ind][0][0])):
ps = passes[ind][0][0][pi]
dat = datas[ind][0][0][pi]
axes[0][0].plot(states[ind][0,ps.inds], states[ind][1,ps.inds], states[ind][2,ps.inds], '-')
if dat is not None:
SNRdB = 10*np.log10(dat['snr'])
det_inds = SNRdB > 10.0
axes[0][1].plot(dat['t']/3600.0, dat['range']*1e-3, '-', label=f'obj{ind}-pass{pi}')
axes[1][0].plot(dat['t']/3600.0, dat['range_rate']*1e-3, '-')
axes[1][1].plot(dat['t']/3600.0, SNRdB, '-')
axes[0][1].plot(dat['t'][det_inds]/3600.0, dat['range'][det_inds]*1e-3, '.r')
axes[1][0].plot(dat['t'][det_inds]/3600.0, dat['range_rate'][det_inds]*1e-3, '.r')
axes[1][1].plot(dat['t'][det_inds]/3600.0, SNRdB[det_inds], '.r')
axes[1][1].set_ylim([0, None])
font_ = 18
axes[0][1].set_xlabel('Time [h]', fontsize=font_)
axes[1][0].set_xlabel('Time [h]', fontsize=font_)
axes[1][1].set_xlabel('Time [h]', fontsize=font_)
axes[0][1].set_ylabel('Two way range [km]', fontsize=font_)
axes[1][0].set_ylabel('Two way range rate [km/s]', fontsize=font_)
axes[1][1].set_ylabel('SNR [dB]', fontsize=font_)
axes[0][1].legend()
dr = 600e3
axes[0][0].set_xlim([eiscat3d.tx[0].ecef[0]-dr, eiscat3d.tx[0].ecef[0]+dr])
axes[0][0].set_ylim([eiscat3d.tx[0].ecef[1]-dr, eiscat3d.tx[0].ecef[1]+dr])
axes[0][0].set_zlim([eiscat3d.tx[0].ecef[2]-dr, eiscat3d.tx[0].ecef[2]+dr])
plt.show()
Total running time of the script: ( 2 minutes 41.576 seconds)