Madagascar - School of Earth Sciences - The Ohio State University

plainspecialSoftware and s/w Development

Dec 14, 2013 (3 years and 8 months ago)

84 views

2D and 3D Seismic Modeling with Madagascar


b
y


Kyle Shalek and Jeff Daniels


The Ohio State University



Preface and Motivation


We needed a 3D seismic modeling package at OSU to compliment our 3D EM
modeling software.

Madagascar is a freeware package tha
t is readily available. The
first challenge was to make the package usable and accessible. Having done this, we
decided to post the instructions for others use. The following is a very terse outline
of a “how to run” Madagascar models.





















Table of Contents


1.)
Introduction

2.)
Xwin


2.1) Installation


2.2) Basic Unix Commands

3.) Running Programs

4.) Programs


4.1) File Operations


4.2) Processing


4.3) Imaging

5.)
Modeling


5.1) awefd2d: 2D Acoustic

time
-
domain

Finite Difference Modeling



5.1.1) Input



5.1.2) Output



5.1.3) Example


5.2) awefd3d: 3D Acoustic

time
-
domain

Finite Difference Modeling



5.
2
.1) Input



5.
2
.2) Output



5.
2
.3) Example


5.3) ewefd2d: 2D Elastic

time
-
domain

Finite Difference Modeling



5.
3
.1) Input



5.
3
.2) Outp
ut



5.
3
.3) Example


5.4) ewefd3d: 3D Elastic

time
-
domain

Finite Difference Modeling



5.
4
.1) Input



5.
4
.2) Output



5.
4
.3) Example


1.) Introduction



Open
-
source software package for data processing and numerical experiments



Large collection of programs (
~300) written in C and C++



Command line environment
.
There is n
o GUI available for modeling. There is a GUI, OpendTect,
available for processing flows

but needs a seismic survey data set to begin
.




Links:


Madagascar Main Page:
http://www.reproducibility.org/wiki/Main_Page


Madagascar Download:
http://www.reproducibility.org/wiki/Download


List of Programs:
http://www.reproducibility.org/RSF/



Currently,
our version of
Madagascar
is stored on a private cluster

that requires Linux or Xwin access
.



Madagascar can be installed on a Windows box utilizing Cygwin.



The examples are best und
erstood by looking at and executing the code yourself.



2.)
Xwin

(X
-
Win32 V9.4
)


2.1) Installation



Ohio State Faculty/Students can download for free at
:
http://osusls.osu.edu/



Follow the on
-
screen instructions. The

license is floating and located in a txt file in the
downloaded folder.




2.2)
Basic Unix Commands



ls: lists all of the files in the current directory



cd

directory_name: changes to the indicated directory



cd .. : goes up one directory



mkdir

directory_nam
e: creates the indicated directory



gedit
: opens a file for editing



cp home2/shalek/test . : copies the directory test to the current directory



ln

s home2/shalek/test . : creates a symbolic link to the test directory in the current directory


3
.) Running
Programs:


Programs can be executed from the command line or the program commands can be piped together
using
SConstruct files (Python scripts)
. The SConstruct file is saved in the dir
ectory you want to use
for output
files. T
he execution commands are li
sted below.



s
cons: executes the SConstruct from command line



scons view: shows the produced images



xtpen image_file: displays a single image of your choice



gedit SConstruct: opens the SConstruct file for editing


4.) Programs


A complete listing of ava
ila
ble programs can be found at:
http://www.reproducibility.org/RSF/

. Find the
program you want to use and the link will take you to a description of what the program does and what
arguments
it accepts. E
xample SConstruct files

and the source code
are

also available from that page.

Useful
and important programs are listed below.


4.
1
)
File Operations
:



< input_file.rsf
: command line format



<
output_file.rsf
: command line format



sfin < file.rsf: displays ba
sic information about rsf files



sfattr < file.rsf: Display dataset attributes



sfdd

