diff --git a/Tutorial/SpecModTutorial.ipynb b/Tutorial/SpecModTutorial.ipynb index 38c4665..f86d11d 100644 --- a/Tutorial/SpecModTutorial.ipynb +++ b/Tutorial/SpecModTutorial.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "39eda767", + "id": "80661db4", "metadata": {}, "source": [ "# SpecMod tutorial\n", @@ -29,7 +29,7 @@ }, { "cell_type": "markdown", - "id": "59f58548", + "id": "01a3ea27", "metadata": {}, "source": [ "## 1. Read and prepare" @@ -38,13 +38,13 @@ { "cell_type": "code", "execution_count": 1, - "id": "e3afbe8f", + "id": "b06e3ef2", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:20.194120Z", - "iopub.status.busy": "2026-08-08T16:05:20.193836Z", - "iopub.status.idle": "2026-08-08T16:05:21.133218Z", - "shell.execute_reply": "2026-08-08T16:05:21.131949Z" + "iopub.execute_input": "2026-08-08T16:49:51.779439Z", + "iopub.status.busy": "2026-08-08T16:49:51.779238Z", + "iopub.status.idle": "2026-08-08T16:49:53.063613Z", + "shell.execute_reply": "2026-08-08T16:49:53.062068Z" } }, "outputs": [], @@ -65,7 +65,7 @@ }, { "cell_type": "markdown", - "id": "4c94213f", + "id": "6606fad9", "metadata": {}, "source": [ "The origin is what distances and theoretical arrivals are measured from.\n", @@ -82,13 +82,13 @@ { "cell_type": "code", "execution_count": 2, - "id": "f1bc897f", + "id": "864322f1", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.135876Z", - "iopub.status.busy": "2026-08-08T16:05:21.135432Z", - "iopub.status.idle": "2026-08-08T16:05:21.139895Z", - "shell.execute_reply": "2026-08-08T16:05:21.138517Z" + "iopub.execute_input": "2026-08-08T16:49:53.066533Z", + "iopub.status.busy": "2026-08-08T16:49:53.066048Z", + "iopub.status.idle": "2026-08-08T16:49:53.071135Z", + "shell.execute_reply": "2026-08-08T16:49:53.069373Z" } }, "outputs": [], @@ -100,13 +100,13 @@ { "cell_type": "code", "execution_count": 3, - "id": "732e9a38", + "id": "b4e74062", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.141881Z", - "iopub.status.busy": "2026-08-08T16:05:21.141655Z", - "iopub.status.idle": "2026-08-08T16:05:21.243229Z", - "shell.execute_reply": "2026-08-08T16:05:21.242244Z" + "iopub.execute_input": "2026-08-08T16:49:53.073283Z", + "iopub.status.busy": "2026-08-08T16:49:53.073067Z", + "iopub.status.idle": "2026-08-08T16:49:53.178931Z", + "shell.execute_reply": "2026-08-08T16:49:53.177544Z" } }, "outputs": [ @@ -131,13 +131,13 @@ { "cell_type": "code", "execution_count": 4, - "id": "86bd514a", + "id": "7f37be38", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.245322Z", - "iopub.status.busy": "2026-08-08T16:05:21.244903Z", - "iopub.status.idle": "2026-08-08T16:05:21.254042Z", - "shell.execute_reply": "2026-08-08T16:05:21.252921Z" + "iopub.execute_input": "2026-08-08T16:49:53.181679Z", + "iopub.status.busy": "2026-08-08T16:49:53.181107Z", + "iopub.status.idle": "2026-08-08T16:49:53.191799Z", + "shell.execute_reply": "2026-08-08T16:49:53.190749Z" } }, "outputs": [ @@ -176,7 +176,7 @@ }, { "cell_type": "markdown", - "id": "fac2ce5b", + "id": "fbef660f", "metadata": {}, "source": [ "**Instrument correction happens here, in the notebook, not inside the\n", @@ -192,13 +192,13 @@ { "cell_type": "code", "execution_count": 5, - "id": "80933dcc", + "id": "df25562f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.256143Z", - "iopub.status.busy": "2026-08-08T16:05:21.255921Z", - "iopub.status.idle": "2026-08-08T16:05:21.765679Z", - "shell.execute_reply": "2026-08-08T16:05:21.764595Z" + "iopub.execute_input": "2026-08-08T16:49:53.194187Z", + "iopub.status.busy": "2026-08-08T16:49:53.193855Z", + "iopub.status.idle": "2026-08-08T16:49:53.782359Z", + "shell.execute_reply": "2026-08-08T16:49:53.781148Z" } }, "outputs": [ @@ -239,13 +239,13 @@ { "cell_type": "code", "execution_count": 6, - "id": "a4af3926", + "id": "396193c1", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.767678Z", - "iopub.status.busy": "2026-08-08T16:05:21.767471Z", - "iopub.status.idle": "2026-08-08T16:05:23.302059Z", - "shell.execute_reply": "2026-08-08T16:05:23.300791Z" + "iopub.execute_input": "2026-08-08T16:49:53.784573Z", + "iopub.status.busy": "2026-08-08T16:49:53.784327Z", + "iopub.status.idle": "2026-08-08T16:49:55.334948Z", + "shell.execute_reply": "2026-08-08T16:49:55.333489Z" } }, "outputs": [], @@ -256,7 +256,7 @@ }, { "cell_type": "markdown", - "id": "23e4c246", + "id": "1286883e", "metadata": {}, "source": [ "## 2. Cut the windows" @@ -264,7 +264,7 @@ }, { "cell_type": "markdown", - "id": "a9667906", + "id": "aeb08889", "metadata": {}, "source": [ "The S-window opens at a fixed fraction of the elapsed P–S time after the P\n", @@ -280,13 +280,13 @@ { "cell_type": "code", "execution_count": 7, - "id": "56e0fd1c", + "id": "45304586", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:23.305017Z", - "iopub.status.busy": "2026-08-08T16:05:23.304814Z", - "iopub.status.idle": "2026-08-08T16:05:23.322224Z", - "shell.execute_reply": "2026-08-08T16:05:23.321012Z" + "iopub.execute_input": "2026-08-08T16:49:55.337432Z", + "iopub.status.busy": "2026-08-08T16:49:55.337141Z", + "iopub.status.idle": "2026-08-08T16:49:55.358539Z", + "shell.execute_reply": "2026-08-08T16:49:55.357428Z" } }, "outputs": [], @@ -300,13 +300,13 @@ { "cell_type": "code", "execution_count": 8, - "id": "a4fe9b32", + "id": "720530dd", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:23.324105Z", - "iopub.status.busy": "2026-08-08T16:05:23.323905Z", - "iopub.status.idle": "2026-08-08T16:05:23.328110Z", - "shell.execute_reply": "2026-08-08T16:05:23.327014Z" + "iopub.execute_input": "2026-08-08T16:49:55.361393Z", + "iopub.status.busy": "2026-08-08T16:49:55.361161Z", + "iopub.status.idle": "2026-08-08T16:49:55.366743Z", + "shell.execute_reply": "2026-08-08T16:49:55.365013Z" } }, "outputs": [ @@ -338,13 +338,13 @@ { "cell_type": "code", "execution_count": 9, - "id": "e9929541", + "id": "9381e25c", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:23.329901Z", - "iopub.status.busy": "2026-08-08T16:05:23.329699Z", - "iopub.status.idle": "2026-08-08T16:05:25.155951Z", - "shell.execute_reply": "2026-08-08T16:05:25.154634Z" + "iopub.execute_input": "2026-08-08T16:49:55.368922Z", + "iopub.status.busy": "2026-08-08T16:49:55.368732Z", + "iopub.status.idle": "2026-08-08T16:49:56.997523Z", + "shell.execute_reply": "2026-08-08T16:49:56.996091Z" } }, "outputs": [], @@ -355,7 +355,7 @@ }, { "cell_type": "markdown", - "id": "a8dd38a4", + "id": "57733e86", "metadata": {}, "source": [ "## 3. Spectra and bandwidth" @@ -363,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "e3fee1df", + "id": "08db6076", "metadata": {}, "source": [ "`spectrum_set_from_streams` transforms both windows, puts the noise on the\n", @@ -378,13 +378,13 @@ { "cell_type": "code", "execution_count": 10, - "id": "87aebed5", + "id": "0c1fa24f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.158367Z", - "iopub.status.busy": "2026-08-08T16:05:25.158112Z", - "iopub.status.idle": "2026-08-08T16:05:25.355124Z", - "shell.execute_reply": "2026-08-08T16:05:25.353938Z" + "iopub.execute_input": "2026-08-08T16:49:57.000245Z", + "iopub.status.busy": "2026-08-08T16:49:56.999975Z", + "iopub.status.idle": "2026-08-08T16:49:57.201338Z", + "shell.execute_reply": "2026-08-08T16:49:57.199880Z" } }, "outputs": [ @@ -408,13 +408,13 @@ { "cell_type": "code", "execution_count": 11, - "id": "b93a32ef", + "id": "1e547105", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.357445Z", - "iopub.status.busy": "2026-08-08T16:05:25.357245Z", - "iopub.status.idle": "2026-08-08T16:05:25.377656Z", - "shell.execute_reply": "2026-08-08T16:05:25.376409Z" + "iopub.execute_input": "2026-08-08T16:49:57.204138Z", + "iopub.status.busy": "2026-08-08T16:49:57.203843Z", + "iopub.status.idle": "2026-08-08T16:49:57.231532Z", + "shell.execute_reply": "2026-08-08T16:49:57.230099Z" } }, "outputs": [], @@ -431,13 +431,13 @@ { "cell_type": "code", "execution_count": 12, - "id": "a93bce39", + "id": "5af0fd70", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.379784Z", - "iopub.status.busy": "2026-08-08T16:05:25.379584Z", - "iopub.status.idle": "2026-08-08T16:05:25.383672Z", - "shell.execute_reply": "2026-08-08T16:05:25.382549Z" + "iopub.execute_input": "2026-08-08T16:49:57.234367Z", + "iopub.status.busy": "2026-08-08T16:49:57.234148Z", + "iopub.status.idle": "2026-08-08T16:49:57.239171Z", + "shell.execute_reply": "2026-08-08T16:49:57.237476Z" } }, "outputs": [ @@ -460,7 +460,7 @@ }, { "cell_type": "markdown", - "id": "dd6078f8", + "id": "cb81de32", "metadata": {}, "source": [ "### Changing ground-motion domain\n", @@ -473,13 +473,13 @@ { "cell_type": "code", "execution_count": 13, - "id": "bc8fef2a", + "id": "508d3f73", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.385510Z", - "iopub.status.busy": "2026-08-08T16:05:25.385339Z", - "iopub.status.idle": "2026-08-08T16:05:25.424633Z", - "shell.execute_reply": "2026-08-08T16:05:25.423478Z" + "iopub.execute_input": "2026-08-08T16:49:57.241350Z", + "iopub.status.busy": "2026-08-08T16:49:57.241166Z", + "iopub.status.idle": "2026-08-08T16:49:57.280015Z", + "shell.execute_reply": "2026-08-08T16:49:57.278838Z" } }, "outputs": [ @@ -507,7 +507,7 @@ }, { "cell_type": "markdown", - "id": "1d8f7fea", + "id": "75de79d6", "metadata": {}, "source": [ "## 4. Fit a source model" @@ -515,7 +515,7 @@ }, { "cell_type": "markdown", - "id": "3dfca7b2", + "id": "b92f4418", "metadata": {}, "source": [ "The model comes from configuration — a Brune source with constant Q by\n", @@ -532,13 +532,13 @@ { "cell_type": "code", "execution_count": 14, - "id": "b5ac93c1", + "id": "000db0ef", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.426632Z", - "iopub.status.busy": "2026-08-08T16:05:25.426437Z", - "iopub.status.idle": "2026-08-08T16:05:25.477264Z", - "shell.execute_reply": "2026-08-08T16:05:25.476242Z" + "iopub.execute_input": "2026-08-08T16:49:57.282340Z", + "iopub.status.busy": "2026-08-08T16:49:57.282151Z", + "iopub.status.idle": "2026-08-08T16:49:57.318987Z", + "shell.execute_reply": "2026-08-08T16:49:57.317867Z" } }, "outputs": [ @@ -563,13 +563,13 @@ { "cell_type": "code", "execution_count": 15, - "id": "f1e68721", + "id": "6abdc313", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.479686Z", - "iopub.status.busy": "2026-08-08T16:05:25.479095Z", - "iopub.status.idle": "2026-08-08T16:05:26.383789Z", - "shell.execute_reply": "2026-08-08T16:05:26.382633Z" + "iopub.execute_input": "2026-08-08T16:49:57.322514Z", + "iopub.status.busy": "2026-08-08T16:49:57.321509Z", + "iopub.status.idle": "2026-08-08T16:49:58.166284Z", + "shell.execute_reply": "2026-08-08T16:49:58.164994Z" } }, "outputs": [ @@ -589,7 +589,7 @@ }, { "cell_type": "markdown", - "id": "6a7eb72f", + "id": "d19e4f9b", "metadata": {}, "source": [ "The guess is only a starting point, and a crude one: it takes the largest\n", @@ -602,13 +602,13 @@ { "cell_type": "code", "execution_count": 16, - "id": "03f89e94", + "id": "f33f5eda", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:26.385827Z", - "iopub.status.busy": "2026-08-08T16:05:26.385633Z", - "iopub.status.idle": "2026-08-08T16:05:26.401731Z", - "shell.execute_reply": "2026-08-08T16:05:26.400676Z" + "iopub.execute_input": "2026-08-08T16:49:58.168880Z", + "iopub.status.busy": "2026-08-08T16:49:58.168643Z", + "iopub.status.idle": "2026-08-08T16:49:58.188051Z", + "shell.execute_reply": "2026-08-08T16:49:58.186792Z" } }, "outputs": [ @@ -725,13 +725,13 @@ { "cell_type": "code", "execution_count": 17, - "id": "15af287a", + "id": "095d9120", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:26.403508Z", - "iopub.status.busy": "2026-08-08T16:05:26.403329Z", - "iopub.status.idle": "2026-08-08T16:05:26.421423Z", - "shell.execute_reply": "2026-08-08T16:05:26.420372Z" + "iopub.execute_input": "2026-08-08T16:49:58.190521Z", + "iopub.status.busy": "2026-08-08T16:49:58.190301Z", + "iopub.status.idle": "2026-08-08T16:49:58.210382Z", + "shell.execute_reply": "2026-08-08T16:49:58.208850Z" } }, "outputs": [], @@ -744,13 +744,13 @@ { "cell_type": "code", "execution_count": 18, - "id": "704b062a", + "id": "514a66be", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:26.423709Z", - "iopub.status.busy": "2026-08-08T16:05:26.423496Z", - "iopub.status.idle": "2026-08-08T16:05:30.421284Z", - "shell.execute_reply": "2026-08-08T16:05:30.419781Z" + "iopub.execute_input": "2026-08-08T16:49:58.212770Z", + "iopub.status.busy": "2026-08-08T16:49:58.212534Z", + "iopub.status.idle": "2026-08-08T16:50:01.491642Z", + "shell.execute_reply": "2026-08-08T16:50:01.490434Z" } }, "outputs": [], @@ -763,13 +763,13 @@ { "cell_type": "code", "execution_count": 19, - "id": "c8d5cdd1", + "id": "54bf2daa", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:30.424820Z", - "iopub.status.busy": "2026-08-08T16:05:30.424540Z", - "iopub.status.idle": "2026-08-08T16:05:30.439861Z", - "shell.execute_reply": "2026-08-08T16:05:30.438685Z" + "iopub.execute_input": "2026-08-08T16:50:01.493730Z", + "iopub.status.busy": "2026-08-08T16:50:01.493496Z", + "iopub.status.idle": "2026-08-08T16:50:01.505023Z", + "shell.execute_reply": "2026-08-08T16:50:01.503985Z" } }, "outputs": [ @@ -922,7 +922,7 @@ }, { "cell_type": "markdown", - "id": "4f76ba73", + "id": "1d7cdb71", "metadata": {}, "source": [ "### The fit is not unique, and that is not a detail\n", @@ -946,13 +946,13 @@ { "cell_type": "code", "execution_count": 20, - "id": "4c83be9e", + "id": "9b806a45", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:30.442104Z", - "iopub.status.busy": "2026-08-08T16:05:30.441839Z", - "iopub.status.idle": "2026-08-08T16:05:31.463489Z", - "shell.execute_reply": "2026-08-08T16:05:31.462097Z" + "iopub.execute_input": "2026-08-08T16:50:01.507171Z", + "iopub.status.busy": "2026-08-08T16:50:01.506985Z", + "iopub.status.idle": "2026-08-08T16:50:02.470569Z", + "shell.execute_reply": "2026-08-08T16:50:02.469397Z" } }, "outputs": [ @@ -1086,7 +1086,7 @@ }, { "cell_type": "markdown", - "id": "0424928a", + "id": "c2c148f3", "metadata": {}, "source": [ "Most stations agree to a fraction of a percent. A few do not, and the top of\n", @@ -1102,13 +1102,13 @@ { "cell_type": "code", "execution_count": 21, - "id": "372771f9", + "id": "493d0dea", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.466011Z", - "iopub.status.busy": "2026-08-08T16:05:31.465712Z", - "iopub.status.idle": "2026-08-08T16:05:31.471262Z", - "shell.execute_reply": "2026-08-08T16:05:31.470165Z" + "iopub.execute_input": "2026-08-08T16:50:02.472588Z", + "iopub.status.busy": "2026-08-08T16:50:02.472382Z", + "iopub.status.idle": "2026-08-08T16:50:02.477907Z", + "shell.execute_reply": "2026-08-08T16:50:02.476853Z" } }, "outputs": [ @@ -1139,13 +1139,13 @@ { "cell_type": "code", "execution_count": 22, - "id": "18e64e8f", + "id": "056e59ac", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.473285Z", - "iopub.status.busy": "2026-08-08T16:05:31.473029Z", - "iopub.status.idle": "2026-08-08T16:05:31.491080Z", - "shell.execute_reply": "2026-08-08T16:05:31.489606Z" + "iopub.execute_input": "2026-08-08T16:50:02.480019Z", + "iopub.status.busy": "2026-08-08T16:50:02.479814Z", + "iopub.status.idle": "2026-08-08T16:50:02.500292Z", + "shell.execute_reply": "2026-08-08T16:50:02.498830Z" } }, "outputs": [], @@ -1161,7 +1161,7 @@ }, { "cell_type": "markdown", - "id": "f52d1081", + "id": "171b6750", "metadata": {}, "source": [ "Look at where they differ: high up the falling limb, where the source corner\n", @@ -1186,13 +1186,13 @@ { "cell_type": "code", "execution_count": 23, - "id": "ea8f5a9b", + "id": "f4f10578", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.493795Z", - "iopub.status.busy": "2026-08-08T16:05:31.493595Z", - "iopub.status.idle": "2026-08-08T16:05:31.498573Z", - "shell.execute_reply": "2026-08-08T16:05:31.497367Z" + "iopub.execute_input": "2026-08-08T16:50:02.503081Z", + "iopub.status.busy": "2026-08-08T16:50:02.502866Z", + "iopub.status.idle": "2026-08-08T16:50:02.508180Z", + "shell.execute_reply": "2026-08-08T16:50:02.507142Z" } }, "outputs": [ @@ -1214,7 +1214,7 @@ }, { "cell_type": "markdown", - "id": "f10f5659", + "id": "4ce09cca", "metadata": {}, "source": [ "### Why the answer comes from many stations, and from two stages\n", @@ -1239,13 +1239,13 @@ { "cell_type": "code", "execution_count": 24, - "id": "82e8c326", + "id": "6b32068b", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.500823Z", - "iopub.status.busy": "2026-08-08T16:05:31.500622Z", - "iopub.status.idle": "2026-08-08T16:05:31.512549Z", - "shell.execute_reply": "2026-08-08T16:05:31.511462Z" + "iopub.execute_input": "2026-08-08T16:50:02.510699Z", + "iopub.status.busy": "2026-08-08T16:50:02.510493Z", + "iopub.status.idle": "2026-08-08T16:50:02.522429Z", + "shell.execute_reply": "2026-08-08T16:50:02.521541Z" } }, "outputs": [ @@ -1253,20 +1253,28 @@ "name": "stdout", "output_type": "stream", "text": [ - "powell event fc = 12.751 Hz\n", - "leastsq event fc = 12.675 Hz\n", + "powell event fc = 11.585 Hz\n", + "leastsq event fc = 11.535 Hz\n", "\n", "worst single station: fc ratio 1.441 -> stress drop 2.99x\n", - "event ensemble: fc ratio 1.006 -> stress drop 1.018x\n" + "event ensemble: fc ratio 1.004 -> stress drop 1.013x\n" ] } ], "source": [ - "# Weighted by inverse hypocentral distance: the nearer station has less path\n", - "# between the source and the sensor, so less of its high-frequency falloff can\n", - "# be attenuation, and its corner is the better constrained of the two.\n", + "# Weighted by inverse distance: the nearer station has less path between the\n", + "# source and the sensor, so less of its high-frequency falloff can be\n", + "# attenuation, and its corner is the better constrained of the two.\n", + "#\n", + "# `repi` — epicentral — because that is what `[windows] distance_metric` says.\n", + "# Which distance you use is not a detail at short range: here the nearest\n", + "# station is 0.89 km epicentral against 2.30 km hypocentral. `rhyp` is built\n", + "# from the source depth and the station *elevation*, so it assumes every\n", + "# sensor is at the surface; where sensor depths are unknown, as they are here,\n", + "# epicentral is the honest choice.\n", + "metric = \"repi\"\n", "weight = pd.Series(\n", - " {id: 1.0 / spectra[id].signal.meta[\"rhyp\"] for id in fits.models}\n", + " {id: 1.0 / spectra[id].signal.meta[metric] for id in fits.models}\n", ")\n", "\n", "event_fc = {}\n", @@ -1287,7 +1295,7 @@ }, { "cell_type": "markdown", - "id": "be9c2b9a", + "id": "c308fef1", "metadata": {}, "source": [ "A factor of three in stress drop at the worst station becomes **under 2%**\n", @@ -1302,13 +1310,13 @@ { "cell_type": "code", "execution_count": 25, - "id": "d0e5ba30", + "id": "d67f1cef", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.514662Z", - "iopub.status.busy": "2026-08-08T16:05:31.514484Z", - "iopub.status.idle": "2026-08-08T16:05:32.299622Z", - "shell.execute_reply": "2026-08-08T16:05:32.298368Z" + "iopub.execute_input": "2026-08-08T16:50:02.524750Z", + "iopub.status.busy": "2026-08-08T16:50:02.524533Z", + "iopub.status.idle": "2026-08-08T16:50:03.295928Z", + "shell.execute_reply": "2026-08-08T16:50:03.294940Z" } }, "outputs": [ @@ -1316,8 +1324,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "t* worst disagreement between minimisers: 0.31%\n", - "omega worst disagreement between minimisers: 2.3e-03 log10 units\n" + "t* worst disagreement between minimisers: 0.23%\n", + "omega worst disagreement between minimisers: 1.8e-03 log10 units\n" ] }, { @@ -1357,45 +1365,45 @@ " \n", " \n", " LV.L001..HHE\n", - " 0.0338\n", - " 0.0337\n", - " -4.9827\n", - " -4.9818\n", + " 0.0326\n", + " 0.0326\n", + " -4.9657\n", + " -4.9649\n", " \n", " \n", " LV.L002..HHE\n", - " 0.0547\n", - " 0.0546\n", - " -5.3057\n", - " -5.3053\n", + " 0.0529\n", + " 0.0529\n", + " -5.2984\n", + " -5.2979\n", " \n", " \n", " LV.L006..HHE\n", - " 0.0834\n", - " 0.0833\n", - " -5.1199\n", - " -5.1193\n", + " 0.0820\n", + " 0.0819\n", + " -5.1084\n", + " -5.1078\n", " \n", " \n", " LV.L007..HHE\n", - " 0.0774\n", - " 0.0772\n", - " -5.6150\n", - " -5.6150\n", + " 0.0752\n", + " 0.0751\n", + " -5.6139\n", + " -5.6138\n", " \n", " \n", " LV.L008..HHE\n", - " 0.0559\n", - " 0.0557\n", - " -5.8243\n", - " -5.8246\n", + " 0.0533\n", + " 0.0532\n", + " -5.8285\n", + " -5.8286\n", " \n", " \n", " LV.L009..HHE\n", - " 0.0439\n", - " 0.0439\n", - " -5.4172\n", - " -5.4162\n", + " 0.0427\n", + " 0.0427\n", + " -5.4009\n", + " -5.4001\n", " \n", " \n", "\n", @@ -1404,12 +1412,12 @@ "text/plain": [ " t* powell t* leastsq log10 omega powell log10 omega leastsq\n", "id \n", - "LV.L001..HHE 0.0338 0.0337 -4.9827 -4.9818\n", - "LV.L002..HHE 0.0547 0.0546 -5.3057 -5.3053\n", - "LV.L006..HHE 0.0834 0.0833 -5.1199 -5.1193\n", - "LV.L007..HHE 0.0774 0.0772 -5.6150 -5.6150\n", - "LV.L008..HHE 0.0559 0.0557 -5.8243 -5.8246\n", - "LV.L009..HHE 0.0439 0.0439 -5.4172 -5.4162" + "LV.L001..HHE 0.0326 0.0326 -4.9657 -4.9649\n", + "LV.L002..HHE 0.0529 0.0529 -5.2984 -5.2979\n", + "LV.L006..HHE 0.0820 0.0819 -5.1084 -5.1078\n", + "LV.L007..HHE 0.0752 0.0751 -5.6139 -5.6138\n", + "LV.L008..HHE 0.0533 0.0532 -5.8285 -5.8286\n", + "LV.L009..HHE 0.0427 0.0427 -5.4009 -5.4001" ] }, "execution_count": 25, @@ -1441,7 +1449,7 @@ }, { "cell_type": "markdown", - "id": "7f80cf74", + "id": "7efdbc23", "metadata": {}, "source": [ "0.3% in $t^*$ and about 0.002 in $\\log_{10}\\Omega$ — which is 0.003 magnitude\n", @@ -1463,7 +1471,7 @@ }, { "cell_type": "markdown", - "id": "4ab068b5", + "id": "9bb2ff9c", "metadata": {}, "source": [ "### The same thing, as one call\n", @@ -1482,13 +1490,13 @@ { "cell_type": "code", "execution_count": 26, - "id": "0082031d", + "id": "34091db5", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:32.301740Z", - "iopub.status.busy": "2026-08-08T16:05:32.301520Z", - "iopub.status.idle": "2026-08-08T16:05:33.379213Z", - "shell.execute_reply": "2026-08-08T16:05:33.377893Z" + "iopub.execute_input": "2026-08-08T16:50:03.298051Z", + "iopub.status.busy": "2026-08-08T16:50:03.297784Z", + "iopub.status.idle": "2026-08-08T16:50:04.380268Z", + "shell.execute_reply": "2026-08-08T16:50:04.378693Z" } }, "outputs": [ @@ -1496,11 +1504,11 @@ "name": "stdout", "output_type": "stream", "text": [ - "fc = 12.75 from 28 channels, weighted by inverse_distance\n", - " stage-1 range 2.835 to 59.64 (445% of the event value)\n", + "fc = 11.59 from 28 channels, weighted by inverse_distance\n", + " stage-1 range 2.835 to 59.64 (490% of the event value)\n", "\n", - "by hand : 12.7512 Hz\n", - "API : 12.7512 Hz\n", + "by hand : 11.5851 Hz\n", + "API : 11.5851 Hz\n", "agree : True\n" ] } @@ -1518,7 +1526,7 @@ }, { "cell_type": "markdown", - "id": "a1c9c1a4", + "id": "a8ef817f", "metadata": {}, "source": [ "`describe()` prints the spread as well as the mean, and that is deliberate. A\n", @@ -1534,21 +1542,26 @@ "wrong, and averaging it in moves the event value for every other station.\n", "\n", "Two things make that lever bigger than it first looks. Inverse-distance\n", - "weighting is concentrated — the nearest two channels carry about a fifth of\n", - "the total weight — and stress drop goes as $f_c^3$, so a modest shift in the\n", - "corner is a large shift in the thing you are reporting." + "weighting is concentrated — on epicentral distance here the nearest two\n", + "channels carry over 40% of the total weight — and stress drop goes as $f_c^3$,\n", + "so a modest shift in the corner is a large shift in the thing you report.\n", + "\n", + "Which distance measure you choose feeds straight into this, and\n", + "`specmod.distance` makes it a registry for that reason: `repi` and `rhyp` are\n", + "implemented, and `rrup`/`rjb` are registered but raise, since a point source\n", + "has no rupture surface to measure from." ] }, { "cell_type": "code", "execution_count": 27, - "id": "23a48dfd", + "id": "4d563fb6", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:33.381166Z", - "iopub.status.busy": "2026-08-08T16:05:33.380995Z", - "iopub.status.idle": "2026-08-08T16:05:33.386019Z", - "shell.execute_reply": "2026-08-08T16:05:33.384855Z" + "iopub.execute_input": "2026-08-08T16:50:04.382465Z", + "iopub.status.busy": "2026-08-08T16:50:04.382271Z", + "iopub.status.idle": "2026-08-08T16:50:04.387143Z", + "shell.execute_reply": "2026-08-08T16:50:04.385991Z" } }, "outputs": [ @@ -1557,12 +1570,12 @@ "output_type": "stream", "text": [ "weight carried by the nearest channels:\n", - " nearest 1: 11.4%\n", - " nearest 2: 22.8%\n", - " nearest 4: 36.1%\n", - " nearest 8: 56.3%\n", + " nearest 1: 20.7%\n", + " nearest 2: 41.4%\n", + " nearest 4: 52.5%\n", + " nearest 8: 68.2%\n", "\n", - "the single nearest is UR.AQ04.00.HHN at 2.30 km\n" + "the single nearest is UR.AQ04.00.HHN at 0.89 km (repi)\n" ] } ], @@ -1570,7 +1583,7 @@ "import numpy as np\n", "\n", "ids = list(staged.contributing)\n", - "distance = np.array([spectra[i].signal.meta[\"rhyp\"] for i in ids])\n", + "distance = np.array([spectra[i].signal.meta[metric] for i in ids])\n", "w = 1 / distance\n", "w = w / w.sum()\n", "order = np.argsort(-w)\n", @@ -1579,19 +1592,19 @@ "for k in (1, 2, 4, 8):\n", " print(f\" nearest {k:2d}: {100 * w[order[:k]].sum():5.1f}%\")\n", "print()\n", - "print(f\"the single nearest is {ids[order[0]]} at {distance[order[0]]:.2f} km\")" + "print(f\"the single nearest is {ids[order[0]]} at {distance[order[0]]:.2f} km ({metric})\")" ] }, { "cell_type": "code", "execution_count": 28, - "id": "4ad2371c", + "id": "add4bb15", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:33.387853Z", - "iopub.status.busy": "2026-08-08T16:05:33.387680Z", - "iopub.status.idle": "2026-08-08T16:05:34.417373Z", - "shell.execute_reply": "2026-08-08T16:05:34.416061Z" + "iopub.execute_input": "2026-08-08T16:50:04.389289Z", + "iopub.status.busy": "2026-08-08T16:50:04.389072Z", + "iopub.status.idle": "2026-08-08T16:50:05.666198Z", + "shell.execute_reply": "2026-08-08T16:50:05.664533Z" } }, "outputs": [ @@ -1599,10 +1612,10 @@ "name": "stdout", "output_type": "stream", "text": [ - "all channels fc = 12.751 Hz (28 channels)\n", - "without AQ04 fc = 14.774 Hz (26 channels)\n", - " change in fc +15.9%\n", - " change in stress drop 1.56x\n", + "all channels fc = 11.585 Hz (28 channels)\n", + "without AQ04 fc = 15.602 Hz (26 channels)\n", + " change in fc +34.7%\n", + " change in stress drop 2.44x\n", "\n", " UR.AQ04.00.HHE: matched exclude='AQ04' at station\n", " UR.AQ04.00.HHN: matched exclude='AQ04' at station\n" @@ -1630,7 +1643,7 @@ }, { "cell_type": "markdown", - "id": "02dddca6", + "id": "3f9b319d", "metadata": {}, "source": [ "One quality-control decision, a factor of 1.5 in stress drop. That is not an\n", @@ -1651,13 +1664,13 @@ { "cell_type": "code", "execution_count": 29, - "id": "79a89ac6", + "id": "ba244a42", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:34.419665Z", - "iopub.status.busy": "2026-08-08T16:05:34.419462Z", - "iopub.status.idle": "2026-08-08T16:05:37.185366Z", - "shell.execute_reply": "2026-08-08T16:05:37.184045Z" + "iopub.execute_input": "2026-08-08T16:50:05.668903Z", + "iopub.status.busy": "2026-08-08T16:50:05.668683Z", + "iopub.status.idle": "2026-08-08T16:50:08.758298Z", + "shell.execute_reply": "2026-08-08T16:50:08.756603Z" } }, "outputs": [ @@ -1665,14 +1678,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "powell 28 channels vote, event fc 12.751 Hz\n" + "powell 28 channels vote, event fc 11.585 Hz\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "leastsq 22 channels vote, event fc 5.228 Hz\n", + "leastsq 22 channels vote, event fc 5.407 Hz\n", "\n", "with the ensemble held fixed at all 28:\n" ] @@ -1681,14 +1694,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "powell 28 channels vote, event fc 12.751 Hz\n" + "powell 28 channels vote, event fc 11.585 Hz\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "leastsq 28 channels vote, event fc 12.675 Hz\n" + "leastsq 28 channels vote, event fc 11.535 Hz\n" ] } ], @@ -1707,7 +1720,7 @@ }, { "cell_type": "markdown", - "id": "5fcd0afb", + "id": "71acbbd4", "metadata": {}, "source": [ "So changing the minimiser changes *which stations vote*, not only how each one\n", @@ -1723,7 +1736,7 @@ }, { "cell_type": "markdown", - "id": "a083b0a6", + "id": "54c400bd", "metadata": {}, "source": [ "## 5. Save the results" @@ -1731,7 +1744,7 @@ }, { "cell_type": "markdown", - "id": "a2aebd59", + "id": "bd03f6be", "metadata": {}, "source": [ "Two formats, because the data is used two different ways.\n", @@ -1750,13 +1763,13 @@ { "cell_type": "code", "execution_count": 30, - "id": "76c8c038", + "id": "9071b6ce", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:37.189039Z", - "iopub.status.busy": "2026-08-08T16:05:37.188751Z", - "iopub.status.idle": "2026-08-08T16:05:37.383084Z", - "shell.execute_reply": "2026-08-08T16:05:37.381766Z" + "iopub.execute_input": "2026-08-08T16:50:08.762125Z", + "iopub.status.busy": "2026-08-08T16:50:08.761787Z", + "iopub.status.idle": "2026-08-08T16:50:08.943258Z", + "shell.execute_reply": "2026-08-08T16:50:08.942038Z" } }, "outputs": [ @@ -1783,13 +1796,13 @@ { "cell_type": "code", "execution_count": 31, - "id": "935c8085", + "id": "9093db2d", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:37.386461Z", - "iopub.status.busy": "2026-08-08T16:05:37.385575Z", - "iopub.status.idle": "2026-08-08T16:05:37.425304Z", - "shell.execute_reply": "2026-08-08T16:05:37.424051Z" + "iopub.execute_input": "2026-08-08T16:50:08.946449Z", + "iopub.status.busy": "2026-08-08T16:50:08.945512Z", + "iopub.status.idle": "2026-08-08T16:50:08.983333Z", + "shell.execute_reply": "2026-08-08T16:50:08.982016Z" } }, "outputs": [ @@ -1812,7 +1825,7 @@ }, { "cell_type": "markdown", - "id": "d93d8ed8", + "id": "38f5c621", "metadata": {}, "source": [ "---\n", diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 629cefc..a95dc23 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1341,6 +1341,245 @@ of each station" is something the code can express, and "rotate to transverse" cannot be written correctly until the code can ask a channel which way it points. +#### A basic Mw, and the spreading exponent + +Worth building, and nearly free once the constants have units: everything it +needs is already produced. `llpsp` from the fit is `log10(Omega)` on +**displacement** — `sources.motion_scaling` converts the displacement model to +whatever motion was recorded, so the fitted plateau is `Omega` directly rather +than a velocity plateau needing correction — and `rhyp` is on every trace. + + M0 = 4 * pi * rho * c**3 * r**n * Omega / (R_c * F) + Mw = (2/3) * (log10(M0) - 9.1) + +**The exponent `n` is 1, not 2, and this is worth stating because it is easy +to say the other one.** Body-wave *amplitude* from a point source decays as +`1/r` in a homogeneous whole space; `1/r**2` is how *energy* decays. The +moment formula corrects an amplitude, so it multiplies by `r`. + +Measured on the 28 PNR windows through the two-stage fit, with `rho = 2500` +kg/m^3, `beta = 2500` m/s, `R_c = 0.63`, `F = 2`: + +| spreading | median M0 | median Mw | +|---|---|---| +| `1/r` | 1.67e13 N m | **+2.75** | +| `1/r**2` | 1.55e17 N m | **+5.39** | +| catalogue Mw 1.6 | 3.16e11 N m | +1.60 | + +`1/r**2` is nearly four magnitude units out, which is the distance term +applied twice; `1/r` lands within reach of the constants. **The exponent is +settled by this; the absolute calibration is not.** + +**The constants and the units are settled — from Holt (2019), the author's +doctoral thesis, Chapter 1 §1.4 and Chapter 2 Eq. 2.7**, which is the +Edwards et al. (2010) spectral method this package implements. Read directly +rather than inferred: + +| symbol | value | what it is | stated in | +|---|---|---|---| +| `rho` | 2600 kg/m^3 (Utah); 2800 generic | density **at the source** | Ch. 2 §2.2; Ch. 1 §1.4 | +| `beta` | 3500 m/s generic; 3520 m/s for the Utah normalisation | S velocity **at the source** | Ch. 1 §1.4; Ch. 2 §2.3 | +| `R_0` | **1000 m** | *reference source distance* — the distance at which the source spectrum is defined | Ch. 2 §2.2 | +| `F` | 2 | free surface, for **vertically incident SH** | Ch. 2 §2.2; the SH qualifier is Ch. 1 §1.4 | +| `Theta-lambda-Phi` | **0.55** | average radiation pattern **of SH propagation** over the focal sphere | Ch. 2 §2.2; Boatwright (1978) is cited in Ch. 1 §1.4 | +| partition factor | (Boore, 2003) | splits energy between vertical and horizontal ground motion | **Ch. 1 §1.4 only — not among the constants Ch. 2 enumerates for Eq. 2.7** | + +Every row above was read from the chapters directly, and the last one is the +correction. Ch. 1 describes the generic constant `C` and lists a partition +factor among its parts; Ch. 2, defining the Utah calculation this package +follows, enumerates the constants of Eq. 2.7 as `Omega_0`, `beta` at the +source, `rho`, `R_0`, `F` and the radiation pattern — and **no partition +factor**. So the implementation should not carry one by default. If a study +wants one it is a `[model]` entry it has to set, not a constant folded into +the formula. + +`beta = 3520 m/s` is worth recording alongside the generic 3500 because it is +where the Utah numbers actually come from: Ch. 2 §2.3 normalises the observed +plateaus to a theoretical Mw 3.5 source at `VS = 3520 m/s`, the value nearest +3500 in the Herrmann et al. (2011) WUS model at source depth. A cube of a +velocity is in the moment expression, so a 0.6% difference in `beta` is a 1.7% +difference in `M0` — negligible against everything else here, and worth +knowing is negligible rather than assuming so. + +Two things this corrects in earlier drafts of this section. + +**`0.55` is not a different average of the same quantity — it is a different +quantity.** It is the **SH** radiation pattern, which pairs with `F = 2` for +vertically incident SH. The textbook 0.63 is the RMS over total S. So the +coefficient is not "not the textbook value", it is the right value for the +phase and component actually being measured — and that couples straight back +to the rotation question above: a formulation using the SH radiation pattern +wants the transverse component, not an arbitrary horizontal. + +**`R_0 = 1000 m` is a reference distance, not a unit bridge.** An earlier +draft of this document guessed it might be bridging kilometres and metres. It +is not; it is where the source spectrum is defined, and the spreading model +carries the observation from there to the site. + +**And that is why both answers about the unit of `R` are right.** There are two +distances. `R_0` in the moment expression is in **metres**, alongside `rho` in +kg/m^3 and `beta` in m/s, giving `M0` in newton-metres. The *geometric +spreading* model `S(R)` is a separate term in Eq. 1.14, and its `R` is in +**kilometres** — Table 2.1's caption says so outright, "R in all cells is +hypocentral distance in kilometres", and tabulates the exponent piecewise. +Both of the thesis's own models, read from the table: + + Holt et al. [O] — original + 0.88 +/- 0.02 1 < R <= 40 km + 2.93 +/- 0.28 40 < R <= 63 km + 0.50 +/- 0.33 63 < R <= 100 km + 1.36 +/- 0.07 100 < R <= 400 km + + Holt et al. [R] — refined, and the one the thesis recommends + 0.90 +/- 0.01 1 < R <= 43 km + 2.57 +/- 0.07 43 < R <= 76 km + 0.44 +/- 0.08 76 < R <= 136 km + 1.54 +/- 0.04 136 < R <= 400 km + +The column heading is `-alpha`, so these are decay exponents: amplitude goes +as `R**-alpha`. **`[R]` is the preferred model**, and §2.7 gives the grounds — +lower uncertainty on every slope, more events resolved (218 against 201), and +the `Mw`-`Mc` relation moving closer to `Mw`-`ML`, which it should. An earlier +draft of this section quoted `[O]` alone; if a Utah spreading table is ever +shipped as study data, `[R]` is the one to ship. + +Worth noting the first segment independently supports the exponent argument +above, and slightly more strongly in the preferred model: the inverted +near-field decay is **0.88** and **0.90**, close to the theoretical body-wave +1, and nowhere near 2. + +It also shows the spreading is not a single power law but a piecewise +empirical function, inverted per region. So `S(R)` should be a registered +model in its own right — theoretical `1/R` as the default a study can start +from, with region-specific spreading as the thing an operator supplies, in the +same shape as `WEIGHT_MODELS` and `NOISE_MODELS`. + +**Bilinear, and why the microseismic regime is not the regional one.** Every +spreading model in the thesis breaks at 40-50 km. The PNR data used throughout +this repository spans **2.3 to 22.9 km** — entirely inside the first segment of +all of them — so the regional tables offer no guidance at these distances, and +a shape fitted for 1-400 km should not be assumed to extend downward. A +bilinear form with a break inside the microseismic range is the right thing to +allow. + +It is not something a single event can determine, though, and it is worth +recording why rather than fitting one and believing it. Measured on the 28 PNR +channels, regressing `log10(Omega)` on `log10(R)` — which for one event is +exactly the spreading, since `Omega_source` is a constant: + +| model | exponents | rms residual | +|---|---|---| +| single segment | 0.73 | 0.265 log10 | +| bilinear, break searched | 0.41 then 1.89, break 13.5 km | 0.253 log10 | + +The bilinear fit buys **5%** in rms for two extra parameters, across **one +decade** of distance. That is not evidence of a break; it is a hinge finding +scatter. And the scatter is the point: 0.265 log10 units is a factor of 1.8, +which for a single event is site response and radiation pattern, neither of +which is separable from spreading when every station contributes exactly one +distance. + +That is the argument for the non-parametric inversion rather than against +bilinear. Spreading becomes separable from site only across a dataset where +each station sees many distances and each distance is sampled by many +stations. One event cannot do it, and the apparent 0.73 here should not be +read as a measurement of anything. + +So: allow bilinear and tabulated forms, **default to theoretical `1/R`**, and +resist fitting a spreading exponent to a single event — the default is +defensible precisely because it is not fitted. + +**The calibration anchor stays Magna, not the thesis's Utah catalogue.** Utah +is tempting — it is a published catalogue of >200 `Mw` with regional constants +and bilinear `Mw`-`ML` relationships attached, which is more validation +material than Magna offers. It is not the right anchor anyway, because the +spreading estimate is the weaker half of it: Utah's `S(R)` is the piecewise +parametric table above, with segment boundaries chosen and an exponent fitted +per segment, while Magna used the **non-parametric G(R) inversion** already +noted in §5.2.5 — no imposed functional form and no chosen breakpoints. + +That ordering matters more than the size of the catalogue, because `S(R)` is +the term this section is least able to check by other means. `rho`, `beta` and +the radiation pattern are all bounded by physics and by the literature; a +spreading function is bounded by the inversion that produced it. Validating +against the dataset with the better-constrained `S(R)` tests the part of the +formula that most needs testing. + +The consequence for scope is worth stating plainly: the non-parametric +inversion is not in this repository (§5.2.5), so `S(R)` for the Magna +comparison has to come in as **data** — a tabulated `G(R)` in +`studies/magna_2020.toml` — rather than being recomputed. A spreading registry +that can be fed a table, not only a functional form, is therefore a +requirement rather than a nicety. + +The calibration is still the deliverable rather than the derivation. Pin the +study's own `rho` and `c`, state the unit of every input, and assert the Mw of +a known event in a test. This section is the case for it: two people can each +be certain about "the unit of R" and both be right, because there are two +different `R`s in the same workflow, and only a test that runs the numbers +end to end distinguishes them. + +*Sourced from Holt (2019), "Addressing Uncertainty in Earthquake Magnitudes +Commonly Used in Modern Seismic Hazard Assessment", University of Liverpool — +Ch. 1 §1.4, Ch. 2 §2.2 (method and constants), §2.3 and Table 2.1 (spreading), +and §2.7 (which spreading model is preferred).* + +**Verification status, checked rather than assumed.** Chapters 1 and 2 have +since been read directly, as separate PDF exports (2.0 MB and 5.5 MB, both +well under the limit that blocked the 44 MB original). An earlier draft of +this section expected that to settle the algebraic forms. **It does not, and +the reason has changed: the equations are not in the file at all.** + +This was checked rather than inferred from a failed extraction. The pages +carrying Eq. 1.12, 1.13, 1.14, 2.6, 2.7 and 2.8 declare **only Times New +Roman fonts and no image XObject** — no math font, no vector drawing, nothing +an extractor could be failing to decode. Rendering those pages produces the +equation *numbers* against blank space. The Word-to-PDF export dropped the +embedded equation objects, and inline math went with them, which is why the +prose reads "where, is the long period spectral displacement plateau at the +source". **So the algebraic form of Eq. 1.12, 1.13, 2.7 and 2.8 remains +unverified, and no re-export of the PDF will fix it** — that closes off the +route this section previously suggested. Settling it needs the original +document, or Edwards et al. (2010), which is the published method anyway. + +What the chapters *do* settle is every constant in the table above, each now +read from the prose of the section that defines it rather than inferred, plus +Table 2.1 in full — it survives as real text and matches the exponents quoted +above exactly. The one substantive correction is the partition factor, noted +under the table. + +`Mw` itself is a near miss. Ch. 1's footnote 10 confirms the thesis quotes +Hanks and Kanamori in **SI**, "the equivalent relation in SI units of Newton +meters (N·m) where 1 N·m = 1x10^7 dyne·cm" — so `M0` in N·m is the right +input and the formula is the SI one. The numeric constant is inside the +dropped equation, so **the 9.1 in `Mw = (2/3)(log10(M0) - 9.1)` is still +unverified from this source**, even though the unit convention around it now +is. + +Two claims in this section *are* independently confirmed by that prose, and +were inferences before: + +- **The spreading model's distance is in kilometres.** Eq. 2.6's description + states it outright — "is a geometrical spreading model and is distance + (km)" — rather than it being read off the units of Table 2.1. +- **The measurement is on the horizontal component of the Sg phase**, "referred + to as SgH". That is the pairing this section infers from `Theta-lambda-Phi` + being the SH average: the phase, the component and the radiation-pattern + constant are one choice, made once, in the published method. + +Two other things that must be true before the number means anything, both +from §4.7 above: `Omega` should be the combined horizontal rather than one +component, and the phase constants must match the phase actually measured. + +So the shape is a `magnitude` module taking a `StagedFit` and a small typed +set of medium constants — density, velocity, radiation pattern, free-surface +factor, spreading exponent — with the exponent **configurable and defaulting +to 1.0**, since a study fitting an empirical spreading term is a legitimate +thing to want and should have to say so. `[model]` holds them, `studies/*.toml` +pins them, and the output carries its unit. A test asserting the Mw of the PNR +event to within a stated tolerance is what stops a factor of 10^9 from a km/s +velocity going unnoticed — the whole point of §4.7's units discussion. + #### Suggested shape - `preprocess.rotate_to_rt(st)` reading the stored back-azimuth, and a @@ -1353,6 +1592,72 @@ points. the way it now demonstrates both minimisers: showing that the choice moves the answer is what stops it being invisible. +#### The published processing recipe, and the four defaults it sources + +The equations did not survive the export, but Ch. 2 §2.2 describes the +*processing* in prose, and that survived intact. It is the most directly +useful thing the chapters give this repository, because it is a published +specification for choices §4.4 and §4.5 currently justify on internal +measurement alone. In order, for each record: + +1. Broadband stations (`HH` and `BH`) within 400 km; instrument response + deconvolved; horizontals rotated to the **transverse** component to isolate + SgH. +2. Pg and Sg velocities of 5.9 and 3.3 km/s predict the arrivals; a **20 s + window** starts at **80% of the predicted Pg-Sg interval**, the 20% margin + being deliberate slack for pick uncertainty. +3. Transform: **zero-pad to 2^N**, de-mean, **five 3-pi prolate tapers** + (Lees and Park, 1995), FFT, **multiply by the sampling period**. +4. Noise window: record start to **75% of the interval between start and the + predicted Pg**; same transform; normalised by the ratio of signal and noise + window lengths, **excluding the zero-padding**. +5. Fit only where spectral **SNR >= 3**, by Powell minimisation, with a common + source corner frequency taken from the station of smallest misfit. +6. Event `Mw` is the mean of station estimates after rejecting beyond + **2.5 sigma**, and requires **at least 3 stations**. + +Four of those land on defaults this repository already has, which is worth +recording because in each case the justification here was internal: + +- **`MultitaperEstimator` defaults to `time_bandwidth=3.0, n_tapers=5`** — + exactly "five 3-pi prolate tapers", `NW = 3` giving `2NW - 1 = 5`. §4.4 + could previously only say the 3 was "the literal 3 passed positionally to + `mtspec`, with no way to configure it". It now has a published source, and + the pair is a considered choice rather than an inherited constant. +- **"Multiply by the sampling period" is `one_sided_fas`'s `|spec| * dt`**, + arrived at there from the requirement that padding change only the frequency + sampling. Same normalisation, reached from the other end. +- **The zero-padding exclusion agrees with that too.** The recipe is careful + that the signal/noise length ratio uses unpadded lengths — padding must not + move an amplitude. That is the same rule as `one_sided_fas`'s refusal to + rescale for padding, and `SnrConfig.scale_parseval`'s + `sqrt(len(signal)/len(noise))` is its amplitude-domain form. Worth checking + when the noise path is next touched that the ratio there is taken on the + cut lengths and not on `n_fft`. +- **SNR >= 3 matches `SnrConfig.tolerance = 3.0`** — though *applied + differently*, and the difference is exactly the distinction the config + already draws. The thesis gates the **fit band** per spectrum: fit where the + ratio holds. The Magna paper gates **spectrum selection**: keep spectra + above 3 in three fixed bands, which is `assert_bandwidths` with `bands`. + Same threshold, two different jobs, and the code can express both. + +One divergence, and it is a knowing one. The recipe pads to `2^N` where +`resolve_n_fft` defaults to `"fast"`, the next 5-smooth length. That is a +speed choice measured in §4.4 — `pow2` overshoots a 65537-sample record to +131072 where 65610 will do — and it is safe precisely because of the point +above: with normalisation keyed to `dt`, padding changes the frequency grid +and nothing else, so `"fast"` and `"pow2"` differ in sampling, not amplitude. +**`n_fft="pow2"` remains available and is what reproduces the published +workflow exactly**, which is the setting a study file should pin when +reproduction is the goal. + +The remaining steps — the 20 s window at 80% of the Pg-Sg interval, the noise +window at 75%, the 2.5-sigma station rejection and the 3-station minimum — +are event-level aggregation this package does not do yet, and they belong with +the `magnitude` module above rather than with the estimators. They are +recorded here so the module has a specification to build to instead of +inventing thresholds. + ### 4.8 Configuration: semantic groups, layered overrides, recorded provenance Scientific parameters are currently scattered across three places with no @@ -1875,13 +2180,31 @@ channel. Every exclusion is recorded with a reason naming the level it matched at, because "excluded" is not actionable and "matched exclude='AQ04' at station" is. -That this matters is measured, not assumed. Under inverse hypocentral distance -weighting the nearest two channels carry 22.8% of the weight and the nearest -four carry 36.1%, so dropping the single nearest station moves the event -corner from 12.751 Hz to 14.774 Hz — 16%, which is 1.56x in stress drop. The -choice of weighting moves it too: 12.751 (inverse hypocentral), 11.585 -(inverse epicentral), 11.048 (uniform). Both are therefore registry choices -with the published one as the default, not constants. +That this matters is measured, not assumed. Under the configured weighting — +inverse *epicentral* distance — the nearest two channels carry **41.4%** of the +weight and the nearest four 52.5%, so dropping the single nearest station moves +the event corner from 11.585 Hz to 15.602 Hz: **35%, which is 2.44x in stress +drop**. The choice of weighting moves it too: 11.585 (epicentral), 12.751 +(hypocentral), 11.048 (uniform). + +**Which distance is itself configured, and `specmod.distance` is the registry +for it.** `[geometry] distance_measure` is read now; it existed from the start as +`[windows] distance_metric` and was read by nothing, in the wrong section — +cutting a window does not depend on how distance is measured. That is not bookkeeping at these ranges: the nearest +PNR station is **0.89 km epicentral against 2.30 km hypocentral**, a factor of +2.57, while the farthest agree to 1.00 — so the two measures disagree most +exactly where the inverse-distance weight is largest. Epicentral is the honest +default here because `rhyp` is built from the source depth and the station +*elevation*, assuming every sensor sits at the surface; where sensor depths are +unknown, as they are for this deployment, that assumption is unverifiable and +the inventory's placeholder channel `depth` of `123456.0` is the tell. + +`rrup` and `rjb` are registered and **raise**. Both need a rupture surface, +and for a point source they degenerate exactly to hypocentral and epicentral — +so a silent fallback would produce plausible numbers that are wrong for any +event large enough to justify asking for them. The error names what they would +need, which puts the requirement where whoever adds finite-fault support will +read it. **One trap, found while testing and now pinned.** `require_pass` drops a station whose stage-1 fit ended against a bound. `pass_fitting` asks whether diff --git a/src/specmod/config/__init__.py b/src/specmod/config/__init__.py index a1f7b01..51705ec 100644 --- a/src/specmod/config/__init__.py +++ b/src/specmod/config/__init__.py @@ -16,6 +16,7 @@ AcquireConfig, Config, FittingConfig, + GeometryConfig, ModelConfig, SmoothingConfig, SnrConfig, @@ -29,6 +30,7 @@ "AcquireConfig", "Config", "FittingConfig", + "GeometryConfig", "ModelConfig", "Provenance", "ResolvedConfig", diff --git a/src/specmod/config/sections.py b/src/specmod/config/sections.py index db79519..4da429c 100644 --- a/src/specmod/config/sections.py +++ b/src/specmod/config/sections.py @@ -66,7 +66,6 @@ class WindowsConfig: #: run used s=3.4; 2.9 is the shipped default and is kept as such. p_velocity: float = 5.9 s_velocity: float = 2.9 - distance_metric: Literal["repi", "rhyp"] = "repi" #: Used when an S pick is missing: s_time = p_time + emergency_ratio * (p - o). emergency_ratio: float = 1.7 @@ -239,7 +238,7 @@ class FittingConfig: #: How stations are weighted into the event value. The published choice is #: inverse hypocentral distance: the nearer station has less path, so less #: of its falloff can be attenuation. See ``specmod.staged.WEIGHT_MODELS``. - event_weighting: str = "inverse_hypocentral_distance" + event_weighting: str = "inverse_distance" #: Which channels contribute to the event value, as shell globs matched #: against the trace id and each of its SEED components — so ``"AQ07"`` @@ -271,12 +270,45 @@ class VizConfig: plot_columns: int = 3 +@dataclass(frozen=True, slots=True) +class GeometryConfig: + """Source-to-site geometry. + + Its own section because more than one stage needs it. Distance feeds the + ensemble weighting of the two-stage fit (:mod:`specmod.staged`) and the + geometric spreading a moment calculation corrects for, and a setting two + consumers each keep their own copy of is how the two come to disagree. + + It lived in ``[windows]`` until there was a second reader, which was the + wrong home even then: cutting a window does not depend on how distance is + measured. + """ + + #: Which distance, resolved through :data:`specmod.distance.DISTANCE_MEASURES`. + #: + #: ``repi`` is the default and is the honest one wherever sensor depths are + #: not known. ``rhyp`` is built from the source depth and the station + #: *elevation*, so it assumes every sensor sits at the surface — for a + #: borehole deployment that is wrong by the burial depth, and nothing in + #: the metadata says so. + #: + #: The choice is not a detail at short range: on the PNR data the nearest + #: station is 0.89 km epicentral against 2.30 km hypocentral, a factor of + #: 2.57, while the farthest agree to 1.00 — so anything weighted by inverse + #: distance is most sensitive to it exactly where it matters most. + #: + #: ``rrup`` and ``rjb`` are registered and raise: both need a rupture + #: surface, and for a point source they degenerate to ``rhyp`` and ``repi``. + distance_measure: str = "repi" + + @dataclass(frozen=True, slots=True) class Config: """The whole resolved configuration.""" acquire: AcquireConfig = field(default_factory=AcquireConfig) windows: WindowsConfig = field(default_factory=WindowsConfig) + geometry: GeometryConfig = field(default_factory=GeometryConfig) transform: TransformConfig = field(default_factory=TransformConfig) smoothing: SmoothingConfig = field(default_factory=SmoothingConfig) snr: SnrConfig = field(default_factory=SnrConfig) diff --git a/src/specmod/distance.py b/src/specmod/distance.py new file mode 100644 index 0000000..99afd2c --- /dev/null +++ b/src/specmod/distance.py @@ -0,0 +1,190 @@ +"""Source-to-site distance, as a registry rather than a stat name. + +Which distance you mean is a modelling choice, and at short range it is not a +small one. On the PNR data the nearest station is **0.89 km epicentral against +2.30 km hypocentral** — a factor of 2.57 — while the farthest agree to 1.00. +Anything weighted by inverse distance, or corrected for geometric spreading, +therefore depends on the choice most strongly at exactly the station that +matters most. + +Two are implemented here because they are the two a point source supports. +Both read a value :func:`specmod.preprocess.set_stream_distance` has already +computed: + +``repi`` + Epicentral. Horizontal distance from the epicentre. +``rhyp`` + Hypocentral. Slant distance from the hypocentre. + +**Epicentral is the honest choice when sensor depths are unknown**, and that +is more often than it sounds. ``rhyp`` is built from the source depth and the +station *elevation*, which silently assumes every sensor sits at the surface. +For a borehole deployment that is wrong by the burial depth, and nothing in +the metadata announces it — the PNR inventory records channel ``depth`` as +``123456.0``, a placeholder, so on that dataset ``rhyp`` is an assumption +wearing a measurement's name. + +Finite-fault measures +--------------------- +``Rrup`` (closest distance to the rupture surface) and ``Rjb`` (Joyner-Boore, +closest horizontal distance to the surface projection of the rupture) are the +measures ground-motion work generally wants, and they are **not implemented** +— deliberately, rather than by omission. + +Both need a rupture *surface*: strike, dip, length, width and a hypocentre +position on it. SpecMod carries a point source, so there is nothing to compute +them from, and a version that quietly degenerated to ``rhyp`` and ``repi`` +would be worse than an error — those are exactly what `Rrup` and `Rjb` reduce +to for a point source, so the substitution would be invisible in the output +and wrong for any event large enough to warrant asking. + +They are registered all the same, raising with what they would need. A name +that resolves to a clear failure is a better extension point than a name that +does not resolve at all, and it puts the requirement where someone adding +finite-fault support will read it. + +The registry is the same shape as :data:`specmod.transforms.ESTIMATORS`, +:data:`specmod.core.noise.NOISE_MODELS` and +:data:`specmod.staged.WEIGHT_MODELS`, so a study names a distance the way it +names anything else and the choice travels with the resolved configuration. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any, Protocol, runtime_checkable + +import numpy as np + +from .config import load_config + +if TYPE_CHECKING: # pragma: no cover + from collections.abc import Sequence + + from numpy.typing import NDArray + +__all__ = [ + "DISTANCE_MEASURES", + "DistanceMeasure", + "Epicentral", + "FiniteFaultDistance", + "Hypocentral", + "get_distance_measure", + "resolve_distance_measure", +] + + +@runtime_checkable +class DistanceMeasure(Protocol): + """One distance per channel, in kilometres.""" + + name: str + + def distances( + self, spectra: Any, ids: Sequence[str] + ) -> NDArray[np.float64]: ... # pragma: no cover + + +@dataclass(frozen=True, slots=True) +class _FromMeta: + """A distance already computed onto the trace metadata.""" + + key: str + name: str + + def distances(self, spectra: Any, ids: Sequence[str]) -> NDArray[np.float64]: + out = np.empty(len(ids), dtype=np.float64) + for i, id in enumerate(ids): + meta = spectra[id].signal.meta + if self.key not in meta: + raise ValueError( + f"{id} carries no {self.key!r}, so its {self.name} distance " + f"is unknown. Set the geometry with " + f"specmod.preprocess.set_stream_distance." + ) + value = float(meta[self.key]) + if value <= 0: + raise ValueError( + f"{id} has {self.key}={value}, which is not a distance" + ) + out[i] = value + return out + + +@dataclass(frozen=True, slots=True) +class Epicentral(_FromMeta): + key: str = "repi" + name: str = "epicentral" + + +@dataclass(frozen=True, slots=True) +class Hypocentral(_FromMeta): + key: str = "rhyp" + name: str = "hypocentral" + + +@dataclass(frozen=True, slots=True) +class FiniteFaultDistance: + """``Rrup`` and ``Rjb``: registered, and not implemented. + + Raising here rather than omitting the name is the point. For a point source + these degenerate exactly to hypocentral and epicentral, so an + implementation that silently fell back would produce plausible numbers that + are wrong for any event big enough to justify asking for them. + """ + + name: str + needs: str + + def distances(self, spectra: Any, ids: Sequence[str]) -> NDArray[np.float64]: + raise NotImplementedError( + f"{self.name} is not implemented. It needs {self.needs}, and " + f"SpecMod carries a point source — there is no rupture surface to " + f"measure from. For a point source {self.name} degenerates to " + f"{'hypocentral' if self.name == 'rrup' else 'epicentral'}; name " + f"that instead if it is what you mean, rather than getting it by " + f"accident." + ) + + +#: Registered distance measures, resolved by name from configuration. +DISTANCE_MEASURES: dict[str, Any] = { + "repi": Epicentral, + "rhyp": Hypocentral, + "rrup": lambda: FiniteFaultDistance( + name="rrup", needs="a rupture surface — strike, dip, length and width" + ), + "rjb": lambda: FiniteFaultDistance( + name="rjb", + needs="the surface projection of a rupture — strike, dip, length and width", + ), +} + + +def get_distance_measure(name: str) -> DistanceMeasure: + """Resolve a registered measure by name.""" + try: + factory = DISTANCE_MEASURES[name] + except KeyError: + raise ValueError( + f"Unknown distance measure {name!r}. " + f"Available: {sorted(DISTANCE_MEASURES)}." + ) from None + measure: DistanceMeasure = factory() + return measure + + +def resolve_distance_measure( + measure: str | DistanceMeasure | None = None, +) -> DistanceMeasure: + """A measure from a name, an instance, or the configuration. + + ``None`` takes ``[geometry] distance_measure``, which is the project-wide + choice. It lived in ``[windows]`` and had no reader at all until this + module; cutting a window does not depend on how distance is measured. + """ + if measure is None: + measure = str(load_config().config.geometry.distance_measure) + if isinstance(measure, str): + return get_distance_measure(measure) + return measure diff --git a/src/specmod/staged.py b/src/specmod/staged.py index a25a7c7..1de6421 100644 --- a/src/specmod/staged.py +++ b/src/specmod/staged.py @@ -92,6 +92,7 @@ import numpy as np from .config import load_config +from .distance import DistanceMeasure, resolve_distance_measure from .fitting import FitSpectra if TYPE_CHECKING: # pragma: no cover @@ -134,35 +135,24 @@ class InverseDistance: modelling choice rather than a derivation, which is why this is a registry and not a hardcoded expression. - ``metric`` names a trace stat — ``rhyp`` for hypocentral, ``repi`` for - epicentral. A station missing it is a hard error rather than a silent - weight of zero: a geometry that was never set is a broken run, not a - station that should quietly stop contributing. + **Which distance is itself a choice**, and at short range not a small one: + see :mod:`specmod.distance`. ``measure=None`` takes the project-wide + setting, so a study that has decided on epicentral does not have to say so + again here. """ - metric: str = "rhyp" + #: ``None`` means "whatever the configuration says". Distance is needed by + #: geometric spreading as well as by weighting, so the choice belongs in + #: one place rather than being restated per consumer. + measure: str | DistanceMeasure | None = None name: str = "inverse_distance" def weights( self, table: Any, spectra: Any, ids: Sequence[str] ) -> NDArray[np.float64]: - out = np.empty(len(ids), dtype=np.float64) - for i, id in enumerate(ids): - meta = spectra[id].signal.meta - if self.metric not in meta: - raise ValueError( - f"{id} has no {self.metric!r}, so it cannot be weighted by " - f"distance. Set the geometry with " - f"specmod.preprocess.set_stream_distance, or choose a " - f"weighting that does not need it." - ) - distance = float(meta[self.metric]) - if distance <= 0: - raise ValueError( - f"{id} has {self.metric}={distance}, which is not a distance" - ) - out[i] = 1.0 / distance - return out + distances = resolve_distance_measure(self.measure).distances(spectra, ids) + weights: NDArray[np.float64] = 1.0 / distances + return weights @dataclass(frozen=True, slots=True) @@ -216,8 +206,13 @@ def weights( #: Registered weightings, resolved by name from ``[fitting] event_weighting``. WEIGHT_MODELS: dict[str, Any] = { - "inverse_hypocentral_distance": lambda: InverseDistance(metric="rhyp"), - "inverse_epicentral_distance": lambda: InverseDistance(metric="repi"), + # Follows the configured distance measure, so a project-wide choice is + # honoured in one place. The shipped default. + "inverse_distance": InverseDistance, + # And explicit spellings, for a study that wants to say which it used + # regardless of what the rest of the configuration says. + "inverse_hypocentral_distance": lambda: InverseDistance(measure="rhyp"), + "inverse_epicentral_distance": lambda: InverseDistance(measure="repi"), "uniform": Uniform, "inverse_variance": InverseVariance, } diff --git a/studies/magna_2020_paper.toml b/studies/magna_2020_paper.toml index 754f11b..0763da7 100644 --- a/studies/magna_2020_paper.toml +++ b/studies/magna_2020_paper.toml @@ -47,7 +47,6 @@ remove_response = false # Group velocities from Pechmann et al. (2007). p_velocity = 5.9 s_velocity = 3.4 -distance_metric = "rhyp" # 20 s window opening at 80% of the elapsed Pg-Sg time, on the transverse. s_start_ratio = 0.8 s_length = 20.0 @@ -129,3 +128,9 @@ motion = "velocity" [fitting] # "We use Powell's minimization technique (Powell, 1964; Press et al., 1997)." method = "powell" + +[geometry] +# The published run used hypocentral distance. Kept as it was: this is the +# record of what was done, not a recommendation. `repi` is the shipped default +# because it does not assume sensor depths that are often unknown. +distance_measure = "rhyp" diff --git a/tests/test_staged.py b/tests/test_staged.py index 12ba95f..c3aa673 100644 --- a/tests/test_staged.py +++ b/tests/test_staged.py @@ -20,6 +20,12 @@ obspy = pytest.importorskip("obspy") from specmod.config import load_config # noqa: E402 +from specmod.distance import ( # noqa: E402 + Epicentral, + Hypocentral, + get_distance_measure, + resolve_distance_measure, +) from specmod.fitting import FitSpectra # noqa: E402 from specmod.pipeline import spectrum_set_from_streams # noqa: E402 from specmod.staged import ( # noqa: E402 @@ -60,8 +66,10 @@ def test_it_reproduces_the_weighted_mean_computed_by_hand( stage1 = FitSpectra(spectra) stage1.fit_spectra() table = stage1.table.set_index("id") + # `repi`, because that is what `[geometry] distance_measure` says and + # the weighting now reads it. It used to hardcode `rhyp`. weights = np.array( - [1.0 / float(spectra[id].signal.meta["rhyp"]) for id in table.index] + [1.0 / float(spectra[id].signal.meta["repi"]) for id in table.index] ) by_hand = float((table["fc"].to_numpy() * weights).sum() / weights.sum()) @@ -91,7 +99,7 @@ def test_the_two_stages_are_both_kept(self, pnr_windows: Any) -> None: def test_it_reads_the_configured_parameter_and_weighting(self) -> None: fitting = load_config().config.fitting assert fitting.event_parameter == "fc" - assert fitting.event_weighting == "inverse_hypocentral_distance" + assert fitting.event_weighting == "inverse_distance" assert fitting.include == () assert fitting.exclude == () assert fitting.require_pass is True @@ -277,6 +285,7 @@ def test_describe_still_says_something_useful(self, pnr_windows: Any) -> None: class TestWeighting: def test_the_registry_resolves_and_rejects_by_name(self) -> None: assert isinstance(get_weight_model("uniform"), Uniform) + assert isinstance(get_weight_model("inverse_distance"), InverseDistance) assert isinstance( get_weight_model("inverse_hypocentral_distance"), InverseDistance ) @@ -334,6 +343,53 @@ class Pair: with pytest.raises(ValueError, match="set_stream_distance"): InverseDistance().weights(None, spectra, ["XX.A..HHZ"]) + def test_the_configured_distance_measure_is_what_gets_used( + self, pnr_windows: Any + ) -> None: + """`[geometry] distance_measure` had no reader before this. + + It matters at short range. Hypocentral and epicentral converge far from + the source and diverge near it — on these windows the nearest station + is 0.89 km epicentral against 2.30 km hypocentral — and since this + weighting is by *inverse* distance the disagreement lands hardest on + the station carrying the most weight. + """ + spectra = _spectra(pnr_windows) + assert load_config().config.geometry.distance_measure == "repi" + + configured = fit_event(spectra).value + epicentral = fit_event(spectra, weighting="inverse_epicentral_distance").value + hypocentral = fit_event(spectra, weighting="inverse_hypocentral_distance").value + + assert configured == pytest.approx(epicentral, rel=1e-12) + assert abs(hypocentral / epicentral - 1) > 0.05, ( + "the two distance measures now agree; if the geometry has changed " + "so that this no longer matters, say so here instead" + ) + + +class TestDistanceMeasures: + def test_the_registry_resolves_and_rejects(self) -> None: + assert isinstance(get_distance_measure("repi"), Epicentral) + assert isinstance(get_distance_measure("rhyp"), Hypocentral) + with pytest.raises(ValueError, match="Unknown distance measure"): + get_distance_measure("nope") + + @pytest.mark.parametrize("name", ["rrup", "rjb"]) + def test_finite_fault_measures_refuse_rather_than_degenerate( + self, name: str + ) -> None: + """For a point source these reduce exactly to hypocentral and + epicentral, so a silent fallback would give plausible numbers that are + wrong for any event large enough to justify asking for them.""" + with pytest.raises(NotImplementedError, match="rupture surface"): + get_distance_measure(name).distances(None, ["XX.A..HHZ"]) + + def test_resolve_takes_a_name_an_instance_or_the_configuration(self) -> None: + assert resolve_distance_measure("rhyp").name == "hypocentral" + assert resolve_distance_measure(Epicentral()).name == "epicentral" + assert resolve_distance_measure().name == "epicentral" # configured + class TestReporting: def test_spread_reports_what_the_mean_hides(self, pnr_windows: Any) -> None: