from rsf.proj import *

# get data
Fetch('cmps-tp.HH','blake')

# CMP (common midpoint) gathers
Flow('cmps','cmps-tp.HH',
     'dd form=native | reverse which=2')
# one CMP
Flow('cmp','cmps',
     '''
     window f3=950 n3=1 max1=6 |
     put o2=0.0 d2=1
     ''')
Plot('cmp',
     '''
     grey title="CMP gather" 
     unit1=s label2=Offset unit2=km
     ''')

# Near-offset section
Result('noff','cmps',
       '''
       window n2=1 n3=1024 |
       grey title="Near Offset Section" 
       unit1=s label2=Distance unit2=km
       ''')

# offset maps
Flow('off1',None,'math n1=24 o1=0.4 d1=0.1  output=x1')
Flow('off2',None,'math n1=24 o1=2.8 d1=0.05 output=x1')
Flow('off','off1 off2','cat axis=1 ${SOURCES[1]}')
Flow('offs','off','spray n=111')

# Velocity analysis
###################

vscan = '''
vscan offset=${SOURCES[1]} 
v0=1.4 nv=61 dv=0.005 half=n semblance=y
'''
pick = 'pick rect1=20 rect2=3 vel0=1.5'

# compute semblance for one CMP
Flow('vscan','cmp off',vscan)
Plot('vscan',
     '''
     grey color=j allpos=y title="Velocity Scan" 
     unit1=s label2=Velocity unit2=km/s pclip=100
     ''')

# pick maximum semblance for one CMP
Flow('pick','vscan',pick)
Plot('pick0','pick',
     '''
     graph transp=y yreverse=y min2=1.4 max2=1.7 
     plotcol=7 plotfat=10 pad=n wanttitle=n wantaxis=n
     ''')
Plot('pick1','pick',
     '''
     graph transp=y yreverse=y min2=1.4 max2=1.7 
     plotcol=0 plotfat=1 pad=n wanttitle=n wantaxis=n
     ''')
Plot('vscan2','vscan pick0 pick1','Overlay')

# compute semblance for every 10th CMP
Flow('vscans','cmps offs','window j3=10 | ' + vscan)

# pick max semblance for every 10th CMP
Flow('picks0','vscans',pick)

# interpolate picks on the original grid
Flow('picks','picks0',
     'transp | remap1 n1=1105 d1=0.05 o1=0 | transp')
Result('picks',
       '''
       grey color=j scalebar=y barreverse=y 
       allpos=y bias=1.4
       title="Stacking Velocity" 
       label1=Time unit1=s label2=Distance unit2=km
       ''')

# Normal moveout and stack
##########################

nmo = '''
nmo offset=${SOURCES[1]} 
velocity=${SOURCES[2]} half=n
'''

Flow('nmo','cmp off pick',nmo)
Plot('nmo',
     '''
     grey title="Normal Moveout" 
     unit1=s label2=Offset unit2=km
     ''')
Result('nmo','cmp vscan2 nmo','SideBySideAniso')

Flow('nmos','cmps off picks',nmo)
Flow('stack','nmos','stack')
Result('stack',
       '''
       window n2=1024 |
       grey title=Stack 
       unit1=s label2=Distance unit2=km
       ''')

# Convert from time to depth
############################

# Dix inversion
Flow('semb','vscans picks0','slice pick=${SOURCES[1]}')
Flow('vel0','picks0 semb',
     'dix rect1=20 rect2=2 weight=${SOURCES[1]}')

# interpolate on the original grid
Flow('vel','vel0',
     'transp | remap1 n1=1105 d1=0.05 o1=0 | transp')
Result('vel',
       '''
       grey color=j scalebar=y barreverse=y 
       allpos=y bias=1.4
       title="Interval Velocity" 
       label1=Time unit1=s label2=Distance unit2=km
       ''')

Flow('image','stack vel',
     '''
     time2depth velocity=${SOURCES[1]} 
     intime=y dz=0.005 nz=1001
     ''')
Result('image',
       '''
       window n2=1024 min1=3 | grey title=Image
       label1=Depth unit1=km label2=Distance unit2=km
       ''')

End()