from rsf.proj import *
from rsf.recipes.beg import server as private

Fetch('cmp.HH','acig',private)
Flow('cmp','cmp.HH','dd form=native | tpow tpow=2 | mutter v0=1.4 half=n')

def grey(title,other=''):
    return '''
    grey title="%s" font=2 labelsz=12 titlesz=12 labelfat=4 titlefat=4
    scalebar=y  label1="Time (s)" label2="Offset (km)" %s
    ''' % (title,other)

Plot('cmp',grey('Data','clip=228'))

# Interpolate near offsets
Flow('cmp0','cmp','window max2=1 | pad beg2=5')
Flow('cmp1','cmp0','reverse which=2 opt=n | cat axis=2 $SOURCE')
Flow('mask0','cmp','math output=1 | window max2=1 | pad beg2=5')
Flow('mask1','mask0','reverse which=2 opt=n | cat axis=2 $SOURCE')
Flow('dip1','cmp1','math output="x2*%g/(x1+0.001)" ' % (0.05/(0.004*1.5*1.5)))
Flow('dip2','cmp1 dip1 mask1',
     'twodip2 eps=100 lam=10 dip1=${SOURCES[1]} mask=${SOURCES[2]} q0=0')
Flow('mis','cmp1 dip2 mask1',
     '''
     planemis2 dip=${SOURCES[1]} mask=${SOURCES[2]} verb=y prec=0 niter=10000
     ''')
Flow('mis2','mis cmp','window min2=0 n2=6 | cat axis=2 ${SOURCES[1]}')
Flow('mis3','mis2','window f2=1 | reverse which=2 opt=n | cat axis=2 $SOURCE')

#Plot('mis2','grey title="Near Offsets Interpolated" ')

# Predict multiples
###################
Flow('ccmp','mis2','pad n1=2048 | fft1 | fft3')
Flow('mult','ccmp','add mode=p $SOURCE | fft3 inv=y | fft1 inv=y | window n1=1000')
Plot('mult','window f2=6 | ' + grey('SRME-predicted Multiples'))

# Mask the important part
#########################
Flow('mask','mult','math output=1 | mutter hyper=y t0=0.7 v0=2 half=n | smooth rect1=5 rect2=5')
Flow('cmp2','mask mis2','add mode=p ${SOURCES[1]}')
Flow('mult2','mask mult','add mode=p ${SOURCES[1]}')

# Estimate dips
###############
Flow('mdip','mult2','dip rect1=20 rect2=10 liter=40 pmin=0')
Flow('vdip','cmp2',
     'math output="%g*x2/(x1+0.004)" ' % (0.05/(2.5*2.5*0.004)))
Flow('mask2','mdip','mask max=5 | dd type=float | smooth rect1=5 rect2=5')

Flow('mdip2','cmp2 mdip mask2 vdip',
     '''
     twodip2 eps=30  dip2=${SOURCES[1]} dip1=${SOURCES[3]} mask=${SOURCES[2]} verb=y 
     ''')

Result('srme','cmp mult','SideBySideAniso')


Flow('vcmp','cmp','vscan semblance=y v0=1 nv=100 dv=0.02 half=n')
Plot('vcmp',grey(' Data','color=j allpos=y scalebar=y maxval=0.75'))


#============== Separate signal/noise [method 2: seisigk] ======================
Flow('cmp2s','cmp2','window n2=64')
Flow('mdip2s','mdip2','window n2=64')
Flow('masks','mask','window n2=64')
Flow('comps','cmp2s mdip2s','seisigk dips=${SOURCES[1]} verb=y niter=1000')

Flow('noissIRLS','comps masks','window f3=1 | add mode=p ${SOURCES[1]}')
Flow('signsIRLS','mis2 noissIRLS','window n2=64 | add scale=1,-1 ${SOURCES[1]}')
Flow('vsignsIRLS','signsIRLS','vscan semblance=y v0=1 nv=100 dv=0.02 half=n')
Flow('vnoissIRLS','noissIRLS','vscan semblance=y v0=1 nv=100 dv=0.02 half=n')

Plot('noissIRLS','window f2=6 | ' + grey('Seislet-IRLS (noise)','clip=228'))
Plot('signsIRLS','window f2=6 | ' + grey('Seislet-IRLS (signal)','clip=228'))
Plot('vsignsIRLS',grey('Seislet-IRLS (signal)','color=j allpos=y maxval=0.75'))
Plot('vnoissIRLS',grey('Seislet-IRLS (noise)','color=j allpos=y maxval=0.75'))

#============== Separate signal/noise [method 2: mcaseislet] ======================
# obtain similar separation effect with only 10-20 iterations, more efficient!
Flow('compsmca','cmp2s mdip2s','mcaseislet dips=${SOURCES[1]} verb=y order=1 type=l eps=0.01 niter=15 pclip=10 pscale=13')

Flow('noissmca','compsmca masks','window f3=1 | add mode=p ${SOURCES[1]}')
Flow('signsmca','mis2 noissmca','window n2=64 | add scale=1,-1 ${SOURCES[1]}')
Flow('vsignsmca','signsmca','vscan semblance=y v0=1 nv=100 dv=0.02 half=n')
Flow('vnoissmca','noissmca','vscan semblance=y v0=1 nv=100 dv=0.02 half=n')

Plot('noissmca','window f2=6 | ' + grey('Seislet-MCA (noise)','clip=228'))
Plot('signsmca','window f2=6 | ' + grey('Seislet-MCA (signal)','clip=228'))
Plot('vsignsmca',grey('Seislet-MCA (signal)','color=j allpos=y maxval=0.75'))
Plot('vnoissmca',grey('Seislet-MCA (noise)','color=j allpos=y maxval=0.75'))

Plot('markA',None,'box font=2 x0=4.5 y0=1.5 label="A" xt=-0.5 yt=0.5 size=0.5')
Plot('markB',None,'box font=2 x0=7.5 y0=1.9 label="B" xt=0.5 yt=0.5 size=0.5')
Plot('markC',None,'box font=2 x0=8.2 y0=4.2 label="C" xt=0.5 yt=0.5 size=0.5')
Plot('markD',None,'box font=2 x0=7.5 y0=5 label="D" xt=0.5 yt=0.5 size=0.5')

Plot('vnoissA','vnoissIRLS markA markB markC markD','Overlay')
Plot('vnoissmcaA','vnoissmca markA markB markC markD','Overlay')
Plot('vcmpA','vcmp markA markB markC markD','Overlay')

Result('signal','signsIRLS signsmca','SideBySideAniso')
Result('nois','noissIRLS noissmca','SideBySideAniso')

Result('vsignal','vsignsIRLS vsignsmca vcmp','SideBySideAniso')
Result('vnois','vnoissA vnoissmcaA vcmpA','SideBySideAniso')


End()