Observing a set of passes

Out:

Space object 1: <Time object: scale='utc' format='mjd' value=53005.0>:
a    : 7.2000e+06   x : -1.2858e+06
e    : 1.0000e-01   y : 2.3487e+06
i    : 7.5000e+01   z : 6.3831e+06
omega: 0.0000e+00   vx: -2.1177e+03
Omega: 7.9000e+01   vy: -6.8611e+03
anom : 6.2000e+01   vz: 2.8722e+03
Parameters: C_D=2.3, m=1.0, C_R=1.0, d=1.0

Space object 1: <Time object: scale='utc' format='mjd' value=53005.0>:
a    : 7.2000e+06   x : -1.2822e+06
e    : 1.0000e-01   y : 2.3603e+06
i    : 7.5000e+01   z : 6.3782e+06
omega: 0.0000e+00   vx: -2.1203e+03
Omega: 7.9000e+01   vy: -6.8563e+03
anom : 6.1900e+01   vz: 2.8852e+03
Parameters: C_D=2.3, m=1.0, C_R=1.0, d=1.0
Temporal points: 177678
/home/danielk/IRF/IRF_GITLAB/pyant/pyant/coordinates.py:89: RuntimeWarning: invalid value encountered in greater
  proj[proj > 1.0] = 1.0
/home/danielk/IRF/IRF_GITLAB/pyant/pyant/coordinates.py:90: RuntimeWarning: invalid value encountered in less
  proj[proj < -1.0] = -1.0
/home/danielk/IRF/IRF_GITLAB/SORTS/sorts/passes.py:341: RuntimeWarning: invalid value encountered in less
  check_st = zenith_st_ang < 90.0-station.min_elevation
     t [s]    rx0 az    rx0 el    rx1 az    rx1 el    rx2 az    rx2 el  Controller    Target
----------  --------  --------  --------  --------  --------  --------  ------------  -------------
  0.121569  197.459    65.9292  234.018    68.8548  207.813    79.3268  Tracker       Cool object 1
 17.3843    176.736    75.5346  246.31     80.3619   93.1725   86.6929  Tracker       Cool object 1
 34.6471    121.568    79.5048   -3.0253   85.1947   52.1712   75.8492  Tracker       Cool object 1
 51.9099     78.9973   73.1485   31.588    74.8055   47.599    65.3952  Tracker       Cool object 1
 69.1727     64.3553   64.2488   37.3099   64.7108   45.9324   56.4559  Tracker       Cool object 1
 86.4354     57.9827   55.9974   39.6176   56.0348   45.0938   49.0171  Tracker       Cool object 1
103.698      54.5242   48.9105   40.8829   48.7827   44.6042   42.8685  Tracker       Cool object 1
120.961      52.3813   42.9503   41.6958   42.7618   44.2945   37.7639  Tracker       Cool object 1
138.224      50.9369   37.9469   42.2721   37.7423   44.0895   33.4827  Tracker       Cool object 1
155.487      49.9059   33.7186   42.7092   33.5171   43.951    29.8467  Tracker       Cool object 1
172.749      49.1393   30.1081   43.0576   29.9175   43.8572   26.7175  Tracker       Cool object 1
--------------------------------------- Performance analysis --------------------------------------
 Name                                                   |   Executions | Mean time     | Total time
--------------------------------------------------------+--------------+---------------+--------------
 Obs.Param.:calculate_observation:get_state             |            3 | 1.96428e-02 s | 26.28 %
 Obs.Param.:calculate_observation:enus,range,range_rate |            3 | 4.16517e-04 s | 0.56 %
 Obs.Param.:calculate_observation:snr-step:gain         |           30 | 4.88706e-03 s | 65.38 %
 Obs.Param.:calculate_observation:snr-step:snr          |           30 | 2.59320e-05 s | 0.35 %
 Obs.Param.:calculate_observation:snr-step              |           30 | 4.93180e-03 s | 65.98 %
 Obs.Param.:calculate_observation:generator             |            3 | 5.45375e-02 s | 72.96 %
 Obs.Param.:calculate_observation                       |            3 | 7.46603e-02 s | 99.89 %
 total                                                  |            1 | 2.24236e-01 s | 100.00 %
---------------------------------------------------------------------------------------------------

import pathlib
from tabulate import tabulate
import numpy as np
import matplotlib.pyplot as plt

import sorts
eiscat3d = sorts.radars.eiscat3d
from sorts.controller import Tracker
from sorts.scheduler import StaticList, ObservedParameters
from sorts import SpaceObject
from sorts.profiling import Profiler

from sorts.propagator import SGP4
Prop_cls = SGP4
Prop_opts = dict(
    settings = dict(
        out_frame='ITRF',
    ),
)
prop = Prop_cls(**Prop_opts)


objs = [
    SpaceObject(
        Prop_cls,
        propagator_options = Prop_opts,
        a = 7200e3,
        e = 0.1,
        i = 75,
        raan = 79,
        aop = 0,
        mu0 = mu0,
        epoch = 53005.0,
        parameters = dict(
            d = 1.0,
        ),
    )
    for mu0 in [62.0, 61.9]
]

