from rsf.proj import *

par = dict(dx=0.005,dz=0.005,
           minx=35,maxx=45,
           minz=0,maxz=6,
           f=15,nt=6000,dt=0.0005,
           sx=40,sz=0.01,rz=0.00625,
           jsnap=100,jdata=10)

par['kt'] = 1.2/(par['dt']*par['f'])
par['nz'] = par['maxz']/par['dz']
par['nx'] = (par['maxx']-par['minx'])/par['dx']
velsegy = 'vel_z6.25m_x12.5m_exact.segy'
densegy = 'density_z6.25m_x12.5m.segy'

Fetch([velsegy+'.gz',densegy+'.gz'],dir='2004_BP_Vel_Benchmark',server='ftp://software.seg.org',top='pub/datasets/2D')

Flow('BPvelocity.segy',velsegy+'.gz','gzip -d')
Flow('BPvelocity','BPvelocity.segy','segyread | put d1=0.00625 label1="Z" unit1="km" d2=0.0125 label2="X" unit2="km"')
Flow('BPdensity.segy',densegy+'.gz','gzip -d')
Flow('BPdensity','BPdensity.segy','segyread | put d1=0.00625 label1="Z" unit1="km" d2=0.0125 label2="X" unit2="km" ' %par)
Plot('BPvelocity','grey color=j mean=y scalebar=y title="BP Velocity Model"')
Flow('velocity','BPvelocity','window max1=%(maxz)d max2=%(maxx)g min2=%(minx)g | add scale=0.001 |transp | sinc n1=%(nx)d d1=%(dx)g o1=%(minx)g | transp' %par)
Flow('density','BPdensity','window max1=%(maxz)d max2=%(maxx)g min2=%(minx)g | transp | sinc n1=%(nx)d o1=%(minx)g d1=%(dx)g | transp'%par)

Flow('wav',None,'spike n1=%(nt)d d1=%(dt)g o1=0 k1=%(kt)d | ricker1 frequency=%(f)d ' %par )
Flow('acousticwavelet','wav','transp plane=12 | pad n1=2')

Flow('sx',None,'math n1=1 output=%(sx)g' %par)
Flow('sz',None,'math n1=1 output=%(sz)g' %par)
Flow('source','sx sz','cat axis=1 ${SOURCES[1]}')
Flow('rx',None,'math n1=%(nx)d d1=%(dx)g o1=%(minx)g output=x1' %par)
Flow('rz',None,'math n1=%(nx)d d1=%(dx)g o1=%(minx)g output=%(rz)g' %par)
Flow('receivers','rx rz','cat axis=2 ${SOURCES[1]} | transp')
Plot('velocity','grey color=j mean=y scalebar=y title="Velocity Model subsection"')
Flow('arecfield awavefield','acousticwavelet density receivers source velocity','''
awefd2d
den=${SOURCES[1]}
rec=${SOURCES[2]}
sou=${SOURCES[3]}
vel=${SOURCES[4]}
wfl=${TARGETS[1]}
verb=y dabc=y snap=y jsnap=%(jsnap)d jdata=%(jdata)d free=y
''' %par)
Plot('awavefield','grey gainpanel=a title="Acoustic wave propagation"',view=1)
Result('awavefield','window n3=1 f3=40 | grey title="Acoustic wave propagation"')
Result('arecfield','transp | grey gainpanel=a title="Acoustic receiver data"')

# . . Elastic modelling
Flow('elasticwavelet','wav','put n2=1 n3=1 | transp plane=13 | pad n2=2')
# Explicit about c13 c33
Flow('c11','velocity density','math vp=${SOURCES[0]} den=${SOURCES[1]} output="den*vp^2" | smooth rect1=4 rect2=4 repeat=3')
Flow('Vs','velocity','math vp=${SOURCES[0]} output="vp/sqrt(3)" ')
Flow('c55','Vs density','math vs=${SOURCES[0]} den=${SOURCES[1]} output="den*vs^2" | smooth rect1=4 rect2=4 repeat=3')
Flow('c33','c11','cp')
Flow('c13','c11 c55','math c11=${SOURCES[0]} c55=${SOURCES[1]} output="c11-2*c55" ')
Flow('ccc','c11 c33 c55 c13','cat axis=3 ${SOURCES[1:4]} ')
Flow('erecfield ewavefield','elasticwavelet density receivers source ccc',
'''
ewefd2d
den=${SOURCES[1]}
rec=${SOURCES[2]}
sou=${SOURCES[3]}
ccc=${SOURCES[4]}
wfl=${TARGETS[1]}
dabc=y snap=y verb=y jsnap=%(jsnap)d jdata=%(jdata)d
ssou=y nb=20 nbell=11
''' %par)

Flow('ewfldZ','ewavefield','window n3=1 f3=0 min2=%(minx)g max2=%(maxx)g max1=%(maxz)g min1=%(minz)g' % par)
Flow('ewfldX','ewavefield','window n3=1 f3=1 min2=%(minx)g max2=%(maxx)g max1=%(maxz)g min1=%(minz)g' % par)

Plot('ewfldZ','window n3=1 f3=40 | grey max1=6 pclip=99.5 title="Elastic wave propagation (Z component)"')
Plot('ewfldX','window n3=1 f3=40 | grey max1=6 pclip=99.5 title="Elastic wave propagation (X component)"')
Result('ewfld','ewfldZ ewfldX','SideBySideAniso')

Flow('erfldZ','erecfield','window n2=1 f2=0')
Flow('erfldX','erecfield','window n2=1 f2=1')

Plot('erfldZ','grey transp=n title="Elastic receiver data (Z component)" ')
Plot('erfldX','grey transp=n title="Elastic receiver data (X component)" ')
Result('erfld','erfldZ erfldX','SideBySideAniso')

End()