Escuchar terremotos: trazar rayos por un planeta entero

Una onda S que alcanza el límite núcleo–manto está acabada. El núcleo externo es líquido, un líquido no tiene resistencia a la cizalla y una onda de cizalla necesita algo que cizallar: así que a 2.891,5 km de profundidad la fase termina. Lo que eso deja al otro lado de un gran terremoto es una zona de sombra sísmica del tamaño de un continente, y es una de las observaciones que nos dijeron que el centro del planeta está fundido, décadas antes de que nadie pudiera modelarlo.

Esa sombra aparece en pantalla en Tremor cada vez que seleccionas un evento grande, y cómo llega hasta ahí es todo el diseño en miniatura. Ningún código dibuja una sombra. No existe en todo el renderizador ningún umbral de distancia angular. Trazamos rayos a través de un modelo de velocidades, le preguntamos a cada uno dónde está en el instante actual de reproducción y escribimos las respuestas en un búfer reservado de antemano.

Los rayos que murieron en el núcleo escriben NaN. El renderizador de líneas corta su serie en un NaN y empieza otra. Los huecos son las zonas de sombra.

El modelo de tiempos de viaje ak135, en 22 nodos

Todo lo que viene después se apoya en una tabla. ak135 es un modelo de referencia unidimensional estándar de la Tierra: para cada profundidad, una velocidad de compresión y otra de cizalla promediadas sobre el planeta entero — lo bastante aproximado como para describirlo en un párrafo, lo bastante bueno como para que los sismólogos sigan localizando terremotos reales con él. Nuestra versión son 22 filas de [depthKm, vp, vs], interpoladas linealmente a trozos en el radio en lugar de en la profundidad, a unos 0,1 km/s del modelo publicado. El propio código lo llama grado visual y lo dice en serio: es una herramienta de dibujo, no una tabla de tiempos de viaje ak135 contra la que cronometrarías la lectura de una fase.

El truco de la tabla está en cómo se codifican las discontinuidades: una profundidad que aparece dos veces es una discontinuidad. El constructor recorre la lista de nodos, la parte en segmentos continuos allí donde una profundidad se repite y anota el salto. Nada más en el motor conoce los límites por su nombre: el Moho, la 410, la 660, el límite núcleo–manto y el límite del núcleo interno son una sola estructura de datos, y añadir otro sería añadir dos filas.

El núcleo externo líquido no es un caso especial dentro del motor. Es un cero en una tabla.

Ese cero hace un trabajo enorme. Cuando un rayo S alcanza el límite núcleo–manto, el paso de refracción consulta la velocidad de cizalla al otro lado, no encuentra nada hacia lo que refractarse y termina el rayo: el trazador nunca llega a enterarse de que ha chocado con un núcleo. Por debajo del límite del núcleo interno la velocidad de cizalla vuelve, porque el núcleo interno es sólido; el trazador propagaría allí tan contento si alguna onda de cizalla pudiera llegar. Ninguna llega, y una aserción comprueba que ninguna lo ha hecho jamás.

Curvar 360 rayos a través de un planeta

Una Tierra con simetría esférica te regala una dimensión: un rayo que sale de una fuente nunca abandona el plano que contiene la fuente y el centro del planeta. Así que el trazador es bidimensional: avanza cada rayo por ese plano de círculo máximo y lo guarda como un radio, una distancia angular y un tiempo acumulado — tres floats por muestra almacenada, en un Float32Array plano, con una muestra conservada cada seis pasos.

La integración es teoría de rayos clásica: con índice de refracción n = 1/v, la ecuación del rayo d/ds(n·d̂) = ∇n dice que un rayo se curva hacia el índice más alto, es decir, hacia el material más lento. En un medio de simetría radial el gradiente apunta a lo largo del radio, así que toda la curvatura se reduce a una docena de líneas de aritmética escalar, ejecutadas una vez por cada paso de ocho kilómetros.

if (r > 0.5) {
  const dndr = -velocityGradientAt(r, phase) / (v * v);
  const k = (dndr * DS_KM) / r;
  const n = 1 / v;
  const wx = n * dx + k * x;
  const wz = n * dz + k * z;
  const wl = Math.hypot(wx, wz);
  if (wl > 0) {
    dx = wx / wl;
    dz = wz / wl;
  }
}

