from rsf.proj import *

Fetch('salt_slow_desp.HH','segeage')
Flow('salt','salt_slow_desp.HH',
     '''
     dd form=native |
     math output=1/input |
     scale dscale=0.001 |
     put d1=0.02 d2=0.02 d3=0.02
     label1=Z label2=X label3=Y
     unit1=km unit2=km unit3=km
     ''')

Result('salt',
       '''
       byte clip=1 bias=2.5 |
       grey3 color=j title="SEG/EAGE Salt Model" 
       frame2=350 frame3=350 frame1=100 point1=0.4
       ''')

Result('salt3','salt',
       '''
       byte clip=1 bias=2.5 |
       grey3 frame1=100 frame2=350 frame3=350 point1=0.4 point2=0.7
       label1=Z label2=X label3=Y flat=n
       color=j title="SEG/EAGE Salt Model" 
       ''')


#remap1 n1=159 d1=%g order=2 |
#     transp plane=12 memsize=1000 | remap1 n1=159 d1=%g order=2 | transp plane=12 memsize=1000 |
#     transp plane=13 memsize=1000 | remap1 n1=159 d1=%g order=2 | transp plane=13 memsize=1000

Flow('win','salt',
     '''
     window n2=210 n3=210 f2=200 f3=200 
     ''')

Result('win',
       '''
       byte clip=1 bias=2.5 bar=bar.rsf barreverse=y |
       grey3 color=j title="Portion of SEG/EAGE Salt Model" flat=n
       frame2=100 frame3=100 frame1=75 scalebar=y bartype=v
       barlabel=Velocity barunit=km/s barreverse=y point1=0.66 point2=0.66
       ''')


dx5 = 209*0.02/78
dx10 = 209*0.02/158

for w in (5,10):
    dx = 209*0.02/(16*w-2)
    sol = 'sol%d' % w
    helm = 'helm%d' % w

#    Flow(helm,sol,
#         '''
#         real |
#         put o1=0 o2=4 o3=4 d1=%g d2=%g d3=%g |
#         remap1 n1=210 d1=0.02 order=4 |
#         transp plane=12 memsize=1000 | remap1 n1=210 d1=0.02 order=2 | transp plane=12 memsize=1000 |
#         transp plane=13 memsize=1000 | remap1 n1=210 d1=0.02 order=2 | transp plane=13 memsize=1000
#         ''' % (dx,dx,dx))

#    Result(helm,
#           '''
#           byte gainpanel=all |
#           grey3 title="Solution at %d Hz" flat=n
#           frame2=100 frame3=100 frame1=75 
#           point1=0.66 point2=0.66
#           label1=Z label2=X label3=Y unit1=km unit2=km unit3=km
#           ''' % w)

Flow('salt2','salt','window min2=2 max2=12 min3=2 max3=12 j3=2')

# Lowrank decomposition

nt=1751
dt=0.001

Flow('fft','salt2','fft1 | fft3 axis=2 pad=1 | fft3 axis=3 pad=1')
Flow('right left','salt2 fft',
     'isolr3 seed=2010 npk=10 eps=0.00001 dt=%g fft=${SOURCES[1]} left=${TARGETS[1]}' % dt)

Flow('source',None,
     '''
     spike n1=%d d1=%g k1=100 | 
     ricker1 frequency=15
     '''%(nt,dt))

Flow('refl',None,
     '''
     spike n1=210 n2=501 n3=251 d1=0.02 d2=0.02 d3=0.04 o2=2 o3=2 unit1=km k1=3 k2=251 k3=126 |
     smooth rect1=2 rect2=2 rect3=3 repeat=3
     ''')

Flow('wave3 snaps','source refl left right',
     '''
     fftwave3 ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y
     snap=25 snaps=${TARGETS[1]}
     ''')

Result('wave3',
       '''
       byte gainpanel=all |
       grey3 frame1=50 frame2=250 frame3=125 point1=0.4 point2=0.7
       wanttitle=n label1=Z label2=X label3=Y flat=n
       ''')

Plot('snaps2','snaps',
     '''
     put n4=71 |
     byte gainpanel=all |
     grey4 frame1=50 frame2=250 frame3=125 point1=0.4 point2=0.7
     wanttitle=n label1=Z label2=X label3=Y flat=n
     ''')

End()