Add automatic first motion polarity inversion

[comoglu]

Summary

• Add a grid-search based focal mechanism inversion library (firstmotion.h/cpp) that determines fault plane solutions from P-wave first motion polarities
• Integrate into scolv with an "Auto" button on the FM plot that triggers inversion using takeoff angles obtained directly from arrival data or computed via the configured travel time table (libtau, LOCSAT, homogeneous)
• Include 16 unit tests with hand-computed reference values

How it works

1. Collects P-wave first motion polarities from arrivals with associated azimuths and takeoff angles
2. Grid searches strike (0-360°), dip (0-90°), rake (-180° to 180°) at configurable spacing (default 5°)
3. Predicts polarity at each station using the P-wave radiation pattern: amplitude = (n·r)(d·r)
4. Selects the mechanism with fewest misfitting stations (configurable max misfit fraction, default 20%)
5. Computes both nodal planes via nd2dc() and displays the result

Test plan

• 16 unit tests pass (polarity prediction, azimuthal gap, full inversion)
• Tested on 3 real events in scolv with correct focal mechanism results

[gempa-lukas]

Hey Mustafa,

I tested both PRs.

When I press the 'autofm' button, the result is coming almost immediately.

I tested several events and I think the results are not correct:

The debug output says just one polarity is wrong, clearly more are wrong with the proposed FM:

12:15:56 [debug] FM inversion: 17 observations, grid=5.0 deg, maxMisfitFrac=0.20, azGap=81.9
12:15:56 [info] FM inversion: NP1=40.0/25.0/-160.0 NP2=291.7/81.7/-66.3 misfit=1/17 (0.06) gap=81.9
12:15:56 [info] Auto FM inversion: NP1=40.0/25.0/-160.0 NP2=291.7/81.7/-66.3 misfit=0.06 (1/17) gap=81.9
12:15:56 [info] Auto FM: 1 misfitting arrivals: 41

12:21:43 [debug] FM inversion: 19 observations, grid=5.0 deg, maxMisfitFrac=0.20, azGap=68.1
12:21:44 [info] FM inversion: NP1=10.0/50.0/-5.0 NP2=103.2/86.2/-139.9 misfit=0/19 (0.00) gap=68.1
12:21:44 [info] Auto FM inversion: NP1=10.0/50.0/-5.0 NP2=103.2/86.2/-139.9 misfit=0.00 (0/19) gap=68.1

[comoglu]

Hey Lukas,
Thanks for reviewing it. Do you mind sharing your full event XMLs? that I can also debug things on my side.

[gempa-lukas]

Attached are the xml files of both events above.

gempa2026clqhm.xml
gempa2026cnlvo.xml

[comoglu]

ga2026cxmhxx.xml
can you please this event? I made some updates/bugfixing.

Looks like you use PAIC which allows you to pick polarities automatically. Which is cool , on ht eother hand small number of polarities might not be perfect.

[gempa-lukas]

Currently, I'm finalizing a complete overhaul of the P-AIC plugin, with more complex and 'better' algorithms for automatic polarity determination.

Using first your test event. In an scolv offline way, I could reproduce the initial FM, but I believe this still no good solution. Also the text in the new pop-up window says '3/39' as misfit, but clearly if you count the misfitting polarities in the beachball it are more.

After dispatching your event in my database, somehow a polarity got lost, but the resulting FM was looking way better. Still the same issue as above, the shown station misfit number is not the same as counting optically.

However I tested again on more events. I attach the xml files for each screenshot for you to test. I think those should be benchmarks in terms of 'easy to compute FMs'. In some cases it looks again as the FM just needs to be mirrored along a specific axis to be valid one in my opinion.

gempa2026bntxz.xml

This one looks (almost) correct, I wonder why the algorithm did not include the most upper one - probably due to the 5° grid limitation?

gempa2026dhypo.xml

[comoglu]

So, there is progress. I have a few additional guesses to improve it. Again, thank you for testing.
It would be great to have access to P-AIC for making these tests easier.