In [1]:

%matplotlib inline

In [2]:

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as apu
from pycraf import atm
from pycraf import conversions as cnv

WARNING: AstropyDeprecationWarning: astropy.utils.compat.funcsigs is now deprecated - use inspect instead [astropy.utils.compat.funcsigs]

In [3]:

obs_alt_m = 100.1
earth_radius = atm.EARTH_RADIUS
freq_grid = [1] * apu.GHz
atm_layers_cache = atm.atm_layers(freq_grid, atm.profile_highlat_winter)

layer_angles = np.linspace(-0.001, 0.02, 361)
radii = atm_layers_cache['radii']
heights = atm_layers_cache['heights']
ref_index = atm_layers_cache['ref_index']


plt.close()
fig = plt.figure(figsize=(14, 8))
for r in earth_radius + heights:
    plt.plot(r * np.sin(layer_angles), r * np.cos(layer_angles) - earth_radius, 'k--', alpha=0.5)

for elevation in np.linspace(-0.04, -0.01, 21):
# for elevation in [0.88]:
    path_params, _, _ = atm.atm._prepare_path(
        elevation, obs_alt_m / 1000., radii, heights, ref_index
        # max_arc_length=0.134,
        # max_path_length=15.,
        )
    plt.plot(path_params.x_n, path_params.y_n - earth_radius, '-')

plt.xlim((-0.1, 15.1))
plt.ylim(((obs_alt_m - 12) / 1000, (obs_alt_m + 1) / 1000))

plt.title('Path propagation through layered atmosphere')
plt.xlabel('Projected distance (km)')
plt.ylabel('Height above ground (km)')
# plt.savefig('/tmp/path_prop.jpg', dpi=300)
plt.show()

In [4]:

np.degrees(path_params.delta_n)[-1]

Out[4]:

8.933501565830346

In [5]:

n1, n2 = 1.0003089446689986, 1.0003089007075814
alpha_n = np.arange(0, 90, 1)
beta_n = np.degrees(np.arcsin(n2 / n1 * np.sin(np.radians(alpha_n))))
plt.plot(alpha_n, alpha_n - beta_n, 'rx')
beta_n = np.degrees(np.arcsin(n1 / n2 * np.sin(np.radians(alpha_n))))
plt.plot(alpha_n, alpha_n - beta_n, 'bx')
plt.show()

In [6]:

# obs_alt_m = 1000.  # m
elevations = np.arange(-0.5, 0.25, 0.0001)
h_ends, refractions = [], []
for elevation in elevations:
    ret = atm.atm._path_endpoint(
        elevation, obs_alt_m / 1000., radii, heights, ref_index,
        # max_path_length=30,
        max_arc_length=1,
        )
    (
        a_n, r_n, h_n, x_n, y_n, alpha_n, beta_n, delta_n, layer_idx,
        refraction,
        is_space_path,
        ) = ret    
    h_ends.append(h_n)
    refractions.append(refraction)

In [24]:

fig = plt.figure(figsize=(14, 8))
plt.plot(elevations, h_ends, '-')
plt.xlim((-0.3, 0.05))
plt.ylim((-0.01, 0.9))
plt.show()

In [8]:

fig = plt.figure(figsize=(14, 8))
plt.plot(elevations, refractions, '-')
plt.show()

In [14]:

obs_alt_m = 100.1
target_alt_m = 810
arc_length = 1.
elev_final, h_final = atm.atm._find_elevation(
    2., obs_alt_m / 1000, target_alt_m / 1000, arc_length,
    radii, heights, ref_index, niter=150, interval=25,
    )

In [15]:

elev_final, h_final

Out[15]:

(-0.019485916559211967, 0.8186206445152493)

In [16]:

fig = plt.figure(figsize=(14, 8))
plt.plot(elevations, h_ends, '-')
plt.plot(elev_final, h_final, 'x')
plt.xlim((-0.05, 0.05))
plt.ylim((0.7, 0.9))
plt.show()

In [ ]:

obs_alt_m = 100.1
target_alt_m = 112.

plt.close()
fig = plt.figure(figsize=(14, 8))
for r in earth_radius + heights:
    plt.plot(r * np.sin(layer_angles), r * np.cos(layer_angles) - earth_radius, 'k--', alpha=0.5)

# for elevation in np.linspace(-0.04, -0.01, 21):
for elevation in [-0.01]:
    path_params, _, _ = atm.atm._prepare_path(
        elevation, obs_alt_m / 1000., radii, heights, ref_index,
        max_arc_length=0.134,
        # max_path_length=15.,
        )
    plt.plot(path_params.x_n, path_params.y_n - earth_radius, '-', lw=5)

    # backpropagation
    elevation_r = 90 - (180. - np.degrees(path_params.alpha_n[-1]))
    target_alt_r = path_params.h_n[-1]
    path_params_r, _, _ = atm.atm._prepare_path(
        elevation_r, target_alt_r, radii, heights, ref_index, 
        max_arc_length=0.134,
        # max_path_length=15.,
        )
    x_n_r, y_n_r = path_params_r.x_n, path_params_r.y_n
    delta_n_r = path_params_r.delta_n
    ddelta_n = -delta_n_r + delta_n_r[-1] - delta_n_r
    x_n_r, y_n_r = (
        np.cos(ddelta_n) * x_n_r + np.sin(ddelta_n) * y_n_r,
        -np.sin(ddelta_n) * x_n_r + np.cos(ddelta_n) * y_n_r,
        )
    plt.plot(x_n_r, y_n_r - earth_radius, '-')
    
plt.xlim((-0.1, 15.1))
plt.ylim(((obs_alt_m - 12) / 1000, (obs_alt_m + 1) / 1000))

plt.title('Path propagation through layered atmosphere')
plt.xlabel('Projected distance (km)')
plt.ylabel('Height above ground (km)')
# plt.savefig('/tmp/path_prop.jpg', dpi=300)
plt.show()

In [ ]:

obs_alt_r

In [ ]:

path_params.alpha_n, path_params.h_n

In [ ]:

elevation = 90 - (180. - np.degrees(path_params.alpha_n[-1]))
elevation, np.degrees(path_params.alpha_n[-1])

In [ ]:

elevation = 9  # deg
obs_alt = 1000  # m
path_params, refraction, _, _ = atm.atm._prepare_path(
    elevation, obs_alt, atm.profile_highlat_winter, max_path_length=50.
    )
ret = atm.atm._path_endpoint(
    elevation, obs_alt, atm.profile_highlat_winter, max_path_length=50.
    )

path_params, refraction, is_space_path, weather = atm.atm._prepare_path(
    elevation, obs_alt, atm.profile_highlat_winter, max_path_length=5.
    )

In [ ]:

%timeit _ = atm.atm._prepare_path(90, 10, atm.profile_highlat_winter, max_path_length=50.)
# 235 µs ± 2.42 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit _ = atm.atm._prepare_path(-90, 50100, atm.profile_highlat_winter, max_path_length=50.)
# 217 µs ± 1.36 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit _ = atm.atm._prepare_path(1, 100, atm.profile_highlat_winter, max_path_length=50.)
# 246 µs ± 5.18 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit _ = atm.atm._prepare_path(-2, 3000, atm.profile_highlat_winter, max_path_length=50.)
# 182 µs ± 2.28 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit _ = atm.atm._path_endpoint(-2, 3000, atm.profile_highlat_winter, max_path_length=50.)
# 40.7 µs ± 1.36 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)