from rsf.proj import *

#########################
# Real Data
#########################
Fetch('hector.HH','ground')

Flow('dat','hector.HH','dd form=native | pad beg2=5 | put o2=0.03 ')

Result('dat',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       grey label2="Offset (km)" label1="Time (s)" title=n title="Raw data"
       screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
       ''')

## Result('wdat','dat',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset (km)" label1="Time (s)"
##        title=n title="Raw data"
##        screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
##        transp=y yreverse=y poly=y
##        ''')

## Flow('trace','dat','window n2=1 f2=10')
## Result('trace',
##        '''
##        graph wanttitle=n screenratio=0.8 crowd1=0.75
##        label2=Amplitude label1=Time unit1=s crowd2=0.3 title=Inverse
##        ''')
## Flow('spec','trace','spectra')
## Result('spec',
##        '''
##        graph title=Spectra
##        ''')

#########################
# Ltf transform
#########################
Flow('ltft','dat',
     '''
     ltft rect=20 verb=n nw=200 dw=0.5 niter=25
     ''',split=[2,93],reduce="cat axis=3")
Result('ltft',
       '''
       window n3=88 f3=5 max2=50 | put o3=0.017 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=23 memsize=1000 | 
       grey3 color=j frame1=375 frame3=15 frame2=18 label1=Time flat=n 
       unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
       screenht=10 screenratio=0.7
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f
       ''')

Flow('thr1','ltft',
     '''
     real | transp plane=23 memsize=1000 | 
     mutter t0=-0.2 v0=1.1 | cut min3=12
     ''' )#mutter inner=y v0=0.3 |
Flow('thr2','ltft',
     '''
     imag | transp plane=23 memsize=1000 |
     mutter t0=-0.2 v0=1.1 | cut min3=12
     ''')#mutter inner=y v0=0.3 |
Flow('complx','thr1 thr2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=23 memsize=1000
     ''')
Plot('complx',
     '''
     window n3=88 f3=5 max2=50 | put o3=0.017 |
     math output="abs(input)" | real |
     byte allpos=y gainpanel=40 pclip=100 |
     grey3 color=j frame1=375 frame2=18 frame3=15 label1=Time flat=y 
     unit1=s label3=Distance label2=Frequency unit3=Hz
     point1=0.8 point2=0.2 wanttitle=n
     ''',view=1)

Flow('maskthr1','ltft',
     '''
     real | transp plane=23 memsize=1000 |
     math output=1. |
     mutter t0=-0.2 v0=1.1 | cut min3=12 |
     math output="1-input"
     ''' )#mutter inner=y v0=0.3 |
Flow('maskthr2','ltft',
     '''
     imag | transp plane=23 memsize=1000 |
     math output=1. |
     mutter t0=-0.2 v0=1.1 | cut min3=12 |
     math output="1-input"
     ''')#mutter inner=y v0=0.3 |
Flow('maskcomplx','maskthr1 maskthr2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=23 memsize=1000
     ''')
Result('mask','maskcomplx',
       '''
       window n3=88 f3=5 max2=50 | put o3=0.017 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=23 memsize=1000 | 
       grey3 frame1=375 frame3=15 frame2=18 label1=Time flat=n 
       unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
       screenht=10 screenratio=0.7 color=g
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f
       ''')

Flow('iltft','complx','ltft inv=y verb=n')
Result('iltft',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       grey label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 title="Difference"
       labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       ''')
## Result('wiltft','iltft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Difference"
##        labelfat=4 font=2 titlefat=4  clip=2.74044
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

Flow('siltft','dat iltft','add scale=1,-1 ${SOURCES[1]}')
Result('siltft',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       grey label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 title="LTF decomposition"
       labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       ''')

## Result('wsiltft','siltft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="LTF decomposition"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

# t-f-k transform
Flow('tfk','ltft','fft3 axis=3')
Result('tfk',
       '''
       window min3=0 max2=50 j1=5 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=12 memsize=1000 | transp plane=23 memsize=1000 |
       grey3 color=j frame1=15 frame3=75 frame2=48 label1="\F5 F \F-1" flat=n 
       unit3=s label3=Time label2="\F5 K \F-1" unit2=1/km unit1=Hz
       screenht=10 screenratio=0.7
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
       ''')

