from rsf.proj import *

par = dict(
    ny=501,
    nx=501,
    nz=501,
    dy=0.004,
    dx=0.004,
    dz=0.004,
    y0=0.0,
    x0=0.0,
    z0=0.0,

    ns=700,
    dt=0.0004,

    vp0_1=2500,
    vs0_1=1200,
    eps_1=0.25,
    del_1=-0.25,
    gam_1=0.3,
    the_1=0.0,
    vp0_2=3600,
    vs0_2=1800,
    eps_2=0.2,
    del_2=0.1,
    gam_2=0.05,
    the_2=0.0,

        bd=30,

    hny=250,
    hnx=250,
    hnz=250,

    itaper=0,

    seed=2012,
    eps=1.e-6,
    npk=30,

    oy=1.0,
    ox=1.0,

    )

# =================================================================================
Flow(['vp0','vs0','epsi','del','the','gam'], None,
     '''
         twolayer3dti
         vs0=${TARGETS[1]}
         epsi=${TARGETS[2]}
         del=${TARGETS[3]}
         the=${TARGETS[4]}
         gam=${TARGETS[5]}
         ny=%d
         nx=%d
         nz=%d
         dy=%g
         dx=%g
         dz=%g
         vp0_1=%g
         vs0_1=%g
         eps_1=%g
         del_1=%g
         gam_1=%g
         the_1=%g
         vp0_2=%g
         vs0_2=%g
         eps_2=%g
         del_2=%g
         gam_2=%g
         the_2=%g
         ''' % (par['ny'],par['nx'],par['nz'],par['dy'],par['dx'],par['dz'],
                par['vp0_1'],par['vs0_1'],par['eps_1'],par['del_1'],par['gam_1'],par['the_1'],
                par['vp0_2'],par['vs0_2'],par['eps_2'],par['del_2'],par['gam_2'],par['the_2'])
    )

Plot('vp0','vp0',
      '''
      byte |
      grey3 polarity=n scalebar=n frame1=251 frame2=221 frame3=201 wanttitle=n flat=y screenratio=1
      ''')

# =================================================================================
# calculate pseudo-pure-mode qSV-wavefield
# =================================================================================
name1='''
PseudoPureSVxy PseudoPureSVz PseudoPureSV
'''

Flow(['PseudoPureSVxy',  'PseudoPureSVz', 'PseudoPureSV'],'vp0 vs0 epsi del',
      '''
         vti3dpseudosv
         vp0=${SOURCES[0]}
         vs0=${SOURCES[1]}
         epsi=${SOURCES[2]}
         del=${SOURCES[3]}
         PseudoPureSVz=${TARGETS[1]}
         PseudoPureSV=${TARGETS[2]}
         ns=%d 
         dt=%g
                 bd=%d
     ''' % (par['ns'],par['dt'],par['bd'])
)

Flow('interface','PseudoPureSV',
     '''
     math output=.83333333333333333334
     '''
    )

Plot('interface',
     '''
     graph3 polarity=n min=0 max=2 plotfat=1 dash=5 polarity=n scalebar=n wanttitle=n wantaxis=n flat=y screenratio=1.
     '''
     )

for qq in Split(name1):
       Plot(qq,
        '''
                byte gainpanel=all |
        grey3 polarity=n color= scalebar=n frame1=251 frame2=221 frame3=201 screenratio=1 wanttitle=n flat=y 
        ''')

name11='''
PseudoPureSV2dy1km PseudoPureSV2dx1km PseudoPureSV2dx2dyDiff
'''

Flow('PseudoPureSV2dy1km','PseudoPureSV',
        '''
         window n3=1 min3=%g
        ''' % (par['oy'])
        )
Flow('PseudoPureSV2dx1km','PseudoPureSV',
        '''
         window n2=1 min2=%g
        ''' % (par['ox'])
        )
Flow('PseudoPureSV2dx2dyDiff','PseudoPureSV2dy1km PseudoPureSV2dx1km',
        '''
         math x=${SOURCES[0]} y=${SOURCES[1]} output='x-y'
         '''
        )

for qq in Split(name11):
       Plot(qq,
        '''
        grey polarity=n color= scalebar=n screenratio=1 wanttitle=n
                ''')

for qq in Split(name1):
    Result(qq+'Interf',[qq,'interface'],'Overlay')

Result('vp0'+'Interf',['vp0','interface'],'Overlay')
# =================================================================================
# calculate pseudo-pure-mode SH-wavefield
# =================================================================================
name2='''
SHx SHy SH
'''

Flow(['SHx',  'SHy', 'SH'],'vp0 vs0 epsi del gam',
     '''
         vti3dpseudosh
         vp0=${SOURCES[0]}
         vs0=${SOURCES[1]}
         epsi=${SOURCES[2]}
         del=${SOURCES[3]}
         gam=${SOURCES[4]}
         SHy=${TARGETS[1]}
         SH=${TARGETS[2]}
         ns=%d 
         dt=%g
         bd=%d
     ''' % (par['ns'],par['dt'],par['bd'])
)

for qq in Split(name2):
       Plot(qq,
        '''
        byte gainpanel=all |
        grey3 polarity=n color= scalebar=n frame1=251 frame2=221 frame3=201 screenratio=1 wanttitle=n flat=y 
        ''')

for qq in Split(name2):
    Result(qq+'Interf',[qq,'interface'],'Overlay')

End()