Avanzar en línea recta funciona hasta que un rayo se encuentra con una discontinuidad, y ahí la aproximación se nota al instante: cruza el límite núcleo–manto unos kilómetros tarde y el ángulo de refracción sale mal y el rayo acaba donde no debería. Así que los cruces se resuelven de forma exacta. stepToRadius() halla la raíz de la cuadrática que sitúa el rayo justo sobre el límite, el paso se trunca ahí y se aplica la ley de Snell en la base tangencial/normal local: sin θ₂ = sin θ₁ · v₂/v₁. Si eso pasa de uno, no hay rayo transmitido y la dirección se refleja — reflexión total interna, que es como aparecen las fases reflejadas en el núcleo sin haberlas pedido. Si no, se reconstruye la dirección refractada y se aparta 0,02 km, para que el paso siguiente no pueda volver a detectar el mismo cruce.

Trescientos sesenta rayos se abren en abanico desde el hipocentro entre 0,25° y 179,75° de ángulo de salida, con la profundidad de la fuente acotada a 1–700 km. Un abanico completo cuesta 10–20 ms por fase — en el hilo principal, sin Web Worker, algo defendible únicamente porque ocurre una vez por selección y no una vez por fotograma.

Las fases que nunca escribimos

La sismología tiene un alfabeto de fases rico — P, S, PKP, PKIKP, PcP, ScS — y la forma tentadora de dibujarlo es tratar cada una como un caso especial: dibuja este arco para una fase de núcleo, oculta la P directa pasado aquel ángulo. Tremor no tiene nada de eso. El trazador conoce la velocidad, los gradientes, los límites y la ley de Snell; todas las fases de la imagen se deducen de esas cuatro cosas.

El abanico S se detiene en el núcleo porque allí la velocidad de cizalla es cero. El abanico P hace algo más interesante: la velocidad de compresión cae de 13,66 a 8,0 km/s al entrar en el núcleo externo, así que los rayos se refractan bruscamente hacia abajo, las llegadas directas a la superficie se agotan pasados los 104° y la energía que atraviesa el núcleo reaparece como PKP mucho más lejos. Los dos bordes de la sombra de la P salen de la refracción, no de un if, que es justo por lo que hace falta comprobarlos.

Así que el motor viene con una autocomprobación ejecutable, que corre en desarrollo, con el mismo espíritu que la aserción de la ISS en el Orrery: un motor de física que nadie puede verificar a ojo necesita una prueba que sí pueda.

• El tiempo de viaje crece de forma monótona a lo largo de todos los rayos. (Un rayo que retrocede en el tiempo significa que el integrador se ha roto.)
• La P directa a 60° llega cerca de los 600 s, con ±60 s de margen; a 90°, cerca de los 780 s con ±80 s.
• La S directa a 60° llega cerca de los 1.090 s, con ±110 s de margen.
• La última llegada de P directa a la superficie cae dentro de 90°–118°: lo bastante ancho para tolerar un modelo de grado visual, lo bastante estrecho para cazar una refracción rota.
• Existen llegadas de PKP más allá de 130°.
• Nunca se encuentra ningún rayo S a más de 5 km por debajo del límite núcleo–manto.

Una isócrona por fotograma, cero reservas de memoria

Un rayo trazado es un camino por el espacio y por el tiempo; un frente de onda es un corte transversal. sampleIsochron() busca por bisección, en las muestras de cada rayo, el intervalo que contiene el instante actual de reproducción e interpola una posición a partir de él — un par radio/distancia por rayo, dentro de un array reservado de antemano. Los rayos que nunca alcanzaron ese instante, o que murieron en el núcleo, escriben (NaN, NaN).

Los frentes se dibujan como LineSegments y no como una polilínea, y esa elección sostiene el resto. Una polilínea uniría la última muestra viva antes de un hueco con la primera de después, soldando una cuerda que cruzaría la sombra de lado a lado. Los segmentos no lo hacen: un NaN termina una serie, el siguiente par válido empieza otra, y la sombra se talla sola en la geometría. La cinta de resplandor se construye con quads independientes por la misma razón — una tira compartida volvería a soldar el hueco que las líneas acaban de dejar.

Aquí no se reconstruye nada por fotograma: 720 huecos de muestra, 2.880 vértices de línea y 1.440 quads de cinta se reservan una sola vez, y setDrawRange hace el trabajo. Cada frente se dibuja dos veces, reflejado a ±delta, porque un frente de onda en una Tierra de simetría esférica es una superficie de revolución en torno al eje de la fuente y el corte expone dos semiplanos de ella. Los frentes aguantan 120 s de tiempo simulado y se apagan al envejecer — la P del casi blanco a la brasa, la S del menta pálido al verde azulado — para que una reproducción movida se lea como una secuencia y no como una mancha.

