1package locate
2
3import (
4 _ "embed"
5 "encoding/json"
6 "ewdetect/calc"
7 "ewdetect/config"
8 "fmt"
9 "image/color"
10 "math"
11 "sort"
12
13 sm "github.com/flopp/go-staticmaps"
14 "github.com/fogleman/gg"
15 "github.com/golang/geo/s2"
16 "github.com/rs/zerolog/log"
17 "gonum.org/v1/gonum/optimize"
18 "gonum.org/v1/plot/palette/moreland"
19)
20
21type ComputedTableStruct struct {
22 Angles []float64
23 Depths []float64
24 PTable [][]float64 `json:"p_table"` // Usage: PTable[depth][angle]
25 STable [][]float64 `json:"s_table"`
26}
27
28//go:embed assets/computed_tables.json
29var computedTableJson []byte
30
31var ComputedTable ComputedTableStruct
32
33func Init() {
34 err := json.Unmarshal(computedTableJson, &ComputedTable)
35 if err != nil {
36 log.Error().Err(err).Msg("locate: failed to unmarshal computed tables")
37 }
38 for i := range ComputedTable.Angles {
39 ComputedTable.Angles[i] = ComputedTable.Angles[i] / 180 * math.Pi // to rad
40 }
41 log.Info().Msg("locate: initialized computed tables")
42}
43
44func Interpolate2D(depth, angle float64, table [][]float64) float64 {
45 // Find bounding indices
46 dIdxLow, dIdxHigh, dFrac := findBoundingIndices(ComputedTable.Depths, depth)
47 aIdxLow, aIdxHigh, aFrac := findBoundingIndices(ComputedTable.Angles, angle)
48
49 // Get the four surrounding table values
50 t11 := table[dIdxLow][aIdxLow]
51 t12 := table[dIdxLow][aIdxHigh]
52 t21 := table[dIdxHigh][aIdxLow]
53 t22 := table[dIdxHigh][aIdxHigh]
54
55 // Perform bilinear interpolation
56 tInterpLow := (1-aFrac)*t11 + aFrac*t12
57 tInterpHigh := (1-aFrac)*t21 + aFrac*t22
58 tInterpolated := (1-dFrac)*tInterpLow + dFrac*tInterpHigh
59
60 /*log.Debug().
61 Float64("depth", depth).
62 Float64("angle", angle).
63 Float64("interpolated_value", tInterpolated).
64 Msg("locate: performed 2D interpolation")*/
65
66 return tInterpolated
67}
68
69// InterpolatePWave performs bilinear interpolation on the P-wave table.
70func InterpolatePWave(depth, angle float64) float64 {
71 return Interpolate2D(depth, angle, ComputedTable.PTable)
72}
73
74// InterpolateSWave performs bilinear interpolation on the S-wave table.
75func InterpolateSWave(depth, angle float64) float64 {
76 return Interpolate2D(depth, angle, ComputedTable.STable)
77}
78
79// findBoundingIndices locates the indices surrounding a value in a sorted slice using binary search.
80func findBoundingIndices(arr []float64, val float64) (lowIdx, highIdx int, frac float64) {
81 n := len(arr)
82 if val <= arr[0] {
83 return 0, 1, 0 // Just in case we are out of bounds
84 }
85 if val >= arr[n-1] {
86 return n - 2, n - 1, 1 // Just in case we are out of bounds (bis)
87 }
88
89 // Binary search: Find the smallest index where arr[i] >= val
90 idx := sort.Search(n, func(i int) bool { return arr[i] >= val })
91
92 // Ensure we get two valid bounding indices
93 lowIdx = idx - 1
94 highIdx = idx
95
96 if idx >= len(arr) {
97 lowIdx = len(arr) - 2
98 highIdx = len(arr) - 1
99 }
100
101 // Compute interpolation fraction
102 frac = (val - arr[lowIdx]) / (arr[highIdx] - arr[lowIdx])
103
104 return lowIdx, highIdx, frac
105}
106
107type Observation struct {
108 StationName string
109 Lat float64 //rad
110 Lon float64 //rad
111 PWaveArrival float64 // seconds from event start epoch
112 SWaveArrival float64 // seconds from event start epoch
113}
114
115type Guess struct {
116 Lat float64 //rad
117 Lon float64 //rad
118 Epoch float64 // seconds from event start epoch (is normally negative, unless seismographs open a wormhole and travel back in time)
119 Depth float64
120}
121
122func Square(f float64) float64 {
123 return f * f
124}
125
126func GreatCircleAngle(lat1, lon1, lat2, lon2 float64) float64 { // rad
127 return math.Acos(math.Sin(lat1)*math.Sin(lat2) + math.Cos(lat1)*math.Cos(lat2)*math.Cos(lon1-lon2))
128}
129
130func ErrorFunction(observations *[]Observation, guess Guess) float64 {
131 cumulErrorSquared := 0.0
132 for _, observation := range *observations {
133 diffAngle := GreatCircleAngle(guess.Lat, guess.Lon, observation.Lat, observation.Lon)
134 pDelay := observation.PWaveArrival - guess.Epoch
135 sDelay := observation.SWaveArrival - guess.Epoch
136 cumulErrorSquared += Square((sDelay-InterpolateSWave(guess.Depth, diffAngle))/sDelay) + Square((pDelay-InterpolatePWave(guess.Depth, diffAngle))/pDelay)
137 }
138 log.Debug().
139 Float64("error_squared", cumulErrorSquared).
140 Msg("locate: calculated error function")
141 return cumulErrorSquared // no unit
142}
143
144// http://epubs.siam.org/doi/pdf/10.1137/S1052623496303470
145// https://pkg.go.dev/gonum.org/v1/gonum/optimize#NelderMead
146func NelderMeadOptimization(eventName string, observations *[]Observation, drawDebugMap bool) Guess {
147 log.Info().
148 Str("event_name", eventName).
149 Int("observation_count", len(*observations)).
150 Msg("locate: starting Nelder-Mead optimization")
151
152 problem := optimize.Problem{
153 Func: func(x []float64) float64 {
154 return ErrorFunction(observations, Guess{
155 Lat: x[0], Lon: x[1], Epoch: x[2], Depth: x[3]})
156 },
157 }
158 sumLatX := 0.0
159 sumLatY := 0.0
160 sumLonX := 0.0
161 sumLonY := 0.0
162 for _, observation := range *observations {
163 sumLatX += math.Cos(observation.Lat)
164 sumLatY += math.Sin(observation.Lat)
165 sumLonX += math.Cos(observation.Lon)
166 sumLonY += math.Sin(observation.Lon)
167 }
168 meanLat := math.Atan2(sumLatY/float64(len(*observations)), sumLatX/float64(len(*observations)))
169 meanLon := math.Atan2(sumLonY/float64(len(*observations)), sumLonX/float64(len(*observations)))
170 result, err := optimize.Minimize(problem, []float64{meanLat, meanLon, 0, 0}, nil, &optimize.NelderMead{})
171 if err != nil {
172 log.Error().Err(err).Msg("locate: optimization failed")
173 }
174
175 if config.Debug && drawDebugMap {
176 log.Debug().Msg("locate: generating debug map")
177 minLat := result.X[0]
178 maxLat := result.X[0]
179 minLon := result.X[1]
180 maxLon := result.X[1]
181 for _, observation := range *observations {
182 if observation.Lat > maxLat {
183 maxLat = observation.Lat
184 }
185 if observation.Lon > maxLon {
186 maxLon = observation.Lon
187 }
188 if observation.Lat < minLat {
189 minLat = observation.Lat
190 }
191 if observation.Lon < minLon {
192 minLon = observation.Lon
193 }
194 }
195 minLat -= config.DebugBoundingBoxPadding
196 minLon -= config.DebugBoundingBoxPadding
197 maxLat += config.DebugBoundingBoxPadding
198 maxLon += config.DebugBoundingBoxPadding
199 var evals [][]float64
200 minEval := math.MaxFloat64
201 maxEval := 0.0
202 for i := 0; i < config.DebugFunctionAreas; i++ {
203 var evalLine []float64
204 for j := 0; j < config.DebugFunctionAreas; j++ {
205 lat := float64(i)/float64(config.DebugFunctionAreas)*minLat + (1-float64(i)/float64(config.DebugFunctionAreas))*maxLat
206 lon := float64(j)/float64(config.DebugFunctionAreas)*minLon + (1-float64(j)/float64(config.DebugFunctionAreas))*maxLon
207 eval := math.Sqrt(ErrorFunction(observations, Guess{
208 Lat: lat, Lon: lon, Epoch: result.X[2], Depth: result.X[3]}))
209 if eval > maxEval {
210 maxEval = eval
211 }
212 if eval < minEval {
213 minEval = eval
214 }
215 evalLine = append(evalLine, eval)
216
217 }
218 evals = append(evals, evalLine)
219 }
220 ctx := sm.NewContext()
221 ctx.SetSize(int(config.DebugResolution*math.Abs(calc.AngleDiff(maxLon, minLon)/calc.AngleDiff(maxLat, minLat))), config.DebugResolution)
222 ctx.SetBoundingBox(s2.RectFromLatLng(s2.LatLngFromDegrees(minLat*180/math.Pi, minLon*180/math.Pi)).AddPoint(s2.LatLngFromDegrees(maxLat*180/math.Pi, maxLon*180/math.Pi)))
223 ctx.AddObject(
224 sm.NewMarker(
225 s2.LatLngFromDegrees(result.X[0]*180/math.Pi, result.X[1]*180/math.Pi),
226 color.RGBA{0xff, 0, 0, 0xff},
227 16.0,
228 ),
229 )
230
231 for _, observation := range *observations {
232 ctx.AddObject(
233 sm.NewMarker(
234 s2.LatLngFromDegrees(observation.Lat*180/math.Pi, observation.Lon*180/math.Pi),
235 color.RGBA{0, 0xff, 0, 0xff},
236 16.0,
237 ),
238 )
239 }
240
241 colormap := moreland.ExtendedBlackBody()
242 colormap.SetMax(1)
243 colormap.SetMin(0)
244 colormap.SetAlpha(0.5)
245
246 for i := 0; i < config.DebugFunctionAreas; i++ {
247 for j := 0; j < config.DebugFunctionAreas; j++ {
248 latLow := float64(i)/float64(config.DebugFunctionAreas)*minLat + (1-float64(i)/float64(config.DebugFunctionAreas))*maxLat
249 lonLow := float64(j)/float64(config.DebugFunctionAreas)*minLon + (1-float64(j)/float64(config.DebugFunctionAreas))*maxLon
250 latHi := float64(i+1)/float64(config.DebugFunctionAreas)*minLat + (1-float64(i+1)/float64(config.DebugFunctionAreas))*maxLat
251 lonHi := float64(j+1)/float64(config.DebugFunctionAreas)*minLon + (1-float64(j+1)/float64(config.DebugFunctionAreas))*maxLon
252 remappedEval := (math.Log(evals[i][j]/maxEval+math.E) - math.Log(minEval/maxEval+math.E)) / (math.Log(1+math.E) - math.Log(minEval/maxEval+math.E))
253 currColor, err := colormap.At(remappedEval)
254 if err != nil {
255 log.Warn().Msg("locate: failed to get color from colormap")
256 currColor = color.RGBA{0, 0, 0, 0}
257 }
258 ctx.AddObject(
259 sm.NewArea(
260 []s2.LatLng{
261 s2.LatLngFromDegrees(latLow*180/math.Pi, lonLow*180/math.Pi),
262 s2.LatLngFromDegrees(latHi*180/math.Pi, lonLow*180/math.Pi),
263 s2.LatLngFromDegrees(latHi*180/math.Pi, lonHi*180/math.Pi),
264 s2.LatLngFromDegrees(latLow*180/math.Pi, lonHi*180/math.Pi),
265 },
266 color.RGBA{0, 0, 0, 0},
267 currColor,
268 1,
269 ),
270 )
271
272 }
273 }
274
275 ctx.OverrideAttribution(fmt.Sprintf("%s - EWDetect - Louis \"OnTake\" Dalibard - 2025.", ctx.Attribution()))
276
277 img, err := ctx.Render()
278 if err != nil {
279 log.Error().Err(err).Msg("locate: failed to render debug map")
280 }
281
282 if err := gg.SavePNG("debug/event-maps/"+eventName+".png", img); err != nil {
283 log.Error().Err(err).Msg("locate: failed to save debug map")
284 }
285 log.Debug().Str("event_name", eventName).Msg("locate: debug map generated")
286 }
287
288 log.Info().
289 Float64("lat", result.X[0]).
290 Float64("lon", result.X[1]).
291 Float64("epoch", result.X[2]).
292 Float64("depth", result.X[3]).
293 Msg("locate: optimization complete")
294
295 return Guess{
296 Lat: result.X[0], Lon: result.X[1], Epoch: result.X[2], Depth: result.X[3]}
297}