Simulate scanning for objects

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
2020-10-15 13:25:53.795 INFO    ; Scanner:init
Temporal points obj 0: 4555
--------------------------------------- Performance analysis --------------------------------------
 Name                                                   |   Executions | Mean time     | Total time
--------------------------------------------------------+--------------+---------------+--------------
 equidistant_sampling                                   |            1 | 6.91414e-05 s | 0.00 %
 get_state                                              |            1 | 5.40323e-02 s | 0.42 %
 find_passes                                            |            1 | 5.85890e-03 s | 0.05 %
 Obs.Param.:calculate_observation:get_state             |            3 | 5.40149e-02 s | 1.27 %
 Obs.Param.:calculate_observation:enus,range,range_rate |            3 | 6.99520e-04 s | 0.02 %
 Scanner:generator:point_radar:_point_station[tx]       |         8361 | 1.10018e-04 s | 7.22 %
 Scanner:generator:point_radar:_point_station[rx]       |        25083 | 1.96978e-04 s | 38.76 %
 Scanner:generator:point_radar                          |         8361 | 8.46843e-04 s | 55.54 %
 Obs.Param.:calculate_observation:snr-step:gain         |         8149 | 6.24864e-04 s | 39.95 %
 Obs.Param.:calculate_observation:snr-step:snr          |         8149 | 2.17805e-05 s | 1.39 %
 Obs.Param.:calculate_observation:snr-step              |         8361 | 6.45645e-04 s | 42.35 %
 Obs.Param.:calculate_observation:generator             |            3 | 4.17384e+00 s | 98.23 %
 Obs.Param.:calculate_observation                       |            3 | 4.22882e+00 s | 99.52 %
 observe_passes                                         |            1 | 1.26873e+01 s | 99.53 %
 total                                                  |            1 | 1.27473e+01 s | 100.00 %
---------------------------------------------------------------------------------------------------

/home/danielk/IRF/IRF_GITLAB/SORTS/examples/simulate_scanning.py:146: 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.propagator import SGP4
Prop_cls = SGP4
Prop_opts = dict(
    settings = dict(
        out_frame='ITRF',
    ),
)

end_t = 600.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,
        ),
    ),
]


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



datas = []
passes = []
states = []
for ind in range(len(objs)):
    p.start('equidistant_sampling')
    t = sorts.equidistant_sampling(
        orbit = objs[ind].state,
        start_t = 0,
        end_t = end_t,
        max_dpos=1e3,
    )
    p.stop('equidistant_sampling')

    print(f'Temporal points obj {ind}: {len(t)}')

    p.start('get_state')
    states += [objs[ind].get_state(t)]
    p.stop('get_state')

    p.start('find_passes')
    #rename cache_data to something more descriptive
    passes += [eiscat3d.find_passes(t, states[ind], cache_data = True)]
    p.stop('find_passes')

    p.start('observe_passes')
    data = scheduler.observe_passes(passes[ind], space_object = objs[ind], snr_limit=False)
    p.stop('observe_passes')

    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()