Умножение np.linalg.norm с сеткой не работает

Question

Я пытаюсь построить следующие функции

#!/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;

integration_time = 1e-3 # seconds
delay_chip =  1/1.023e6 # seconds

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 reflection_vector(r):
    return (r_r - r)/np.linalg.norm(r_r - r)

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

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_shift(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 = 0 # Hz
doppler_increment = 500 # Hz
doppler_end = 4000 # 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_delay = 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()

Каждая функция работает индивидуально, и контурный график для функции time_delay также работает. Но когда я пытаюсь построить функцию doppler_shift, деление матрицы на вывод linalg.norm обрывается:

ValueError: операнды не могут передаваться вместе с фигурами (3,) (500,500)

Как мне написать так, чтобы он работал с сеткой?