up [pdf]
from rsf.proj import *
from rsf.recipes.uncert import uncert

Fetch('beinew.HH','midpts')
Flow('bei','beinew.HH','dd form=native | transp plane=23 | transp plane=34')

uncert('bei',
       nv=96,       # continuation steps
       v0=1.3,      # initial velocity
       dv=0.0125,   # velocity step
       nx=250,      # lateral dimension
       nh=48,        # number of offsets
       padt=1024,   # time padding
       padt2=2048,  # extra time padding
       vslope=0.67, #
       vx0=1.5,     #
       padx=521,    # lateral padding
       n1=876,      # vertical extent
       v1=1.8,      # other velocity
       dx=0.0335,   # lateral sampling
       x0=7.705,    # lateral origin
       rect1=15,    # vertical smoothing
       rect2=15)    # lateral  smoothing

Flow('left0','bei-unc','math output=x1-0.5*input')
Flow('rite0','bei-unc','math output=x1+0.5*input')
Flow('left1','bei-unc2','math output=x2-0.5*input')
Flow('rite1','bei-unc2','math output=x2+0.5*input')
Flow('left','left0 left1','cmplx ${SOURCES[:2]}',stdin=0)
Flow('rite','rite0 rite1','cmplx ${SOURCES[:2]}',stdin=0)

Flow('arr','left rite',
     'cat axis=3 ${SOURCES[1]} | transp plane=13')
Plot('arr',
     '''
     window j2=4 j3=10 | put n3=1 n2=4788 |
     graph title="Structural Uncertainty"
     min1=0 max1=3 min2=7.705 max2=16.0465 pad=n
     transp=y yreverse=y plotcol=6
     wheretitle=b wherexlabel=t
     label1=Time unit1=s label2=Lateral unit2=km       
     ''')

Plot('agc','bei-agc',
     '''
     window max1=3 | grey wanttitle=n wantaxis=n
     ''')

Result('arr','agc arr','Overlay')

Result('vlf','bei-vlf',
       '''
       byte gainpanel=all |
       transp plane=23 |
       grey3 frame1=500 frame2=99 frame3=55 flat=n
       title="Velocity Continuation" point1=0.75 point2=0.75
       label3=Velocity unit3=km/s label2=Distance unit2=km 
       label1=Time unit1=s
       ''')

Result('npk','bei-npk','Overlay')

Flow('scan','bei-vlf1','window n3=1 min3=10 | mutter half=n v0=0.67 x0=1.5')
Plot('scan','grey color=j allpos=y title="Semblance Scan" label1=Time unit1=s label2=Velocity unit2=km/s')

Flow('pick','scan','pick rect1=15') # 'window n2=1 min2=10')

def graph(col,fat):
    return '''
    graph transp=y yreverse=y min2=1.3125 max2=2.5 pad=n plotcol=%d plotfat=%d
    wantaxis=n wanttitle=n
    ''' % (col,fat)

Plot('pick0','pick',graph(0,10))
Plot('pick1','pick',graph(7,1))

Result('scan','scan pick0 pick1','Overlay')

Flow('slice','bei-vlf','window n3=1 min3=10 | mutter half=n v0=0.67 x0=1.5')
Plot('slice','grey title="Midpoint Slice" label1=Time unit1=s label2=Velocity unit2=km/s ')

Result('slice','slice pick0 pick1','Overlay')

Flow('tslice','bei-vlf','window n1=1 min1=2 | transp')
Plot('tslice','grey title="Time Slice" label1=Distance unit1=km label2=Velocity unit2=km/s ')

Flow('tpick','bei-npk','window n1=1 min1=2') # 'window n2=1 min2=10')

Plot('tpick0','tpick',graph(0,10))
Plot('tpick1','tpick',graph(7,1))

Result('tslice','tslice tpick0 tpick1','Overlay')

Flow('ddv','bei-ddv','window n2=1 min2=10')
Flow('ppick','pick ddv','add ${SOURCES[1]} scale=1,0.25')
Flow('mpick','pick ddv','add ${SOURCES[1]} scale=1,-0.25')

