RF calculation using PREM model

Hi all,
at the moment I’m trying to calculate iterative RFs by using the PREM velocity model.
For some distances I get ‘empty’ receiver functions, it is happening when doing: rfstats(stream, phase=phase, dist_range=dist_range, tt_model=tt_model). My assumption is that there is a problem with the dist_range and the phase. I set P phase with a distance range of 30-95°, Pdiff for 95-105°, PP for 105-140° and PKP for 140-165°, which worked perfectly fine for the RF calculation using the iasp91 model.
Now using the prem model, I get problems with for example distances of 134° and 139°, 134° needs to be associated with a PKP, but it is not working to associate 139° with the PKP phase. I checked the arrivals using the taupy model, where I get PKP arrivals and I also checked the ray paths using cake. I also tried the phases P, PP, PKP and Pdiff (just to try, when it is working), but the rfstream is always getting empty after using rfstats.
So, do you have any idea how to solve this problem? Did I miss anything? How can I see, which distance range rf wants me to use for which phase, if it is not conform with the taupy model?

This is how my script looks like at the moment:

tt_model='prem'
for i in range(0,len(epis)): ####epis is epicentral distances
    if 30<= epis[i] and epis[i] <= 95:
        phase='P'
        dist_range=(30,95)
    elif 105<= epis[i] and epis[i] <= 140:
        phase='PP'
        dist_range=(105,140)
    elif 140<= epis[i] and epis[i] <= 165:
        phase='PKP'
        dist_range=(140,165)
    elif 95<= epis[i] and epis[i] <= 105:
        phase='Pdiff'
        dist_range=(95,105)

    zrt=obspy.core.read('...........')
    zrt.detrend("demean")
    for s in range(0,3):
        zrt[s].stats.station_elevation = stel
        zrt[s].stats.station_latitude = stla
        zrt[s].stats.station_longitude =stlo
        zrt[s].stats.event_latitude=elats[i]
        zrt[s].stats.event_longitude =elons[i]
        zrt[s].stats.event_depth =edeps[i]
        zrt[s].stats.event_time =UTCDateTime(edates[i])
        
    stream = RFStream(zrt)
    rfstats(stream, phase=phase, dist_range=dist_range, tt_model=tt_model)

Any help is very appreciated! Many thanks! :slight_smile:

Hi Thea,

and the same script is working for tt_model='iasp91'?

The dist_range parameter automatically excludes all traces outside this distance range, but for all phases apart from P and S it has to be set by the user. You are checking the dist parameter explicitly beforehand. Then you can also set dist_range=None to not use this check. I formatted your code example with ``` ticks.

Hi Tom,
I solved the problem. It was totally my mistake…I used the script for two different stations, but forgot to actualize the epicentral distances. So, there were small differences between the epicentral distance I loaded & used for the dist_range determination and the one that was probably calculated during the RF process.

changing from iasp91 to prem, I needed to change the distance range for the PP from 110 - 146 instead of 110-145, which is understandable.

Sorry for the ‘false alarm’…
Thanks for your fast reply to my issue!

Glad that you solved it!