from rsf.proj import *

par = dict(
    nx=401,
    nz=401,
    dx=0.005,
    dz=0.005,
    x0=0.0,
    z0=0.0,

    ns=601,
    dt=0.0005,

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

    seed=2012,
    eps=1.e-6,
    npk=30,
    ireconstruct=1,
    xrec1=1000.0,
    zrec1=500.0,
    xrec2=1000.0,
    zrec2=1500.0,
    )

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

# =================================================================================
# calculate separated wave-mode:         yes
# =================================================================================
name1='''
Elasticx Elasticz ElasticSepP ElasticSepSV
'''

name2='''
Polxp1 Polzp1 Errpolxp1 Errpolzp1
Polxp2 Polzp2 Errpolxp2 Errpolzp2
'''

if(par['ireconstruct']==1):
    Flow(['Elasticx',  'Elasticz', 'ElasticSepP', 'ElasticSepSV',
          'Errpolxp1', 'Errpolxp2',  'Polxp1', 'Polxp2', 'Errpolzp1', 'Errpolzp2', 'Polzp1', 'Polzp2'],
          'vp0  vs0  epsi del the',
         '''
         tti2delrsep2p
         vp0=${SOURCES[0]}
         vs0=${SOURCES[1]}
         epsi=${SOURCES[2]} 
         del=${SOURCES[3]}
         the=${SOURCES[4]}
         Elasticz=${TARGETS[1]}
         ElasticSepP=${TARGETS[2]}
         ElasticSepSV=${TARGETS[3]}
         Errpolxp1=${TARGETS[4]}
         Errpolxp2=${TARGETS[5]}
         Polxp1=${TARGETS[6]}
         Polxp2=${TARGETS[7]}
         Errpolzp1=${TARGETS[8]}
         Errpolzp2=${TARGETS[9]}
         Polzp1=${TARGETS[10]}
         Polzp2=${TARGETS[11]}
         ns=%d 
         dt=%g
         eps=%g
         seed=%d
         npk=%d
         ireconstruct=%d
         xrec1=%g
         zrec1=%g
         xrec2=%g
         zrec2=%g
         ''' % (par['ns'],par['dt'],par['eps'],par['seed'],par['npk'],par['ireconstruct'],par['xrec1'],par['zrec1'],par['xrec2'],par['zrec2'])
    )
else:
    Flow(['Elasticx',  'Elasticz',  'ElasticSepP', 'ElasticSepSV'],
         'vp0  vs0  epsi del the',
         '''
         tti2delrsep2p
         vp0=${SOURCES[0]}
         vs0=${SOURCES[1]}
         epsi=${SOURCES[2]} 
         del=${SOURCES[3]}
         the=${SOURCES[4]}
         Elasticz=${TARGETS[1]}
         ElasticSepP=${TARGETS[2]}
         ElasticSepSV=${TARGETS[3]}
         ns=%d 
         dt=%g
         eps=%g
         seed=%d
         npk=%d
         ireconstruct=%d
         xrec1=%g
         zrec1=%g
         xrec2=%g
         zrec2=%g
         ''' % (par['ns'],par['dt'],par['eps'],par['seed'],par['npk'],par['ireconstruct'],par['xrec1'],par['zrec1'],par['xrec2'],par['zrec2'])
    )

Result('vp0','vp0',
      '''
      grey color=t scalebar=y bias=1.5 allpos=n barreverse=y wanttitle=n screenratio=1.
      ''')

Flow('interface','Elasticx',
     '''
     math output="(%d-1)*%g*7/12"
     ''' % (par['nz'],par['dz'])
    )

Plot('interface',
     '''
     graph yreverse=y color= transp=n pad=n plotfat=1
     min1=0 max1=2 
     min2=0 max2=2
     wanttitle=n wantaxis=n screenratio=1.
     '''
     )

for qq in Split(name1):
        Plot(qq,
        '''
        grey polarity=n scalebar=n screenratio=1. wanttitle=n
        axisfat=5 axiscol=7 labelsz=10
        ''')

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

if(par['ireconstruct']==1):
   for tt in Split(name2):
        Result(tt,
        '''
        grey color=t polarity=n scalebar=y screenratio=0.87 wanttitle=n pclip=100
        ''')
# =================================================================================
# Vector decomposition
# =================================================================================
name3='''
ElasticX ElasticZ ElasticPx ElasticPz ElasticSVx ElasticSVz
'''

name4='''
Decompxp1 Decompzp1 Decompxzp1 Errdecxp1 Errdeczp1 Errdecxzp1
Decompxp2 Decompzp2 Decompxzp2 Errdecxp2 Errdeczp2 Errdecxzp2
'''

Flow(['ElasticX',  'ElasticZ',  'ElasticPx', 'ElasticPz', 'ElasticSVx', 'ElasticSVz',
      'Decompxp1', 'Decompzp1', 'Decompxzp1', 'Errdecxp1', 'Errdeczp1', 'Errdecxzp1',
      'Decompxp2', 'Decompzp2', 'Decompxzp2', 'Errdecxp2', 'Errdeczp2', 'Errdecxzp2'],
      'vp0  vs0  epsi del the',
         '''
         tti2delrdecomp2p
         vp0=${SOURCES[0]}
         vs0=${SOURCES[1]}
         epsi=${SOURCES[2]} 
         del=${SOURCES[3]}
         the=${SOURCES[4]}
         ElasticZ=${TARGETS[1]}
         ElasticPx=${TARGETS[2]}
         ElasticPz=${TARGETS[3]}
         ElasticSVx=${TARGETS[4]}
         ElasticSVz=${TARGETS[5]}
         Decompxp1=${TARGETS[6]}
         Decompzp1=${TARGETS[7]}
         Decompxzp1=${TARGETS[8]}
         Errdecxp1=${TARGETS[9]}
         Errdeczp1=${TARGETS[10]}
         Errdecxzp1=${TARGETS[11]}
         Decompxp2=${TARGETS[12]}
         Decompzp2=${TARGETS[13]}
         Decompxzp2=${TARGETS[14]}
         Errdecxp2=${TARGETS[15]}
         Errdeczp2=${TARGETS[16]}
         Errdecxzp2=${TARGETS[17]}
         ns=%d 
         dt=%g
         eps=%g
         seed=%d
         npk=%d
         ireconstruct=%d
         xrec1=%g
         zrec1=%g
         xrec2=%g
         zrec2=%g
         ''' % (par['ns'],par['dt'],par['eps'],par['seed'],par['npk'],par['ireconstruct'],par['xrec1'],par['zrec1'],par['xrec2'],par['zrec2'])
    )

for qq in Split(name3):
        Plot(qq,
        '''
        grey polarity=n scalebar=n screenratio=1. wanttitle=n
        axisfat=5 axiscol=7 labelsz=10
        ''')

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

if(par['ireconstruct']==1):
   for tt in Split(name4):
        Result(tt,
        '''
        grey color=t polarity=n scalebar=y screenratio=0.87 wanttitle=n pclip=100
        ''')
End()