from rsf.proj import *

Fetch('dune3D.H','mideast')

Flow('dat','dune3D.H',
     'dd form=native | window n3=1 f3=2 n1=500 f1=100 | scale dscale=100.')

Result('dune-dat','dat',
       'grey pclip=95 label1="Time (s)" label2="Offset (m)" ')

Flow('noiz','dat',
     '''
     bandpass fhi=20. |
     mutter half=n v0=1500 t0=0.8 hyper=y tp=0.12 |
     window f1=90 | pad beg1=90
     ''')
Flow('mask','noiz',
     'math output="input*input" | smooth rect1=5 rect2=5 | mask min=0.01')
Flow('mask2','mask','dd type=float')

Flow('dat2','mask2 dat','add mode=p ${SOURCES[1]}')

Flow('freq','dat2','twofreq2 verb=y')

for case in range(3):
    eps = (0.1,5,0.05)[case]
    Flow('exp%d' % case,'dat2 freq',
         'expsignoi eps=%g freq=${SOURCES[1]} niter=100 verb=y' % eps)

Flow('ex','exp2 mask2 dat dat2',
     '''
     window n3=1 | add mode=p ${SOURCES[1]} |
     add scale=1,1,-1 ${SOURCES[2:4]}
     ''')
Flow('expsn','exp2 ex',
     'window f3=1 | cat ${SOURCES[1]} | byte pclip=95 gainpanel=all')

Plot('exps','expsn',
     '''
     window f3=1 | grey label1="Time (s)" label2="Offset (m)"
     title="Estimated Signal"
     ''')
Plot('expn','expsn',
     '''
     window n3=1 | grey label1="Time (s)" label2="Offset (m)"
     title="Estimated Noise"
     ''')
Result('dune-exp','exps expn','OverUnderAniso')

Flow('sdip','exp1 mask',
     'window n3=1 | dip verb=y order=2 rect1=10 rect2=10 mask=${SOURCES[1]}')
Flow('ndip','exp2 mask',
     'window f3=1 | dip verb=y order=2 rect1=10 rect2=10 mask=${SOURCES[1]}')

Flow('sn','dat2 ndip sdip freq',
     '''
     explanesignoi niter=2000 verb=y ndip=${SOURCES[1]} sdip=${SOURCES[2]}
     freq=${SOURCES[3]} eps=0.01
     ''')
Flow('sig','sn mask2 dat dat2',
     '''
     window n3=1 | add mode=p ${SOURCES[1]} |
     add scale=1,1,-1 ${SOURCES[2:4]}
     ''')
Flow('signoi','sn sig',
     'window f3=1 | cat ${SOURCES[1]} | byte pclip=95 gainpanel=all')

Plot('sig','signoi',
     '''
     window f3=1 | grey label1="Time (s)" label2="Offset (m)"
     title="Estimated Signal"
     ''')
Plot('noi','signoi',
     '''
     window n3=1 | grey label1="Time (s)" label2="Offset (m)"
     title="Estimated Noise"
     ''')

Result('dune-sn','sig noi','OverUnderAniso')
Result('dune-sn2','sig noi','SideBySideIso')

End()