Flow('tfkmask1','tfk',
     '''
     real | transp plane=12 memsize=1000 | transp plane=23 memsize=1000 | 
     mutter v0=3. inner=y | cut max3=0.1
     ''' )#mutter inner=y v0=0.3 |
Flow('tfkmask2','tfk',
     '''
     imag | transp plane=12 memsize=1000 | transp plane=23 memsize=1000 |
     mutter v0=3. inner=y | cut max3=0.1
     ''')#mutter inner=y v0=0.3 |
Flow('tfkcomplx','tfkmask1 tfkmask2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=23 memsize=1000 | transp plane=12 memsize=1000
     ''')
## Result('tfkcomplx',
##        '''
##        window min3=0 max2=50 j1=5 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=12 memsize=1000 | transp plane=23 memsize=1000 |
##        grey3 color=j frame1=10 frame3=75 frame2=1 label1="\F5 F \F-1" flat=n 
##        unit3=s label3=Time label2="\F5 K \F-1" unit2=1/km unit1=Hz
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f o2num=0 d2num=15 n2tic=6 movie=2
##        ''')

Flow('masktfk1','tfk',
     '''
     real | transp plane=12 memsize=1000 | transp plane=23 memsize=1000 |
     math output=1. |
     mutter v0=3. inner=y | cut max3=0.1 |
     math output="1-input"
     ''' )#mutter inner=y v0=0.3 |
Flow('masktfk2','tfk',
     '''
     imag | transp plane=12 memsize=1000 | transp plane=23 memsize=1000 | 
     math output=1. |
     mutter v0=3. inner=y | cut max3=0.1 |
     math output="1-input"
     ''')#mutter inner=y v0=0.3 |
Flow('tfkmask','masktfk1 masktfk2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=23 memsize=1000 | transp plane=12 memsize=1000
     ''')
Result('masktfk','tfkmask',
       '''
       window min3=0 max2=50 j1=5 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=12 memsize=1000 | transp plane=23 memsize=1000 |
       grey3 frame1=15 frame3=75 frame2=48 label1="\F5 F \F-1" flat=n 
       unit3=s label3=Time label2="\F5 K \F-1" unit2=1/km unit1=Hz
       screenht=10 screenratio=0.7 color=g
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
       ''')

Flow('itfk','tfkcomplx',
     '''
     fft3 axis=3 inv=y |
     ltft inv=y verb=n |
     mutter v0=3.5 t0=-0.2 inner=n
     ''')
## Result('itfk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
##        title="LTFK transform"
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
##        ''')
Flow('sitfk','dat itfk','add scale=1,-1 ${SOURCES[1]}')
Result('sitfk',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       grey label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 title="LTFK decomposition"
       labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       ''')
## Result('wsitfk','sitfk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="LTFK Transform"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

# x-f-k transform
Flow('tdat','dat','transp plane=12 memsize=1000')
Flow('lxkt','tdat',
     '''
     ltft rect=5 verb=n niter=50 nw=192 dw=0.1531865 w0=0
     ''',split=[2,750],reduce="cat axis=3")
Flow('tlxkt','lxkt',
     '''
     transp plane=23 memsize=1000 |
     transp plane=12 memsize=1000
     ''')
Flow('fxk','tlxkt','fft3 axis=1')
Result('fxk',
       '''
       window min1=0 max1=50 | window n2=88 f2=5 | put o2=0.017 d1=0.166662 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=23 memsize=1000 | 
       grey3 color=j frame1=45 frame3=18 frame2=96 label1="\F5 F \F-1" flat=n 
       unit1=Hz label3=Offset label2="\F5 K \F-1" unit3=km unit2=1/km
       screenht=10 screenratio=0.7
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
       ''')
Flow('fxkmask1','fxk',
     '''
     real | transp plane=23 memsize=1000 | transp plane=12 memsize=1000 |
     mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5
     ''' )#mutter inner=y v0=0.3 |
Flow('fxkmask2','fxk',
     '''
     imag | transp plane=23 memsize=1000 | transp plane=12 memsize=1000 |
     mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5
     ''')#mutter inner=y v0=0.3 |
Flow('fxkcomplx','fxkmask1 fxkmask2',
     '''
     cmplx ${SOURCES[1:2]} |
     transp plane=12 memsize=1000 |
     transp plane=23 memsize=1000
     ''')
## Result('fxkcomplx',
##        '''
##        window min1=0 max1=50 | window n2=88 f2=5 | put o2=0.017 d1=0.166662 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 color=j frame1=45 frame3=18 frame2=96 label1="\F5 F \F-1" flat=n 
##        unit1=Hz label3=Offset label2="\F5 K \F-1" unit3=km unit2=1/km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
##        ''')

Flow('maskfxk1','fxk',
     '''
     real | transp plane=23 memsize=1000 | transp plane=12 memsize=1000 |
     math output=1. |
     mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5 |
     math output="1-input"
     ''' )#mutter inner=y v0=0.3 |
Flow('maskfxk2','fxk',
     '''
     imag | transp plane=23 memsize=1000 |  transp plane=12 memsize=1000 |
     math output=1. |
     mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5 |
     math output="1-input"
     ''')#mutter inner=y v0=0.3 |
Flow('fxkmask','maskfxk1 maskfxk2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=12 memsize=1000 |
     transp plane=23 memsize=1000
     ''')
Result('maskfxk','fxkmask',
       '''
       window min1=0 max1=50 | window n2=88 f2=5 | put o2=0.017 d1=0.166662 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=23 memsize=1000 | 
       grey3 frame1=45 frame3=18 frame2=96 label1="\F5 F \F-1" flat=n 
       unit1=Hz label3=Offset label2="\F5 K \F-1" unit3=km unit2=1/km
       screenht=10 screenratio=0.7 color=g
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
       ''')

Flow('ifxk','fxkcomplx',
     '''
     fft3 axis=1 inv=y |
     transp plane=12 memsize=1000 |
     transp plane=23 memsize=1000 |
     ltft inv=y verb=n |
     transp plane=12 memsize=1000 |
     mutter v0=3.5 t0=-0.2 inner=n
     ''')

## Result('ifxk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
##        title="LXFK Transform"
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
##        ''')
Flow('sifxk','dat ifxk','add scale=1,-1 ${SOURCES[1]}')
Result('sifxk',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       grey label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 title="LXFK decomposition"
       labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       ''')
## Result('wsifxk','sifxk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="LXFK Transform"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

# Cascading LTF transform
# Ltf transform
Flow('cltft','sifxk',
     '''
     ltft rect=20 verb=n nw=200 dw=0.5 niter=25
     ''',split=[2,93],reduce="cat axis=3")
## Result('cltft',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 color=j frame1=375 frame3=15 frame2=18 label1=Time flat=n 
##        unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f
##        ''')

Flow('cthr1','cltft',
     '''
     real | transp plane=23 memsize=1000 | 
     mutter t0=-0.2 v0=1.1 | cut min3=12
     ''' )#mutter -0.5 v0=0.55 |
Flow('cthr2','cltft',
     '''
     imag | transp plane=23 memsize=1000 |
     mutter t0=-0.2 v0=1.1 | cut min3=12
     ''')#mutter -0.5 v0=0.55 |
Flow('ccomplx','cthr1 cthr2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=23 memsize=1000
     ''')
## Result('ccomplx',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        grey3 color=j frame1=375 frame2=18 frame3=15 label1=Time flat=y 
##        unit1=s label3=Distance label2=Frequency unit3=Hz
##        point1=0.8 point2=0.2 wanttitle=n
##        ''')

Flow('cmaskthr1','cltft',
     '''
     real | transp plane=23 memsize=1000 |
     math output=1. |
     mutter t0=-0.5 v0=0.55 | cut min3=12 |
     math output="1-input"
     ''' )#mutter inner=y v0=0.3 |
Flow('cmaskthr2','cltft',
     '''
     imag | transp plane=23 memsize=1000 |
     math output=1. |
     mutter t0=-0.5 v0=0.55 | cut min3=12 |
     math output="1-input"
     ''')#mutter inner=y v0=0.3 |
Flow('cmaskcomplx','cmaskthr1 cmaskthr2',
     '''
     cmplx ${SOURCES[1:2]} | 
     transp plane=23 memsize=1000
     ''')
## Result('cmask','cmaskcomplx',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 frame1=375 frame3=15 frame2=18 label1=Time flat=n 
##        unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f
##        ''')

Flow('ciltft','ccomplx','ltft inv=y verb=n')
## Result('ciltft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 | 
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Noise"
##        ''')
Flow('csiltft','sifxk ciltft','add scale=1,-1 ${SOURCES[1]}')
Result('csiltft',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       grey label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7
       title="Cascading LXFK and LTF decompositions"
       labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       ''')
## Result('wcsiltft','csiltft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Cascading LXFK and LTF Transforms"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

## Flow('noise','dat csiltft','add scale=1,-1 ${SOURCES[1]}')
## Result('noise',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Difference"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        ''')
## Result('wnoise','noise',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Difference"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')
## #########################
## # S transform
## #########################
## Flow('st','dat','st',split=[2,93],reduce="cat axis=3")
## Result('st',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 color=j frame1=375 frame3=22 frame2=18 label1=Time flat=n 
##        unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f
##        ''')

## # x-f-k transform
## Flow('tstdat','dat','transp plane=12 memsize=1000')
## Flow('stxkt','tstdat','st',split=[2,750],reduce="cat axis=3")
## Flow('tstxkt','stxkt',
##      '''
##      transp plane=23 memsize=1000 |
##      transp plane=12 memsize=1000
##      ''')
## Flow('stfxk','tstxkt','fft3 axis=1')
## Result('stfxk',
##        '''
##        window min1=0 max1=50 | window n2=88 f2=5 | put o2=0.017 d1=0.166662 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 color=j frame1=45 frame3=18 frame2=96 label1="\F5 F \F-1" flat=n 
##        unit1=Hz label3=Offset label2="\F5 K \F-1" unit3=km unit2=1/km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
##        ''')
## Flow('stfxkmask1','stfxk',
##      '''
##      real | transp plane=23 memsize=1000 | transp plane=12 memsize=1000 |
##      mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5
##      ''' )#mutter inner=y v0=0.3 |
## Flow('stfxkmask2','stfxk',
##      '''
##      imag | transp plane=23 memsize=1000 | transp plane=12 memsize=1000 |
##      mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5
##      ''')#mutter inner=y v0=0.3 |
## Flow('stfxkcomplx','stfxkmask1 stfxkmask2',
##      '''
##      cmplx ${SOURCES[1:2]} |
##      transp plane=12 memsize=1000 |
##      transp plane=23 memsize=1000
##      ''')
## Result('stfxkcomplx',
##        '''
##        window min1=0 max1=50 | window n2=88 f2=5 | put o2=0.017 d1=0.166662 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 color=j frame1=45 frame3=18 frame2=96 label1="\F5 F \F-1" flat=n 
##        unit1=Hz label3=Offset label2="\F5 K \F-1" unit3=km unit2=1/km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
##        ''')

## Flow('stmaskfxk1','stfxk',
##      '''
##      real | transp plane=23 memsize=1000 | transp plane=12 memsize=1000 |
##      math output=1. |
##      mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5 |
##      math output="1-input"
##      ''' )#mutter inner=y v0=0.3 |
## Flow('stmaskfxk2','stfxk',
##      '''
##      imag | transp plane=23 memsize=1000 |  transp plane=12 memsize=1000 |
##      math output=1. |
##      mutter v0=1.6 inner=n | cut min3=1.5 | cut n3=5 |
##      math output="1-input"
##      ''')#mutter inner=y v0=0.3 |
## Flow('stfxkmask','stmaskfxk1 stmaskfxk2',
##      '''
##      cmplx ${SOURCES[1:2]} | 
##      transp plane=12 memsize=1000 |
##      transp plane=23 memsize=1000
##      ''')
## Result('stmaskfxk','stfxkmask',
##        '''
##        window min1=0 max1=50 | window n2=88 f2=5 | put o2=0.017 d1=0.166662 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 frame1=45 frame3=18 frame2=96 label1="\F5 F \F-1" flat=n 
##        unit1=Hz label3=Offset label2="\F5 K \F-1" unit3=km unit2=1/km
##        screenht=10 screenratio=0.7 color=g
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f o2num=0 d2num=10 n2tic=6
##        ''')

## Flow('istfxk','stfxkcomplx',
##      '''
##      fft3 axis=1 inv=y |
##      transp plane=12 memsize=1000 |
##      transp plane=23 memsize=1000 |
##      st inv=y |
##      transp plane=12 memsize=1000 |
##      mutter v0=3.5 t0=-0.2 inner=n
##      ''')

## Result('istfxk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
##        title="LXFK Transform"
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
##        ''')
## Flow('stifxk','dat istfxk','add scale=1,-1 ${SOURCES[1]}')
## Result('stifxk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="LXFK Transform"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        ''')
## Result('wstifxk','stifxk',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="LXFK Transform"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

## # Cascading STTF transform
## # sttf transform
## Flow('csttft','stifxk','st',split=[2,93],reduce="cat axis=3")
## Result('csttft',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 color=j frame1=375 frame3=15 frame2=18 label1=Time flat=n 
##        unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f
##        ''')

## Flow('cstthr1','csttft',
##      '''
##      real | transp plane=23 memsize=1000 | 
##      mutter t0=-0.2 v0=1.1 | cut min3=12
##      ''' )#mutter -0.5 v0=0.55 |
## Flow('cstthr2','csttft',
##      '''
##      imag | transp plane=23 memsize=1000 |
##      mutter t0=-0.2 v0=1.1 | cut min3=12
##      ''')#mutter -0.5 v0=0.55 |
## Flow('cstcomplx','cstthr1 cstthr2',
##      '''
##      cmplx ${SOURCES[1:2]} | 
##      transp plane=23 memsize=1000
##      ''')
## Result('cstcomplx',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        grey3 color=j frame1=375 frame2=18 frame3=15 label1=Time flat=y 
##        unit1=s label3=Distance label2=Frequency unit3=Hz
##        point1=0.8 point2=0.2 wanttitle=n
##        ''')

## Flow('cstmaskthr1','csttft',
##      '''
##      real | transp plane=23 memsize=1000 |
##      math output=1. |
##      mutter t0=-0.5 v0=0.55 | cut min3=12 |
##      math output="1-input"
##      ''' )#mutter inner=y v0=0.3 |
## Flow('cstmaskthr2','csttft',
##      '''
##      imag | transp plane=23 memsize=1000 |
##      math output=1. |
##      mutter t0=-0.5 v0=0.55 | cut min3=12 |
##      math output="1-input"
##      ''')#mutter inner=y v0=0.3 |
## Flow('cstmaskcomplx','cstmaskthr1 cstmaskthr2',
##      '''
##      cmplx ${SOURCES[1:2]} | 
##      transp plane=23 memsize=1000
##      ''')
## Result('cstmask','cstmaskcomplx',
##        '''
##        window n3=88 f3=5 max2=50 | put o3=0.017 |
##        math output="abs(input)" | real |
##        byte allpos=y gainpanel=40 pclip=100 |
##        transp plane=23 memsize=1000 | 
##        grey3 frame1=375 frame3=15 frame2=18 label1=Time flat=n 
##        unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
##        screenht=10 screenratio=0.7
##        point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f
##        ''')

## Flow('cisttft','cstcomplx','st inv=y')
## Result('cisttft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 | 
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Noise"
##        ''')
## Flow('csisttft','stifxk cisttft','add scale=1,-1 ${SOURCES[1]}')
## Result('csisttft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Cascading LXFK and LTF Transforms"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        ''')
## Result('wcsisttft','csisttft',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        agc rect1=50 repeat=5 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Cascading LXFK and LTF Transforms"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

## Flow('stnoise','dat csisttft','add scale=1,-1 ${SOURCES[1]}')
## Result('stnoise',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Difference"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        ''')
## Result('wstnoise','stnoise',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |
##        wiggle label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 title="Difference"
##        labelfat=4 font=2 titlefat=4
##        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
##        transp=y yreverse=y poly=y
##        ''')

#########################
# RT transform
#########################
# Comparison with RTT
for nw in (2,4,8):
    rad = 'rad%d' %nw
    Flow(rad,'dat',
         'transp | radial vmin=0 vmax=2.0 nv=1020 tp=0. nw=%d | transp' %
         nw)
##     Result(rad,
##            '''
##            window n2=1000 f2=20 |
##            grey label2="Radial velocity (km/s)" label1="Time (s)"
##            screenht=10 screenratio=1.4 clip=300 wanttitle=n
##            labelfat=4 font=2 titlefat=4
##            ''')
    band = 'band%d' % nw
    Flow(band,rad,'bandpass flo=3')
##     Result(band,
##            '''
##            window n2=1000 f2=20 |
##            grey label2="Radial velocity (km/s)" label1="time (s)"
##            screenht=10 screenratio=1.4 clip=300
##            wanttitle=n title="B-%d bandpass"
##            ''' % nw)    
    result = 'result%d' % nw
    Flow(result,band,
         '''
         transp | radial tp=0. nw=%d inv=y | transp
         ''' % nw)
##     Result(result,
##            '''
##            window n2=88 f2=5 | put o2=0.017 | 
##            grey label2="Offset (km)" label1="Time (s)" 
##            screenht=10 screenratio=1.4 clip=300
##            wanttitle=n labelfat=4 font=2 titlefat=4
##            ''' )

    err = 'err%d' % nw
    Flow(err,[rad,'dat'],
         '''
         transp | radial tp=0. nw=%d inv=y | transp |
         add scale=1,-1 ${SOURCES[1]}
         ''' % nw)
##     Result(err,
##            '''
##            window n2=88 f2=5 | put o2=0.017 | 
##            grey label2="Offset (km)" label1="time (s)" 
##            screenht=10 screenratio=1.4 clip=300
##            wanttitle=n title="B-%d error"
##            ''' % nw)

##     noi = 'noi%d' % nw
##     Flow(noi,['dat',result],
##          '''
##          add scale=1,-1 ${SOURCES[1]}
##          ''')
##     Result(noi,
##            '''
##            window n2=88 f2=5 | put o2=0.017 | 
##            grey label2="Offset (km)" label1="time (s)" 
##            screenht=10 screenratio=1.4 clip=300
##            wanttitle=n title="B-%d noise model"
##            ''' % nw)

Result('resu8','result8',
       '''
       window n2=88 f2=5 | put o2=0.017 | 
       agc rect1=50 repeat=5 |
       grey label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 
       title="RTT + high-pass filter" labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
       ''')
Result('wresu8','result8',
       '''
       window n2=88 f2=5 | put o2=0.017 | 
       agc rect1=50 repeat=5 |
       wiggle label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 
       title="RTT + high-pass filter" labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       transp=y yreverse=y poly=y
       ''')

## Result('noi8',
##        '''
##        window n2=88 f2=5 | put o2=0.017 | 
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 clip=300 wanttitle=n title="model noise"
##        ''')

## # Adaptive subtraction
## nt = 750
## ns = 11
## 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')

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

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

## # Adaptive Subtraction using regularized nonstationary regression
## Flow('aflt apre','shifts dat pef',
##      '''
##      lpf match=${SOURCES[1]} pred=${TARGETS[1]} pef=${SOURCES[2]}
##      rect1=10 rect2=10 niter=400
##      ''')
## Flow('asig','dat apre','add scale=1,-1 ${SOURCES[1]}')
## Plot('apre',
##      '''
##      window n2=88 f2=5 | put o2=0.017 | 
##      grey label2="Offset (km)" label1="Time (s)" 
##      screenht=10 screenratio=1.4 clip=300 labelfat=4 font=2 titlefat=4
##      ''')
## Plot('asig',
##      '''
##      window n2=88 f2=5 | put o2=0.017 | 
##      grey label2="Offset (km)" label1="Time (s)" 
##      screenht=10 screenratio=1.4 clip=300 labelfat=4 font=2 titlefat=4
##      ''')

## Plot('azero','aflt',
##      '''
##      window n3=1 f3=10 |
##      grey bias=0.1 scalebar=y label2="Offset (km)" label1="Time (s)" 
##      screenht=10 screenratio=1.4 wanttitle=n labelfat=4 font=2 titlefat=4
##      ''')
## Plot('acsum','aflt',
##      '''
##      stack axis=3 norm=n |
##      grey bias=0.1 scalebar=y label2="Offset (km)" label1="Time (s)" 
##      screenht=10 screenratio=1.4 wanttitle=n labelfat=4 font=2 titlefat=4
##      ''')

## Plot('cube','aflt', 
##      '''
##      window n2=88 f2=5 | put o2=0.017 | 
##      byte gainpanel=all allpos=n bias=0.1 |
##      grey3  wanttitle=n title="3D filter volume"
##      point2=0.5 point1=0.8 
##      frame1=300 frame2=44 frame3=11 flat=n color=j
##      ''')

## # Adaptive Subtraction with stationary matching filter
## Flow('flt pre','shifts dat',
##      '''
##      lpf match=${SOURCES[1]} pred=${TARGETS[1]}
##      rect1=1000 rect2=1000 niter=400
##      ''')
## Flow('sig','dat pre','add scale=1,-1 ${SOURCES[1]}')
## Plot('pre',
##      '''
##      window n2=88 f2=5 | put o2=0.017 | 
##      grey label2="Offset (km)" label1="time (s)" 
##      screenht=10 screenratio=1.4 clip=300
##      wanttitle=n title="Predicted noise (Stationary)"
##      ''')
## Plot('sig',
##      '''
##      window n2=88 f2=5 | put o2=0.017 | 
##      grey label2="Offset (km)" label1="time (s)" 
##      screenht=10 screenratio=1.4 clip=300
##      title="Predicted signal (Stationary)"
##      ''')

## Plot('zero','flt',
##      '''
##      window n3=1 f3=10 |
##      grey color=j bias=0.1 scalebar=y
##      formatbar="%3.3f" screenht=10 screenratio=1.4 wanttitle=n 
##      ''')
## Plot('csum','flt',
##      '''
##      stack axis=3 norm=n |
##      grey color=j bias=0.0 scalebar=y
##      formatbar="%3.3f" screenht=10 screenratio=1.4 wanttitle=n 
##      ''')

## Flow('rad','asig',
##      'transp | radial vmin=0 vmax=2.0 nv=1020 tp=0. nw=8 | transp')
## Plot('rad',
##      '''
##      window n2=1000 f2=20 |
##      grey label2="Radial velocity (km/s)" label1="time (s)"
##      screenht=10 screenratio=1.4 clip=300 wanttitle=n title="Adaptive"
##      ''')

#########################
# FK filter
#########################
# Comparison with FK filtering
Flow('fk','dat','pad end1=750 | pad end2=93 | fft1 | fft3 axis=2')
## Result('fks','fk',
##        '''
##        real |
##        grey color=j screenht=10 screenratio=1.4
##        title="F-K spectra"
##        ''')
Flow('rmutter','fk','real | mutter v0=3. ')
Flow('imutter','fk','imag | mutter v0=3. ')
## Result('rfks','rmutter',
##        '''
##        grey color=j screenht=10 screenratio=1.4
##        title="F-K spectra"
##        ''')

Flow('mutter','rmutter imutter','cmplx ${SOURCES[:2]}')
Flow('ifk','mutter',
     '''
     fft3 axis=2 inv=y |
     fft1 inv=y |
     window n2=93 | window n1=750 |
     mutter v0=3.5 t0=-0.2 inner=n
     ''')
## Result('ifks','ifk',
##        '''
##        spectra2 |
##        grey color=j crowd1=0.3 wanttitle=n
##        labelfat=4 font=2 titlefat=4
##        ''')
Result('ifk',
        '''
        window n2=88 f2=5 | put o2=0.017 |
        agc rect1=50 repeat=5 |
        grey label2="Offset" unit2=km label1="Time" unit1=s
        screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
        title="FK filter"
        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
        ''')
Result('wifk','ifk',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       wiggle label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
       title="FK filter"
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5
       transp=y yreverse=y poly=y
       ''')

## Flow('inoi','dat ifk',
##      '''
##      add scale=1,-1 ${SOURCES[1]}
##      ''')
## Result('inoi',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |        
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 wanttitle=n labelfat=4 font=2 titlefat=4
##        ''')

#########################
# Bandpass filter
#########################
# Comparison with highpass filtering
Flow('brshot','dat','bandpass flo=10')
Result('brshot',
        '''
        window n2=88 f2=5 | put o2=0.017 |
        agc rect1=50 repeat=5 |
        grey label2="Offset" unit2=km label1="Time" unit1=s
        screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
        title="High-pass filter"
        parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
        ''' )
Result('wbrshot','brshot',
       '''
       window n2=88 f2=5 | put o2=0.017 |
       agc rect1=50 repeat=5 |
       wiggle label2="Offset" unit2=km label1="Time" unit1=s
       screenht=10 screenratio=0.7 labelfat=4 font=2 titlefat=4
       title="High-pass filter"
       parallel2=n format2=%3.1f n2tic=6 o2num=0. d2num=0.5.
       transp=y yreverse=y poly=y
       ''' )

## Flow('brnoi','dat brshot',
##      '''
##      add scale=1,-1 ${SOURCES[1]}
##      ''')
## Result('brnoi',
##        '''
##        window n2=88 f2=5 | put o2=0.017 |        
##        agc rect1=50 repeat=5 |
##        grey label2="Offset" unit2=km label1="Time" unit1=s
##        screenht=10 screenratio=0.7 wanttitle=n labelfat=4 font=2 titlefat=4
##        ''')

## #########################
## # RTT + Median filter
## #########################
## # Comparison with Median filter in RTT domain
## Flow('median','rad8','mf boundary=y nfw=31')
## Result('median',
##        '''
##        window n2=1000 f2=20 |
##        grey label2="Radial velocity (km/s)" label1="time (s)"
##        screenht=10 screenratio=1.4 clip=300
##        wanttitle=n title="MF filtering in RTT"
##        ''')

## Flow('res','rad8 median',
##      '''
##      add scale=1,-1 ${SOURCES[1]}
##      ''')
## Result('res',
##        '''
##        window n2=1000 f2=20 |
##        grey label2="Radial velocity (km/s)" label1="time (s)"
##        screenht=10 screenratio=1.4 clip=300 wanttitle=n title="Residual"
##        ''')

## Flow('rmedian','res',
##      '''
##      transp | radial tp=0. nw=8 inv=y | transp
##      ''')
## Result('rmedian',
##        '''
##        window n2=88 f2=5 | put o2=0.017 | 
##        grey label2="Offset (km)" label1="time (s)" 
##        screenht=10 screenratio=1.4 clip=300  wanttitle=n title="MF result"
##        ''' )

## #########################
## # LRMF
## #########################
## # Comparison with LRMF
## Flow('dat1','hector.HH','dd form=native')

## Result('dat1',
##        '''
##        grey label2="Offset (km)" label1="Time (s)" title="Data"
##        screenht=10 screenratio=1.4 clip=300
##        ''')
## Flow('lrmf','dat1',
##      'lrmf nfw=31 tmax=0.5 tmin=0. vmax=0.6. vmin=0.2 boundary=n')
## Result('lrmf',
##        '''
##        grey label2="Offset (km)" label1="Time (s)" title="Result"
##        screenht=10 screenratio=1.4 clip=300
##        ''')
## Flow('signal','dat1 lrmf','add scale=1,-1 ${SOURCES[1]} ')
## Result('signal',
##        '''
##        grey label2="Offset (km)" label1="Time (s)" title="Result"
##        screenht=10 screenratio=1.4 clip=300
##        ''')
## #########################
## # STFT
## #########################
Flow('stft','dat',
     '''
     stft ntw=100
     ''',split=[2,93],reduce="cat axis=3")
Result('stft',
       '''
       window n3=88 f3=5 max2=50 | put o3=0.017 |
       math output="abs(input)" | real |
       byte allpos=y gainpanel=40 pclip=100 |
       transp plane=23 memsize=1000 | 
       grey3 color=j frame1=375 frame3=3 frame2=18 label1=Time flat=n 
       unit1=s label2=Offset label3="\F5 F \F-1" unit2=km
       screenht=10 screenratio=0.7
       point1=0.85 point2=0.8 wanttitle=n labelfat=4 font=2 titlefat=4
       parallel2=n format2=%3.1f
       ''')

End()