## 
 # elastic modeling; wavefield separation
 ## testing the sensitivity of anisotropy parameters 
# to the separation of VTI elastic wavefields

from rsf.proj import *
import sys
sys.path.append('Python')
import fdmod,pot,fdd,spk,stiff
# ------------------------------------------------------------
par = {
    'nx':600, 'ox':0, 'dx':0.002,  'lx':'x', 'ux':'km',
    'nz':600, 'oz':0, 'dz':0.002,  'lz':'z', 'uz':'km',
    'nt':1601,'ot':0, 'dt':0.0002, 'lt':'t', 'ut':'s',
    'kt':150,
    'jsnap':200,
    'height':10,
    'nb':0,
    'frq':75,
    'ratio':1
    }
fdmod.param(par)
par['labelattr']=par['labelattr']+''' labelsz=10
n2tic=7 o2num=0 d2num=0.2
n1tic=7 o1num=0 d1num=0.2
'''
nframe=8
order=8
# ------------------------------------------------------------
# source/receiver coordinates
# ------------------------------------------------------------
fdmod.point('ss',
            par['ox']+(par['nx']/2*par['dx']),
            par['oz']+(par['nz']/2*par['dz']),par)
fdmod.horizontal('rr',0,par)

Plot('rr',fdmod.rrplot('',par))
Plot('ss',fdmod.ssplot('',par))

# ------------------------------------------------------------
# model parameters
# ------------------------------------------------------------
Flow('zero',None,
     '''
     spike nsp=1 mag=0.0
     n1=%(nz)d o1=%(oz)g d1=%(dz)g 
     n2=%(nx)d o2=%(ox)g d2=%(dx)g |
     put label1=%(lz)s label2=%(lx)s unit1=%(uz)s unit2=%(ux)s
     ''' % par)

Flow('vp',     'zero','add add=3.00')
Flow('vs',     'zero','add add=1.50')
Flow('ro',     'zero','add add=1.00')
Flow('epsilon','zero','add add=+0.25')
Flow('delta',  'zero','add add=-0.29')
Flow('nu',     'zero','add add=0.00')

Plot('vp',fdmod.cgrey('bias=2 allpos=y',par))
Plot('vs',fdmod.cgrey('bias=1 allpos=y',par))
Plot('ro',fdmod.cgrey('bias=1 allpos=y',par))
Plot('epsilon',fdmod.cgrey('',par))
Plot('delta',fdmod.cgrey('bias=-0.29',par))

for k in (['vp','vs','ro','epsilon','delta']):
    Result(k,[k,'ss'],'Overlay')

stiff.iso2d('cI','vp','vs','ro',par)
stiff.tti2d('cA','vp','vs','ro','epsilon','delta','nu',par)

# ------------------------------------------------------------
# elastic source
# ------------------------------------------------------------
fdmod.wavelet('wav_',par['frq'],par)
Flow('ver','wav_','math output="+1*input"')
Flow('hor','wav_','math output="0*input"')
Flow('wave',['ver','hor'],
     '''
     cat axis=2 space=n ${SOURCES[1:2]} |
     transp plane=12 |
     transp plane=23 |
     transp plane=12
     ''')
fdmod.ewavelet('wave','',par)

# EWE modeling: output displacements

fdmod.ewefd2('dAu','uA','wave','cA','ro','ss','rr','ssou=n opot=n  anitype=v nbell=5',par)
pot.displacementsC('uA','uA','uAz','uAx',4,'',par)
pot.cliptogether( 'uA','uAz','uAx','"u\_z"','"u\_x"',1,'',par)
# ------------------------------------------------------------
# derivative operators
fdd.derivatives(par)
spk.delt('spk',64,64,par)
#M stationary
fdd.separatorD('dzM','dxM','spk','cA','y','x','sine',1.0,order,25,25,par)
fdd.oneirST('mop','dzM','dxM',7,7,'color=F pclip=100',par)
#--------------------------------------------------
epsilon=[0.25, 0.4, 0, 0]
delta=[-0.29, -0.1, -0.3, 0]

for ani in range(0,4,1):
    tag=str(ani)
    print epsilon[ani]

    Flow('epsilon'+tag,'zero','add add=%f'%epsilon[ani])
    Flow('delta'+tag,  'zero','add add=%f'%delta[ani])
    Flow('nu'+tag,     'zero','add add=0.00')


    stiff.tti2d('cA'+tag,'vp','vs','ro','epsilon'+tag,'delta'+tag,'nu'+tag,par)
    fdd.separatorD('dzM'+tag,'dxM'+tag,'spk','cA'+tag,'y','x','sine',1.0,order,25,25,par)
    fdd.oneirST('mop'+tag,'dzM'+tag,'dxM'+tag,25,25,'color=F',par)

    pot.potentials(   'pA'+tag,'uAz','uAx','dzM'+tag,'dxM'+tag,'y','','q',par)
    pot.cliptogether( 'pA'+tag,'pA'+tag+'p','pA'+tag+'s',"qP","qS",1,'pclip=99.5',par)

End()