Editing Guide to programming with madagascar (section)

==Python program==

<syntaxhighlight lang="python">
#!/usr/bin/env python

import sys
import numpy
import m8r

c0=-30./12.
c1=+16./12.
c2=- 1./12.

par = m8r.Par()
verb = par.bool("verb",False) # verbosity

# setup I/O files
Fw=m8r.Input()
Fv=m8r.Input ("vel")
Fr=m8r.Input ("ref")
Fo=m8r.Output()

# Read/Write axes
at = Fw.axis(1); nt = at['n']; dt = at['d']
az = Fv.axis(1); nz = az['n']; dz = az['d']
ax = Fv.axis(2); nx = ax['n']; dx = ax['d']

Fo.putaxis(az,1)
Fo.putaxis(ax,2)
Fo.putaxis(at,3)

dt2 =    dt*dt
idz = 1/(dz*dz)
idx = 1/(dx*dx) 

# read wavelet, velocity & reflectivity
ww = numpy.zeros(nt,'f');      Fw.read(ww)
vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)
rr = numpy.zeros([nz,nx],'f'); Fr.read(rr)

# allocate temporary arrays
um = numpy.zeros([nz,nx],'f')
uo = numpy.zeros([nz,nx],'f')
up = numpy.zeros([nz,nx],'f')
ud = numpy.zeros([nz,nx],'f')

# MAIN LOOP
for it in range(nt):
    if verb:
        sys.stderr.write("\b\b\b\b\b %d" % it)

    ud[2:-2,2:-2] = \
    c0* uo[2:-2,2:-2] * (idx + idz)        + \
    c1*(uo[2:-2,1:-3] + uo[2:-2,3:-1])*idx + \
    c2*(uo[2:-2, :-4] + uo[2:-2,4:  ])*idx + \
    c1*(uo[1:-3,2:-2] + uo[3:-1,2:-2])*idz + \
    c2*(uo[ :-4,2:-2] + uo[4:  ,2:-2])*idz

    # inject wavelet
    ud = ud - ww[it] * rr

    # scale by velocity
    ud= ud *vv*vv

    # time step
    up = 2*uo - um + ud * dt2
    um =   uo
    uo =   up

if verb:
    sys.stderr.write("\n")
    Fo.write(uo)

sys.exit(0)
</syntaxhighlight>

[[Image:wavepython.png|frame|center|Wave propagation snapshot.]]

#Open files for input/output. <syntaxhighlight lang="python">
Fw=m8r.Input()
Fv=m8r.Input ("vel")
Fr=m8r.Input ("ref")
Fo=m8r.Output()
</syntaxhighlight>  
#Read axes from input files; write axes to output file. <syntaxhighlight lang="python">
at = Fw.axis(1); nt = at['n']; dt = at['d']
az = Fv.axis(1); nz = az['n']; dz = az['d']
ax = Fv.axis(2); nx = ax['n']; dx = ax['d']

Fo.putaxis(az,1)
Fo.putaxis(ax,2)
Fo.putaxis(at,3)
</syntaxhighlight>   
#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="python">
ww = numpy.zeros(nt,'f');      Fw.read(ww)
vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)
rr = numpy.zeros([nz,nx],'f'); Fr.read(rr)
</syntaxhighlight>   
#Allocate temporary arrays. <syntaxhighlight lang="python">
um = numpy.zeros([nz,nx],'f')
uo = numpy.zeros([nz,nx],'f')
up = numpy.zeros([nz,nx],'f')
ud = numpy.zeros([nz,nx],'f')
</syntaxhighlight>   
#Loop over time. <syntaxhighlight lang="python">
for it in range(nt):
</syntaxhighlight>   
#Compute Laplacian: <math>\Delta U</math>. <syntaxhighlight lang="python">
    ud[2:-2,2:-2] = \
    c0* uo[2:-2,2:-2] * (idx + idz)        + \
    c1*(uo[2:-2,1:-3] + uo[2:-2,3:-1])*idx + \
    c2*(uo[2:-2, :-4] + uo[2:-2,4:  ])*idx + \
    c1*(uo[1:-3,2:-2] + uo[3:-1,2:-2])*idz + \
    c2*(uo[ :-4,2:-2] + uo[4:  ,2:-2])*idz
</syntaxhighlight>   
#Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="python">
    ud = ud - ww[it] * rr
</syntaxhighlight>   
#Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="python">
    ud= ud *vv*vv
</syntaxhighlight>   
#Time step:  <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>  <syntaxhighlight lang="python">
    up = 2*uo - um + ud * dt2
    um =   uo
    uo =   up
</syntaxhighlight>