Métodos numéricos de propagación submarina

La pérdida de propagación en forma cerrada de la propagación submarina del sonido no sabe nada del perfil de velocidad del sonido, del fondo ni de la superficie. Cuando la refracción y los contornos deciden la respuesta, el campo hay que calcularlo: esta guía cubre los cuatro métodos numéricos del módulo underwater, la física que cada uno discretiza, y cómo elegir entre ellos y las formas cerradas. Los cuatro toman el mismo perfil $c(z)$ independiente de la distancia (estratificado horizontalmente) y una superficie de presión nula, y siguen a Jensen, Kuperman, Porter y Schmidt, Computational Ocean Acoustics; los dos métodos basados en rayos pueden además inclinar el fondo con un perfil de profundidad lineal a trozos (sección 3), la única dependencia de la distancia del módulo.

• normal_modes resuelve el problema de autovalores de Sturm-Liouville separado en profundidad por diferencias finitas y suma los modos propagantes en la pérdida de propagación. Validado contra los modos exactos de la guía de ondas ideal (de presión nula).
• ray_trace integra las ecuaciones de trayectoria de rayos (Runge-Kutta, vectorizadas sobre todos los rayos a la vez) a través de un perfil de velocidad del sonido, reflejando en la superficie y el fondo, y arrastra el tiempo de propagación de cada rayo como un estado más de la misma integración. Validado contra los arcos de círculo de un gradiente lineal y contra la forma cerrada del tiempo de propagación a lo largo de ellos. eigenrays busca después en ese abanico trazado los caminos que unen la fuente con un receptor y lista sus llegadas, validado contra la red de imágenes de la guía de ondas ideal, donde cada llegada es forma cerrada.
• gaussian_beams cuelga un haz gaussiano de cada uno de esos rayos y los suma en un campo de pérdida de propagación, que se mantiene finito en las cáusticas donde la amplitud clásica de un rayo es infinita y decae en las zonas de sombra donde no está definida. Validado contra la divergencia esférica en campo libre, el campo de dos rayos del espejo de Lloyd y la suma de fuentes imagen de la guía de ondas ideal.
• parabolic_equation hace avanzar en distancia la PE estándar (Tappert) con el algoritmo split-step de Fourier. Validado contra la divergencia esférica en campo libre; coincide en tendencia con la pérdida de propagación de los modos normales.

Tres paneles, tres preguntas distintas sobre la misma agua. A la izquierda, rayos en un perfil de Munk cuyo mínimo es de 1500,0 m/s a 1271 m frente a los 1548,5 m/s de la superficie: 48 m/s de diferencia bastan para plegar sobre sí mismo el abanico de ±12° y formar zonas de convergencia, y el rayo de ±12° todavía barre de 143 m a 3997 m de profundidad. En el centro, el mismo entorno como campo de la PE, donde las zonas aparecen como interferencia y no como líneas. A la derecha, los dos métodos sobre una guía de ondas distinta, somera, de 200 m: siguen la misma tendencia con un desvío medio de 1,6 dB, pero la diferencia distancia a distancia llega a 25 dB en los mínimos de interferencia, que es la manera honesta de decir qué significa «coincidir en tendencia».

import matplotlib.pyplot as plt
import numpy as np
from phonometry import underwater

# Un perfil de Munk de velocidad del sonido en aguas profundas, para los
# paneles de rayos y de PE.
z = np.linspace(0.0, 5000.0, 60)
eta = 2.0 * (z - 1300.0) / 1300.0
c = 1500.0 * (1.0 + 0.00737 * (eta - 1.0 + np.exp(-eta)))

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) rayos a través del perfil de Munk, hasta 100 km.
underwater.ray_trace(z, c, source_depth=1000.0,
                     launch_angles_deg=np.linspace(-12.0, 12.0, 13),
                     max_range=100e3, n_steps=6000).plot(ax=axes[0], language="es")

# (b) el campo de la PE split-step de Fourier sobre el mismo entorno.
underwater.parabolic_equation(50.0, z, c, source_depth=1000.0,
                              max_range=100e3, range_step=25.0,
                              n_depth_points=1024).plot(ax=axes[1], language="es")

# (c) modos frente a PE en una guía de ondas *distinta*, somera, de 200 m.
r = np.linspace(100.0, 20_000.0, 400)
nm = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1530.0],
                             source_depth=30.0, receiver_depth=120.0,
                             ranges_m=r, n_depth_points=800)
pe = underwater.parabolic_equation(50.0, [0.0, 200.0], [1500.0, 1530.0],
                                   source_depth=30.0, max_range=20e3,
                                   range_step=20.0, n_depth_points=512)
zi = int(np.argmin(np.abs(pe.depths - 120.0)))
axes[2].plot(nm.ranges / 1000.0, nm.propagation_loss, label="Modos normales")
axes[2].plot(pe.ranges / 1000.0, pe.propagation_loss[zi], label="PE")
axes[2].invert_yaxis()
axes[2].set(xlabel="Distancia [km]", ylabel="Pérdida de propagación [dB]")
axes[2].legend()
plt.show()

import numpy as np
from phonometry import underwater

# Un perfil de Munk de aguas profundas.
z = np.linspace(0.0, 5000.0, 60)
eta = 2.0 * (z - 1300.0) / 1300.0
c = 1500.0 * (1.0 + 0.00737 * (eta - 1.0 + np.exp(-eta)))

rays = underwater.ray_trace(z, c, source_depth=1000.0,
                    launch_angles_deg=np.linspace(-12.0, 12.0, 21), max_range=100e3)
rays.plot()   # trayectorias de rayos / zonas de convergencia

# Guía de ondas somera isovelocidad: modos y PE.
modes = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1500.0],
                        source_depth=50.0, receiver_depth=100.0)
print(modes.wavenumbers.size, "modos propagantes")
field = underwater.parabolic_equation(50.0, [0.0, 200.0], [1500.0, 1500.0],
                              source_depth=50.0, max_range=20e3)
field.plot()  # campo PL sobre distancia x profundidad

normal_modes devuelve un NormalModeResult (wavenumbers, mode_functions, propagation_loss); ray_trace un RayTraceResult (ranges, depths, travel_times, arc_lengths y los recuentos de reflexiones por contorno, por rayo); gaussian_beams un GaussianBeamResult (el campo propagation_loss, más el rayo central y la anchura de cada haz); parabolic_equation un ParabolicEquationResult (el campo propagation_loss). Los cuatro difieren en lo que aproximan, no en lo que modelan, así que elegir entre ellos es una cuestión de régimen: la tabla de la sección 6 es la respuesta corta.

En un océano estratificado horizontalmente la ecuación de Helmholtz se separa en coordenadas cilíndricas, $p(r,z)=\Phi (r)\Psi (z)$, y el factor de profundidad obedece un problema de autovalores de Sturm-Liouville (Jensen, Ec. 5.3):

con superficie de presión nula en $z=0$ y, en el fondo $z=D$, $\Psi (D)=0$ para un lecho de presión nula o $d\Psi /dz{|}_D=0$ para uno rígido. Cada autofunción ${\Psi }_m(z)$ es una onda estacionaria en profundidad que viaja en distancia como $e^{i{k}_{rm}r}$ con su propio número de onda horizontal $k_{rm}$; solo propagan los modos con $k_{rm}$ real, el resto son evanescentes y se extinguen al cabo de unas pocas profundidades de agua. El campo es la suma modal (Ec. 5.14),

cada modo ponderado por su excitación en la profundidad de la fuente ${\Psi }_m(z_{\mathrm{s}})$ y su amplitud en la del receptor ${\Psi }_m(z)$, y la pérdida de propagación coherente sale como $PL=-20{\log }_{10}|p(r,z)/p_0(1\mathrm{m})|$ (Ec. 5.15). normal_modes discretiza la ecuación en profundidad por diferencias finitas (un problema de autovalores tridiagonal simétrico) sobre una malla lo bastante fina para mantener honestos los autovalores cercanos al corte, y avisa cuando un modo retenido queda demasiado cerca de su banda de error de discretización.

El número de modos es la física: un canal isovelocidad de profundidad $D$ transporta $M\approx kD/\pi$ modos propagantes, así que un canal de 200 m a 50 Hz transporta 13 y el mismo canal a 5 kHz transportaría más de 1300, que es la frecuencia a la que los rayos pasan a ser la descripción económica. Baja frecuencia y aguas someras significan pocos modos y una descripción compacta y esencialmente exacta. La guía de ondas ideal de presión nula es el oráculo de validación, con $k_{rm}=\sqrt{k^2-(m\pi /D{)}^2}$ en forma cerrada:

import numpy as np
from phonometry import underwater

