obspy.taup arrival and rayp calculation of conversion waves like Sp

[JouSF-1409]

Avoid duplicates

• I searched existing issues

Bug Summary

Hi, I'm calculating different ray parameters for conversion waves like S40p, S100p for ccp stacking.
( I'm full aware of notes on https://docs.obspy.org/packages/obspy.taup.html about Phase naming.
It first works very on Ps wave, but goes wrong on Sp.

In summary, it returns multiple result and no one matches the result cals by TauP toolkit on java.

It's the first time I write Issues on such a fancy interface, so please forgive my poor management. I have no idea where this problem should be post to.

from obspy.taup import TauPyModel
mod = TauPyModel(model = 'iasp91')
arrivals_45 =  mod.get_travel_times(source_depth_in_km=45, distance_in_degree=55,phase_list=['S45.0p'])
arrivals_35 =  mod.get_travel_times(source_depth_in_km=45, distance_in_degree=55,phase_list=['S35.0p'])

but it returns 2 arrivals and none of which matches result from TauP toolkit in Java. So I'm quite confused about which one to believe.

>>> print(arrivals_45)
2 arrivals
        S45.0p phase arrival at 1019.377 seconds
        S45.0p phase arrival at 1019.619 seconds
>>> print(arrivals_35)
1 arrivals
        S35p phase arrival at 1020.996 seconds

taup time -deg 55 -h 45 -ph S45p
Model: iasp91
Distance   Depth   Phase   Travel    Ray Param  Takeoff  Incident  Purist    Purist
  (deg)     (km)   Name    Time (s)  p (s/deg)   (deg)    (deg)   Distance   Name
-----------------------------------------------------------------------------------
   55.00    45.0   S45p    1021.00    13.423     32.95    44.44    55.00   * S35p

Besides, if I use phase list for different layers conversion wave, the output seems to be a random order.

My question is:

1. Is there any details I didn't read from manual about conversion waves
2. which arrivals should I believe.
3. Is there any convenient method to output rayp or arrivals by a given order.

Code to Reproduce

No response

Error Traceback

No response

ObsPy Version?

1.3.1

Operating System?

Windows, Ubuntu

Python Version?

3.10.12

Installation Method?

conda

[ThomasLecocq]

I don't know how the code(s) work for conversions, but for reflections, you NEED to have a "discontinuity" in the model, see point 4. in this part of the doc: https://docs.obspy.org/packages/obspy.taup.html#phase-naming-in-obspy-taup in order to have accurate timings (otherwise the closest-by discontinuity will be used)

[JouSF-1409]

I once thought obspy.taup would act like taup java toolkit as I mentioned or as it was wrote on the doc, which means,
when I hand it 'S45p', It will acutually cal S35p(S_Moho p) since discontinuity in 35km is closest to 45km in model iasp91.

but it returns a ambiguous answer. In the future, I may interp a vel model and save it in tvel format instead of using default iasp91.

Now, I'm just curious about how it works and why.

[ThomasLecocq]

i also wonder about the fact that your source is AT the interface, what happens if you put it 1m below or above?

[JouSF-1409]

it seems if my interface is close enough, obspy.taup will take them as the same.
but 45km is clearly not close enough, since they do not share the same arrivals.

here I change source to 55km and try again.

so I check class member for purist_name, I guess this may relate to real phase cals in obspy.taup as taup toolkit do.
It seems we have a discontinuity just at 45km

[megies]

Maybe @johnrudge has some input on this

[johnrudge]

I need to look at this more closely, but it's something to do with asking for a converted phase at a depth where we don't have a discontinuity specified.

print(TauPyModel(model = 'iasp91').model.s_mod.v_mod.get_discontinuity_depths())

[   0.    20.    35.   210.   410.   660.  2889.  5153.9 6371. ]

The isap91 model has a discontinuity at 35 km, so S35p is ok. We shouldn't be able to ask for S45p really -- it looks like the Java toolkit maps that to S35p, but obspy may be mapping it instead to S210p. If I get chance I'll look in detail what it's doing. Until then I recommend only asking for a converted phase where the model has a specified discontinuity!