< file.rsf col=1 format=data_type
: converts the data file into

any type




sfdisfil < file.rsf col=1: prints the file to the command window with indicated number of
columns



sf
segyread: Convert a SEG
-
Y or SU dataset to RSF



sfsegy2rsf: Convert a SEG
-
Y dataset to RSF



sfsegywrite: Convert an RSF dataset to SEGY or SU

4.2) Processing



sfmath: mathematical operations on data files



sfhwt2d: Computes rays and wavefronts/2
-
D Huygens wave
front tracing traveltimes



sfwindow: windows the data set



sftime2depth: Time
-
to
-
depth conversion in V(z)



sfdepth2time: Conversion from depth to time in a V(z) medium



sfvscan: Velocity analysis



sfnmo: Normal moveout



sffinstack: DMO and stack by finite
-
differ
ence offset continuation



sfagc: automatic gain control



sfbandpass: Bandpass filtering



sfboxsmooth: Multi
-
dimensional smoothing with boxes



sfsmooth: Multi
-
dimensional triangle smoothing.



sfclip: clip the data



sffft1: Fast Fourier Transform along the first
axis



sffft1: FFT transform on extra axis



sfnoise: Add random noise



sfswtdenoise: De
-
noising using stationary wavelet transform



sfspectra: Frequency spectra



sfspectra2: Frequency spectra in 2D



sftransp: Transpose two axes in a dataset



sfstack: Stack a datas
et over one of the dimensions



Sfmig45: Finite
-
difference modeling/migration: 15
-

and 45
-
degree approximation



Sfmigsteep: 3
-
D Kirchhoff time migration for antialiased steep dips.



Sfkirmod: Kirchhoff 2.5
-
D modeling with analytical Green's functions



Sfkirmod3
:

Kirchhoff 3
-
D modeling with analytical Green's functions.



Sfafmod: Time
-
domain acoustic modeling


4.3) Imaging



sfwiggle: plot data with wiggly traces



sfplotrays: Plot rays



sfplotrays3: Plot rays in 3D with OpenGL



sfgraph: 1D graph plot



sfgrey: 2D graysca
le raster plots



sfgrey3: cube plots



sfcontour: contour plots



sfcontour3: 3D contour plots

5.) Modeling


5.1)
awefd2d:

2D Acoustic

time
-
domain

Finite Difference Modeling



5.1.1) Input
: See

the example code below for how the input and output is used.



vel => Vp



den => density



sou => source coordinates



rec => receiver coordinates



wfl => wavelet



5
.1.2) Output



data => traces



wave



5.1.3) Example



############################################
##############

# Created by Kyle Shalek #

# Last Modified: 9
-
4
-
09






#

# Creates a 2D model with velocity and density inputs. #

# Acoustic FD modeling with awefd2d (Psava) #

# Single trace output






#

####
######################################################


from rsfproj import *

import string

########################################################

#1) create velocity model

# Depth in km

xmax = 4.0

zmax = 2.0


layers = ((0.1,0.1,0.1,0.1,0.1),



(1.1,1.1
,1.1,1.1,1.1),



(1.3,1.3,1.3,1.3,1.3))


layer1=0.1

layer2=1.1
-
0.1

layer3=1.3
-
1.1

layer4=2
-
1.3


# Velocity (Vp) in km/s

vp = (0.0,


6,


4.6657,


6)

#Velocity (Vs) in km/s

vs = (0.0,


2,


4,


2)


# Density in g/cc, converted t
o kg/ckm

densities = (0.1*1000000000000,



2.8*1000000000000,



2.5*1000000000000,



2.8*1000000000000)


def arr2str(array,sep=' '):


return string.join(map(str,array),sep)


n1 = len(layers[0])

n2 = len(layers)


Flow('layers',None,


'''


echo %s


n1=%d n2=%d o1=0 d1=%g


data_format=ascii_float in=$TARGET


''' % (string.join(map(arr2str,layers),' '),


n1,n2,xmax/(n1
-
1)))