# Un canal isovelocidad de 200 m a 50 Hz: kD/pi = 2 f D / c = 13,3.
modes = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1500.0],
                                source_depth=50.0, receiver_depth=100.0)
print(modes.wavenumbers.size)                   # 13  modos propagantes
k = 2 * np.pi * 50.0 / 1500.0
print(round(float(modes.wavenumbers[0]), 5))    # 0.20885  kr1 calculado
print(round(np.sqrt(k**2 - (np.pi / 200.0) ** 2), 5))   # 0.20885  exacto

El modo $m$ tiene $m-1$ ceros interiores, y la profundidad de la fuente decide qué modos existen siquiera en el campo: el cuarto modo tiene un cero justo a 50 m en un canal de 200 m, así que una fuente ahí no lo excita. La escalera del centro es el número de modos entrando en propagación de uno en uno conforme sube la frecuencia, siguiendo $M\approx kD/\pi$; el panel de la derecha es la interferencia de los 13 modos que sobreviven a 50 Hz. Un resultado calculado dibuja los tres con modes.plot().

import matplotlib.pyplot as plt
import numpy as np

fig, axes = plt.subplots(1, 3, figsize=(16, 5.4))
r = np.linspace(100.0, 20_000.0, 500)
res = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1500.0],
                              source_depth=50.0, receiver_depth=100.0,
                              ranges_m=r, n_depth_points=800)
for m in range(4):
    axes[0].plot(res.mode_functions[m], res.mode_depths, label=f"m = {m + 1}")
axes[0].axhline(50.0, linestyle="--")
axes[0].invert_yaxis()
axes[0].set(xlabel="Función modal Psi_m(z)", ylabel="Profundidad [m]")
axes[0].legend()

freqs = np.arange(10.0, 201.0, 1.0)
counts = [underwater.normal_modes(f, [0.0, 200.0], [1500.0, 1500.0],
                                  source_depth=50.0, receiver_depth=100.0,
                                  n_depth_points=1200).wavenumbers.size
          for f in freqs]
axes[1].step(freqs, counts, where="post", label="modos devueltos")
axes[1].plot(freqs, 2.0 * freqs * 200.0 / 1500.0, "--", label="M = kD/pi")
axes[1].set(xlabel="Frecuencia [Hz]", ylabel="Modos propagantes M")
axes[1].legend()

res.plot(ax=axes[2], language="es")
plt.show()

Esa guía de ondas ideal cumple dos papeles, y no conviene confundirlos. Es el oráculo exacto contra el que se comprueba el cálculo, que es para lo que se usa aquí, y es además un fondo físicamente extremo: un lecho de presión nula refleja con módulo unidad e inversión de fase, así que no transmite nada al sedimento y no absorbe nada. Una pérdida de propagación modal o PE calculada sobre él subestima, por tanto, la pérdida de cualquier fondo marino real, y la diferencia se ensancha con el alcance porque se acumula rebote a rebote. La opción rígida es el extremo contrario, ofrecida por la misma razón de validación y no como modelo de sedimento. Para un fondo real en aguas someras, tasa el lecho con los regímenes de Weston de propagación submarina del sonido, cuyo ángulo crítico y gradiente de pérdida por reflexión codifican exactamente lo que estas condiciones de contorno omiten, y lee la respuesta modal como una cota superior del nivel recibido. O dale a gaussian_beams el propio lecho fluido (sección 4), que tasa cada rebote en el fondo con el coeficiente de Rayleigh que estos dos extremos idealizan.

En el límite de alta frecuencia la ecuación de Helmholtz colapsa en la ecuación eikonal, y sus características son los rayos: trayectorias integradas del sistema de primer orden (Jensen, Ecs. 3.23-3.24)

con $s$ la longitud de arco y $(\xi ,\zeta )$ la lentitud del rayo. En un perfil independiente de la distancia la lentitud horizontal se conserva a lo largo de cada rayo, que es la ley de Snell en forma continua, $\cos \theta (z)/c(z)=\cos {\theta }_0/c(z_0)$: un rayo se curva hacia la velocidad del sonido más baja, se aplana según crece $c$ y gira donde $c(z_{\mathrm{t}})=c(z_{\mathrm{s}})/\cos {\theta }_0$. En un gradiente lineal $c(z)=c_0+gz$ los arcos son exactamente circulares con radio $R=c_0/(g\cos {\theta }_0)$, la forma cerrada contra la que se valida el trazado de rayos; en un perfil de aguas profundas la familia de rayos se reenfoca periódicamente en las zonas de convergencia de la figura de la sección 1. ray_trace integra todos los ángulos de lanzamiento a la vez con un esquema Runge-Kutta de cuarto orden y paso fijo, reflejando en la superficie y en el fondo, y arrastra con ellos el tiempo de propagación como un tercer estado del mismo paso ($dt/dr=1/(\xi c^2)$, con $\xi =\cos {\theta }_0/c(z_{\mathrm{s}})$ el invariante de Snell):

import numpy as np
from phonometry import underwater

# Una capa profunda isoterma: c sube 0,017 (m/s)/m con la presión, así que
# un rayo lanzado 6 grados hacia abajo desde 100 m gira de vuelta donde
# Snell da c(z_t) = c(z_s)/cos(6 grados).
z = [0.0, 1000.0]
c = [1490.0, 1507.0]                        # gradiente lineal, g = 0,017 1/s
rays = underwater.ray_trace(z, c, source_depth=100.0,
                            launch_angles_deg=[6.0], max_range=40e3,
                            n_steps=20000)
z_turn = (1491.7 / np.cos(np.radians(6.0)) - 1490.0) / 0.017
print(round(z_turn, 1))                     # 583.3  profundidad de giro analítica
print(round(float(rays.depths.max()), 1))   # 583.3  trazada

# El arco fija también el tiempo hasta ese giro: la velocidad del sonido se
# cancela y solo quedan el ángulo de lanzamiento y el gradiente,
# t = (1/g) ln[(1 + sin theta_0) / cos theta_0].
th = np.radians(6.0)
r_turn = np.sin(th) / (np.cos(th) / 1491.7 * 0.017)   # 9222,6 m
t_turn = np.log((1.0 + np.sin(th)) / np.cos(th)) / 0.017
print(round(t_turn, 3))                                        # 6.171  analítico
traced = np.interp(r_turn, rays.ranges[0], rays.travel_times[0])
print(round(float(traced), 3))                                 # 6.171  trazado

Cuanto más empinado es el lanzamiento, más hondo gira el rayo: en la misma agua el rayo de 2° gira a 154 m y el de 10° a 1454 m. El rayo de 6° lleva la validación de la propia página, los 583,3 m de profundidad de giro analítica, que el trazado reproduce hasta la décima de metro que imprime, y la superposición discontinua es el arco de círculo exacto de radio $R=c_0/(g\cos {\theta }_0)$ contra el que se comprueba el cálculo. Un resultado dibuja su propio abanico con rays.plot().

import matplotlib.pyplot as plt
import numpy as np

c_top, grad, z_source = 1490.0, 0.017, 100.0
c_source = c_top + grad * z_source
fig, (ax_c, ax_r) = plt.subplots(1, 2, figsize=(13, 6),
                                 gridspec_kw={"width_ratios": [1, 3.4]})
zz = np.linspace(0.0, 2000.0, 200)
ax_c.plot(c_top + grad * zz, zz)
ax_c.invert_yaxis()
ax_c.set(xlabel="c(z) [m/s]", ylabel="Profundidad [m]")

for angle in (2.0, 4.0, 6.0, 8.0, 10.0):
    rays = underwater.ray_trace([0.0, 2000.0], [c_top, c_top + grad * 2000.0],
                                source_depth=z_source,
                                launch_angles_deg=[angle], max_range=32e3,
                                n_steps=20000)
    ax_r.plot(rays.ranges[0] / 1000.0, rays.depths[0], label=f"{angle:.0f} grados")
    z_turn = (c_source / np.cos(np.radians(angle)) - c_top) / grad
    print(f"{angle:4.0f} grados -> analitica {z_turn:7.1f} m, "
          f"trazada {float(rays.depths.max()):7.1f} m")
ax_r.set_ylim(1760.0, -70.0)
ax_r.set(xlabel="Distancia [km]", ylabel="Profundidad [m]")
ax_r.legend()
plt.show()

Los rayos compran geometría y tiempos: caminos, profundidades de giro, distancias de las zonas de convergencia y el tiempo de propagación a lo largo de cada uno, a un coste independiente de la frecuencia. Como el tiempo viaja en las mismas cuatro etapas Runge-Kutta que la trayectoria, y toma su velocidad del sonido de la interpolación que esas etapas ya hacen, describe el camino que de verdad se devuelve y no una segunda lectura de él, y coincide con la forma cerrada de un gradiente constante (Medwin y Clay 1998, Ec. (3.3.20)) hasta unos $10^{-14}$ s. RayTraceResult sigue sin informar de ningún nivel propio: la amplitud clásica por camino pertenece a los eigenrayos de justo debajo, que son los caminos que un receptor dado de verdad recoge, y un campo que se mantenga finito donde la amplitud clásica diverge es el de la sección 4, sobre los mismos rayos y con el mismo integrador.

