Calculation of takeoff angles for layered velocity model

[krisvanneste]

Claudio,

As mentioned during the WCEE, I looked into the calculation of travel times and takeoff angles when using a layered velocity model for use in sourcespec. Some years ago, I wrote a python class for doing that, based on a Fortran module in hypo71. You can find it here Except for the plot functions, this class does not depend on other custom modules, only numpy and scipy. It can be initialized as follows:

vmodel = CrustalVelocityModel(layer_top_depths, VP, VS, name='sourcespec velocity model')

I added a new calc_tt_and_angles method which calculates both travel times and takeoff angles (similar to the _wave_arrival_* functions in sourcespec), as well as incidence angles in the most efficient way, which works as follows:

result = vmodel.calc_tt_and_angles(Zf, Zs, Repi, wave, wave_type)
(travel_time, takeoff_angle, incidence_angle, wave_type) = result

with Zf = focal depth in km, Zs = station depth (in km), Repi = epicentral distance in km, wave = "P" or "S", wave_type = "REF", "DIR", "REFL" or None (= fastest arrival). It should not be too difficult for you to test this.

However, this method depends on quite a number of other methods. I don't think it would be practical to combine all of them in a single function (it would be similar to the original Fortran code, which is quite ugly...). Maybe we can copy this class in sourcespec and strip all unnecessary methods (e.g., plot methods, methods related to reflected waves, ...)?

Or is there a way to do the same calculations with layered velocity models in obspy??

[claudiodsf]

Thanks Kris!

It should not be too difficult for you to test this.

How that? Are you going to propose a PR?

Or is there a way to do the same calculations with layered velocity models in obspy??

It is possible to define a custom model in obspy.taup, but I never tryed that. Maybe you could make a comparison between your approach and ObsPy's?

[krisvanneste]

Thanks Kris!

It should not be too difficult for you to test this.

How that? Are you going to propose a PR?

I can do that if you like (but it would be better to wait until you refactor the processing submodule), but actually I meant that you can download the velocity_model.py file, import CrustalVelocityModel, instantiate it with a velocity model and run the calc_tt_and_angles method, so you can see what it looks like before incorporating it in sourcespec.

Or is there a way to do the same calculations with layered velocity models in obspy??

It is possible to define a custom model in obspy.taup, but I never tryed that. Maybe you could make a comparison between your approach and ObsPy's?

I will have another look, but I think taup makes spherical calculations rather than planar.

[claudiodsf]

I can do that if you like (but it would be better to wait until you refactor the processing submodule),

I should work on that today 😉

But actually I meant that you can download the velocity_model.py file, import CrustalVelocityModel, instantiate it with a velocity model and run the calc_tt_and_angles method, so you can see what it looks like before incorporating it in sourcespec.

Ah, ok. Sorry, I just noticed the download link 👀

[krisvanneste]

I'm looking into building a custom taup model in obspy, but it's already clear this is not the solution. I had to set the max_radius of the model to the Earth radius to avoid an error when initializing the model. Now generation of the model is running, but it's already been running for 10 minutes. Furthermore, I can see in the code the TauP model needs to be written to a file before it can be used...

[claudiodsf]

Ok, so let's try with your solution!

[krisvanneste]

I finally implemented this in PR https://github.com/SeismicSource/sourcespec/pull/78