d = 0.01 # sampling steps


Flow('mod1','layers',


'''


spline o1=0 d1=%g n1=%d
|


unif2 d1=%g n1=%d v00=%s |


dd form=native


''' % (d,int(1.5+xmax/d),


d,int(1.5+zmax/d),


arr2str(vp,','),))

Flow('rho','layers',


'''


spline o1=0 d1=%g n1=%d |


unif2 d1=%g n1=%d v00=%s |


dd form=nat
ive


''' % (d,int(1.5+xmax/d),


d,int(1.5+zmax/d),


arr2str(densities,','),))


Result('mod1',


'''


grey color=j title="Velocity Model"


allpos=y titlesz=8 labelsz=6 screenratio=0.5


scalebar=y barlabel='v
elocity (km/s)' barlabelsz=6


label1="Depth (km)"


label2="Distance (km)"


''')


# Contour the layers

Plot('modline','mod1','contour title=" " label1='' label2='' wantaxis=n screenratio=0.5')

############################################
##################


#2) 2D Acoustic wave FD Modeling

# Source location

sx = 2

sz = 0.1

Flow('source',None,

'spike n1=3 nsp=3 k1=1,2,3 mag=%g,%g,1 o1=0 o2=0' % (sx,sz))

# Receiver location

rx=2

rz=0.1

Flow ('receiver',None,

'spike n1=3 nsp=3 k1=1,2,3 mag=%g
,%g,1 o1=0 o2=0' % (rx,rz))


########################

time=0.5

timesteps=1500 # timesteps*0.001=sec

frequency=20


##########################



#2)a) Source wavelet

Flow('wavelet',None,


'''


spike nsp=1 n1=%g d1=%g k1=%d |


ricker1 frequency=%
g |


transp


''' % (timesteps,0.001,200,frequency))


#2)b) Awefd modeling

Flow('data wave','wavelet mod1 source receiver rho',


'''


awefd2d verb=y free=y expl=y snap=y jsnap=10 dabc=y


db=0 dbx=0 o1=0 o2=0


vel=${SOURCES[1]}


sou=${SOURCES[2]}


rec=${SOURCES[3]}


den=${SOURCES[4]}


wfl=${TARGETS[1]}


''')


#################################################################


#2)c) Movie of wave

Plot('wave',


'''


grey gainpanel=all title=Wave label
1=Depth unit1=km label2=Lateral unit2=km


''',view=1)

Result('wave1','wave',


'''


grey color=j scalebar=y barlabel='Amplitude' barlabelsz=6


gainpanel=all title=Wave


label1=Depth unit1=km label2=Lateral unit2=km screenratio=0.5


''')


# Wavefield s
napshot at 'time' value

Plot('snap','wave',


'''


window n1=201 n2=401 n3=1 min3=%g |


grey color=j gainpanel=all title="Wave Snapshot"


label1=Depth unit1=km label2=Lateral unit2=km screenratio=0.5


scalebar=y barlabel='Amplitude' barl
abelsz=6


'''%(time))

##################################################################


#2)d) Overlay wave and layer contours

Result('snap1',['snap','modline'],'Overlay')

#2)e) Single Trace

Result('trace','data',


'''


transp |


spline o1=0 d1=0.001
n1=%g |


wiggle label1=Time unit1=s label2=Amplitude unit2='' title='Vertical Trace with Source'


'''% (timesteps))

Result('trace_window','data',


'''


transp |


spline o1=0 d1=0.001 n1=%g |


window min1=0.3 |


wiggle label1=Time unit1=s label2=Amplitude
unit2='' title='Vertical Trace without Source'


'''% (timesteps))

End()

################################################################################





5.2)

awefd3d:

3D Acoustic

ti
me
-
domain

Finite Difference Modeling



5.
2
.1) Input
: See the example code below for how the input and output is used.



vel

=> Vp