Pasar un par bathymetry de nodos, (ranges_m, depths_m), sustituye el fondo horizontal por un perfil de profundidad lineal a trozos — la frontera facetada de la Fig. 3.20 de Jensen — mientras el perfil de velocidad del sonido sigue siendo independiente de la distancia. El integrador localiza cada cruce contra la polilínea interpolada y refleja el rayo especularmente sobre la faceta local (Ec. 3.121), así que un rebote en una pendiente de ángulo $\beta$ gira el rayo $2\beta$: los rebotes pendiente arriba empinan el rayo, los rebotes pendiente abajo lo aplanan, y esa es toda la física en una línea de la propagación en cuña. El invariante de Snell pasa entonces a ser una constante de cada rayo solo entre rebotes de fondo, y una consecuencia se afronta en vez de taparse: un rayo empinado más allá de la vertical correría hacia atrás en distancia, cosa que un integrador cuya variable independiente es la distancia no puede llevar, así que se termina en ese rebote y sus muestras son NaN de ahí en adelante — un dibujo simplemente se acaba donde el rayo dio la vuelta. En una cuña isovelocidad todo rayo son líneas rectas y espejos, así que cada ángulo de abajo es exacto:

import numpy as np
from phonometry import underwater

# Una cuña isovelocidad: el fondo sube de 1000 m a 650 m en 4 km, una
# pendiente arriba de 5 grados. Un rebote en la superficie solo invierte el
# rayo; un rebote en el fondo lo empina el doble de la pendiente.
beta = np.radians(5.0)
wedge = underwater.ray_trace(
    [0.0, 1000.0], [1500.0, 1500.0], source_depth=400.0,
    launch_angles_deg=[20.0], max_range=4000.0, n_steps=1601,
    bathymetry=([0.0, 4000.0], [1000.0, 1000.0 - 4000.0 * np.tan(beta)]))
r, z = wedge.ranges[0], wedge.depths[0]
events = wedge.bottom_reflections[0] + wedge.surface_reflections[0]
incline = np.degrees(np.arctan(np.abs(np.diff(z)) / np.diff(r)))
for leg in range(int(events[-1]) + 1):
    seg = np.flatnonzero(events[:-1] == leg)[1:-1]
    print(leg, round(float(np.median(incline[seg])), 1))
# 0 20.0   como se lanzó
# 1 30.0   primer rebote de fondo: empinado 2 beta
# 2 30.0   la superficie solo invierte el signo
# 3 40.0   segundo rebote de fondo: 2 beta otra vez

La polilínea empieza en $r=0$, de modo que la columna de agua en la fuente queda declarada y no extrapolada, y continúa horizontal más allá de su último nodo, la misma prolongación con la que vive el perfil de velocidad del sonido; un vértice refleja especularmente sobre una de sus dos facetas, y un accidente batimétrico más estrecho que un paso de distancia puede esconderse entre dos muestras de la búsqueda de cruces, así que n_steps tiene que resolver la batimetría además de los rayos. eigenrays rechaza un trazado en pendiente, y dice por qué: cada rebote en la pendiente rota el invariante de Snell, así que los toques de fondo de una llegada ya no comparten un único ángulo rasante y la convención de amplitud de un coeficiente por camino de abajo deja de sostenerse.

El abanico dibuja todos los caminos que el perfil admite; a un receptor solo llegan unos pocos. Esos son los eigenrayos, «los rayos que pasan por ese punto» (Jensen §3.3.5.2), y su lista — un retardo, un ángulo de salida y otro de llegada, los contadores de toques en los contornos y una amplitud compleja por camino — es la estructura de llegadas que la ecuación del sonar consume como multitrayecto, el esqueleto de una respuesta al impulso del canal, y aquello contra lo que ecualiza el trabajo de comunicaciones. eigenrays toma el abanico trazado y un receptor y los encuentra como raíces: hay un eigenrayo allí donde la profundidad de un rayo a la distancia del receptor cruza la profundidad del receptor, así que un par de rayos adyacentes del abanico que lo encajonan acota uno, y cada acotación se cierra por bisección sobre trazados nuevos a través del mismo perfil, nunca interpolando entre los rayos trazados — los dos rayos de una acotación pueden llevar historias de rebotes distintas, y una mezcla de ambos es un camino que el agua no contiene (el riesgo de interpolación de Jensen §3.7.5.1). Cada llegada es un rayo real, con su tiempo de propagación, sus ángulos y sus contadores de rebote leídos de su propia trayectoria, y con la divergencia $q$ de la sección siguiente montada en la misma marcha en su forma real de fuente puntual (Ec. 3.63), de modo que la amplitud clásica de la Ec. (3.65) pertenece exactamente al rayo que acertó:

# El canal isóvelo de 100 m se despliega en imágenes especulares, así que
# cada llegada se comprueba a mano: la n-ésima vuela en línea recta hasta una
# imagen del receptor, a sqrt(r^2 + z_n^2) metros.
fan = underwater.ray_trace([0.0, 100.0], [1500.0, 1500.0], source_depth=36.0,
                           launch_angles_deg=np.arange(-48.0, 48.5, 0.5),
                           max_range=600.0, n_steps=201)
arr = underwater.eigenrays(fan, receiver_range=500.0, receiver_depth=46.0)
print(arr.travel_times.size)          # 11 llegadas dentro del abanico de 48 grados

# El camino directo: 500 m en distancia, 10 m hacia abajo, sin toques, 1/R.
print(round(float(arr.travel_times[0]) * 1e3, 2))       # 333.4 ms
print(round(float(abs(arr.amplitudes[0])) * 500.1, 6))  # 1.0: eso es 1/R

# El primer eco salió de la fuente hacia arriba y llega desde arriba, con el
# signo invertido una vez por la superficie de presión nula.
print(int(arr.surface_reflections[1]), int(arr.bottom_reflections[1]))  # 1 0
print(round(float(arr.launch_angles[1]), 2))            # -9.31 grados: arriba
arr.plot(language="es")   # el esqueleto de la respuesta al impulso

El canal del propio fragmento, dibujado: una barra por eigenrayo, con la cabeza en la pérdida de ese único camino. El camino directo llega el primero y el más fuerte; el primer eco le sigue 4,4 ms por detrás, con el signo invertido por la superficie; y cada toque adicional despliega una imagen más del receptor, más tardía, más empinada y algo más débil escalera abajo. Sobre reflectores perfectos nada salvo $1/R$ atenúa un camino, así que 147 ms de dispersión de retardos cuestan apenas 3 dB: esta lista es el multitrayecto que consume la ecuación del sonar y el esqueleto sobre el que se construye una respuesta al impulso del canal.

Una sola lista sirve para todas las frecuencias. Las amplitudes están normalizadas a presión unidad a 1 m, en el mismo convenio $e^{-i\omega t}$ que el resto del módulo, y la búsqueda no pregunta nunca por una frecuencia: el tono que oye un receptor es $p(\omega )={\sum }_j{a}_j{e}^{i\omega {\tau }_j}$, su pérdida $-20\lg |p(\omega )|$ en la misma escala que los otros tres métodos, y la respuesta al impulso de una banda es esa suma transformada. Cada amplitud lleva un $(-1)$ por toque en la superficie y el coeficiente del fondo por toque en el fondo — el $\mp 1$ de los reflectores perfectos o, con el mismo FluidSeabed que toman los haces, la $\mathcal{R}$ de Rayleigh del ángulo rasante propio de cada rayo, módulo y fase — y el $(-i{)}^m$ de la Ec. (3.79) por cada cáustica cruzada, la forma discreta del $-\pi /2$ que la rama de la raíz cuadrada del método de haces gasta de forma continua. Validado contra el único entorno cuya estructura de llegadas es forma cerrada entera, la red de imágenes de la guía de ondas ideal: el censo coincide imagen por imagen, y tiempos, ángulos y amplitudes concuerdan hasta $10^{-13}$ s, $10^{-11}$ grados y $10^{-13}$ relativo; y contra el rayo parabólico exacto del perfil $n^2$-lineal (Ec. 3.195), cuyo tiempo de propagación en forma cerrada, integrado por cuadratura hasta la precisión de máquina, el eigenrayo refractado hacia arriba reproduce hasta $2\times 10^{-7}$ s, siendo el residuo el muestreo del perfil y no la búsqueda.

