Several users received an email on how to unlock their account. Due to a brute force login attack, accounts are automatically locked after a number of failed attempts. We recommend to regularly change passwords and use 2FA whenever possible. If you would like to activate 2-factor authentication on Gitlab, you can do so by logging in and then to Preferences -> Account. You will need an "authenticator" app and a compatible phone (we do not support SMS codes on Gitlab). If you have further questions regarding online safety, please check the website https://www.tudelft.nl/en/privacy-security.

Commit 71132192 authored by Andreas Theodosiou's avatar Andreas Theodosiou 馃惂

馃帀 Initial code commit

parent 8ec63fec
"""
====================================
DRAMA Documentation (:mod:`drama`)
====================================
.. currentmodule:: DRAMA
DRAMA: **D**\elft **RA**\dar **M**\odeling and **P**\erformance **A**\ nalysis
"""
"""
===========================================
Constants Package (:mod:`drama.constants`)
===========================================
.. currentmodule:: drama.constants
Constants
=========
.. autosummary::
:toctree: generated/
"""
from drama.constants.constants import *
"""
**Custom Added Constants**
:r_earth: Earth radius [m]
:m_earth: Earth mass [Kg]
:c: speed of light [m/s]
:g_star: gravitational constant g_star [m^3/km/s^2]
:gm_earth: g_star * mass of earth [m^3/s^2]
:au: astronomical unit [m]
:gm_sun: g_star * mass of sun [m^3/s^2]
:r_sun: radius of the sun [m]
:p_sun: radiation pressure of the sun [N/m^2]
:gm_moon: g_star * mass of moon [m^3/s^2]
:d_moon_earth: mean distance moon to earth [m]
:r_moon: radius of moon [m]
:h: Planck's constant [Js]
:omega_earth: angular velocity of earth [1/s] WGS84
:pi: pi to 30 digits
:j2: zonal coefficient
:j3: zonal coefficient
:j4: zonal coefficient
:j5: zonal coefficient
"""
import numpy as np
from scipy.constants import *
r_earth = 6378.137e3 # Earth equat. radius [m]
m_earth = 5.9742e24 # Earth mass [Kg]
c = 2.99792458e8 # speed of light [m/s]
g_star = 6.673e-11 # gravitational constant g_star [m^3/km/s^2]
gm_earth = 398600.4415e9 # g_star * mass of earth [m^3/s^2]
gm_sun = 1.32712440018e20 # g_star * mass of sun [m^3/s^2]
au = 149597870.691e3 # astronomical unit [m]
r_sun = 6.96e8 # radius of the sun [m]
p_sun = 4.560e-6 # radiation pressure of the sun [N/m^2]
w_sun = 1.99099299e-07 # mean apparent revolution speed of the sun around the earth [rad/s]
gm_moon = 4902.801e9 # g_star * mass of moon [m^3/s^2]
d_moon_earth = 384400e3 # mean distance moon to earth [m]
r_moon = 1738e3 # radius of moon [m]
h = 6.626068765e-34 # Planck's constant [Js]
omega_earth = 0.729211585530e-4 # angular velocity of earth [1/s] WGS84
pi = 3.141592653589793238462643383279 # pi to 30 digits
j2 = 0.00108263 # zonal coefficient
j3 = -0.00254e-3 # zonal coefficient
j4 = -0.00161e-3 # zonal coefficient
j5 = -0.246e-6 # zonal coefficient
# Sea water rel.dielectric constant
epsilon_sw = np.complex(73, 18)
# Water refractive index (approximate)
n_w = 4.0 / 3.0
# Impedance of vacuum [Ohm]
etha_0 = physical_constants["characteristic impedance of vacuum"][0]
# equatorial and polar radii according to different standards in meters
r_equatorial = {"wgs84": 6378.137e3, "iau": 6378.160e3, "krasovsky": 6378.245e3}
r_polar = {"wgs84": 6356.75231425e3, "iau": 6356.775e3, "'krasovsky": 6356.86301877e3}
"""
====================================
Geometry Package (:mod:`drama.geo`)
====================================
.. currentmodule:: drama.geo
"""
from drama.geo.geometry import (
inc_to_sr,
inc_to_gr,
inc_to_look,
look_to_inc,
look_to_sr,
max_look_angle,
gr_to_geo,
sr_to_geo,
ecef_to_geodetic,
pt_get_intersection_ellipsoid,
eci_to_ecef,
create_LoS,
QuickRadarGeometry,
)
from drama.geo.interferometry import Baseline
from drama.geo.timing import (
blind_ranges,
timing_diagram,
mode_on_timing,
conf_to_timing,
)
from drama.geo.swath_geo import (
SingleSwathBistatic,
SingleSwath,
line_of_sight,
)
This diff is collapsed.
This diff is collapsed.
import numpy as np
import drama.constants as cnst
from drama.geo import geometry as sargeo
__author__ = "Paco Lopez Dekker"
__email__ = "F.LopezDekker@tudeft.nl"
class Baseline:
"""A class to convert between baselines
Parameters
----------
b :
Baseline: if it is a two element list then first element is
the horizontal baseline, and the second is the vertical
(radial) one. Else, it is the perpendicular baseline, defaults
to False
theta_i :
incident angle or array of incident angles, in degrees
h_orb :
height of sensor, in m.
f0 :
reasonant frequency, defaults to 5.4 GHz
Returns
-------
"""
def __init__(self, b, theta_i, h_orb, bistatic=False, f0=5.4e9):
"""Initialise Baseline.
"""
# Look angle
self.sargeo = sargeo.QuickRadarGeometry(h_orb)
self.theta_i = np.radians(theta_i)
self.theta_l = self.sargeo.inc_to_look(self.theta_i)
self.h_orb = h_orb
self.f0 = f0
if type(b) is tuple:
self.b_h = b[0]
self.b_r = b[1]
self.b_p = self.b_h * np.cos(self.theta_l) + self.b_r * np.sin(
self.b_r
)
else:
self.b_p = b
self.b_h = None
self.b_r = None
self.sr = self.sargeo.look_to_sr(self.theta_l)
self.bistfact = 1 if bistatic else 2
def kz(self, delta_grg=0):
"""
Parameters
----------
delta_grg :
Shift of target in ground range, useful for ice (Default value = 0)
Returns
-------
"""
# phase_int = phase_m - phase_s = -2 k0 (R_m - R_s) = 2 k0 (R_s-R_m)
# positive phase means R_s > R_m -> lower height
k0 = 2 * np.pi * self.f0 / cnst.c
# Modify the effective baseline. Here we take a coordinate system with the vertical direction
# perpendicular to the surface
d_b_p = -delta_grg * np.cos(self.theta_i)
b_p = self.b_p + d_b_p
kz = -self.bistfact * k0 * b_p / self.sr / np.sin(self.theta_i)
return kz
def dkz(self, delta_grg=0):
"""
Parameters
----------
delta_grg :
Shift of target in ground range, useful for ice (Default value = 0)
Returns
-------
"""
# phase_int = phase_m - phase_s = -2 k0 (R_m - R_s) = 2 k0 (R_s-R_m)
# positive phase means R_s > R_m -> lower height
k0 = 2 * np.pi * self.f0 / cnst.c
# Modify the effective baseline. Here we take a coordinate system with the vertical direction
# perpendicular to the surface
d_b_p = -delta_grg * np.cos(self.theta_i)
dkz = -self.bistfact * k0 * d_b_p / self.sr / np.sin(self.theta_i)
return dkz
def set_kz(self, kz):
"""
Parameters
----------
kz :
Returns
-------
"""
k0 = 2 * np.pi * self.f0 / cnst.c
self.b_p = -kz * self.sr * np.sin(self.theta_i) / self.bistfact / k0
This diff is collapsed.
This diff is collapsed.
import os
import unittest
import numpy as np
import xarray as xr
from drama.io import cfg
import drama.geo as sargeo
class TestSwathGeo(unittest.TestCase):
def setUp(self):
self.dau = np.linspace(-400e3, 400e3, 17)
self.inc_m = np.linspace(20, 50, 31)
par_dir = os.path.join(os.path.expanduser("~/Code/stereoid"), "PAR")
runid = "2019_1"
self.par_file = os.path.join(par_dir, ("Hrmny_%s.cfg" % runid))
self.par_file_second = os.path.join(
par_dir, ("Hrmny_%s_orbit_later.cfg" % runid)
)
conf = cfg.ConfigFile(self.par_file)
form_id = conf.formation.id
drama_dir = os.path.expanduser("~/Code/drama/drama")
save_dir = os.path.join(
os.path.join(os.path.join(drama_dir, "geo"), "test_data"), runid
)
self.save_dir = os.path.join(save_dir, form_id)
def test_single_swath_bistatic(self):
# Attention! This dataset was generated using the harmony_ati
# formation config, only compare it to that formation
sample_ds = xr.open_dataset(os.path.join(self.save_dir, "bist_geo.nc"))
inc_s = np.zeros((17, 31))
bist_ang_az = np.zeros_like(inc_s)
for dau_i in range(17):
swth_bst = sargeo.SingleSwathBistatic(
par_file=self.par_file, dau=self.dau[dau_i]
)
inc_s[dau_i] = np.degrees(
swth_bst.inc2slave_inc(np.radians(self.inc_m))
)
bist_ang_az[dau_i] = np.degrees(
swth_bst.inc2bist_ang_az(np.radians(self.inc_m))
)
inc_m_xr = xr.DataArray(
self.inc_m,
dims=["inc"],
name="inc_m",
coords={"inc": self.inc_m.astype(np.int)},
)
inc_s_xr = xr.DataArray(
inc_s,
dims=["dau", "inc"],
name="inc_s",
coords={
"inc": self.inc_m.astype(np.int),
"dau": (self.dau / 1e3).astype(np.int),
},
)
bist_ang_az_xr = xr.DataArray(
bist_ang_az,
dims=["dau", "inc"],
name="bist_ang_az",
coords={
"inc": self.inc_m.astype(np.int),
"dau": (self.dau / 1e3).astype(np.int),
},
)
dau_xr = xr.DataArray(
self.dau,
dims=["dau"],
name="dau",
coords={"dau": (self.dau / 1e3).astype(np.int)},
)
geo = xr.Dataset(
{
"inc_s1": inc_m_xr,
"inc_hrmny": inc_s_xr,
"bist_ang_az": bist_ang_az_xr,
"along_track_distance": dau_xr,
}
)
self.assertIsNone(
xr.testing.assert_allclose(geo, sample_ds), "datasets do not match"
)
def test_single_swath_bistatic_iterator(self):
swth_bst = sargeo.SingleSwathBistatic(
par_file=self.par_file, n_orbits=2
)
swt_bst_later = sargeo.SingleSwathBistatic(
par_file=self.par_file_second
)
swath_list = []
for swath in swth_bst:
swath_list.append(swath)
self.assertIsNone(
np.testing.assert_allclose(
swath_list[1].master_swath.incident, swt_bst_later.master_inc
),
"incident angles are the same",
)
def test_single_swath_bistatic_iterator_first(self):
swth_bst = sargeo.SingleSwathBistatic(
par_file=self.par_file, n_orbits=2
)
swt_bst_later = sargeo.SingleSwathBistatic(
par_file=self.par_file_second
)
swath_list = []
for swath in swth_bst:
swath_list.append(swath)
self.assertIsNone(
np.testing.assert_allclose(
swath_list[0].master_swath.incident, swth_bst.master_inc
),
"incident angles are the same",
)
def test_single_swath_bistatic_iterator_different(self):
swth_bst = sargeo.SingleSwathBistatic(
par_file=self.par_file, n_orbits=2
)
swt_bst_later = sargeo.SingleSwathBistatic(
par_file=self.par_file_second
)
swath_list = []
for swath in swth_bst:
swath_list.append(swath)
self.assertTrue(
np.any(
np.not_equal(
swath_list[0].master_swath.incident,
swt_bst_later.master_inc,
)
),
"incident angles are the same",
)
"""timing groups together functions to compute a suitable timing configuration for a given swath.
It takes into account the tilt of the antenna and its size."""
from collections import namedtuple
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
from drama import constants as const
from drama.geo import geometry as sargeo
blind_out = namedtuple(
"BlindRanges",
["inc_s", "inc_e", "look_s", "look_e", "slant_range_s", "slant_range_e"],
)
def blind_ranges(
angular_range, orbit_h, PRF, duty_cycle, angle_is_look=False, radian=False
):
"""Code to find the positions of the blind ranges
:author: Paco Lopez-Dekker
Parameters
----------
angular_range :
tuple with start and end angle (typically,
indicent angle, in degree)
orbit_h :
orbit height
PRF :
PRF
duty_cycle :
a number between 0 and 0.5
angle_is_look :
set to True if angular range are look angles.
Defaults to False
radian :
set to True if input and output angles are in radians (Default value = False)
Returns
-------
type
a named tuple with the start and end of the blind incident
and look angles, and the start of the blind slant ranges
"""
if radian:
ang_range = np.array(angular_range)
else:
ang_range = np.array(angular_range) * np.pi / 180
if angle_is_look:
inc_range = sargeo.look_to_inc(ang_range, orbit_h)
else:
inc_range = ang_range
inc_v = np.arange(
inc_range[0], inc_range[1], (inc_range[1] - inc_range[0]) / 1000
)
sr_v = sargeo.inc_to_sr(inc_v, orbit_h)
# Find out blind ranges
delta_r_unamb = const.c / 2 / PRF
blind_ranges = np.arange(0, np.max(sr_v), delta_r_unamb)
good_indices = np.where(blind_ranges > np.min(sr_v))
blind_ranges = blind_ranges[good_indices]
if blind_ranges.size == 0:
print("No blind ranges found")
return 0
delta_r_duty = const.c / 2 * (2 * duty_cycle / PRF)
r2inc = interpolate.InterpolatedUnivariateSpline(sr_v, inc_v, k=2)
blind_inc_1 = r2inc(blind_ranges)
blind_inc_2 = r2inc(blind_ranges + delta_r_duty)
res = blind_out(
blind_inc_1 * 180 / np.pi,
blind_inc_2 * 180 / np.pi,
sargeo.inc_to_look(blind_inc_1, orbit_h) * 180 / np.pi,
sargeo.inc_to_look(blind_inc_2, orbit_h) * 180 / np.pi,
blind_ranges,
blind_ranges + delta_r_duty,
)
return res
def timing_diagram(
angular_range,
orbit_h,
PRFs,
duty_cycle,
angle_is_look=False,
radian=False,
nadir_echo_len=10e-6,
plot_tx=True,
plot_echo_incomplete=True,
plot_nadir=True,
):
"""Code to find the positions of the blind ranges
:author: Paco Lopez-Dekker
Parameters
----------
angular_range :
tuple with start and end angle (typically,
indicent angle, in degree)
orbit_h :
orbit height
PRFs :
a vector of PRF values
duty_cycle :
a number between 0 and 0.5
angle_is_look :
set to True if angular range are look angles.
Defaults to False
radian :
set to True if input and output angles are in radians (Default value = False)
nadir_echo_len :
(Default value = 10e-6)
plot_tx :
(Default value = True)
plot_echo_incomplete :
(Default value = True)
plot_nadir :
(Default value = True)
Returns
-------
type
a named tuple with the start and end of the blind incident
and look angles, and the start of the blind slant ranges
"""
if radian:
ang_range = np.array(angular_range)
else:
ang_range = np.array(angular_range) * np.pi / 180
if angle_is_look:
inc_range = sargeo.look_to_inc(ang_range, orbit_h)
else:
inc_range = ang_range
inc_v = np.linspace(inc_range[0], inc_range[1], 1000)
inc_v = inc_v.reshape((inc_v.size, 1))
sr_v = sargeo.inc_to_sr(inc_v, orbit_h)
# Find out blind ranges
PRF = np.array(PRFs)
PRF = PRF.reshape((1, PRF.size))
delta_r_duty = const.c / 2 * (duty_cycle / PRF)
delta_r = delta_r_duty
delta_r_unamb = const.c / 2 / PRF
rgwin_min = np.floor(sr_v.min() / delta_r_unamb.max()).astype(np.int)
rgwin_max = np.floor(sr_v.max() / delta_r_unamb.min()).astype(np.int)
nwin = rgwin_max - rgwin_min + 1
sr_win_start = delta_r_unamb * np.arange(rgwin_min, rgwin_max + 1).reshape(
(nwin, 1)
)
sr_win_end = sr_win_start + delta_r
sr_echowin_end = sr_win_start - delta_r
(gr_win_start, theta_i_win_start, theta_l_win_start, kk) = sargeo.sr_to_geo(
sr_win_start, orbit_h
)
(gr_win_end, theta_i_win_end, theta_l_win_end, kk) = sargeo.sr_to_geo(
sr_win_end, orbit_h
)
(
gr_echowin_end,
theta_i_echowin_end,
theta_l_echowin_end,
kk,
) = sargeo.sr_to_geo(sr_echowin_end, orbit_h)
# sr_v_amb = np.mod(sr_v, delta_r_unamb)
# blind_mask = np.less(sr_v_amb, delta_r_duty)
# Nadir return
sr_nadir = orbit_h
nnadir = np.floor((sr_v.max() - sr_nadir) / delta_r_unamb.min() + 1).astype(
np.int
)
sr_nadir_start = sr_nadir + delta_r_unamb * np.arange(nnadir).reshape(
(nnadir, 1)
)
sr_nadir_end = sr_nadir_start + nadir_echo_len
(
gr_nadir_start,
theta_i_nadir_start,
theta_l_nadir_start,
kk,
) = sargeo.sr_to_geo(sr_nadir_start, orbit_h)
(gr_nadir_end, theta_i_nadir_end, theta_l_nadir_end, kk) = sargeo.sr_to_geo(
sr_nadir_end, orbit_h
)
if angle_is_look:
theta_win_start = theta_l_win_start
theta_win_end = theta_l_win_end
theta_echowin_end = theta_l_echowin_end
theta_nadir_start = theta_l_nadir_start
theta_nadir_end = theta_l_nadir_end
ylabel = "Look angle [deg]"
else:
theta_win_start = theta_i_win_start
theta_win_end = theta_i_win_end
theta_echowin_end = theta_i_echowin_end
theta_nadir_start = theta_i_nadir_start
theta_nadir_end = theta_i_nadir_end
ylabel = "Incident angle [deg]"
if plot_tx or plot_echo_incomplete or plot_nadir:
plt.figure()
for ind in range(nwin):
if plot_tx:
plt.fill_between(
PRF.flatten(),
np.degrees(theta_win_start[ind]),
np.degrees(theta_win_end[ind]),
facecolor="red",