from rsf.proj import *

# Critical parameters
###########################
cut = 12 # cutoff frequency
nf  = 11 # filter length
###########################

# Download data
Fetch('dune3D.H','mideast')

# Plotting macro
def grey(title):
    return '''
    window n1=490 |
    grey clip=2.5 title="%s"
    label1=Time unit1=s label2=Offset unit2=m
    ''' % title

# Select one 2-D slice
Flow('data','dune3D.H',
     '''
     dd form=native |
     window n3=1 f3=2 n1=500 f1=100 |
     scale dscale=100 
     ''')
Result('data',grey('Data'))

# Create noise model by low-pass filtering
Flow('noise0','data',
     '''
     bandpass fhi=%g |
     mutter half=n v0=1500 t0=0.8 hyper=y tp=0.12 |
     cut n1=90
     ''' % cut)
Plot('noise0',grey('Noise Model'))

# Signal = Data - Noise
Flow('signal0','data noise0','add scale=1,-1 ${SOURCES[1]}')
Plot('signal0',grey('Signal Model'))

Result('noise0','noise0 signal0','SideBySideIso')

# Plot spectrum
Plot('spec','data',
     '''
     spectra all=y | 
     graph title="Average Spectrum" max2=3 label2=
     ''')
Plot('nspec0','noise0',
     '''
     spectra all=y |
     graph wanttitle=n wantaxis=n max2=3 plotcol=5 dash=1
     ''')
Plot('sspec0','signal0',
     '''
     spectra all=y |
     graph wanttitle=n wantaxis=n max2=3 plotcol=4 dash=2
     ''')
Result('spec0','spec nspec0 sspec0','Overlay')

# Matching filter program
# program = Program('match.c')

# COMMENT ABOVE AND UNCOMMENT BELOW IF YOU WANT TO USE PYTHON 
program = Command('match.exe','match.py','cp $SOURCE $TARGET')
AddPostAction(program,Chmod(program,0o755))

match = program[0]

# Dot product test 
Flow('filt0',None,'spike n1=%d' % nf)
Flow('dot.test','%s data noise0 filt0' % match,
     '''
     dottest ./${SOURCES[0]} nf=%d
     dat=${SOURCES[1]} other=${SOURCES[2]} 
     mod=${SOURCES[3]}
     ''' % nf,stdin=0,stdout=-1)

# Conjugate-gradient optimization
Flow('filt','data %s noise0 filt0' % match,
     '''
     conjgrad ./${SOURCES[1]} nf=%d niter=%d
     other=${SOURCES[2]} mod=${SOURCES[3]} 
     ''' % (nf,2*nf))

# Extract new noise and signal
Flow('noise','filt %s noise0' % match,
     './${SOURCES[1]} other=${SOURCES[2]}')
Flow('signal','data noise','add scale=1,-1 ${SOURCES[1]}')

End()