from rsf.proj import *

# Download data

Fetch('Nshots.su','nankai')

# Convert to RSF

Flow('shots tshots','Nshots.su',
     'suread suxdr=y tfile=${TARGETS[1]} | window min1=5')

# Extract shots

Flow('shots2 smask','shots',
     'intbin mask=${TARGETS[1]} xk=tracf yk=fldr')

Result('smask',
       '''
       dd type=float |
       stack axis=1 norm=n |
       graph symbol=x title="Nankai Trough"
       label2="Traces per Shot Gather" 
       ''')

Flow('offset','tshots',
     '''
     window n1=1 f1=11 squeeze=n | dd type=float |
     intbin head=$SOURCE xk=tracf yk=fldr
     ''')

# Display individual shots

Flow('mask1707','smask','window n2=1 min2=1707')

Flow('shot1707','shots2 mask1707',
     '''
     window n3=1 min3=1707 squeeze=n |
     headerwindow mask=${SOURCES[1]}
     ''')
Flow('offset1707','offset mask1707',
     '''
     window n3=1 min3=1707 squeeze=n |
     headerwindow mask=${SOURCES[1]}
     ''')

Result('shot1707','shot1707 offset1707',
       '''
       window min1=5.8 max1=7.8 |
       wiggle xpos=${SOURCES[1]} yreverse=y transp=y poly=y
       title="fldr 1707" label2=Offset
       ''')

Flow('fx1707','shot1707','spectra')
Result('fx1707','grey color=j allpos=y title="FX spectra"')

# Resample to 4 ms

Flow('subsamples','shots',
     '''
     bandpass fhi=125 | window j1=2
     ''')
Result('spectra','subsamples',
       'spectra all=y | graph title="Average Spectra"')

# Extract CMPs
Flow('cmps cmask','subsamples',
     '''
     intbin mask=${TARGETS[1]} xk=tracf yk=cdp
     ''')

Flow('offs','tshots',
     '''
     window n1=1 f1=11 squeeze=n | dd type=float |
     intbin xk=tracf yk=cdp head=$SOURCE
     ''')

Result('cmask',
       '''
       dd type=float |
       stack axis=1 norm=n |
       graph symbol=x title="Nankai Trough"
       label2="Traces per CMP Gather" 
       ''')

# Apply t^2 gain
Flow('cmpt','cmps','pow pow1=2')
Flow('mask','cmask','window n2=1 min2=1280')

Flow('cmp1','cmps mask',
     '''
     window n3=1 min3=1280 | headerwindow mask=${SOURCES[1]}
     ''')

Flow('cmpt1','cmpt mask',
     '''
     window n3=1 min3=1280 | headerwindow mask=${SOURCES[1]}
     ''')
Flow('off','offs mask',
     '''
     window n3=1 min3=1280 squeeze=n | headerwindow mask=${SOURCES[1]}
     ''')
Plot('cmp1','cmp1 off',
     '''
     wiggle xpos=${SOURCES[1]} title="CMP 1280"
     yreverse=y transp=y poly=y label2=Offset unit2=m
     wherexlabel=t wheretitle=b
     ''')
Plot('cmpt1','cmpt1 off',
     '''
     wiggle xpos=${SOURCES[1]} title="CMP 1280 (Gained)"
     yreverse=y transp=y poly=y label2=Offset unit2=m
     wherexlabel=t wheretitle=b
     ''')
Result('cmp','cmp1 cmpt1','SideBySideAniso')


# Velocity analysis and NMO
Flow('vscan','cmpt1 off mask',
     '''
     vscan half=n offset=${SOURCES[1]} mask=${SOURCES[2]}
     v0=1400 nv=101 dv=10 semblance=y 
     ''')
Plot('vscan','grey color=j allpos=y title="Velocity Scan" unit2=m/s')

Flow('pick','vscan',
     'mutter inner=y half=n t0=5 x0=1400 v0=75 | pick v0=1500 rect1=25')
Plot('pick',
     '''
     graph transp=y yreverse=y plotcol=7 plotfat=3
     pad=n min2=1400 max2=2400 wanttitle=n wantaxis=n
     ''')

Plot('vscanp','vscan pick','Overlay')

Flow('nmo','cmpt1 off mask pick',
     '''
     nmo half=n offset=${SOURCES[1]} mask=${SOURCES[2]}
     velocity=${SOURCES[3]}
     ''')