Tres límites, todos deliberados. El abanico es el espacio de búsqueda: un eigenrayo más empinado que la apertura trazada no existe para él, y un par de llegadas que queden entre los mismos dos rayos del abanico se funde en ninguna acotación, así que la densidad y el semiángulo del abanico son las palancas de completitud, y quedan donde se ven, en la llamada a ray_trace. La escalera de multitrayecto de una guía de fondo duro no tiene fin, así que max_arrivals (64 por defecto) conserva las llegadas más tempranas, las llanas y con menos rebotes que llevan la energía, y avisa cuando trunca. Y la amplitud es clásica a propósito: a un receptor plantado exactamente sobre una cáustica se le responde con el infinito que la Ec. (3.65) de verdad tiene ahí, porque una lista de caminos discretos es teoría de rayos y lo dice; el campo que se mantiene finito sobre esa cáustica es el de la sección siguiente, construido ensanchando estos mismos rayos en haces.

La sección 3 se queda en la geometría a propósito. La amplitud clásica de un rayo (Jensen Ec. 3.65) divide por la divergencia del tubo de rayos $q$, y $q$ se anula allí donde la familia de rayos se pliega sobre sí misma. Ese pliegue es una cáustica: la teoría de rayos responde infinito mientras el campo verdadero solo es intenso, y un rayo que la atraviesa se lleva un $-\pi /2$ del que depende todo el patrón de interferencia posterior (Jensen §3.4.1, Figs. 3.13-3.14). Más allá del último rayo de una familia queda una zona de sombra, donde la teoría de rayos no devuelve un número pequeño sino ningún número.

El trazado de haces gaussianos elimina ambas cosas a la vez, ensanchando cada rayo hasta convertirlo en un haz. La divergencia obedece a las ecuaciones dinámicas de rayos (Jensen Ec. 3.58),

que el integrador de rayos resuelve junto a la trayectoria. Partiendo de condiciones iniciales complejas, $p(0)=1$ y $q(0)=i\omega W_0^2/2$ (Ec. 3.91), cada rayo pasa a ser el eje de un haz de semianchura inicial $W_0$ y frente de onda plano, y el campo de ese haz es (Ec. 3.88)

con $n$ la distancia al rayo central y $\tau$ el tiempo de propagación a lo largo de él. El campo total es esa expresión sumada sobre el abanico de lanzamiento con los pesos de la Ec. (3.92).

Por qué se mantiene finito. La Ec. (3.58) es lineal con coeficientes reales, así que las partes real e imaginaria de $(q,p)$ son dos soluciones reales de ella, y su wronskiano $q_R{p}_I-q_I{p}_R$ se conserva. Los impulsos que la divergencia recibe en un quiebro del perfil y en una reflexión son cizallas de determinante unidad, de modo que tampoco pueden cambiarlo. Arranca en $-\omega W_0^2/2$ y ahí se queda, y esa única constante carga con todo el método: $q$ no puede anularse nunca, así que no queda ninguna singularidad de cáustica que parchear, ningún índice KMAH que contar ni ningún suelo de anchura mínima que imponer; $\mathrm{Im}[p/q]<0$ siempre, así que el haz decae siempre al alejarse de su eje; y la semianchura del haz de la Ec. (3.89) se reduce a $W(s)=2|q(s)|/(\omega W_0)$, una hipérbola en espacio libre con su cintura en la fuente y la distancia de Rayleigh $k{W}_0^2/2$ por escala.

import numpy as np
from phonometry import underwater

# El perfil n^2-lineal de Jensen, c(z) = c0/sqrt(1 + 2,4 z/c0) (Ec. 3.77). Con
# la fuente cerca del fondo los rayos ascendentes giran, su envolvente es una
# cáustica, y por encima el abanico se agota en una zona de sombra.
c_bottom_ref = 1550.0
z_beam = np.linspace(0.0, 1000.0, 201)
c_beam = c_bottom_ref / np.sqrt(1.0 + 2.4 * z_beam / c_bottom_ref)

beams = underwater.gaussian_beams(600.0, z_beam, c_beam, source_depth=992.5,
                                  max_range=2500.0, range_step=25.0,
                                  fan=underwater.BeamFan(max_angle_deg=45.0),
                                  n_depth_points=80)
print(beams.propagation_loss.shape)            # (80, 101)  profundidad x distancia
print(beams.launch_angles.size)                # 439 haces entre +-45 grados
print(round(float(beams.initial_beam_widths[0]), 1))  # 35,9 m: W0, Ec. (3.86)
beams.plot(language="es")   # el campo de PL, en el marco de parabolic_equation

El agua del propio fragmento, sobre una malla más fina: 600 Hz sobre el perfil $n^2$-lineal, la fuente a 7,5 m del fondo, 439 haces entre ±45°. El abanico ascendente gira dentro de la columna, y allí donde se pliega sobre sí mismo la suma de haces responde 59 dB, no infinito — ese pliegue es el arco luminoso que cruza lo alto de la cúpula, y es la cáustica. Por encima del arco no llega ningún rayo, y el campo no se corta en ese borde: sube 88 dB en los 100 m que hay por encima, que es la penumbra graduada que tiene la solución exacta y que la teoría geométrica de rayos no tiene (Jensen Figs. 3.11, 3.17). Las líneas grises finas son los propios rayos, trazados sobre el mismo perfil con ray_trace; de cada uno cuelga un haz. Dos cosas que la figura no puede enseñar: las tres primeras anchuras de haz, unos 108 m de distancia, donde la ponderación en campo lejano de la Ec. (3.92) no tiene a qué converger, y el 18 % de celdas a las que no llega ningún haz, que son exactamente infinitas y se dibujan en el extremo silencioso de la escala.

import matplotlib.pyplot as plt
import numpy as np
from phonometry import underwater

c_ref = 1550.0
z_caustic = np.linspace(0.0, 1000.0, 201)
c_caustic = c_ref / np.sqrt(1.0 + 2.4 * z_caustic / c_ref)

caustic = underwater.gaussian_beams(600.0, z_caustic, c_caustic,
                                    source_depth=992.5, max_range=2500.0,
                                    range_step=12.5,
                                    fan=underwater.BeamFan(max_angle_deg=45.0),
                                    n_depth_points=400)
caustic.plot(language="es")   # el mismo campo; la figura lo dibuja por bandas
plt.show()

El único parámetro libre. $W_0$ lo es, y el libro reconoce sin rodeos que «la elección óptima de estas condiciones iniciales es materia de investigación actual», recomendando de 10 a 50 longitudes de onda. Aquí el valor por defecto es más afilado que una regla empírica, y es una anchura por cada ángulo de lanzamiento en vez de una por corrida. En aguas abiertas el §3.5.1 hace la optimización de forma explícita: derivar la anchura en espacio libre de la Ec. (3.86) respecto del desplazamiento complejo da $W_0=\sqrt{\lambda r_{\max }/\pi }$, la anchura que mejor resuelve el campo en el extremo lejano del recorrido, que es donde peor se resuelve. Es también donde la integral en ángulo de lanzamiento que hay detrás de la Ec. (3.92) es una gaussiana de verdad y no una integral de Fresnel, y por eso no es solo una elección elegante: frente al campo libre a 100 Hz a 2, 5 y 8 km el error en $|p|$ es de $7,5\times 10^{-5}$ con esa anchura, de $2,7\times 10^{-2}$ con la quinta parte de ella y de $4,1\times 10^{-2}$ con quince veces. Un canal somero añade encima su propia exigencia: su campo atrapado es un conjunto discreto de modos separados $\lambda /(2D)$ en el seno del ángulo de lanzamiento, y un haz mezcla ángulos de lanzamiento sobre su divergencia de campo lejano $\lambda /(\pi W_0)$, así que resolver los modos vecinos hasta la mitad de su separación pide $W_0\ge 4D\cos {\theta }_0/\pi$: una anchura por cada ángulo de lanzamiento, máxima para los haces tendidos, cuyos modos se agolpan, y relajada por el coseno para los inclinados, con la huella vertical de todos los haces en el mismo $4D/\pi$. El valor por defecto toma el mayor de los dos, dentro de la banda de 10 a 50 longitudes de onda del libro, y entrega a todo el abanico el óptimo de espacio libre cuando el canal es demasiado profundo en longitudes de onda para que la anchura de guía quepa en la banda. Una versión anterior, en cambio, acotaba $W_0$ a un cuarto de la profundidad de la columna, leyendo el «no grandes en comparación con la profundidad del agua» del libro como un techo sobre la anchura y no sobre la huella; lo que ese tope decía proteger, la contabilidad que pliega un rayo reflejado de vuelta a la columna, ya lo restituye la escalera de imágenes del receptor, y en aguas someras costaba en silencio decibelios de nivel, siempre hacia lo demasiado flojo. Las medidas están dos párrafos más abajo, tope contra valor por defecto sobre los mismos oráculos exactos.