Plot('ppick0','ppick',graph(0,10))
Plot('ppick1','ppick',graph(7,1))

Plot('mpick0','mpick',graph(0,10))
Plot('mpick1','mpick',graph(7,1))

Result('scan2','scan ppick0 mpick0 ppick1 mpick1','Overlay')

Result('slice2','slice ppick0 mpick0 ppick1 mpick1','Overlay')

Flow('smb','bei-vlf1','window n3=1 | mutter half=n v0=0.67 x0=1.5 | put unit2=km/s')
Plot('smb','grey title=Semblance wheretitle=t allpos=y wherexlabel=b')

eps = (10,50)
for rect in eps:
    pick = 'vpick%d' % rect
    Flow(pick,'smb','pick rect1=%d' % rect)
    Plot(pick,'graph pad=n transp=y yreverse=y min2=1.325 max2=2.5 title="epsilon=%g" ' % (0.01*rect))

Result('velpick','smb vpick%d vpick%d' % eps,
       'SideBySideAniso',vppen='txscale=1.5')

Flow('dip','bei-agc','dip rect1=40 rect2=10')

Result('dip','grey color=j scalebar=y title=Dip')

Flow('time','dip','pwpaint i0=125 eps=0.1')
Flow('flat','bei-agc time','iwarp warp=${SOURCES[1]}')

Result('flat','grey title=Flattened')

Flow('time2','time','math output=x1 | iwarp warp=$SOURCE')

Flow('cont','bei-agc',
     '''
     window n2=1 f2=125 max1=2.8 | envelope | 
     max1 | real | window n1=15 | sort | spray axis=2 n=250
     ''')

Flow('hors','time2 cont','inttest1 coord=${SOURCES[1]} interp=lag nw=2')
Flow('hdt','bei-unc  cont','inttest1 coord=${SOURCES[1]} interp=lag nw=2')
Flow('hdx','bei-unc2 cont','inttest1 coord=${SOURCES[1]} interp=lag nw=2')

Flow('xm','hors hdx','math dx=${SOURCES[1]} output="x2-0.5*dx" ')
Flow('tm','hors hdt','math dt=${SOURCES[1]} output="input-0.5*dt" ')

Flow('xp','hors hdx','math dx=${SOURCES[1]} output="x2+0.5*dx" ')
Flow('tp','hors hdt','math dt=${SOURCES[1]} output="input+0.5*dt" ')

Plot('hors',
     '''
     transp |
     graph title="Horizon Uncertainty"
     min2=0 max2=3 pad=n yreverse=y 
     wheretitle=b wherexlabel=t
     label2=Time unit2=s label1=Lateral unit1=km       
     ''')

Plot('hm','xm tm',
     '''
     cmplx ${SOURCES[1]} |
     transp |
     graph title="Horizon Uncertainty"
     min2=0 max2=3 pad=n yreverse=y 
     wheretitle=b wherexlabel=t
     label2=Time unit2=s label1=Lateral unit1=km    
     ''')

Plot('hp','xp tp',
     '''
     cmplx ${SOURCES[1]} |
     transp |
     graph wanttitle=n wantaxis=n
     min2=0 max2=3 pad=n yreverse=y 
     ''')

Result('hors','agc hm hp','Overlay')


End()

sfdd
sftransp
sfhalfint
sfpreconstkirch
sfwindow
sfpad
sfcosft
sfput
sffourvc
sfmath
sfclip2
sfmul
sfdivn
sffourvc2
sfstack
sfgrey
sfmutter
sfscale
sfpick
sfslice
sfagc
sfrefer
sfdip
sfadd
sfpwd
sfspray
sffocus
sfcmplx
sfcat
sfgraph
sfbyte
sfgrey3
sfpwpaint
sfiwarp
sfenvelope
sfmax1
sfreal
sfsort
sfinttest1

data/midpts/beinew.HH