#!/usr/bin/env python
# coding: utf-8

# # Calculate receiver functions from minimal example
# This notebook demonstrates the most simple use case of the [rf](https://github.com/trichter/rf) package.
# 
# If you want to run this notebook locally, have a look at the repository [readme](https://github.com/trichter/notebooks). Dependencies of this notebook are ObsPy and rf.
# 
# We load example data into a RFStream instance by simply calling `read_rf()`. The stream consists of 9 traces with three component data from three different earthquakes. The header values of the first trace in the stream is printed and we can observe that the SAC headers in `stats.sac` are correctly mapped to corresponding entries in `stats`. Raw waveforms of the first event are plotted with ObsPy.

# In[1]:


import matplotlib.pyplot as plt
from rf import read_rf, rfstats

stream = read_rf()
print(stream)
print('\nStats:\n', stream[0].stats)
stream[:3].plot()


# Now, we fill the stats dictionary with entries which we need for receiver function calculation with the `rfstats` function. The print statement automatically displays some important header values. Data is filtered, trimmed and the preprocessed data of the first event is again plotted with ObsPy.

# In[2]:


rfstats(stream)
stream.filter('bandpass', freqmin=0.4, freqmax=1)
stream.trim2(5, 95, 'starttime')
print(stream)
stream[:3].plot(type='relative', reftime=stream[0].stats.onset)


# Finally, we perform the receiver function calculation, i.e. rotation and deconvolution, by calling `RFStream.rf()`. The L components of the receiver functions have the expected peak at 0s. The Q component stack shows a distinct phase at around 9s. That's it.

# In[3]:


stream.rf()
stream.moveout()
stream.trim2(-5, 22, 'onset')
print(stream)
stream.select(component='L').plot_rf()    
stream.select(component='Q').plot_rf()
plt.show()