from rsf.proj import *

for case in Split('data mult ddata dmult'):
    segy = case + '.segy'
    Fetch(segy,'ray')
    Flow(case,segy,'segyread what=d | put label1=Time unit1=s label2=Trace')
    Plot(case,
         '''
         grey title="\F2 %s" clip=0.12 label1="\F2 Time" label2="\F2 Trace"
         parallel2=n format2="%%3.1f"
         ''' % ('Data','Multiple Model')[case[-4:]=='mult'])

nt=251
ns=7

Flow('lag.asc',None,
     'echo %d n1=1 n=%d,100 data_format=ascii_int in=$TARGET' % (nt,nt))
Flow('lag','lag.asc','dd form=native')

Flow('pef.asc','lag',
     'echo -1 a0=1 n1=1 data_format=ascii_float in=$TARGET lag=$SOURCE',
     stdin=0)
Flow('pef','pef.asc','dd form=native')

for c in ('','d'):
    mult = c+'mult'
    shifts = [mult]
    for s in range(1,ns):
        shift = '%sshift-%d' % (c,s)
        Flow(shift,mult,'window f1=%d | pad end1=%d' % (s,s))
        shifts.append(shift)

        shift = '%sshift+%d' % (c,s)
        Flow(shift,mult,'window n1=%d | pad beg1=%d' % (nt-s,s))
        shifts.append(shift)
    Flow(c+'shifts',shifts,'cat ${SOURCES[1:%d]}' % len(shifts))

    data = c+'data'
    pred = c+'pred'
    prim = c+'prim'
    filt = c+'filt'
    Flow([filt,pred],[c+'shifts',data,'pef'],
         '''
         lpf match=${SOURCES[1]} pred=${TARGETS[1]} pef=${SOURCES[2]}
         rect1=5 rect2=3 niter=400
         ''')
    Flow(prim,[data,pred],'add scale=1,-1 ${SOURCES[1]}')

    Plot(prim,'grey title="\F2 Estimated signal" clip=0.12 parallel2=n format2="%3.1f" label1="\F2 Time" label2="\F2 Trace"')
    Plot(pred,'grey title="\F2 Matched multiples" clip=0.12 parallel2=n format2="%3.1f" label1="\F2 Time" label2="\F2 Trace"')

    Flow(filt+'0',filt,'window n3=1     ')
    Plot(filt+'0','grey title="\F2 Zero coefficient" scalebar=y color=j bias=0.1 parallel2=n format2="%3.1f" label1="\F2 Time" label2="\F2 Trace" formatbar="%3.2f" ')

    Flow(filt+'s',filt,'stack axis=3')
    Plot(filt+'s','grey title="\F2 Mean coefficient" scalebar=y color=j bias=0.1 parallel2=n format2="%3.1f" label1="\F2 Time" label2="\F2 Trace" formatbar="%3.2f" ')

    Result(filt,[filt+'0',filt+'s'],'SideBySideIso')

Result('planes','data mult prim pred','TwoRows')
Result('curves','ddata dmult dprim dpred','TwoRows')

End()