for obj in objs: print(obj)

t = sorts.equidistant_sampling(
    orbit = objs[0].state,
    start_t = 0,
    end_t = 3600*6,
    max_dpos=1e3,
)

print(f'Temporal points: {len(t)}')
states0 = objs[0].get_state(t)
states1 = objs[1].get_state(t)

#set cache_data = True to save the data in local coordinates
#for each pass inside the Pass instance, setting to false saves RAM
passes0 = eiscat3d.find_passes(t, states0, cache_data = False)
passes1 = eiscat3d.find_passes(t, states1, cache_data = False)


#just create a controller for observing 10 points of the first pass
ps = passes0[0][0][0]
use_inds = np.arange(0,len(ps.inds),len(ps.inds)//10)
e3d_tracker = Tracker(radar = eiscat3d, t=t[ps.inds[use_inds]], ecefs=states0[:3,ps.inds[use_inds]])
e3d_tracker.meta['target'] = 'Cool object 1'

class MyStaticList(StaticList, ObservedParameters):

    def __init__(self, radar, controllers, profiler=None, logger=None):
        super().__init__(
            radar=radar,
            controllers=controllers,
            profiler=profiler,
            logger=logger,
        )

    def generate_schedule(self, t, generator):
        data = np.empty((len(t),len(self.radar.rx)*2+1), dtype=np.float64)
        data[:,0] = t
        names = []
        targets = []
        for ind,mrad in enumerate(generator):
            radar, meta = mrad
            names.append(meta['controller_type'].__name__)
            targets.append(meta['target'])
            for ri, rx in enumerate(radar.rx):
                data[ind,1+ri*2] = rx.beam.azimuth
                data[ind,2+ri*2] = rx.beam.elevation
        data = data.T.tolist() + [names, targets]
        data = list(map(list, zip(*data)))
        return data


p = Profiler()

scheduler = MyStaticList(radar = eiscat3d, controllers=[e3d_tracker], profiler=p)

sched_data = scheduler.schedule()

rx_head = [f'rx{i} {co}' for i in range(len(scheduler.radar.rx)) for co in ['az', 'el']]
sched_tab = tabulate(sched_data, headers=["t [s]"] + rx_head + ['Controller', 'Target'])

print(sched_tab)

p.start('total')
data0 = scheduler.observe_passes(passes0, space_object = objs[0], snr_limit=False)
p.stop('total')
print(p.fmt(normalize='total'))

data1 = scheduler.observe_passes(passes1, space_object = objs[1], snr_limit=False)

#create a tdm file example
# pth = pathlib.Path(__file__).parent / 'data' / 'test_tdm.tdm'
# print(f'Writing TDM data to: {pth}')

# dat = data0[0][0][0]
# sorts.io.write_tdm(
#     pth,
#     dat['t'],
#     dat['range'],
#     dat['range_rate'],
#     np.ones(dat['range'].shape),
#     np.ones(dat['range_rate'].shape),
#     freq=eiscat3d.tx[0].beam.frequency,
#     tx_ecef=eiscat3d.tx[0].ecef,
#     rx_ecef=eiscat3d.rx[0].ecef,
#     tx_name="EISCAT 3D Skiboten",
#     rx_name="EISCAT 3D Skiboten",
#     oid="Some cool space object",
#     tdm_type="track",
# )


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),
    ],
]

for tx in scheduler.radar.tx:
    axes[0][0].plot([tx.ecef[0]],[tx.ecef[1]],[tx.ecef[2]], 'or')
for rx in scheduler.radar.rx:
    axes[0][0].plot([rx.ecef[0]],[rx.ecef[1]],[rx.ecef[2]], 'og')

for pi in range(len(passes0[0][0])):
    dat = data0[0][0][pi]
    dat2 = data1[0][0][pi]
    if dat is not None:
        axes[0][0].plot(states0[0,passes0[0][0][pi].inds], states0[1,passes0[0][0][pi].inds], states0[2,passes0[0][0][pi].inds], '-', label=f'pass-{pi}')
        axes[0][1].plot(dat['t']/3600.0, dat['range'], '-', label=f'pass-{pi}')
        axes[1][0].plot(dat['t']/3600.0, dat['range_rate'], '-', label=f'pass-{pi}')
        axes[1][1].plot(dat['t']/3600.0, 10*np.log10(dat['snr']), '-', label=f'pass-{pi}')
    if dat2 is not None:
        axes[0][0].plot(states1[0,passes1[0][0][pi].inds], states1[1,passes1[0][0][pi].inds], states1[2,passes1[0][0][pi].inds], '-', label=f'obj2 pass-{pi}')
        axes[0][1].plot(dat2['t']/3600.0, dat2['range'], '-', label=f'obj2 pass-{pi}')
        axes[1][0].plot(dat2['t']/3600.0, dat2['range_rate'], '-', label=f'obj2 pass-{pi}')
        axes[1][1].plot(dat2['t']/3600.0, 10*np.log10(dat2['snr']), '-', label=f'obj2 pass-{pi}')

axes[0][1].legend()
plt.show()