from rsf.proj import *

Flow('spik',None,'spike n1=40 n2=40 k1=20 k2=20 d1=1 d2=1')

Flow('lag.asc',None,
     '''
     echo          1   2
      98  99 100 101 102
     198 199 200 201 202
     n=100,100
     n1=12
     data_format=ascii_int
     in=$TARGET
     ''')

filter = {
    'lap': '''   -288 0
    16 -112 -228 -112 16
    -1   16   0    16 -1''',
    'thin': '''-176 14
    24 -56 -176 -56 24
    -1  24   14  24 -1'''}

for f in filter.keys():
    Flow(f+'.asc','lag',
     '''
     echo  %s lag=$SOURCE a0=342
     n1=12 data_format=ascii_float
     in=$TARGET
     ''' % filter[f],stdin=0)
    Flow(f+'1',['spik',f],'helicon filt=${SOURCES[1]} adj=n')
    Flow(f+'2',['spik',f],'helicon filt=${SOURCES[1]} adj=y')
    Flow(f+'imp',[f+'1',f+'2'],'add ${SOURCES[1]}')

lag = range(1,13)+range(85,113)+range(185,203)
Flow('lag2.asc',None,
     '''
     echo %s n1=%d n=100,100
     data_format=ascii_int
     in=$TARGET
     ''' % (string.join(map(str,lag),' '),len(lag)))

for asc in ['lag','lag2']+filter.keys():
    Flow(asc,asc+'.asc','dd form=native')

Flow('fact flag','thin lag2','wilson lagin=${SOURCES[1]} lagout=${TARGETS[1]}')

Flow('fact2','spik fact','helicon filt=${SOURCES[1]}')
Flow('spik2','thinimp fact',
     '''
     helicon filt=${SOURCES[1]} div=y adj=n |
     helicon filt=${SOURCES[1]} div=y adj=y
     ''')

plot = 'grey pclip=100 wantaxis=n crowd=0.88 gpow=0.7 title="%s" '
winplot = 'window n1=10 n2=10 f1=15 f2=15 | ' + plot

Plot('filt','lapimp',winplot % 'laplacian')
Plot('auto','thinimp',winplot % 'autocorrelation')
Plot('fact','fact2',winplot % 'Wilson factor')
Plot('spik','spik2',winplot % 'Wilson inversion')

Result('laplac','filt auto fact spik','TwoRows')

Flow('spike1',None,'spike n1=40 n2=40 nsp=2 k1=11,16 k2=8,3   mag=1,-1')
Flow('spike2',None,'spike n1=40 n2=40 nsp=2 k1=31,28 k2=24,16 mag=1,-1')
Flow('inp','spike1 fact spike2',
     'helicon filt=${SOURCES[1]} | add ${SOURCES[2]}')
Flow('div','inp fact', 'helicon filt=${SOURCES[1]} div=y')
Flow('div2','div fact','helicon filt=${SOURCES[1]} div=y adj=y')

Plot('inp', plot % 'input')
Plot('div', plot % 'input/filter')
Plot('div2',plot % "(input/filter)/filter'")

Result('thin42','inp div div2','SideBySideAniso',vppen='txscale=2.4')

End()