from rsf.proj import *

Fetch('seismic.segy','viking')

# Convert to RSF
Flow('seismic tseismic','seismic.segy',
     'segyread tfile=${TARGETS[1]}')

for key in ('cdp','offset'):
    Flow(key,'tseismic',
         'dd type=float | headermath output=%s' % key)

wiggle = '''
wiggle transp=y yreverse=y pclip=99 wanttitle=n
xpos=${SOURCES[1]} label2=Offset unit2=m poly=y
'''

# Capture a single CMP
Flow('cdpmask','cdp','mask min=265 max=265')
Flow('cdp265','seismic cdpmask',
     'headerwindow mask=${SOURCES[1]} | pow pow1=2')
Flow('offset265','offset cdpmask',
     'headerwindow mask=${SOURCES[1]}')
Result('cdp','cdp265 offset265',wiggle)

# NMO with water velocity
Flow('vel','cdp265','window n2=1 | math output=1500')
Flow('nmo','cdp265 offset265 vel',
     'nmo half=n offset=${SOURCES[1]} velocity=${SOURCES[2]}')
Result('nmo','nmo offset265',wiggle)

# Tau-p transform
Flow('taup','nmo offset265 ',
     'radon parab=y x0=1000 np=251 dp=0.002 p0=-0.25 offset=${SOURCES[1]}')
Result('taup','grey wanttitle=n label2=p label1="\F10 t\F3 " ')

# Remove multiples
Flow('nmo2','taup','radon parab=y x0=1000 adj=n nx=60 ox=-3237 dx=50')
Flow('sign2','taup nmo2',
     '''
     mutter v0=0.02 |
     radon parab=y x0=1000 adj=n nx=60 ox=-3237 dx=50 |
     add scale=-1,1 ${SOURCES[1]}
     ''')
Flow('sign','sign2 vel','inmo half=n velocity=${SOURCES[1]}')
Result('sign',
       '''
       wiggle transp=y yreverse=y pclip=99 wanttitle=n poly=y
       ''')

# Velocity analysis
Flow('scan','cdp265 offset265',
     'vscan semblance=y half=n offset=${SOURCES[1]} dv=15 v0=1450 nv=200')
Plot('scan','grey color=j allpos=y wanttitle=n')

Flow('scan2','sign',
     'vscan semblance=y half=n dv=15 v0=1450 nv=200')
Plot('scan2','grey color=j allpos=y wanttitle=n')

Result('scan','scan scan2','SideBySideAniso')

End()