Contra qué está validado. La divergencia esférica en campo libre, hasta $10^{-3}$ dB, que es la única comparación que fija a la vez la normalización de la amplitud y el convenio de fase; el campo de dos rayos del espejo de Lloyd con una reflexión en superficie, hasta 0,01 dB; y la suma de fuentes imagen de la guía de ondas ideal de presión nula, hasta 0,0004 dB con el abanico abierto a 88 grados. Sobre un lecho disipativo la misma red de imágenes lleva $\mathcal{R}(\theta {)}^n$ al ángulo en forma cerrada de cada imagen, y los haces la siguen hasta 0,07 dB en el peor caso a ambos lados del ángulo crítico. Esa misma guía desarrollada sobre sus modos en vez de sobre sus imágenes (Ec. 5.13) es una segunda forma cerrada de un mismo campo exacto, y concuerda hasta 0,03 dB y $8\times 10^{-4}$ rad: es la comparación que alcanza la fase absoluta, y pone a los haces en el mismo plano en el que ya se sostienen normal_modes y parabolic_equation. La última dobla los rayos, cosa que ninguna de las otras hace, porque en un canal isovelocidad $c^{\prime \prime }$ se anula y se lleva consigo el coeficiente de acoplamiento de las ecuaciones dinámicas de rayos: el perfil $n^2$-lineal de la Ec. (3.77) tiene $k^2(z)$ exactamente lineal, así que sus modos son funciones de Airy con autovalores en forma cerrada, y los haces, con su anchura por ángulo por defecto, los siguen hasta +0,19 dB en media en 200 m de agua a 200 Hz y hasta +0,39 dB en una segunda guía independiente de 100 m a 250 Hz. Esos dos cortes son además la medida que retiró el viejo tope de un cuarto de la profundidad, que sobre esos mismos oráculos salía +3,08 y +4,12 dB demasiado flojo. Las dos condiciones de contorno salen de la propia suma de haces en vez de imponerse: el campo en una superficie o un fondo de presión nula es 3 partes en $10^5$ de su valor a media columna, y un fondo rígido lo duplica.

La absorción de volumen del agua está desactivada por defecto, así que todos los números de arriba se midieron sin ella y siguen siendo reproducibles tal cual están impresos; el precio es que el nivel más allá de unos kilómetros a frecuencias de sonar resulta optimista. Pasar absorption ("francois-garrison", "ainslie-mccolm" o "thorp", o un VolumeAbsorption que nombra su propia agua; los mismos nombres, argumentos y valores por defecto que seawater_absorption) multiplica cada haz por $e^{-\alpha s}$ con $s$ la longitud de arco a lo largo de su rayo central, que es la §3.6.2 de Jensen tal como está impresa: perturbar la eikonal con la velocidad del sonido compleja que implica una pérdida de volumen deja en pie los rayos reales y cuelga de cada uno el factor $e^{-{\int }_0^s\alpha d{s}^{\prime }}$ (Ec. 3.116). No es $\alpha r$ sobre la distancia horizontal, el atajo que la misma sección señala que «se usa en muchos modelos de rayos»: un trayecto a 60° es el doble de largo que la distancia que cubre, y los rebotes múltiples más inclinados de una guía de ondas son exactamente las llegadas que la absorción tiene que drenar. El integrador de rayos acumula $s$ con las mismas etapas Runge-Kutta que colocan el rayo ($ds/dr=1/(\xi c)$), ray_trace la expone por rayo como arc_lengths, y el coeficiente, un solo $\alpha$ por ejecución evaluado a la frecuencia y la profundidad de la fuente, queda registrado en el resultado como absorption_coefficient.

Un lecho marino disipativo es la otra pérdida que el valor por defecto deja fuera, y en aguas someras la dominante: un fondo perfecto devuelve intactos todos los múltiples inclinados, cuando lo que el lecho se queda es lo que más pesa en la pérdida de propagación real en aguas someras. Pasar un FluidSeabed como bottom (la misma descripción fluida que acepta seabed_reflection, con water_density para el agua de encima) sustituye el reflector perfecto por el coeficiente de Rayleigh al ángulo rasante propio de cada haz: cada toque en el fondo multiplica la amplitud compleja del haz por $\mathcal{R}$, módulo y fase a la vez (Jensen §3.6.3, Ecs. 3.125-3.126). La fase no es un refinamiento prescindible: por debajo del ángulo crítico $|\mathcal{R}|=1$ y solo la fase distingue el lecho disipativo de uno perfecto, y desplaza las franjas de interferencia de todos los trayectos que tocan el fondo. El ángulo al que se tasa cada haz es el que fija el invariante de Snell en el fondo, $\cos \theta =\xi c(D)$, el mismo en todos sus toques haga lo que haga el perfil de encima entre medias, así que su coeficiente se evalúa una sola vez, exacto, y se eleva al recuento de toques en el fondo del integrador; ray_trace expone esos recuentos por contorno junto a sus longitudes de arco, que es todo lo que una amplitud necesita de la geometría. Validado contra la red de imágenes con $\mathcal{R}(\theta {)}^n$ al ángulo en forma cerrada de cada imagen: 0,07 dB en el peor caso en distancias cuyas imágenes dominantes caen por encima y por debajo del ángulo crítico, donde despojar a la $\mathcal{R}$ del oráculo de su fase lo mueve 5,7 dB.

El fondo también puede inclinarse aquí: el mismo par bathymetry que toma ray_trace (sección 3) pone el rayo central de cada haz sobre la frontera facetada, y el par dinámico cruza cada rebote en pendiente con el impulso de reflexión de las Ecs. (3.122)-(3.123) evaluado en la faceta local. Dos de los dispositivos de la guía horizontal se generalizan con él, cada uno con su coste declarado. La escalera de imágenes del receptor pliega cada columna de receptores sobre su propia faceta — profundidad e inclinación, porque un plano de plegado equivocado en la inclinación de la faceta desplaza una imagen del primer pliegue metros frente a una longitud de onda — de modo que la pila vertical de espejos se convierte en el abanico diédrico en torno al vértice local, exacto para una sola faceta; a pendiente cero es la escalera horizontal bit a bit, cosa que un test deja clavada. Y los peldaños envueltos del abanico dejan que las colas analíticas de los haces lleven las llegadas que ningún eje integrado puede: un camino que subió la pendiente, se empinó más allá de la vertical y volvió es un peldaño envuelto, analítico en vez de integrado. El propio haz, una vez empinado más allá de la vertical, se termina en ese rebote con su peso a cero — una suma que avanza en distancia conserva solo lo que aún viaja hacia delante, exactamente igual que la ecuación parabólica de un solo sentido no conserva retrodispersión — y una cola no se acepta más allá de dos extensiones de la marcha desde su muestra, porque más allá estaría poniendo precio a geometría de cuña que la polilínea nunca describió (medido: hasta 11 dB de campo inventado sobre un fondo inclinado al que ningún haz llega siquiera, 0,001 dB con el presupuesto en su sitio). Validado contra la cuña ideal, que tiene solución exacta por imágenes — un abanico de $2\pi /\beta$ fuentes imagen con signo sobre una circunferencia en torno al vértice, construido en el test con geometría pura y probado contra las dos condiciones de contorno antes de juzgar: un solo rebote de faceta a 0,05 dB, y a 1,85 dB en el peor caso / 0,67 dB en media del campo completo a través de una cuña fina de 2,8° en la que cada celda es multitrayecto denso, la mayor parte llegadas cerca de su propio punto de giro, números que se mueven en el cuarto decimal al partir en dos el paso de distancia. El lecho fluido disipativo de arriba no puede combinarse con una pendiente — un ángulo rasante por haz es un hecho de fondo horizontal, y el método se niega a fingir lo contrario — así que una pasada en pendiente toma los reflectores perfectos de bottom.

Dónde se acaba. Cuatro límites, en el orden en que aprietan.

