up [pdf] 
from rsf.proj import *
from rsf.prog import RSFROOT

def Grey(data,other):
        Result(data,'grey label2=Offset unit2="km" label1=Time unit1="s" title="" labelsz=10 labelfat=4 font=2 titlesz=10 titlefat=4 screenht=10.24 screenratio=1.3 wherexlabel=t wheretitle=b color=g bartype=v clip=10113000 %s'%other)

def Graph(data,other):
        Result(data,'graph label1="Frequency" label2="Amplitude" unit2= unit1="Hz" labelsz=10 labelfat=4 font=2 titlesz=10 titlefat=4 title="" wherexlabel=b wheretitle=t %s' %other)

# Download data 
Fetch('wz.25.H','wz')

# Convert and window
Flow('data','wz.25.H',
     '''
     dd form=native | window min2=-2 max2=2 | 
     put label1=Time label2=Offset unit1=s unit2=km
     ''')
Flow('field','data','pow pow1=2 | cut n2=2 f2=20')

gamma = 2
Flow('med1','data','window n1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)
Flow('med2','data','window f1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)

Flow('dg','med1 med2','math m2=${SOURCES[1]} output="log(input/m2)/log(3)" ')

gamma = 2.3316
Flow('mmed1','data','window n1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)
Flow('mmed2','data','window f1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)

Flow('dg2','mmed1 mmed2','math m2=${SOURCES[1]} output="log(input/m2)/log(3)" ')

gamma = 2.40025

Flow('mmmed1','data','window n1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)
Flow('mmmed2','data','window f1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)

Flow('dg3','mmmed1 mmmed2','math m2=${SOURCES[1]} output="log(input/m2)/log(3)" ')

gamma = 2.41203

Flow('mmmmed1','data','window n1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)
Flow('mmmmed2','data','window f1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)

Flow('dg4','mmmmed1 mmmmed2','math m2=${SOURCES[1]} output="log(input/m2)/log(3)" ')

gamma = 2.409273

Flow('mmmmmed1','data','window n1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)
Flow('mmmmmed2','data','window f1=1000 | math output="x1^%g*abs(input)" | median | median' % gamma)

Flow('dg5','mmmmmed1 mmmmmed2','math m2=${SOURCES[1]} output="log(input/m2)/log(3)" ')

Flow('field-1','field','bandpass flo=25')
Flow('field-2','field','bandpass flo=10')
Plot('field','grey title=raw cloi=2.40113e+07')
Plot('field-1','grey title="Bandpass" cloi=2.40113e+07')
Plot('field-2','grey title="Bandpass" cloi=2.40113e+07')

Flow('dif-1','field field-1','add scale=1,-1 ${SOURCES[1]} ')
Flow('dif-2','field field-2','add scale=1,-1 ${SOURCES[1]} ') #Ground roll 1

Plot('dif-1','grey title="Dif 1" cloi=2.40113e+07')
Plot('dif-2','grey title="Dif 1" cloi=2.40113e+07')



Flow('field-1-00','field','bandpass flo=23')
Flow('dif-1-00','field field-1-00','add scale=1,-1 ${SOURCES[1]} ')


Flow('field-11','dif-1-00 field','mutter x0=0 v0=3 | add scale=-1,1 ${SOURCES[1]}')  # Target signal
Plot('field-11','grey title="Target" cloi=2.40113e+07')
Flow('dif-11','field field-11','add scale=1,-1 ${SOURCES[1]}')  # Target signal
Plot('dif-11','grey title="Dif 11" cloi=2.40113e+07')

Flow('field-11-u','field-11','window n1=750')
Flow('field-11-d','field-11','window f1=750')
Flow('dif-11-u','dif-11','window n1=750')
Flow('dif-11-d','dif-11','window f1=750')

Flow('dif-111-u field-111-u','dif-11-u field-11-u','ortho rect1=10 rect2=10 sig=${SOURCES[1]} sig2=${TARGETS[1]}')
Flow('dif-111-d field-111-d','dif-11-d field-11-d','ortho rect1=3 rect2=3 sig=${SOURCES[1]} sig2=${TARGETS[1]}')
Flow('field-ortho','field-111-u field-111-d','cat axis=1 ${SOURCES[1]}')
Flow('dif-ortho','dif-111-u dif-111-d','cat axis=1 ${SOURCES[1]}')

Plot('field-ortho','grey title="Ortho" cloi=2.40113e+07')
Plot('dif-ortho','grey title="Dif 4" cloi=2.40113e+07')

Flow('simi1','field-1 dif-1','similarity other=${SOURCES[1]} rect1=10 rect2=10')
Flow('simi4','field-ortho dif-ortho','similarity other=${SOURCES[1]} rect1=10 rect2=10')
Plot('simi1','grey color=j title="Simi1" scalebar=y')
Plot('simi4','grey color=j title="Simi2" scalebar=y')
Result('compsimi1','simi1 simi4','SideBySideAniso')

Result('comp1','field field-1 dif-1 field-ortho dif-ortho','SideBySideAniso')
Result('comp11','field field-1 dif-1 field-11 dif-11','SideBySideAniso')
Result('comp2','field field-2 dif-2','SideBySideAniso')






#######################################################################
# Adaptive matching filtering
#######################################################################
# Matching filter program
match = Program('match.c')[0]
nf  = 5 # filter length
# Dot product test 
Flow('filt0',None,'spike n1=%d' % nf)
Flow('dot.test','%s field dif-2 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','field %s dif-2 filt0' % match,
     '''
     conjgrad ./${SOURCES[1]} nf=%d niter=%d
     other=${SOURCES[2]} mod=${SOURCES[3]} 
     ''' % (nf,100))

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

Grey('field','title="Raw data"')
Grey('field-1','title="fl=25 Hz"')
Grey('field-2','title="fl=10 Hz"')
Grey('field-3','title="Adaptive subtraction"')
Grey('dif-1','title="fl=25 Hz"')
Grey('dif-2','title="fl=10 Hz"')
Grey('dif-3','title="Adaptive subtraction"')
Grey('field-ortho','title="Orthogonalized"')
Grey('dif-ortho','title="Orthogonalized"')

Flow('zooma-1','field-1','window f1=875 n1=500 f2=10 n2=20')
Flow('zooma-2','field-2','window f1=875 n1=500 f2=10 n2=20')
Flow('zooma-3','field-3','window f1=875 n1=500 f2=10 n2=20')
Flow('zooma-ortho','field-ortho','window f1=875 n1=500 f2=10 n2=20')

Flow('zoomb-1','dif-1','window f1=950 n1=500 f2=57 n2=20')
Flow('zoomb-2','dif-2','window f1=950 n1=500 f2=57 n2=20')
Flow('zoomb-3','dif-3','window f1=950 n1=500 f2=57 n2=20')
Flow('zoomb-ortho','dif-ortho','window f1=950 n1=500 f2=57 n2=20')

Grey('zooma-1','title="Zoomed A (fl=25 Hz)"')
Grey('zooma-2','title="Zoomed A (fl=10 Hz)"')
Grey('zooma-3','title="Zoomed A (Adaptive)"')
Grey('zooma-ortho','title="Zoomed A (Ortho)"')
Grey('zoomb-1','title="Zoomed B (fl=25 Hz)"')
Grey('zoomb-2','title="Zoomed B (fl=10 Hz)"')
Grey('zoomb-3','title="Zoomed B (Adaptive)"')
Grey('zoomb-ortho','title="Zoomed B (Ortho)"')

## Creating framebox1
x=0.5
y=-0.2
w=1.0
w1=1

Flow('frame1.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame1','frame1.asc',
        '''
        dd type=complex form=native |
        graph min1=0 max1=4 min2=-2 max2=2 pad=n plotfat=15 plotcol=4  screenht=10.24 screenratio=1.3
        wantaxis=n wanttitle=n yreverse=y 
        ''')

## Creating framebox2
x=2.9
y=-0.1
w=1.0
w1=1

Flow('frame2.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        string.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame2','frame2.asc',
        '''
        dd type=complex form=native |
        graph min1=0 max1=4 min2=-2 max2=2 pad=n plotfat=15 plotcol=2 screenht=10.24 screenratio=1.3
        wantaxis=n wanttitle=n yreverse=y 
        ''')

## Create label A
Plot('labela',None,
        '''
        box x0=3.2 y0=5.55 label="A" xt=0.5 yt=0.5 length=0.75 
        ''')

## Create label B
Plot('labelb',None,
        '''
        box x0=5.5 y0=5.5 label="B" xt=-0.5 yt=0.5 length=0.75
        ''')

Result('field-1-0','Fig/field-1.vpl frame1 labela','Overlay')
Result('field-2-0','Fig/field-2.vpl frame1 labela','Overlay')
Result('field-3-0','Fig/field-3.vpl frame1 labela','Overlay')
Result('field-ortho-0','Fig/field-ortho.vpl frame1 labela','Overlay')
Result('dif-1-0','Fig/dif-1.vpl frame2 labelb','Overlay')
Result('dif-2-0','Fig/dif-2.vpl frame2 labelb','Overlay')
Result('dif-3-0','Fig/dif-3.vpl frame2 labelb','Overlay')
Result('dif-ortho-0','Fig/dif-ortho.vpl frame2 labelb','Overlay')

Flow('field-f','field','spectra all=y')
Flow('field-1-f','field-1','spectra all=y')
Flow('field-2-f','field-2','spectra all=y')
Flow('field-3-f','field-3','spectra all=y')
Flow('field-ortho-f','field-ortho','spectra all=y')

Flow('field-fs','field-f field-ortho-f field-1-f field-2-f field-3-f','cat axis=2 ${SOURCES[1:5]} | window max1=50')
Graph('field-fs','plotfat=10 plotcol="7,3,5,4,6"')



End()

sfdd
sfwindow
sfput
sfpow
sfcut
sfmath
sfmedian
sfbandpass
sfgrey
sfadd
sfmutter
sfortho
sfcat
sfsimilarity
sfspike
sfdottest
sfconjgrad
sfgraph
sfbox
sfspectra

data/wz/wz.25.H