## 
 # elastic modeling; wavefield separation
 ## Constant VTI media
# study operator trucation

from rsf.proj import *
import sys
import fdmod,fdd,spk,stiff,pot


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']+'''
titlesz=15  labelsz=9
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.tti2d('cI','vp','vs','ro','zero','zero','zero',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)

# ------------------------------------------------------------
# modeling
# ------------------------------------------------------------
# EWE modeling: output displacements for isotroic and VTI models
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',9,'pclip=98',2,par)
fdmod.emovie('uAmovie','uA',9,'pclip=98',2,par)

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

fdd.derivatives(par)
spk.delt('spk',64,64,par)


fdd.separatorD('dzI','dxI','spk','cI','y','x','sine',1.0,8,4,1,par)
#M stationary
fdd.separatorD('dzM','dxM','spk','cA','y','x','sine',1.0,8,25,25,par)
#A non-sationary
fdd.separatorD('dzA','dxA','spk','cA','n','x','sine',1.0,8,25,25,par)




# impulse responses
fdd.oneirST('cop','dzC','dxC',7,7,'color=e pclip=100',par)
fdd.oneirST('iop','dzI','dxI',7,7,'color=e pclip=100',par)
fdd.oneirST('mop','dzM','dxM',7,7,'color=e pclip=100',par)

# ------------------------------------------------------------
# isotropic displacements/potentials
pot.displacementsC('uI','uI','uIz','uIx',4,'',par)
pot.cliptogether( 'uI','uIz','uIx','"u\_z"','"u\_x"',1,'',par)

pot.potentials(   'qI','uIz','uIx','dzC','dxC','y','','',par)
pot.cliptogether( 'qI','qIp','qIs',"P","S",1,'',par)

pot.potentials(   'pI','uIz','uIx','dzI','dxI','y','','',par)
pot.cliptogether( 'pI','pIp','pIs',"P","S",1,'',par)
# ------------------------------------------------------------
# anisotropic displacements/potentials
pot.displacementsC('uA','uA','uAz','uAx',4,'',par)
pot.cliptogether( 'uA','uAz','uAx','"u\_z"','"u\_x"',1,'',par)

pot.potentials(   'qA','uAz','uAx','dzC','dxC','y','','q',par)
pot.cliptogether( 'qA','qAp','qAs',"qP","qS",1,'',par)

pot.potentials(   'pA','uAz','uAx','dzM','dxM','y','','q',par)
pot.cliptogether( 'pA','pAp','pAs',"qP","qS",1,'',par)


fdmod.box('box1',0.41,0.8,0.41,0.8, par)
Plot('box1',fdmod.qqplot('screenht=10.25 plotcol=1',par))
fdmod.box('box2',0.21,1.0,0.21,1.0, par)
Plot('box2',fdmod.qqplot('screenht=10.25 plotcol=3',par))
fdmod.box('box3',0.01,1.19,0.01,1.19, par)
Plot('box3',fdmod.qqplot('screenht=10.25 plotcol=6',par))
Plot('box',['box1','box2','box3'],'Overlay')

# ------------------------------------------------------------
# different operator sizes
for i in range(1,6,2):
    n=i*5
    tag=str(i)

    fdd.separatorD('dzM'+tag,'dxM'+tag,'spk','cA','y','x','sine',1.0,8,n,n,par)

    fdd.oneirST_box('mop'+tag,'dzM'+tag,'dxM'+tag,25,25,'box',
                    '''pclip=99 format1=%3.0f format2=%3.0f''',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=100',par)
End()