Plot('tpowg','cmpt1','grey title="CMP 1280" labelsz=12 titlesz=18')
Plot('nmog','nmo','grey title="Normal Moveout" labelsz=12 titlesz=18')
Plot('nmo','tpowg nmog','SideBySideAniso')

Result('vscan','cmpt1 vscanp nmo','SideBySideAniso')

# Apply to all CMPs

Flow('vscans','cmpt offs cmask',
     '''
     vscan half=n offset=${SOURCES[1]} mask=${SOURCES[2]}
     v0=1400 nv=101 dv=10 semblance=y nb=5
     ''')
Flow('picks','vscans',
     '''
     mutter inner=y half=n t0=5 x0=1400 v0=75 |
     pick v0=1500 rect1=25 rect2=10
     ''')

Result('picks',
       '''
       window | grey color=j mean=y scalebar=y title="NMO Velocity"
       label2=CMP barreverse=y barlabel=Velocity barunit=m/s
       ''')

Flow('nmos','cmpt offs cmask picks',
     '''
     nmo half=n offset=${SOURCES[1]} mask=${SOURCES[2]}
     velocity=${SOURCES[3]}
     ''')

Flow('stack','nmos','stack')

Result('stack','window min1=5.5 max1=8.0 | grey title="NMO Stack" ')

# Velocity continuation
Flow('first','stack',
     '''
     put o2=0 d2=16.667 label2=Distance unit2=m o3=0 |
     cosft sign2=1 | pad beg1=1250 |
     stolt vel=1400 
     ''')

Flow('mig0','first',
     '''
     window min1=5.5 max1=8.0 | 
     cosft sign2=-1
     ''')

Result('first','mig0',
       '''
       grey title="Migration with 1400 m/s"
       ''')

Flow('velcon','first',
     '''
     vczo nv=101 dv=10 v0=1400 verb=y |
     window min1=5.5 max1=8.0 | 
     transp plane=23 memsize=1000 |
     cosft sign2=-1
     ''')
Result('velcon',
       '''
       byte gainpanel=a | 
       grey3 frame1=300 frame2=200 frame3=50 
       point1=0.8 point2=0.8 flat=n unit3=m/s
       title="Velocity Continuation"
       ''')
Plot('movie','velcon','grey title="velocity continuation"',view=1)

# Migration with the stacking velocity
Flow('mpick','picks','window min1=5.5 max1=8.0')
Flow('mstack','velcon mpick',
     '''
     transp plane=23 memsize=1000 |
     slice pick=${SOURCES[1]}
     ''')
Result('mstack',
       'grey title="Migration with Stacking Velocity"')

# Separating diffractions

Flow('dip','mig0','dip rect1=20 rect2=10 order=2')
Result('dip',
       '''
       grey color=j scalebar=y wanttitle=n barlabel=Slope barunit=samples
       ''')

Flow('refl','mig0 dip','pwspray dip=${SOURCES[1]} ns=5 order=2 reduce=stack')
Result('refl','window min1=5.5 max1=8 | grey title="Reflections" ')

Flow('diff','mig0 refl','add scale=1,-1 ${SOURCES[1]}')
Result('diff','window min1=5.5 max1=8 | grey title="Diffractions" ')

Flow('stoltd','diff',
     '''
     cosft sign2=1 |
     pad beg1=1250
     ''')

Flow('velcond','stoltd',
     '''
     vczo nv=101 dv=10 v0=1400 verb=y pad2=3000 |
     window min1=5.5 max1=8 | 
     transp plane=23 memsize=1000 |
     cosft sign2=-1
     ''')
Plot('velcond','grey title="Velocity Continuation"',view=1)

Flow('focus','velcond',
     '''
     transp plane=23 | 
     pad beg1=1250
     focus rect1=40 rect3=20 |
     window min1=5.5 |
     math output="1/abs(input)" |
     clip clip=8 |
     scale axis=2
     ''')

Flow('fpik','focus','pick v0=1500 rect1=25 rect2=10')

Result('fpik',
       '''
       window | grey color=j allpos=y bias=1500 scalebar=y title="Focusing Velocity"
       label2=CMP barreverse=y barlabel=Velocity barunit=m/s
       ''')

Flow('mdif','velcond fpik',
     '''
     transp plane=23 memsize=1000 |
     slice pick=${SOURCES[1]}
     ''')
Result('mdif',
       'grey title="Focused Diffractions"')

Flow('mfpik','velcon fpik',
     '''
     transp plane=23 memsize=1000 |
     slice pick=${SOURCES[1]}
     ''')
Result('mfpik',
       'grey title="Migration with Focusing Velocity"')


End()