den => density



sou => source coordinates



rec => receiver coordinates



wfl => wavelet




5.
2
.2) Output



data => traces



wave


5.
2
.
3) Example
: This code uses the python files pplot.py and fdmod.py which are located in
/book/packages directory. The files have commonly used modeling routines piped together
allowing the SConstruct files

to be

cleaner and easier to read.






##########
##############################################################

from rsfproj import *

import fdmod


#
------------------------------------------------------------

par = {


'nx':151, 'ox':0, 'dx':0.002, 'lx':'x', 'ux':'km',


'ny':151, 'oy':0, 'dy':0.0
02, 'ly':'y', 'uy':'km',


'nz':151, 'oz':0, 'dz':0.002, 'lz':'z', 'uz':'km',


'nt':601, 'ot':0, 'dt':0.0002, 'lt':'t', 'ut':'s',


'kt':100,


'jsnap':100,


'height':10,


'nb':20,


'frq':45,


'nbell':5


}

fdmod.param(par)

par['nf
rame']=5

par['iframe']=4

par['dabc']='y'


#
------------------------------------------------------------

# Thomsen parameters

par['vp']=2.0

par['vs']=1.0

par['ro']=2.0


#
------------------------------------------------------------

par['kz']=2./3.*par['nz'
]


#
------------------------------------------------------------

fdmod.wavelet('wav_',par['frq'],par)

#
------------------------------------------------------------

# acoustic source

Flow( 'wava','wav_','transp')

Result('wava','transp |' + fdmod.waveplot
('',par))


#
------------------------------------------------------------

# source/receiver coordinates

#
------------------------------------------------------------


xsou=par['ox']+(par['nx']/2*par['dx']);

ysou=par['oy']+(par['ny']/2*par['dy']);

zsou=par
['oz']+(par['nz']/2*par['dz']);


center=fdmod.center3d(1.1*xsou,1.1*ysou,1.1*zsou,par)


fdmod.point3d('ss
-
3d',xsou,ysou,zsou,par)

fdmod.horizontal3d('rr
-
3d',0,par)


#
------------------------------------------------------------

Flow('zero
-
3d',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


n3=%(ny)d o3=%(oy)g d3=%(dy)g |


put label1=%(lz)s label2=%(lx)s label3=%(ly)s unit1=%(uz)s unit2=%(ux)s unit3=%(uy)s


''' % par)


#
-----------
-------------------------------------------------


Flow('vp
-
3d','zero
-
3d','math output="%(vp)g"' %par)

Flow('vs
-
3d','zero
-
3d','math output="%(vs)g"' %par)

Flow('ro
-
3d','zero
-
3d','math output="%(ro)g"' %par)

Result('vp
-
3d',fdmod.cgrey3d('allpos=y bias=2.0'+
center,par))

Result('vs
-
3d',fdmod.cgrey3d('allpos=y bias=1.0'+center,par))

Result('ro
-
3d',fdmod.cgrey3d('allpos=y bias=2.0'+center,par))


#
-------------------------------------------------------------

# 3D acoustic modeling

fdmod.awefd3d('da
-
3d','wa
-
3d','
wava','vp
-
3d','ro
-
3d','ss
-
3d','rr
-
3d','',par)

Result('wa
-
3d','window n4=1 f4=%(iframe)d |'%par


+ fdmod.cgrey3d('pclip=99.9'+center,par))

Result('da
-
3d',


'''


put


n1=%(nx)d o1=%(ox)g d1=%(dx)g


n2=%(ny)d o2=%(oy)g d2=%(dy)g


n3=%(nt)d o3=%(ot)g d3=%(dt)g |


transp plane=23 |


transp plane=12 |


''' % par


+ fdmod.dgrey3d('pclip=99.9'+center+' frame1=%d' % (0.65*par['nt']) ,par))


End()





5.3)

ewefd2d
: 2D Elastic

time
-
domain

Finite Difference Modeling



5.
3
.1) Input
: See the example code below for how the input and output is used.



ccc

=> Vp, Vs, den used to calculate stiffness



den

=> density



sou

=> source coordinates



rec

=> rece
iver coordinates



wfl

=> wavelet



5.
3
.2) Output



data

=> traces



wave




5.
3
.3) Example
: This code uses the python files pplot.py, stiffness, and fdmod.py which
are located in /book/packages directory. The files have commonly used modeling routines piped
t
ogether allowing the SConstruct files to be cleaner and easier to read.

###############################################################

# Created by Kyle Shalek

# Oct. 2009

# 2D elastice finite difference modeling

##########################################
#####################

from rsfproj import *

import fdmod,stiffness


#
------------------------------------------------------------

par = {


'nx':151, 'ox':0, 'dx':0.002, 'lx':'x', 'ux':'km',


'ny':151, 'oy':0, 'dy':0.002, 'ly':'y', 'uy':'km',


'
nz':151, 'oz':0, 'dz':0.002, 'lz':'z', 'uz':'km',


'nt':2001, 'ot':0, 'dt':0.0002, 'lt':'t', 'ut':'s',


'kt':100,


'jsnap':100,


'height':10,


'nb':20,


'frq':45,


'nbell':5


}

fdmod.param(par)

par['nframe']=5

par['iframe']=4

par['
dabc']='y'

#
-------------------------------------------------------------

def arr2str(array,sep=' '):


return string.join(map(str,array),sep)

#
-------------------------------------------------------------

# Create Layers

layers = ((0.01,0.01,0.01,0.01,0
.01),



(0.15,0.15,0.15,0.15,0.15),



(0.25,0.25,0.25,0.25,0.25))


n1 = len(layers[0])

n2 = len(layers)


layers1 = string.join(map(arr2str,layers),' ')


Flow('layers',None,


'''


echo %s


n1=%g n2=%g o1=%g d1=%g


data_format=ascii_float i
n=$TARGET


''' % (layers1,n1,n2,par['ox'],par['dx'] ))


#
------------------------------------------------------------

# Thomsen parameters


# Velocity (Vp) in km/s

vp =( 2,


2,


6,


2)

vp1= arr2str(vp,',')


#Velocity (Vs) in km/s

vs = (1,



1,


4,


1)

vs1=arr2str(vs,',')


# Density in g/cc, converted to kg/ckm

rho = (2*1000000000000,


2*1000000000000,


6*1000000000000,


2*1000000000000)

rho1=arr2str(rho,',')


#
----------------------------------------------------
-------

Flow('vp
-
2d','layers',


'''


spline o1=%g d1=%g n1=%d |


unif2 o1=%g d1=%g n1=%d v00=%s |


dd form=native |


put label1=%s label2=%s unit1=%s unit2=%s


''' % (par['ox'],par['dx'],par['nx'],



par['oz'],par['dz'],par['nz']
,vp1,



par['lz'],par['lx'],par['uz'],par['ux']))


Flow('vs
-
2d','layers',


'''


spline o1=%g d1=%g n1=%d |


unif2 o1=%g d1=%g n1=%d v00=%s |


dd form=native |


put label1=%s label2=%s unit1=%s unit2=%s


''' % (par['ox'],par['dx']
,par['nx'],



par['oz'],par['dz'],par['nz'],vs1,



par['lz'],par['lx'],par['uz'],par['ux']))


Flow('ro
-
2d','layers',


'''


spline o1=%g d1=%g n1=%d |


unif2 o1=%g d1=%g n1=%d v00=%s |


dd form=native |


put label1=%s label2=%s uni
t1=%s unit2=%s


''' % (par['ox'],par['dx'],par['nx'],



par['oz'],par['dz'],par['nz'],rho1,



par['lz'],par['lx'],par['uz'],par['ux']))


#
------------------------------------------------------------

par['kz']=2./3.*par['nz']


#
-----------------
-------------------------------------------

fdmod.wavelet('wav_',par['frq'],par)


#
------------------------------------------------------------

# elastic source

#
------------------------------------------------------------

Flow('souz','wav_','math output
=input*1')

Flow('soux','wav_','math output=input*0')


Flow('wave
-
2d',['souz','soux'],


'''


cat axis=2 space=n ${SOURCES[1:2]} |


transp plane=12 |


transp plane=23 |


transp plane=12


''')

fdmod.ewavelet('wave
-
2d','',par)


#
------
------------------------------------------------------

# source/receiver coordinates

#
------------------------------------------------------------


xsou=par['ox']+(par['nx']/2*par['dx']);

#zsou=par['oz']+(par['nx']/8*par['dx']);

zsou=0.03


fdmod.point('ss
-
2d',xsou,zsou,par)

fdmod.horizontal('rr
-
2d',zsou,par)

#fdmod.point('rr
-
2d',xsou,zsou,par)

Plot('rr
-
2d',fdmod.rrplot('title="Receiver Line"',par))

Plot('ss
-
2d',fdmod.ssplot('title="Source"',par))


#
---------------------------------------------------------
---

Plot('vp
-
2d',fdmod.cgrey('allpos=y bias=2.0 title="Vp"',par))

Plot('vs
-
2d',fdmod.cgrey('allpos=y bias=1.0 title="Vs"',par))

Plot('ro
-
2d',fdmod.cgrey('allpos=y bias=2.0 title="Rho"',par))

Result('vp
-
2d
-
s','vp
-
2d ss
-
2d','Overlay')

Result('vp
-
2d
-
r','vp
-
2d

rr
-
2d','Overlay')

#Result('vs
-
2d','vs
-
2d ss
-
2d','Overlay')

#Result('ro
-
2d','ro
-
2d ss
-
2d','Overlay')


stiffness.iso2d('Ic
-
2d','vp
-
2d','vs
-
2d','ro
-
2d',par)


#
------------------------------------------------------------

# 2D elastic modeling

stiffness.cplot
2d('Ic
-
2d',1,1,par)


fdmod.ewefd2d('Id
-
2d','Iw
-
2d','wave
-
2d','Ic
-
2d','ro
-
2d','ss
-
2d','rr
-
2d','ssou=n opot=n',par)

fdmod.emovie( 'Iw
-
2d
-
movie','Iw
-
2d',par['nframe'],'pclip=99.9 color=j',2,par)


Result('Id
-
2d','window n2=1 | transp |' + fdmod.dgrey('title="M
odel Response" ',par))



"""

#
--------------------------------------------------------------

# Ewefd modeling

Flow('Id
-
2d Iw
-
2d','wave
-
2d Ic
-
2d ro
-
2d ss
-
2d rr
-
2d ',


'''


ewefd2d nbell=5 jsnap=10 dabc=y ssou=n opot=n nb=20 ompchunk=1 ompnth=0 verb=
y free=n snap=y


ccc=${SOURCES[1]}


den=${SOURCES[2]}


sou=${SOURCES[3]}


rec=${SOURCES[4]}


wfl=${TARGETS[1]}


''')

"""

# Movie of wave

Plot('Iw
-
2d',


'''


window n3=1 |


grey gainpanel=all title=Wave labe
l1=Depth unit1=km label2=Lateral unit2=km


''',view=1)

Result('Iw
-
2d
-
1','Iw
-
2d',


'''


window n3=1 |


grey color=j scalebar=y barlabel='Amplitude' barlabelsz=6


gainpanel=all title=Wave


label1=Depth unit1=km label2=Lateral unit2=km


''')


Result('I
d
-
2d
-
1','Id
-
2d',


'''


window n2=1 |


transp |



grey color=j


gainpanel=all title=data


label1=time unit1=s label2=Lateral unit2=km


''')


"""


Flow('Id
-
2d
-
1','Id
-
2d','dd form=native line=2 format=float')



Result('Id
-
2d
-
trace','Id
-
2d
-
1',


'''


spline
o1=0 d1=0.0002 n1=2001 |


spline o2=0 d2=0.0002 n2=2001 |


wiggle label1=Time unit1=s label2=Amplitude unit2='' title='Vertical Trace with Source'


''')

"""


End()

#############################################################################


Wavelet


V
elocity Model



Stiffness Diagram
















A B


Wave
-
field
s with different polarizations




Trace from the ‘A’ wave
-
field shown above


‘A’ trace zoomed in



Trace from the ‘B’ wave
-
field shown ab
ove


Model Response

cross
-
section

from a
receiver line




5.4)

ewefd3d:

3D Elastic

time
-
domain

Finite Difference Modeling



5.
4
.1) Input
: See the example code below for how the input and
output is used.



ccc => Vp, Vs, den used to calculate stiffness



den



sou



rec



wfl




5.
4
.2) Output



data => three lines of output



wave




5.
4
.3) Example
: This code uses the python files pplot.py, stiffness, and fdmod.py which
are located in /book/packages di
rectory. The files have commonly used modeling routines piped
together allowing the SConstruct files to be cleaner and easier to read.

##############################################################

# Created by Kyle Shalek

# Nov. 2009

# 3d elastic finite
difference modeling

##############################################################

from rsfproj import *

import fdmod,stiffness


#
------------------------------------------------------------

par = {


'nx':151, 'ox':0, 'dx':0.002, 'lx':'x', 'ux':'km',


'ny':151, 'oy':0, 'dy':0.002, 'ly':'y', 'uy':'km',


'nz':151, 'oz':0, 'dz':0.002, 'lz':'z', 'uz':'km',


'nt':601, 'ot':0, 'dt':0.0002, 'lt':'t', 'ut':'s',


'kt':100,


'jsnap':100,


'height':10,


'nb':20,


'frq':45,


'nbell':5


}

fdmod.param(par)

par['nframe']=5

par['iframe']=4

par['dabc']='y'


#
------------------------------------------------------------

#3D model: layers

def layer(output):


return 'math n1=%g d1=%g n2=%g d2=%g output="%s" ' % (par['nx'],par['dx'],par['n
y'],par['dy'],output)


Flow('one',None,layer('0.01'))

Flow('two',None,layer('0.15'))

Flow('three',None,layer('0.25'))

Flow('struct','one two three','cat ${SOURCES[1:3]}')


# Parameters
-------------------------------------------------

Flow('vp
-
3d','struct'
,'unif3 n1=%g d1=%g v00=2,2,6,2' %(par['nz'],par['dz'] ))

Result('vp
-
3d','vp
-
3d','''


byte allpos=y |


grey3 flat=n color=j wanttitle=n


label1="Depth z (km)" label2="X (km)" label3="Y (km)"


''')

Flow('vs
-
3d','struct','unif3 n1=%g d1=%g v00=1,1,4,1' %(par
['nz'],par['dz'] ))

#Result('vs
-
3d','vs
-
3d','''

#

byte allpos=y |

#

grey3 flat=n color=j wanttitle=n

#

label1="Depth z (km)" label2="X (km)" label3="Y (km)"

#

''')

Flow('ro
-
3d','struct','unif3 n1=%g d1=%g v00=2,2,6,2' %(par['nz'],par['dz'] ))

#Result('ro
-
3
d','ro
-
3d','''

#

byte allpos=y |

#

grey3 flat=n color=j wanttitle=n

#

label1="Depth z (km)" label2="X (km)" label3="Y (km)"

#

''')


#
------------------------------------------------------------

par['kz']=2./3.*par['nz']


#
--------------------------------
----------------------------

fdmod.wavelet('wav_',par['frq'],par)

#
------------------------------------------------------------


#
------------------------------------------------------------

# elastic source

#
--------------------------------------------
----------------

Flow('souz','wav_','math output=input*1')

Flow('soux','wav_','math output=input*0')

Flow('souy','wav_','math output=input*0')


Flow('wave
-
3d',['souz','soux','souy'],


'''


cat axis=2 space=n ${SOURCES[1:3]} |


transp plane=12 |


transp plane=23 |


transp plane=12


''')

fdmod.ewavelet3d('wave
-
3d','',par)


#
------------------------------------------------------------

# source/receiver coordinates

#
------------------------------------------------------------


xsou=par
['ox']+(par['nx']/2*par['dx']);

ysou=par['oy']+(par['ny']/2*par['dy']);

#zsou=par['oz']+(par['nz']/2*par['dz']);

zsou= 0.03


center=fdmod.center3d(1.1*xsou,1.1*ysou,1.1*zsou,par)


fdmod.point3d('ss
-
3d',xsou,ysou,zsou,par)

fdmod.horizontal3d('rr
-
3d',0.03,pa
r)

#fdmod.point3d('rr
-
3d',xsou,ysou,zsou,par)


"""

#
------------------------------------------------------------

Flow('zero
-
3d',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


n3=%(ny)d
o3=%(oy)g d3=%(dy)g |


put label1=%(lz)s label2=%(lx)s label3=%(ly)s unit1=%(uz)s unit2=%(ux)s unit3=%(uy)s


''' % par)


#
------------------------------------------------------------


Flow('vp
-
3d','zero
-
3d','math output="%(vp)g"' %par)

Flow('vs
-
3d
','zero
-
3d','math output="%(vs)g"' %par)

Flow('ro
-
3d','zero
-
3d','math output="%(ro)g"' %par)

Result('vp
-
3d',fdmod.cgrey3d('allpos=y bias=2.0'+center,par))

Result('vs
-
3d',fdmod.cgrey3d('allpos=y bias=1.0'+center,par))

Result('ro
-
3d',fdmod.cgrey3d('allpos=y
bias=2.0'+center,par))

"""

stiffness.iso3d('Ic
-
3d','vp
-
3d','vs
-
3d','ro
-
3d',par)


#
------------------------------------------------------------

# 3D elastic modeling

stiffness.cplot3d('Ic
-
3d',1,1,1,par)


fdmod.ewefd3d('Id
-
3d','Iw
-
3d','wave
-
3d','Ic
-
3d','ro
-
3d','ss
-
3d','rr
-
3d','ssou=n opot=n',par)


Result('Iw
-
3d','window n5=1 f5=%(iframe)d n4=1 |' %par


+ fdmod.cgrey3d('pclip=99.9 flat=n screenlg=5 color=j'+center,par))


Result('Id
-
3d',


'''


window n2=1 |


put


n1=%(
nx)d o1=%(ox)g d1=%(dx)g label1=%(lx)s unit1=%(ux)s


n2=%(ny)d o2=%(oy)g d2=%(dy)g label2=%(ly)s unit1=%(uy)s


n3=%(nt)d o3=%(ot)g d3=%(dt)g label3=%(lt)s unit1=%(ut)s |


transp plane=23 |


transp plane=12 |



''' % par


+ fdmod.dgrey3d('pclip=99.9 flat=n color=j'+center+' frame1=%d' % (0.85*par['nt']) ,par))



for i in range(3):


tag = "%d"%i


Result('Iw
-
3d
-
'+tag,'Iw
-
3d',' window n5=1 f5=4 n4=1 f4=%d |'%i


+ fdmod.cgrey3d('pclip=99.9 f
lat=n color=j'+center+' movie=1 frame1=1',par))



End()

###########################################################################


Wavelet



Vp velocity model



Stiffness



Wave
-
field



Zoomed in
trace
from above

wave
-
field



Different polariza
tion wave
-
field


Single Trace


Another different polarization wave
-
field



Single Trace


Model response for receiver line