• No hay campo cercano, y este es el mayor error de los cuatro. La Ec. (3.92) pesa el abanico ajustándolo a una fuente puntual en campo lejano, y la Ec. (3.88) divide por una distancia cilíndrica que se anula sobre el eje del que sale cada rayo, así que de cerca la suma no tiene a qué converger. La escala a la que se recupera es $W_0$, no una distancia fija: en tres ajustes cuyo $W_0$ abarca de 150 a 437 m, el peor error frente a $20\lg R$ en un medio no acotado es de 17, 13 y 4,1 dB con un cuarto de $W_0$, de unos 0,6 dB con $W_0$, de una centésima de decibelio con $2,5W_0$ y de una milésima a partir de $3W_0$. No leas nada dentro de unas tres anchuras de haz de la fuente, y ten en cuenta que, como la anchura de espacio libre del valor por defecto crece como $\sqrt{r_{\max }}$, un recorrido más largo aleja aún más esa frontera. De cerca, usa parabolic_equation.
• El régimen propio de la teoría de rayos (Jensen §3.4.2): «la longitud de onda debería ser sustancialmente menor que cualquier escala física del problema». Es el límite que más aprieta y el que mejor esconde una respuesta de aspecto verosímil, y una versión anterior de este punto le achacó un error que en realidad era del tope de anchura: a 20 Hz en 100 m de agua, donde la profundidad son 1,3 longitudes de onda y se propagan dos modos, el haz acotado quedaba en un tercio de longitud de onda de ancho y la pérdida salía decibelios alta frente a la suma de fuentes imagen. Retirado el tope, la misma guía queda a 0,03 dB de esa suma con la anchura por defecto. El aprobado es más estrecho de lo que parece: una columna isovelocidad entre reflectores perfectos es pura geometría, que las imágenes plegadas del receptor reproducen exactamente a cualquier frecuencia, así que no dice nada de un canal por el que el campo de baja frecuencia de verdad refracta. Ahí normal_modes sigue siendo el método en el que confiar, exacto al precio de un puñado de modos.
• El abanico está truncado en el max_angle_deg del abanico, y una guía de ondas con dos contornos perfectamente reflectantes es el peor caso posible para eso, porque nada salvo $1/R$ atenúa los rebotes múltiples más inclinados. En la guía ideal de 1000 m a 300 Hz, frente a la suma de fuentes imagen a 2, 5 y 10 km: 0,27, 4,06 y 2,52 dB con los 80 grados por defecto, que bajan a 0,0002, 0,0003 y 0,0004 dB al abrir el abanico a 88 grados. Recortar el oráculo al mismo semiángulo lo mueve 0,25, 3,95 y 2,31 dB, así que esto es el abanico y no el método. Un fondo real y disipativo (el FluidSeabed de más arriba) absorbe esos rebotes y entonces el valor por defecto sobra. Abrir el abanico obliga a recortar range_step con él, porque un paso tiene que resolver $\tan {\theta }_{\max }$ unidades de profundidad de ascenso por unidad de distancia; el método avisa cuando ese emparejamiento está mal.
• La sombra lejana tiene suelo. Cada haz se suma hasta cuatro semianchuras, 140 dB por debajo de su propio eje, así que un receptor al que ningún haz del abanico se acerque tanto recibe exactamente cero y una pérdida infinita. Esa es la cuña sin iluminar fuera de la apertura trazada, no la penumbra graduada justo detrás del rayo límite, que es donde está lo interesante de una zona de sombra y que los haces sí resuelven.

El otro sabor, los haces geométricos (Jensen §3.3.5.5, Ecs. 3.72-3.76), no está implementado. Mantiene $q$ real y toma la anchura del propio tubo de rayos, $W(s)=|q(s)\delta {\theta }_0|$, así que la anchura se anula en una cáustica y hay que apuntalarla con el suelo $\pi \lambda$ de Weinberg-Keenan y el índice KMAH de la Ec. (3.79). Ese es el enfoque parcheado que aquí se evita a propósito, y de todos modos vale la pena citar el veredicto del propio libro: los haces geométricos «han demostrado en general ser más satisfactorios» a baja frecuencia, donde la física hace que el haz sea grande comparado con el canal.

La ecuación parabólica cambia el problema de contorno de Helmholtz por un problema de valores iniciales en distancia. Factoriza la oscilación saliente rápida, $p(r,z)\simeq \psi (r,z)e^{i{k}_{0}r}/\sqrt{k_{0}r}$ con un número de onda de referencia $k_0=\omega /c_0$, y para la energía que viaja dentro de un ángulo pequeño respecto a la horizontal la envolvente obedece la PE estándar (Tappert) (Jensen §6.2):

El algoritmo split-step de Fourier la hace avanzar por descomposición de operadores, alternando dos semipasos físicos, exacto cada uno por separado: la difracción es una multiplicación por $e^{-i{k}_z^2\Delta r/2k_0}$ en el dominio del número de onda vertical, y la refracción una pantalla de fase $e^{i{k}_0(n^2-1)\Delta r/2}$ de vuelta en profundidad, con un par de transformadas por paso de distancia. parabolic_equation arranca de un campo gaussiano adaptado a la fuente puntual y usa una transformada discreta de senos en profundidad, que impone por construcción la superficie y el fondo de presión nula. El precio es la aproximación paraxial: la PE estándar es precisa dentro de ±15-20° de la horizontal, aproximadamente, y la energía más empinada arrastra un error de fase que se nota a corta distancia en guías de ondas someras (Jensen §6.2). La calibración en campo libre es el oráculo: sin gradiente alguno, el campo así avanzado debe reproducir la divergencia esférica,

import numpy as np
from phonometry import underwater

# Campo libre: en una columna isovelocidad de 5000 m, antes de sentir
# ningún contorno, la PE debe reproducir la divergencia esférica,
# PL = 20 lg R.
field = underwater.parabolic_equation(50.0, [0.0, 5000.0], [1500.0, 1500.0],
                                      source_depth=2500.0, max_range=2000.0,
                                      range_step=10.0)
iz = np.argmin(np.abs(field.depths - 2500.0))
ir = np.argmin(np.abs(field.ranges - 1000.0))
print(round(float(field.propagation_loss[iz, ir]), 2))   # 60.0 = 20 lg 1000

y lo hace hasta unos $10^{-4}$ dB con el paso de distancia por defecto. Esa es la comprobación que la PE no puede fallar: campo libre, sin contornos, sin energía empinada. El límite paraxial solo se ve donde faltan esas tres cosas.

En una guía de ondas de 100 m a 50 Hz, cuatro de los seis modos propagantes salen de la fuente a 27°, 37°, 49° y 64°: fuera de la banda para la que está escrita la PE estándar. La consecuencia está a la izquierda: la pérdida de la PE promediada en distancia queda 1,6 dB por encima de la referencia modal, a 300 m y todavía a 5,7 km. El sesgo no se diluye con la distancia, porque una guía de ondas ideal no tiene pérdidas y no descama ningún modo; en un fondo marino real los modos empinados se atenúan antes y la diferencia sí se cierra. Este es el modo de fallo que hay que buscar en la propia salida: un nivel uniformemente un poco bajo, no una forma equivocada.

import matplotlib.pyplot as plt
import numpy as np

r = np.linspace(50.0, 6000.0, 700)
modes = underwater.normal_modes(50.0, [0.0, 100.0], [1500.0, 1500.0],
                                source_depth=30.0, receiver_depth=60.0,
                                ranges_m=r, n_depth_points=800)
pe_field = underwater.parabolic_equation(50.0, [0.0, 100.0], [1500.0, 1500.0],
                                         source_depth=30.0, max_range=6000.0,
                                         range_step=5.0, n_depth_points=1024)
iz2 = int(np.argmin(np.abs(pe_field.depths - 60.0)))
pe_pl = np.interp(r, pe_field.ranges, pe_field.propagation_loss[iz2])

def range_average(pl, width=61):        # incoherente, en el dominio de la energía
    return -10.0 * np.log10(np.convolve(10.0 ** (-pl / 10.0),
                                        np.ones(width) / width, mode="same"))

inner = (r > 300.0) & (r < 5700.0)
sm_nm = range_average(np.asarray(modes.propagation_loss))
sm_pe = range_average(pe_pl)
print(round(float(np.mean(sm_pe[inner] - sm_nm[inner])), 2))   # 1.64 dB

fig, (ax_pl, ax_ang) = plt.subplots(1, 2, figsize=(13.5, 5.6),
                                    gridspec_kw={"width_ratios": [2, 1]})
ax_pl.plot(r / 1000.0, modes.propagation_loss, alpha=0.35)
ax_pl.plot(r / 1000.0, pe_pl, alpha=0.35)
ax_pl.plot(r[inner] / 1000.0, sm_nm[inner], label="Modos normales, promediado")
ax_pl.plot(r[inner] / 1000.0, sm_pe[inner], label="PE, promediada")
ax_pl.set_ylim(30.0, 92.0)
ax_pl.invert_yaxis()
ax_pl.set(xlabel="Distancia [km]", ylabel="Pérdida de propagación [dB]")
ax_pl.legend()

k0 = 2.0 * np.pi * 50.0 / 1500.0
grazing = np.degrees(np.arccos(np.clip(modes.wavenumbers / k0, -1.0, 1.0)))
ax_ang.bar(np.arange(1, grazing.size + 1), grazing)
ax_ang.axhline(20.0, linestyle="--")
ax_ang.set(xlabel="Índice de modo m",
           ylabel="Ángulo rasante arccos(k_rm/k) [grados]")