Una Tierra en corte hecha de geometría, no de recorte

Para ver un frente de onda cruzar el interior hay que abrir el planeta, y la vía obvia — los planos de recorte — es la herramienta equivocada para las capas: trabajo de fragmento en cada píxel de cada capa, y caras huecas al final del proceso. En vez de eso, las capas se construyen incompletas. Una SphereGeometry con phiStart = π y phiLength = 1.5π es una capa de 270° a la que le falta un cuadrante, orientada hacia la cámara por defecto. El agujero está en la malla, así que no cuesta nada.

Hay una capa que no puede usar ese truco. Las costas y la retícula tienen que girar — eso es justamente para lo que están —, así que son lo único de la escena que sí se recorta de verdad en la GPU, mediante dos planos con clipIntersection activado. Los datos de costas se decodifican del módulo que empaquetamos para el Orrery y se reinterpolan esféricamente a lo largo de círculos máximos, con unos 2° por subdivisión, para que las cuerdas largas no se hundan bajo la superficie y desaparezcan dentro de la esfera.

La escena funciona entonces con dos sistemas de referencia, y esta es la pieza de la que estamos calladamente más orgullosos. El sistema del corte no se mueve nunca: la cuña, los estratos, las líneas finas de los límites y el eje polar son estáticos, así que los planos de recorte en coordenadas de mundo son constantes y el dibujo de la cara cortada se hornea en una sola malla indexada y una sola llamada de dibujo. El sistema geográfico — costas, retícula, el catálogo entero — se interpola esféricamente para que el hipocentro seleccionado cabalgue sobre el eje polar +Y, la arista que comparten ambas caras expuestas. Cualquier plano que pase por el eje de la fuente es una sección válida de círculo máximo en una Tierra de simetría esférica, y aquí resulta que dos de ellos son paredes que puedes ver. La sección transversal es honesta con el evento que elijas — e instantánea, sin ninguna animación, para quienes han pedido movimiento reducido.

Cuatro mil terremotos en una sola malla

El catálogo llega de dos feeds resumen del USGS — el último día a cualquier magnitud, más M4,5 y superiores del último mes — a través de un proxy que adelgaza cada evento a una tupla posicional y gestiona su propia caché. Una revalidación a nivel de ruta cachearía un 503 con la misma alegría que los datos, algo que en el Orrery llegó a congelar un estado de «sin conexión» durante seis horas. Aquí el TTL de 90 segundos se ancla a la marca de tiempo del más antiguo de los dos feeds, así que un feed caído se sigue reintentando mientras el sano sigue sirviendo. Más vale rancio que en blanco.

En pantalla, el catálogo entero — hasta 4.000 eventos — es un único InstancedMesh de octaedros. Cada instancia se sitúa a un radio de (6371 − depth) / 1000 unidades de escena: a su profundidad hipocentral real en vez de clavada en la superficie, así que una zona de subducción se lee como una losa inclinada hacia el manto. La escala es logarítmica en la magnitud, 0.012 + 0.02 · 1.9^(mag − 4), acotada por ambos extremos: con tope para que un M8 siga siendo una cuenta y no una luna, y con suelo para que los microterremotos sigan existiendo. El color es la edad: naranja intenso hacia brasa durante las primeras 24 horas, y luego de brasa a pizarra a lo largo de una cola de 30 días. La historia de la última semana se ve antes de hacer clic en nada.

Sonificación de terremotos, con un decodificador miniSEED escrito a mano

Los sismogramas son públicos. Cualquiera puede descargar movimiento del suelo en bruto de los servicios dataselect de la FDSN, y las estaciones IU de la Global Seismographic Network llevan décadas grabando sin interrupción. Así que lo último que Tremor hace con un evento es ofrecerte reproducirlo: sonificación de terremotos real, no un sintetizador imitando uno.

El navegador no puede hacerlo por sí solo, por dos razones independientes. Los datos llegan como miniSEED comprimido con Steim, un formato binario de codificación diferencial del que ningún navegador ha oído hablar, y el servicio no envía cabeceras CORS, así que una página no podría ni siquiera descargar los bytes aunque supiera leerlos. Ambas mitades viven en nuestro servidor: un proxy descarga, un decodificador escrito a mano desempaqueta y un búfer Int16 normalizado llega al cliente.

