Tukukkk

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:

where for a given time t_0, m is the reflected unit vector, n is incident unit vector, v_t is the velocity of the transmitter, 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

elevation = 60*np.pi/180 # rad



# 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



# 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:

Please feel free to provide recommendations about the implementation and style.

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)