# 
 # elastic modeling; wavefield separation
 ## VTI sigsbee model

from rsf.proj import *
import sys
#sys.path.append('Python')
import fdmod,stiff,pot,fdd,spk,sigs,model

par=sigs.param()             # init model parameters
sigs.fdmpar(par)             # add FDM parameters
sigs.rtmpar(par)             # add RTM parameters
sigs.velocity('vall',par)    # read velocity
sigs.win3('vwin','vall',par) # select model window
fdmod.param(par)


fdmod.param(par)

par['frq']=15
par['nt']=4001
par['dt']=0.0005
par['jsnap']=100

par['xs']=0.5
par['ys']=0.5
par['xc']=-12

par['labelattr']=par['labelattr']+'''
titlesz=15 labelsz=8
n2tic=8 o2num=2 d2num=1
wantaxis=y
'''



# ------------------------------------------------------------
# source/receiver coordinates
# ------------------------------------------------------------
fdmod.point('ss',
            par['ox']+(par['nx']/2*par['dx']),
            par['oz']+(par['nz']/2*par['dz']),par)
fdmod.horizontal('rr',par['oz'],par)

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

model.inimodel('ss','rr',par)

stiff.iso2d('cI','vp','vs','ro',par)
stiff.tti2d('cA','vp','vs','ro','epsilon','delta','nu',par)
Flow('zero','vp','math output=0 ')

# ------------------------------------------------------------
# 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)

# ------------------------------------------------------------
# modeling
# ------------------------------------------------------------
# EWE modeling: output displacements
fdmod.ewefd2('dIu','uI','wave','cI','ro','ss','rr','ssou=n opot=n nbell=5 anitype=i',par)
fdmod.ewefd2('dAu','uA','wave','cA','ro','ss','rr','ssou=n opot=n nbell=5 anitype=v',par)

fdmod.emovie('uImovie','uI',20,'pclip=99',2,par,par['xs'],par['ys'],-12)
fdmod.emovie('uAmovie','uA',20,'pclip=99',2,par,par['xs'],par['ys'],-12)

# ------------------------------------------------------------
# derivative operators

fdd.derivatives(par)
spk.delt('spk',64,64,par)
order=8

fdd.separatorD('dzI','dxI','spk','cI','y','x','sine',1.0,order,4,1,par)
#M stationary
Flow('ep','epsilon','math output=.3')
Flow('de','epsilon','math output=.1')
stiff.tti2d('cM','vp','vs','ro','ep','de','nu',par)
fdd.separatorD('dzM','dxM','spk','cM','y','x','sine',1.0,order,25,25,par)
#A non-sationary
fdd.separatorD('dzA','dxA','spk','cA','n','x','sine',1.0,order,25,25,par)



# impulse responses
fdd.oneirST('cop','dzC','dxC',7,7,'',par)
fdd.oneirST('iop','dzI','dxI',7,7,'',par)
fdd.oneirST('mop','dzM','dxM',7,7,'',par)
fdd.oneirNS('aop','dzA','dxA',7,7,'',par)

oframe=24
nframe=28
jframe=4
#print range(oframe,nframe,jframe)
for iframe in range(oframe,nframe,jframe):
    tag = "-f%02d"%iframe
#    tag=''
    # ------------------------------------------------------------
    # anisotropic displacements/potentials
    pot.displacementsC('uA'+tag,'uA','uAz'+tag,'uAx'+tag,iframe,'pclip=99.5',par)
    pot.potentials(    'qA'+tag,     'uAz'+tag,'uAx'+tag,'dzC','dxC','y','pclip=99.','q',par)
    pot.potentials(    'pA'+tag,     'uAz'+tag,'uAx'+tag,'dzA','dxA','n','pclip=99.','q',par)
    pot.potentials(    'mA'+tag,     'uAz'+tag,'uAx'+tag,'dzM','dxM','y','pclip=99.','q',par)
##############################################################
# wavefield on model
    model.wfom('uA'+tag+'-wom','uAz'+tag,'uAx'+tag,'vp',3.15,'"u\_z"','"u\_x"',1,'',par)
    model.wfom('qA'+tag+'-wom','qA'+tag+'p','qA'+tag+'s','vp',3.15,'qP'    ,'qS'    ,1,'',par)
    model.wfom('pA'+tag+'-wom','pA'+tag+'p','pA'+tag+'s','vp',3.15,'qP'    ,'qS'    ,1,'',par)
    model.wfom('mA'+tag+'-wom','mA'+tag+'p','mA'+tag+'s','vp',3.15,'qP'    ,'qS'    ,1,'',par)


End()