from rsf.proj import *

Fetch('Gulfnear.HH','mig')

Flow('data','Gulfnear.HH',
     '''
     dd form=native |
     window f2=200 n2=480 n1=1000 |
     bandpass flo=10 fhi=50 | put o3=0
     ''')

Plot('data',
     '''
     window f2=60 n2=360 max1=2.0 |
     pow pow1=0.5 | 
     grey label1=Time unit1=s label2=Midpoint unit2=m 
     labelfat=2 labelsz=4 titlefat=2 titlesz=5 
     xll=3.0 yll=5.28  xur=11.0 yur=9.77  title="(a)" 
     ''')

vwater=1500
alpha=460
beta=0.5
twater=0.20

Flow('water','data','window n2=1 | spike | cut max1=%g' % twater)

Flow('vel','water',
     '''
     math output="%g+input*(%g * abs(x1 - %g)^%g)" 
     ''' % (vwater,alpha,twater,beta))
Flow('lowvel','water',
     '''
     math output="%g+input*(%g * abs(x1 - %g)^%g)" 
     ''' % (vwater,0.0,twater,beta))

Flow('vel2',      'vel','spray axis=2 n=480 d=25 o=5000')
Flow('lowvel2','lowvel','spray axis=2 n=480 d=25 o=5000')

Flow('migaa','data vel','kirchnew velocity=${SOURCES[1]}')

Flow('migaa2','data vel2','mig2 vel=${SOURCES[1]}')

plotb = '''
window f2=60 n2=360 max1=2.0 | pow pow1=0.5 |
grey label1=Time unit1=s wantaxis2=n
labelfat=2 labelsz=4 titlefat=2 titlesz=5 
xll=3.0 yll=0.32  xur=11.0 yur=4.81  title="(b)" 
'''

plota = '''
window f2=60 n2=360 max1=2.0 | pow pow1=0.5 |
grey label1=Time unit1=s label2=Midpoint unit2=m 
labelfat=2 labelsz=4 titlefat=2 titlesz=5 
xll=3.0 yll=5.28  xur=11.0 yur=9.77  title="(a)" 
'''

Plot('migaa2',plotb)

Result('Gulfaa','data migaa2','Overlay')

Flow('mig90','data vel2','mig2 vel=${SOURCES[1]} antialias=0')

Flow('mig45','data vel2','mig2 vel=${SOURCES[1]} antialias=0 angle=45')

Plot('mig90',plota)
Plot('mig45',plotb)

Result('Gulfaper','mig90 mig45','Overlay')

Plot('migaa-zoom','migaa2',
     '''
     window max1=0.8 min2=9000 max2=14000 | pow pow1=0.5 |
     grey xll=1.5 yll=5.30 xur=6.5  yur=9.77 title="(a)" 
     label1=Time unit1=s label2=Midpoint unit2=m labelsz=4 titlesz=5
     ''')
Plot('mig90-zoom','mig90',
     '''
     window max1=0.8 min2=9000 max2=14000 | pow pow1=0.5 |
     grey xll=1.5 yll=0.32 xur=6.5 yur=4.79 title="(b)" 
     label1=Time unit1=s wantaxis2=n wantaxis2=n labelsz=4 titlesz=5
     ''')
Plot('arti1',None,'box x0=3.2 y0=4.2 xt=1. yt=0.3 size=.15 label="artifacts" ')
Plot('arti2',None,'box x0=5.7 y0=1.8 xt=.4 yt=1.0 size=.15 label="artifacts" ')

Result('Zoom2','migaa-zoom mig90-zoom arti1 arti2','Overlay')

Flow('lowmigaa2','data lowvel2','mig2 vel=${SOURCES[1]}')
Flow('lowmig90', 'data lowvel2','mig2 vel=${SOURCES[1]} antialias=0')

Plot('lowmigaa2',plota)
Plot('lowmig90',plotb)

Result('Miglow','lowmigaa2 lowmig90','Overlay')

End()