import string

from rsf.proj import *

Flow('vrms',None,'math d1=0.004 n1=1001 o1=0 output="x1*125+2000" ')
Flow('vint','vrms','math output="125*sqrt((16+x1)*(16+3*x1))" ')

for vel in ('vrms','vint'):
    Plot(vel,
         '''
         graph transp=y yreverse=y min2=1995 max2=3005
         pad=n wantaxis=n wanttitle=n plotcol=5
         ''')

Flow('synt','vrms',
     '''
     spike d1=0.004 n1=1001 o1=0 nsp=17 n2=128 d2=20 o2=0
     label2=Offset unit2=m
     k1=%s
     mag=1,1,1,1,-1,1,1,-1,1,1,1,-1,1,1,-1,1,1 |
     bandpass flo=4 fhi=20 |
     inmo velocity=$SOURCE half=n
     ''' % string.join(map(str,range (100,916,48)),','), stdin=0)
Plot('synt','grey title="(a) Input" ')

Flow('sdip0','synt',
     '''
     math output="(0.00125*(x2+10)/x1)" |
     mutter half=n v0=10000 t0=0.004 tp=0.1
     ''')
Flow('sdip','synt sdip0',
     '''
     dip order=3 liter=100 pmin=0
     idip=${SOURCES[1]} niter=10 rect1=40 rect2=5
     ''')
Plot('sdip',
     'grey title="(b) Slope" allpos=y scalebar=y clip=2 maxval=2 color=j')

Flow('dipx','sdip','transp | deriv order=3 | transp')
Flow('dipt','sdip','deriv order=3')
Flow('ssdip','sdip dipx dipt',
     'math dx=${SOURCES[1]} dt=${SOURCES[2]} output="input*dt+dx" ')
Plot('ssdip','grey title=Curvature scalebar=y allpos=y minval=0 pclip=90')

Flow('test1','sdip',
     'scale dscale=%g | math output="input*(x1-input*(x2+10))" ' % (0.004/20))
Plot('test1','grey color=j allpos=y scalebar=y clip=0.0004')

Flow('test2','ssdip',
     'scale dscale=%g | math output="input*x1*(x2+10)" ' % (0.004/(20*20)))
Plot('test2','grey color=j allpos=y scalebar=y clip=0.0004')

Flow('spnmo svel','synt sdip',
     'pnmo half=n dip=${SOURCES[1]} vel=${TARGETS[1]}')
Plot('spnmo','grey title="(c) Oriented NMO" ')

Flow('sspnmo ssvel sseta','synt sdip ssdip',
     '''
     pnmo dip=${SOURCES[1]} crv=${SOURCES[2]}
     vel=${TARGETS[1]} eta=${TARGETS[2]} half=n
     ''')
Plot('sspnmo','grey title="Nonhyperbolic Slope Moveout" ')

Result('synt','synt sdip spnmo','SideBySideAniso',vppen='txscale=1.5')

Result('ssynt','synt ssdip sspnmo','SideBySideAniso')

Flow('sdip100','sdip','window min2=150')
Flow('svtr','synt sdip100',
     '''
     window min2=150 |
     pveltran v0=2000 dv=10 nv=101 half=n dip=${SOURCES[1]} niter=20
     ''')
Plot('svtr',
     '''
     envelope |
     grey allpos=y label2=Velocity unit2=m/s
     title="(b) Oriented Velocity Map"
     ''')
Plot('svtr2','svtr vrms','Overlay')

Flow('sdipt','sdip100','deriv')
Flow('svin','synt sdip100 sdipt',
     '''
     window min2=150 |
     pveltran v0=2000 dv=10 nv=101 half=n dip=${SOURCES[1]} niter=20
     interval=y dipt=${SOURCES[2]}
     ''')
Plot('svin',
     '''
     envelope |
     grey allpos=y label2="Interval Velocity" unit2=m/s
     wanttitle=n min2=2000 max2=3000
     ''')
Result('svin','svin vint','Overlay')

Flow('svsc','synt','vscan v0=2000 dv=10 nv=101 half=n')
Plot('svsc',
     '''
     envelope |
     grey allpos=y label2=Velocity unit2=m/s
     title="(a) Velocity Scan"
     ''')
Plot('svsc2','svsc vrms','Overlay')

Result('svsc','svsc2 svtr2','SideBySideAniso',vppen='txscale=1.5')

Plot('svel',
     '''
     window min1=0.5 max1=3.5 min2=150 |
     grey allpos=y bias=2000 scalebar=y title='Oriented Vrms'
     barlabel=Velocity barunit=m/s minval=2500 maxval=2000
     ''')

Plot('svel0','svel vrms',
     '''
     window min2=150 | stack |
     cat axis=2 ${SOURCES[1]} |
     window min1=0.5 max1=3.5 |
     graph dash=1,0 title="Estimated/True Vrms" yreverse=y transp=y
     wherexlabel=t wheretitle=b label2=Velocity unit2=m/s
     ''')

Result('svel','svel svel0','SideBySideAniso')

Plot('ssvel',
     '''
     window min1=0.5 max1=3.5 min2=150 |
     grey allpos=y bias=2000 scalebar=y
     title='Oriented Nonhyperbolic Vrms'
     barlabel=Velocity barunit=m/s minval=2500 maxval=2000
     ''')

Result('ssvel','svel ssvel','SideBySideAniso')

Plot('ssvel0','ssvel vrms',
     '''
     window min2=500 max2=2000 | stack |
     cat axis=2 ${SOURCES[1]} |
     window min1=0.5 max1=3.5 |
     graph dash=1,0 title="Estimated/True Vrms" yreverse=y transp=y
     wherexlabel=t wheretitle=b label2=Velocity unit2=m/s
     ''')

Flow('sseta0','sseta','window min2=500 max2=2000 | stack')

Plot('sseta0',
     '''
     math output=0 | cat axis=2 $SOURCE |
     window min1=0.5 max1=3.5 |
     graph dash=0,1 title="Estimated/True Eta" yreverse=y transp=y
     wherexlabel=t wheretitle=b label2=" " min2=-0.2 max2=0.2
     plotcol=5,6
     ''')

Result('sspar','ssvel0 sseta0','SideBySideAniso')

End()