plt.show()

Los cuatro métodos los configura quien llama, y todos los fragmentos anteriores fijan una malla sin decir por qué. Las reglas son cortas.

La malla en profundidad. La transformada de senos de la PE, y el autoproblema en diferencias finitas de los modos, producen ambos aliasing de un espectro que es empinado en profundidad, así que la malla se dimensiona en longitudes de onda. A 50 Hz, $\lambda =30$ m, y una malla cercana a $\lambda /10$ sobre 5000 m de agua son unos cuantos miles de puntos interiores, que es por lo que el n_depth_points=512 del fragmento de la sección 1 está etiquetado como el ajuste grueso y rápido y no como el correcto. Refina y vuelve a calcular antes de fiarte de un número.

El paso en distancia del split-step. El semipaso de refracción es una pantalla de fase $e^{i{k}_0(n^2-1)\Delta r/2}$, así que debe girar la fase bastante menos de un radián por paso, y hay que apretarlo allí donde el perfil tenga gradientes fuertes. La prueba honesta es la misma: divide $\Delta r$ entre dos y mira si el campo se mueve.

Cuando normal_modes avisa. El aviso dice que un modo retenido queda dentro de la banda de error de las diferencias finitas. El remedio es subir n_depth_points: el error de autovalor es $O(\Delta z^2)$, y los modos cercanos al corte son sus primeras víctimas, que son justo los modos de ángulo pequeño que llevan el campo a larga distancia, así que ignorar el aviso cuesta en silencio la respuesta a las distancias que interesan.

El abanico de rayos y el paso de rayo. Los ángulos de lanzamiento tampoco son una elección libre: se toman de la condición de giro de Snell para que el abanico abarque los rayos que el canal atrapa, en vez de los ±12° redondos del fragmento de arriba. El número de pasos tiene que mantener la longitud de arco por paso pequeña frente a la escala vertical del perfil, o un punto de giro se salta de una zancada.

El abanico de haces y la anchura del haz. gaussian_beams toma ese mismo abanico y le añade dos mandos. n_beams decide si los haces vecinos siguen solapándose allí donde más se han ensanchado: dos haces contiguos distan $s\delta {\theta }_0$ a la longitud de arco $s$ mientras cada uno se ha ensanchado hasta $s\lambda /(\pi W_0)$, así que la condición de solape $\delta {\theta }_0\lesssim \lambda /(\pi W_0)$ es independiente de la distancia y el valor por defecto toma cuatro veces ese margen. Un abanico demasiado basto no se lee como un nivel equivocado sino como un rizado periódico en distancia al paso de los haces, fácil de confundir con interferencia física. beam_width es el $W_0$ de la sección 4, cuyo valor por defecto es el óptimo solo donde el canal no lo recorta, así que en aguas someras pásalo explícitamente. Los tres viajan juntos como el BeamFan que toma el método. Y abrir max_angle_deg obliga a recortar range_step con él, porque un paso tiene que resolver $\tan {\theta }_{\max }$ unidades de profundidad de ascenso por unidad de distancia.

La disciplina de convergencia. Refinar, recalcular, comparar; y solo entonces comparar métodos. Que los modos y la PE coincidan, que es lo que la sección 6 llama la prueba práctica de convergencia, no vale nada si ninguno de los dos ha convergido en su propia malla.

Todas las funciones de propagación del módulo underwater responden a la misma pregunta, «cuánto nivel sobrevive al camino», a distinto precio en física. La terminología sigue en toda la página la ISO 18405:2017 (pérdida por propagación, nivel de fuente, niveles re 1 µPa). Los modelos de velocidad del sonido y de absorción nombrados abajo están implementados y referenciados en la propagación submarina del sonido.

Velocidad del sonido. Las cuatro ecuaciones coinciden a menos de 1 m/s dentro del dominio que comparten: 0,98 m/s a 25 °C, 35 ppt y 1000 m, y 0,52 m/s en el kilómetro superior de un perfil completo, como muestra la comparación de modelos de la página de propagación, así que la elección va de rango de validez, no de exactitud. La forma UNESCO / Chen-Millero por defecto (en la reformulación de Wong y Zhu de 1995) cubre de 0 a 40 °C, de 0 a 40 ppt y de 0 a 1000 bar, la envolvente más amplia, y es el estándar internacional. Del Grosso (1974) queda restringida a 0-30 °C y 30-40 ppt, pero algunos autores la prefieren para trabajo de océano profundo dentro de ese dominio (buena parte de la literatura del canal SOFAR la usa). Mackenzie (1981) cambia presión por profundidad directamente (2-30 °C, 25-40 ppt, 0-8000 m), lo que la convierte en la opción cómoda cuando se dispone de una profundidad de ecosonda y no de una presión de CTD; las otras dos convierten profundidad a presión internamente mediante Leroy y Parthiot (1998). Medwin (1975) es una simplificación de seis términos válida hasta unos 1000 m en aguas someras y cálidas; es la más basta de las cuatro y la que hay detrás de las reglas prácticas $\partial c/\partial T\approx 4,6-0,110T$ m/s por °C y $\partial c/\partial z\approx 0,016$ m/s por metro, así que resérvala para el cálculo mental y las comprobaciones a mano, no para un presupuesto.

Absorción. Francois-Garrison (1982) es la referencia y el valor por defecto: incorpora las relajaciones del ácido bórico, el sulfato de magnesio y el agua pura con sus dependencias completas de temperatura, salinidad, profundidad y pH implícito, y es fiable desde unos 100 Hz hasta 1 MHz. Ainslie-McColm (1998) es una simplificación deliberada de la misma física que se mantiene a menos de un 10 % de ella en ese rango; úsese cuando una fórmula legible importe más que el último punto porcentual. Thorp (1967) depende solo de la frecuencia (lleva incorporada agua a 4 °C cerca de 1000 m) y es anterior a ambas; resérvese para estimaciones rápidas de baja frecuencia por debajo de unas decenas de kHz y para comparar con literatura antigua que la usó.

Ley de divergencia. La divergencia esférica ($20{\log }_{10}R$) describe un frente de onda que se expande libremente en tres dimensiones, antes de que ningún contorno lo confine; el cilíndrico ($10{\log }_{10}R$) describe energía atrapada entre la superficie y el fondo (o en el canal SOFAR) que solo puede expandirse en distancia. La ley "practical" empalma ambos en una distancia de transición $R_0$, físicamente del orden de la profundidad del agua (o del canal): esférico mientras el frente no ha llenado el conducto, cilíndrico cuando ya lo ha hecho. En el ejemplo de 10 kHz de la sección de pérdida de propagación la elección no es cosmética: frente a la misma figura de mérito de 87 dB, la divergencia solo esférica predice detección hasta unos 8,7 km, mientras que la ley práctica con $R_0=1000$ m la estira hasta unos 15,8 km. Cuando la ley de divergencia es la mayor incertidumbre del presupuesto, esa es la señal de dejar la forma cerrada y calcular el campo.

¿Forma cerrada o cálculo numérico? La pérdida de propagación en forma cerrada no sabe nada del perfil de velocidad del sonido, del fondo ni de la superficie; es honesta para caminos cortos, directos y sin contornos, y para presupuestos de sonar de primera pasada. Cuando la refracción y los contornos deciden la respuesta, elíjase el método por frecuencia y geometría (Jensen et al. 2011, cap. 1):

En la práctica las fronteras se difuminan: los rayos siguen siendo útiles a frecuencias sorprendentemente bajas para trabajo de tiempos de propagación, y la PE sigue siendo el caballo de batalla bastante por encima de su régimen formal de ángulos pequeños. Cuando dos de los cuatro coinciden en un caso, como hacen los modos y la PE en la figura de la sección 1, esa coincidencia es la prueba práctica de convergencia. Los dos extremos se reparten el trabajo con limpieza: normal_modes y parabolic_equation se encarecen al subir la frecuencia, porque ambos tienen que resolver la longitud de onda sobre una malla, mientras que el núcleo de rayos no, así que gaussian_beams es el que sigue siendo asequible justo donde los otros dejan de serlo, y es también aquel cuya aproximación está mejor justificada allí.

Encadenando las piezas de punta a punta: una fuente de 140 dB re 1 µPa²/Hz a 10 kHz, un nivel espectral de ruido ambiental de 60 dB, una ganancia de array de 15 dB y un umbral de detección de 8 dB dan la figura de mérito $FOM=140-(60-15)-8=87$ dB, que calcula passive_sonar_equation en la sección de la ecuación del sonar de la guía de propagación. Cada término mueve después el cruce de la curva de pérdida de propagación de la sección de pérdida de propagación de esa guía (10 °C, 35 ppt, 100 m), y ese cruce es la distancia de detección predicha.

