search

paudetseis / PyRaysum

Teleseismic body wave modeling through stacks of (dipping/anisotropic) layers
https://paudetseis.github.io/PyRaysum/
MIT License
41 stars 15 forks source link

PyRaysum does not handle evanescent waves #19

Open pasansherath opened 11 months ago

pasansherath commented 11 months ago

hi @wsja,

I have been trying to use rather thin layers (~2 km) in the velocity model for pyraysum and when I do, the following lines (#593-595) in raysum.f are triggered. The maxvals I get are <0.99 and >0.9. The resultant receiver functions (and seismograms??) are nans.

          write (*,*) 'WARNING in evec_check, maxval: ',maxval,
     &                ' (should be 1). Ignoring phase.'

Is there any workaround for this? Or is this a limitation?

wasjabloch commented 10 months ago

The warning is not specifically linked to the thickness of the layer, but rather indicates that some energy has been lost at an interfaces. This is most often the case when evanescent waves emerge (i.e. reflections at the bottom of a layer, headwaves). Such waves cannot be handled with PyRaysum. However, because their amplitude is usually small, I found it safe to ignore this warning.

I am a bit puzzled that you find that the warning is related to thickness. I usually get it with too steep dips (in relation to ray incidence), strong anisotropy, or when introducing a low velocity layer.

If you really need to, you may get around this warning by ensuring that for all your rays, all phases, and all interfaces, the critical angle of total reflection is not exceeded.

pasansherath commented 10 months ago

Thanks very much for the explanation. The thin layer is infact a low velocity layer, and has a dip of about 20 degrees.

Do you think it is safe to set max_val to a lower value (e.g. 0.90) rather than the default 0.99?

wasjabloch commented 10 months ago

The original behavior was only to warn the user and not bail out the phases. This resulted, however, in high amplitude artifacts that spoiled my receiver function. maxval is set somewhat arbitrary, so you can try to give it another value (e.g., 0.9, don't forget to re-compile your code) and observe what's happening. If the results look satisfying to you, you can proceed.