The module can produce points from a point to another:

D = project(origin=(-50,10,0), end_pt=(-10,30,-100), step=10, km=true)

But I wish it could:

D = project(origin=(-50,10,0), end_pt=(-10,30,-100 [this depth in km]), step=10, km=true).

The problem is that I am trying to find a good way to interpolate a Ray inside the Earth, but normal interpolators do not understand that the Earth has limits and it is not a Cartesian system HAHAHA.

But rays are curved because velocity increases with depth. Interpolating along those rays is a kind of creating flow lines when velocity changes. Otherwise it’s just a linear interpolation between 2 3D points.

I have a few problems in my hands. The project output is great as an initial ray for ray tracing. But as the ray traces I need to interpolate the nodes I have in the ray (say 17 nodes and interpolate 10 times along the path of the ray). This is complicated because the interpolation libraries do not understand the “roundness of the Earth”.

By this I mean, if I interpolate from

-67,10 to 10,20

It works

But from

-175, 0 to 175, 0 (via the minor arc) … I don’t know how to get this done well.

I’m still not sure that I fully understand the problem. Earth is not “round” on the vertical direction but due to velocity increase the rays are curved and any interpolation along the rays paths must take that into account.
You can try also the GMT forum, there are seismologist there as well and may someone will have an idea.

Interpolating along a small circle of a sphere should be simple. Just use normal range with the geographical coordinates and then compute the arc length for each interval knowing the radius of the circle at that latitude. For ellipsoidal computations I think the -G option of mapproject can be useful too. But again, this is not going on depth too.