Elegir estación es un pequeño problema de ordenación con una restricción física dentro. Tremor conoce 31 estaciones IU de banda ancha y quiere cuarenta minutos de registro en la componente vertical, pero la más cercana no es la mejor: por debajo de 2° el instrumento satura con un evento grande, y pasados los 90° el núcleo empieza a ensombrecer la P. Así que la banda telesísmica encabeza la ordenación, la distancia desempata, se prueban tres candidatas por orden y un registro de menos de 45 minutos se rechaza de plano: los archivos necesitan tiempo, y el cliente se adelanta a ese rechazo en vez de dejar que lo descubras tú.

El decodificador miniSEED es prácticamente lo más de bajo nivel que puede ser el código web: un DataView, lecturas big-endian y una cabecera blockette-1000 que dice qué codificación tiene delante. Steim empaqueta las muestras en tramas de 64 bytes con dieciséis palabras de 32 bits, la primera de las cuales es un mapa de nibbles — dieciséis códigos de dos bits que describen cómo están empaquetadas las otras quince. Steim-1 ofrece una, dos o cuatro diferencias por palabra; Steim-2 añade subcódigos que meten cinco diferencias de 6 bits, seis de 5 bits o siete de 4 bits en una palabra, más una ranura de 30 bits para valores atípicos. Las palabras uno y dos de la trama cero no son diferencias en absoluto, sino las constantes de integración X0 y Xn, lo que convierte al flujo en autocomprobable. Ambas codificaciones están implementadas; cuando la sincronización entre registros se desvía más de 1,5 muestras contamos el hueco y concatenamos igualmente, porque para la sonificación la continuidad importa más que la exactitud.

Lo que nos lleva a la parte que no tiene código. La energía sísmica vive en torno a 0,05–10 Hz; la audición empieza cerca de los 20 Hz. La solución de manual es un desplazador de tono, y no escribimos ni una línea de uno, porque la Web Audio API te deja declarar la frecuencia de muestreo de un búfer. Cuarenta muestras por segundo, declaradas como 40 × 150 = 6.000 Hz, se reproducen 150 veces más rápido y 150 veces más agudo — dentro de la ventana legal de 3.000 a 768.000 Hz — y dejan la banda sísmica en 7,5–1500 Hz. Cuarenta minutos se vuelven dieciséis segundos. No es una transformación de los datos: son los datos, leídos deprisa.

El rango dinámico es el problema difícil. La llegada de una P telesísmica puede ser cien veces la coda que la sigue, y reproducida linealmente el visitante oye un tic y luego silencio. Así que el cliente aplica compansión con un exponente cercano a 0,6 — en el navegador, no en el servidor, porque el trabajo del servidor son cuentas honestas y el modelado corresponde al sitio donde el visitante puede oír cómo se aplica. Antes de transportarlo, al registro se le resta la media (SNZO, por ejemplo, está unos 1.500 cuentas por encima de cero) y se normaliza por pico. La minigráfica dibuja el pico por intervalo en vez del RMS, por la misma razón: la llegada es un pico, y promediar es la única operación capaz de esconderlo.

El reproductor en sí es deliberadamente pequeño: un GainNode, rampas de 30 ms programadas con una de antelación, un contexto de audio creado dentro del manejador del clic para que la política de reproducción automática no tenga nada que objetar, y una reproducción que se detiene cuando la pestaña desaparece. El silencio no es propio de Tremor: lee el silenciado global del sintetizador que comparte con todos los juegos del Laboratorio. Silencia un juego y el planeta se calla también.

Una etiqueta en una esquina de la pantalla dice tiempos de viaje de grado visual, por la misma razón que la del Orrery dice propagación de grado visual. Un modelo unidimensional sin correcciones de estación ni estructura tridimensional pone el frente de onda donde el ojo lo espera, no donde lo registraría un instrumento. Fiel al ojo, no al instrumento; no lo uses nunca para nada que no sea asombrarte. Decirlo con claridad no cuesta nada y le compra credibilidad a todo lo demás.

Las lagunas honestas: el trazado de rayos corre en el hilo principal, y la página todavía no tiene cobertura automatizada de extremo a extremo — solo las aserciones de desarrollo, que comprueban la física y no los botones. Ambas cosas están en la lista.