import numpy as np

env = {"temperature": 10.0, "salinity": 35.0, "depth": 100.0}
se = underwater.passive_sonar_equation(
    source_level=140.0, propagation_loss=np.linspace(40.0, 120.0, 400),
    noise_level=60.0, directivity_index=15.0, detection_threshold=8.0)
print(se.figure_of_merit)                                    # 87.0 dB

for law, r0 in (("spherical", None), ("practical", 1000.0)):
    res = underwater.detection_range(se.figure_of_merit, 10e3, law=law,
                                     transition_range=r0, **env)
    print(law, round(float(res.detection_range) / 1000.0, 1))  # 8.7 / 15.8 km

# Dos maneras de perderla: la mitad de ganancia de array, o el doble de frecuencia.
trim = underwater.passive_sonar_equation(
    source_level=140.0, propagation_loss=np.linspace(40.0, 120.0, 400),
    noise_level=60.0, directivity_index=7.5, detection_threshold=8.0)
print(trim.figure_of_merit)                                  # 79.5 dB
print(round(float(underwater.detection_range(
    trim.figure_of_merit, 10e3, law="practical",
    transition_range=1000.0, **env).detection_range) / 1000.0, 1))       # 10.0 km
print(round(float(underwater.detection_range(
    se.figure_of_merit, 20e3, law="practical",
    transition_range=1000.0, **env).detection_range) / 1000.0, 1))       # 5.9 km

La misma fuente y el mismo mar, seis respuestas entre 4,1 y 15,8 km. La ley de divergencia es la palanca mayor: elegir la esférica sola en lugar de la práctica con $R_0=1$ km cuesta 7,1 km a 10 kHz. Reducir a la mitad la ganancia de array cuesta 5,8 km; doblar la frecuencia cuesta 9,9 km, porque $\alpha$ pasa de 0,95 a 3,30 dB/km. Cuando la ley de divergencia es la mayor incertidumbre de un presupuesto, esa es la señal de dejar la forma cerrada y calcular el campo.

import matplotlib.pyplot as plt

ranges = np.linspace(50.0, 30_000.0, 800)
fig, ax = plt.subplots(figsize=(11, 6.4))
cases = (("10 kHz, solo esférica", 10e3, "spherical", None),
         ("10 kHz, práctica R0 = 1 km", 10e3, "practical", 1000.0),
         ("20 kHz, práctica R0 = 1 km", 20e3, "practical", 1000.0))
for label, freq, law, r0 in cases:
    ax.plot(ranges / 1000.0,
            underwater.propagation_loss(ranges, freq, law=law,
                                         transition_range=r0, **env).pl,
            label=label)
for fom in (se.figure_of_merit, trim.figure_of_merit):
    ax.axhline(fom, linewidth=1.4, label=f"FOM = {fom:.1f} dB")
    for _label, freq, law, r0 in cases:
        r50 = float(underwater.detection_range(fom, freq, law=law,
                                               transition_range=r0,
                                               **env).detection_range)
        ax.plot([r50 / 1000.0], [fom], "o")
ax.set_ylim(40.0, 126.0)
ax.invert_yaxis()
ax.set(xlabel="Distancia [km]", ylabel="Pérdida de propagación [dB]")
ax.legend()
plt.show()

Un número de ese presupuesto merece una segunda mirada. Un nivel espectral ambiental de 60 dB a 10 kHz está bastante por encima de lo que produce el mar por sí solo: ocean_ambient_noise da ahí 42 dB re 1 µPa²/Hz con un viento de 15 nudos. Sesenta decibelios describen por tanto un emplazamiento ruidoso o un receptor limitado por su propio ruido, no un océano tranquilo, y darse cuenta de eso es la diferencia entre un presupuesto y un ejercicio de aritmética. Este acoplamiento entre el modelo de absorción, la ley de divergencia y la ecuación del sonar es la razón de que todos vivan en un mismo módulo.

• Cubierto

  normal_modes, ray_trace, gaussian_beams y parabolic_equation implementan los métodos numéricos de Jensen et al. (2011) para un océano estratificado independiente de la distancia: la suma modal de Sturm-Liouville por diferencias finitas (cap. 5, Ecs. 5.3-5.17), la integración Runge-Kutta de las trayectorias de rayos (cap. 3, Ecs. 3.23-3.24), la suma de haces gaussianos construida sobre ella (cap. 3, §3.5, Ecs. 3.88-3.92) y la PE estándar split-step de Fourier (cap. 6), con la terminología de la ISO 18405:2017 en toda la página. ray_trace devuelve además el tiempo de propagación de cada rayo, integrado como un tercer estado del mismo paso Runge-Kutta, y eigenrays busca en un abanico trazado las llegadas que unen una fuente con un receptor, cada una refinada por bisección sobre trazados nuevos y cargada con la amplitud clásica de rayo de la Ec. (3.65) con las fases de contorno y de cáustica en el registro (Ecs. 3.79, 3.125-3.126). ray_trace y gaussian_beams aceptan además un perfil de fondo lineal a trozos, la frontera facetada del cap. 3, que refleja cada rayo sobre la faceta local (Ec. 3.121) con el par dinámico cruzando el rebote según las Ecs. (3.122)-(3.123). Cada uno está validado contra una forma cerrada exacta: los modos de la guía de ondas ideal de presión nula y su suma de fuentes imagen (para los eigenrayos, la lista completa de llegadas de esa red: tiempo, ángulo, contadores y amplitud por imagen), los rayos en arco de círculo de un gradiente lineal junto con el tiempo de propagación publicado a lo largo de ellos (Medwin y Clay 1998, Ec. (3.3.20)), la divergencia esférica en campo libre y, para el fondo en pendiente, el abanico cerrado exacto de imágenes de la cuña ideal, más la coincidencia mutua de la pérdida de propagación modal y de la PE.

• No cubierto

  Los cuatro métodos suponen una columna de agua independiente de la distancia con superficie de presión nula, y el fondo es un reflector perfecto por defecto, de presión nula o (para los modos y los haces) rígido. Los haces pueden cambiarlo por el lecho fluido disipativo del modelo de Rayleigh (solo con fondo horizontal), y los dos métodos basados en rayos pueden inclinar el fondo con un perfil de profundidad lineal a trozos, que es la única dependencia de la distancia que hay: una columna de agua dependiente de la distancia $c(r,z)$ queda excluida deliberadamente, porque la frontera en pendiente tiene un oráculo exacto (el abanico cerrado de imágenes de la cuña ideal) y $c(r,z)$ no tiene ninguno. No hay fondo elástico ni atenuación del sedimento. El trazado de rayos devuelve caminos y tiempos de propagación, no amplitudes; el nivel por camino pertenece a eigenrays, cuya lista de llegadas es teoría de rayos clásica y diverge honestamente para un receptor plantado sobre una cáustica, mientras que gaussian_beams es lo que da un campo finito, sumando haces en vez de buscando rayos. Su alternativa de haces geométricos (cap. 3, §3.3.5.5) con el índice KMAH y el suelo de anchura de Weinberg-Keenan no está implementada. La PE es la forma estándar de ángulos pequeños (Tappert), no una variante de gran angular de Padé. Para la física del lecho elástico que estos métodos fluidos dejan fuera, véase Ondas elásticas y acoplamiento fluido-sólido.

Elige por frecuencia y geometría (Jensen et al. 2011, cap. 1): rayos para alta frecuencia y océano profundo (geometría de los caminos de rayos, tiempos de propagación y zonas de convergencia a un coste independiente de la frecuencia, sobre fondo horizontal o en pendiente, más, a través de eigenrays, la lista de llegadas por receptor), haces gaussianos cuando de esa misma geometría de alta frecuencia haga falta un nivel (el campo PL($z$,$r$), finito en las cáusticas y graduado en las zonas de sombra), modos normales para baja frecuencia en aguas someras (pocos modos propagantes, $M\approx kD/\pi$, la solución de referencia en su régimen) y la ecuación parabólica para caminos largos de un solo sentido a baja frecuencia (el campo completo PL($z$,$r$) con refracción). Cuando dos métodos coinciden en un caso, esa coincidencia es la prueba práctica de convergencia.

Cuando la refracción o los contornos deciden la respuesta: un mínimo de velocidad del sonido que atrapa energía (el canal SOFAR), reflexiones en la superficie y el fondo en aguas someras, o una distancia de detección que oscila con la elección de la ley de divergencia. La forma cerrada $PL=\text{divergencia}+\alpha R$ es honesta para caminos cortos, directos y sin contornos, y para presupuestos de sonar de primera pasada; más allá, calcula el campo.