# Computing Doppler delay on a meshgrid

Objective

Plot the contour of the iso-doppler and iso-delay lines for a transmitter-receiver reflection on a specular plane.

Implementation

This Doppler shift can be expressed as follows: $$f_{D,0}(\vec{r},t_0) = [\vec{V_t} \cdot \vec{m}(\vec{r},t_0) - \vec{V_r} \cdot \vec{n}(\vec{r},t_0)]/\lambda$$

where for a given time $$\t_0\$$, $$\\vec{m}\$$ is the reflected unit vector, $$\\vec{n}\$$ is incident unit vector, $$\\vec{V_t}\$$ is the velocity of the transmitter, $$\\vec{V_r}\$$ is the velocity of the receiver, and $$\\lambda\$$ is the wavelength of the transmitted electromagnetic wave.

The time delay of the electromagnetic wave is just the path it travels divided by the speed of light, assuming vacuum propagation.

#!/usr/bin/env python

import scipy.integrate as integrate
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

h_t = 20000e3 # meters
h_r = 500e3 # meters

# Coordinate Frame as defined in Figure 2
#      J. F. Marchan-Hernandez, A. Camps, N. Rodriguez-Alvarez, E. Valencia, X.
#      Bosch-Lluis, and I. Ramos-Perez, “An Efficient Algorithm to the Simulation of
#      Delay–Doppler Maps of Reflected Global Navigation Satellite System Signals,”
#      IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 8, pp.
#      2733–2740, Aug. 2009.
r_t = np.array([0,h_t/np.tan(elevation),h_t])
r_r = np.array([0,-h_r/np.tan(elevation),h_r])

# Velocity
v_t = np.array([2121, 2121, 5]) # m/s
v_r = np.array([2210, 7299, 199]) # m/s

light_speed = 299792458 # m/s
delay_chip =  1/1.023e6 # seconds

# GPS L1 center frequency is defined in relation to a reference frequency
# f_0 = 10.23e6, so that f_carrier = 154*f_0 = 1575.42e6 # Hz
# Explained in section 'DESCRIPTION OF THE EMITTED GPS SIGNAL' in Zarotny
# and Voronovich 2000
f_0 = 10.23e6 # Hz
f_carrier = 154*f_0;

def doppler_shift(r):
'''
Doppler shift as a contribution of the relative motion of transmitter and
receiver as well as the reflection point.

Implements Equation 14
V. U. Zavorotny and A. G. Voronovich, “Scattering of GPS signals from
the ocean with wind remote sensing application,” IEEE Transactions on
Geoscience and Remote Sensing, vol. 38, no. 2, pp. 951–964, Mar. 2000.
'''
wavelength = light_speed/f_carrier
f_D_0 = (1/wavelength)*(
np.inner(v_t, incident_vector(r)) \
-np.inner(v_r, reflection_vector(r))
)
#f_surface = scattering_vector(r)*v_surface(r)/2*pi
f_surface = 0
return f_D_0 + f_surface

def doppler_increment(r):
return doppler_shift(r) - doppler_shift(np.array([0,0,0]))

def reflection_vector(r):
reflection_vector = (r_r - r)
reflection_vector_norm = np.linalg.norm(r_r - r)
reflection_vector[0] /= reflection_vector_norm
reflection_vector[1] /= reflection_vector_norm
reflection_vector[2] /= reflection_vector_norm
return reflection_vector

def incident_vector(r):
incident_vector = (r - r_t)
incident_vector_norm = np.linalg.norm(r - r_t)
incident_vector[0] /= incident_vector_norm
incident_vector[1] /= incident_vector_norm
incident_vector[2] /= incident_vector_norm
return  incident_vector

def time_delay(r):
path_r = np.linalg.norm(r-r_t) + np.linalg.norm(r_r-r)
path_specular = np.linalg.norm(r_t) + np.linalg.norm(r_r)
return (1/light_speed)*(path_r - path_specular)

# Plotting Area

x_0 =  -100e3 # meters
x_1 =  100e3 # meters
n_x = 500

y_0 =  -100e3 # meters
y_1 =  100e3 # meters
n_y = 500

x_grid, y_grid = np.meshgrid(
np.linspace(x_0, x_1, n_x),
np.linspace(y_0, y_1, n_y)
)

r = [x_grid, y_grid, 0]
z_grid_delay = time_delay(r)/delay_chip
z_grid_doppler = doppler_increment(r)

delay_start = 0 # C/A chips
delay_increment = 0.5 # C/A chips
delay_end = 15 # C/A chips
iso_delay_values = list(np.arange(delay_start, delay_end, delay_increment))

doppler_start = -3000 # Hz
doppler_increment = 500 # Hz
doppler_end = 3000 # Hz
iso_doppler_values = list(np.arange(doppler_start, doppler_end, doppler_increment))

fig_lines, ax_lines = plt.subplots(1,figsize=(10, 4))
contour_delay = ax_lines.contour(x_grid, y_grid, z_grid_delay, iso_delay_values, cmap='winter')
fig_lines.colorbar(contour_delay, label='C/A chips', )

contour_doppler = ax_lines.contour(x_grid, y_grid, z_grid_doppler, iso_doppler_values, cmap='winter')
fig_lines.colorbar(contour_doppler, label='Hz', )

ticks_y = ticker.FuncFormatter(lambda y, pos: '{0:g}'.format(y/1000))
ticks_x = ticker.FuncFormatter(lambda x, pos: '{0:g}'.format(x/1000))
ax_lines.xaxis.set_major_formatter(ticks_x)
ax_lines.yaxis.set_major_formatter(ticks_y)
plt.xlabel('[km]')
plt.ylabel('[km]')

plt.show()


Which produces this presumably right output:

Questions

In order to compute the incident vector from a point $$\r_t\$$ I've implemented the following code:

def incident_vector(r):
incident_vector = (r - r_t)
incident_vector_norm = np.linalg.norm(r - r_t)
incident_vector[0] /= incident_vector_norm
incident_vector[1] /= incident_vector_norm
incident_vector[2] /= incident_vector_norm
return  incident_vector


This works perfectly fine, but I think there must be a cleaner way to write this. I would like to write something like this:

def incident_vector(r):
return (r - r_t)/np.linalg.norm(r - r_t)


But unfortunately it doesn't work with the meshgrid, as it doesn't know how to multiply the scalar grid with the vector grid:

ValueError: operands could not be broadcast together with shapes (3,) (500,500)