Guide to programming with madagascar

This page was created from the LaTeX source in book/rsf/rsf/demo.tex using latex2wiki

This guide demonstrates a simple time-domain finite-differences modeling code in RSF.

Introduction

This section presents time-domain finite-difference modeling written with the RSF library. The C, C++, and Fortran 90 interfaces demonstrate the program. The acoustic wave-equation

$\Delta U-{\frac {1}{v^{2}}}{\frac {\partial ^{2}U}{\partial t^{2}}}=f(t)$

can be written as

$\left[\Delta U-f(t)\right]v^{2}={\frac {\partial ^{2}U}{\partial t^{2}}}\;.$

$\Delta$ is the Laplacian symbol, $f(t)$ is the source wavelet, $v$ is the velocity, and $U$ is a scalar wavefield. A discrete-time step involves the following computations:

$U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}\;,$

where $U_{i-1}$, $U_{i}$ and $U_{i+1}$ represent the propagating wavefield at various time steps.

C program

Wave propagation snapshot.

C program

<syntaxhighlight lang="c"> /* time-domain acoustic FD modeling */

1. include <rsf.h>

int main(int argc, char* argv[]) {

   /* Laplacian coefficients */
   float c0=-30./12.,c1=+16./12.,c2=- 1./12.;

   bool verb;           /* verbose flag */
   sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
   sf_axis at,az,ax;    /* cube axes */
   int it,iz,ix;        /* index variables */
   int nt,nz,nx;
   float dt,dz,dx,idx,idz,dt2;

   float  *ww,**vv,**rr;     /* I/O arrays*/
   float **um,**uo,**up,**ud;/* tmp arrays */

   sf_init(argc,argv);
   if(! sf_getbool("verb",&verb)) verb=0; /* verbose flag */

   /* setup I/O files */
   Fw = sf_input ("in" );
   Fo = sf_output("out");
   Fv = sf_input ("vel");
   Fr = sf_input ("ref");

   /* Read/Write axes */
   at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);
   az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
   ax = sf_iaxa(Fv,2); nx = sf_n(ax); dx = sf_d(ax);

   sf_oaxa(Fo,az,1); 
   sf_oaxa(Fo,ax,2); 
   sf_oaxa(Fo,at,3);

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

   /* read wavelet, velocity & reflectivity */
   ww=sf_floatalloc(nt);     sf_floatread(ww   ,nt   ,Fw);
   vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
   rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);

   /* allocate temporary arrays */
   um=sf_floatalloc2(nz,nx);
   uo=sf_floatalloc2(nz,nx);
   up=sf_floatalloc2(nz,nx);
   ud=sf_floatalloc2(nz,nx);

   for (ix=0; ix<nx; ix++) {

for (iz=0; iz<nz; iz++) { um[ix][iz]=0; uo[ix][iz]=0; up[ix][iz]=0; ud[ix][iz]=0; }

   }

   /* MAIN LOOP */
   if(verb) fprintf(stderr,"

");

   for (it=0; it<nt; it++) {

if(verb) fprintf(stderr,"\b\b\b\b\b %d",it);

/* 4th order laplacian */ for (ix=2; ix<nx-2; ix++) { for (iz=2; iz<nz-2; iz++) { ud[ix][iz] = c0* uo[ix ][iz ] * (idx+idz) + c1*(uo[ix-1][iz ] + uo[ix+1][iz ])*idx + c2*(uo[ix-2][iz ] + uo[ix+2][iz ])*idx + c1*(uo[ix ][iz-1] + uo[ix ][iz+1])*idz + c2*(uo[ix ][iz-2] + uo[ix ][iz+2])*idz; } }

/* inject wavelet */ for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { ud[ix][iz] -= ww[it] * rr[ix][iz]; } }

/* scale by velocity */ for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; } }

/* time step */ for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { up[ix][iz] = 2*uo[ix][iz] - um[ix][iz] + ud[ix][iz] * dt2;

um[ix][iz] = uo[ix][iz]; uo[ix][iz] = up[ix][iz]; } }

/* write wavefield to output */ sf_floatwrite(uo[0],nz*nx,Fo);

   }
   if(verb) fprintf(stderr,"\n");

   exit (0);

} </syntaxhighlight>

1. Declare input, output and auxiliary file tags: Fw for input wavelet, Fv for velocity, Fr for reflectivity, and Fo for output wavefield. <syntaxhighlight lang="c">

   sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */

</syntaxhighlight>

1. Declare RSF cube axes: at time axis, ax space axis, az depth axis. \tiny <syntaxhighlight lang="c">

   sf_axis at,az,ax;    /* cube axes */

</syntaxhighlight>

1. Declare multi-dimensional arrays for input, output, and computations. <syntaxhighlight lang="c">

   float  *ww,**vv,**rr;     /* I/O arrays*/

</syntaxhighlight>

1. Open files for input/output. <syntaxhighlight lang="c">

   Fw = sf_input ("in" );
   Fo = sf_output("out");
   Fv = sf_input ("vel");
   Fr = sf_input ("ref");

</syntaxhighlight>

1. Read axes from input files; write axes to output file. <syntaxhighlight lang="c">

   at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);
   az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
   ax = sf_iaxa(Fv,2); nx = sf_n(ax); dx = sf_d(ax);

   sf_oaxa(Fo,az,1); 
   sf_oaxa(Fo,ax,2); 
   sf_oaxa(Fo,at,3);

</syntaxhighlight>

1. Allocate arrays and read wavelet, velocity, and reflectivity. <syntaxhighlight lang="c">

   ww=sf_floatalloc(nt);     sf_floatread(ww   ,nt   ,Fw);
   vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
   rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);

</syntaxhighlight>

1. Allocate temporary arrays. <syntaxhighlight lang="c">

   um=sf_floatalloc2(nz,nx);
   uo=sf_floatalloc2(nz,nx);
   up=sf_floatalloc2(nz,nx);
   ud=sf_floatalloc2(nz,nx);

</syntaxhighlight>

1. Loop over time. <syntaxhighlight lang="c">

   for (it=0; it<nt; it++) {

</syntaxhighlight>

1. Compute Laplacian: $\Delta U$. <syntaxhighlight lang="c">

for (ix=2; ix<nx-2; ix++) { for (iz=2; iz<nz-2; iz++) { ud[ix][iz] = c0* uo[ix ][iz ] * (idx+idz) + c1*(uo[ix-1][iz ] + uo[ix+1][iz ])*idx + c2*(uo[ix-2][iz ] + uo[ix+2][iz ])*idx + c1*(uo[ix ][iz-1] + uo[ix ][iz+1])*idz + c2*(uo[ix ][iz-2] + uo[ix ][iz+2])*idz; } } </syntaxhighlight>

1. Inject source wavelet: $\left[\Delta U-f(t)\right]$ <syntaxhighlight lang="c">

for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { ud[ix][iz] -= ww[it] * rr[ix][iz]; } } </syntaxhighlight>

1. Scale by velocity: $\left[\Delta U-f(t)\right]v^{2}$ <syntaxhighlight lang="c">

for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { ud[ix][iz] *= vv[ix][iz]*vv[ix][iz]; } } </syntaxhighlight>

1. Time step: $U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}$ <syntaxhighlight lang="c">

for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { up[ix][iz] = 2*uo[ix][iz] - um[ix][iz] + ud[ix][iz] * dt2;

um[ix][iz] = uo[ix][iz]; uo[ix][iz] = up[ix][iz]; } } </syntaxhighlight>

\newpage

C++ program

<syntaxhighlight lang="cpp"> // time-domain acoustic FD modeling

1. include <valarray>
2. include <iostream>
3. include <rsf.hh>
4. include <cub.hh>

1. include "vai.hh"

using namespace std;

int main(int argc, char* argv[]) {

   // Laplacian coefficients
   float c0=-30./12.,c1=+16./12.,c2=- 1./12.;

   sf_init(argc,argv);// init RSF
   bool verb;         // vebose flag
   if(! sf_getbool("verb",&verb)) verb=0;

   // setup I/O files
   CUB Fw("in", "i"); Fw.headin(); //Fw.report();
   CUB Fv("vel","i"); Fv.headin(); //Fv.report();
   CUB Fr("ref","i"); Fr.headin(); //Fr.report();
   CUB Fo("out","o"); Fo.setup(3); 

   // Read/Write axes
   sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);
   sf_axis az = Fv.getax(0); int nz = sf_n(az); float dz = sf_d(az);
   sf_axis ax = Fv.getax(1); int nx = sf_n(ax); float dx = sf_d(ax);

   Fo.putax(0,az); 
   Fo.putax(1,ax); 
   Fo.putax(2,at);
   Fo.headou();

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

   // read wavelet, velocity and reflectivity
   valarray<float> ww( nt    ); ww=0; Fw >> ww;
   valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
   valarray<float> rr( nz*nx ); rr=0; Fr >> rr;
  
   // allocate temporary arrays
   valarray<float> um(nz*nx); um=0;
   valarray<float> uo(nz*nx); uo=0;
   valarray<float> up(nz*nx); up=0;
   valarray<float> ud(nz*nx); ud=0;

   // init ValArray Index counter
   VAI k(nz,nx);

   // MAIN LOOP
   if(verb) cerr << endl;
   for (int it=0; it<nt; it++) {

if(verb) cerr << "\b\b\b\b\b" << it;

// 4th order laplacian for (int ix=2; ix<nx-2; ix++) { for (int iz=2; iz<nz-2; iz++) { ud[k(iz,ix)] = c0* uo[ k(iz ,ix )] * (idx+idz) + c1*(uo[ k(iz ,ix-1)]+uo[ k(iz ,ix+1)]) * idx + c1*(uo[ k(iz-1,ix )]+uo[ k(iz+1,ix )]) * idz + c2*(uo[ k(iz ,ix-2)]+uo[ k(iz ,ix+2)]) * idx + c2*(uo[ k(iz-2,ix )]+uo[ k(iz+2,ix )]) * idz; } }

// inject wavelet ud -= ww[it] * rr;

// scale by velocity ud *= vv*vv;

// time step up=(float)2 * uo - um + ud * dt2; um = uo; uo = up;

// write wavefield to output output Fo << uo;

   }
   if(verb) cerr << endl;

   exit(0);

}

</syntaxhighlight>

1. Declare input, output, and auxiliary file cubes (of type CUB). <syntaxhighlight lang="cpp">

   CUB Fw("in", "i"); Fw.headin(); //Fw.report();
   CUB Fv("vel","i"); Fv.headin(); //Fv.report();
   CUB Fr("ref","i"); Fr.headin(); //Fr.report();
   CUB Fo("out","o"); Fo.setup(3); 

</syntaxhighlight>

1. Declare, read and write RSF cube axes: at time axis, ax space axis, az depth axis. <syntaxhighlight lang="cpp">

   sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);
   sf_axis az = Fv.getax(0); int nz = sf_n(az); float dz = sf_d(az);
   sf_axis ax = Fv.getax(1); int nx = sf_n(ax); float dx = sf_d(ax);

</syntaxhighlight>

1. Declare multi-dimensional valarrays for input, output and read data. <syntaxhighlight lang="cpp">

   valarray<float> ww( nt    ); ww=0; Fw >> ww;
   valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
   valarray<float> rr( nz*nx ); rr=0; Fr >> rr;

</syntaxhighlight>

1. Declare multi-dimensional valarrays for temporary storage. <syntaxhighlight lang="cpp">

   valarray<float> um(nz*nx); um=0;
   valarray<float> uo(nz*nx); uo=0;
   valarray<float> up(nz*nx); up=0;
   valarray<float> ud(nz*nx); ud=0;

</syntaxhighlight>

1. Initialize multidimensional valarray index counter (of type VAI). <syntaxhighlight lang="cpp">

   VAI k(nz,nx);

</syntaxhighlight>

1. Loop over time. <syntaxhighlight lang="cpp">

   for (int it=0; it<nt; it++) {

</syntaxhighlight>

1. Compute Laplacian: $\Delta U$. <syntaxhighlight lang="cpp">

for (int ix=2; ix<nx-2; ix++) { for (int iz=2; iz<nz-2; iz++) { ud[k(iz,ix)] = c0* uo[ k(iz ,ix )] * (idx+idz) + c1*(uo[ k(iz ,ix-1)]+uo[ k(iz ,ix+1)]) * idx + c1*(uo[ k(iz-1,ix )]+uo[ k(iz+1,ix )]) * idz + c2*(uo[ k(iz ,ix-2)]+uo[ k(iz ,ix+2)]) * idx + c2*(uo[ k(iz-2,ix )]+uo[ k(iz+2,ix )]) * idz; } } </syntaxhighlight>

1. Inject source wavelet: $\left[\Delta U-f(t)\right]$ <syntaxhighlight lang="cpp">

ud -= ww[it] * rr; </syntaxhighlight>

1. Scale by velocity: $\left[\Delta U-f(t)\right]v^{2}$ <syntaxhighlight lang="cpp">

ud *= vv*vv; </syntaxhighlight>

1. Time step: $U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}$ <syntaxhighlight lang="cpp">

up=(float)2 * uo - um + ud * dt2; um = uo; uo = up; </syntaxhighlight>

Fortran 90 program

<syntaxhighlight lang="fortran"> ! time-domain acoustic FD modeling program AFDMf90

 use rsf

 implicit none

 ! Laplacian coefficients
 real :: c0=-30./12.,c1=+16./12.,c2=- 1./12.

 logical    :: verb         ! verbose flag
 type(file) :: Fw,Fv,Fr,Fo  ! I/O files
 type(axa)  :: at,az,ax     ! cube axes
 integer    :: it,iz,ix     ! index variables
 integer    :: nt,nz,nx
 real       :: dt,dz,dx
 real       :: idx,idz,dt2

 real, allocatable :: vv(:,:),rr(:,:),ww(:)           ! I/O arrays
 real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays

 call sf_init() ! init RSF
 call from_par("verb",verb,.false.)

 ! Setup I/O files
 Fw=rsf_input ("in")
 Fv=rsf_input ("vel")
 Fr=rsf_input ("ref")
 Fo=rsf_output("out")

 ! Read/Write axes
 call iaxa(Fw,at,1); nt = at%n; dt = at%d
 call iaxa(Fv,az,1); nz = az%n; dz = az%d
 call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d
 
 call oaxa(Fo,az,1)
 call oaxa(Fo,ax,2)
 call oaxa(Fo,at,3)

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

 ! read wavelet, velocity & reflectivity
 allocate(ww(nt));    call rsf_read(Fw,ww)
 allocate(vv(nz,nx)); call rsf_read(Fv,vv)
 allocate(rr(nz,nx)); call rsf_read(Fr,rr)

 ! allocate temporary arrays
 allocate(um(nz,nx)); um=0.
 allocate(uo(nz,nx)); uo=0.
 allocate(up(nz,nx)); up=0.
 allocate(ud(nz,nx)); ud=0.

 ! MAIN LOOP
 do it=1,nt
    if(verb) write (0,*) it

    ud(3:nz-2,3:nx-2) = &
         c0* uo(3:nz-2,3:nx-2) * (idx + idz)            + &
         c1*(uo(3:nz-2,2:nx-3) + uo(3:nz-2,4:nx-1))*idx + &
         c2*(uo(3:nz-2,1:nx-4) + uo(3:nz-2,5:nx  ))*idx + &
         c1*(uo(2:nz-3,3:nx-2) + uo(4:nz-1,3:nx-2))*idz + &
         c2*(uo(1:nz-4,3:nx-2) + uo(5:nz  ,3:nx-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

    ! write wavefield to output
    call rsf_write(Fo,uo)
 end do

 call exit(0)

end program AFDMf90 </syntaxhighlight>

1. Declare input, output, and auxiliary file tags. <syntaxhighlight lang="fortran">

 type(file) :: Fw,Fv,Fr,Fo  ! I/O files

</syntaxhighlight>

1. Declare RSF cube axes: at time axis, ax space axis, az depth axis. <syntaxhighlight lang="fortran">

 type(axa)  :: at,az,ax     ! cube axes

</syntaxhighlight>

1. Declare multi-dimensional arrays for input, output, and computations. <syntaxhighlight lang="fortran">

 real, allocatable :: vv(:,:),rr(:,:),ww(:)           ! I/O arrays
 real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays

</syntaxhighlight>

1. Open files for input/output. <syntaxhighlight lang="fortran">

 Fw=rsf_input ("in")
 Fv=rsf_input ("vel")
 Fr=rsf_input ("ref")
 Fo=rsf_output("out")

</syntaxhighlight>

1. Read axes from input files; write axes to output file. <syntaxhighlight lang="fortran">

 call iaxa(Fw,at,1); nt = at%n; dt = at%d
 call iaxa(Fv,az,1); nz = az%n; dz = az%d
 call iaxa(Fv,ax,2); nx = ax%n; dx = ax%d
 
 call oaxa(Fo,az,1)
 call oaxa(Fo,ax,2)
 call oaxa(Fo,at,3)

</syntaxhighlight>

1. Allocate arrays and read wavelet, velocity, and reflectivity. <syntaxhighlight lang="fortran">

 allocate(ww(nt));    call rsf_read(Fw,ww)
 allocate(vv(nz,nx)); call rsf_read(Fv,vv)
 allocate(rr(nz,nx)); call rsf_read(Fr,rr)

</syntaxhighlight>

1. Allocate temporary arrays. <syntaxhighlight lang="fortran">

 allocate(um(nz,nx)); um=0.
 allocate(uo(nz,nx)); uo=0.
 allocate(up(nz,nx)); up=0.
 allocate(ud(nz,nx)); ud=0.

</syntaxhighlight>

1. Loop over time. <syntaxhighlight lang="fortran">

 do it=1,nt

</syntaxhighlight>

1. Compute Laplacian: $\Delta U$. <syntaxhighlight lang="fortran">

    ud(3:nz-2,3:nx-2) = &
         c0* uo(3:nz-2,3:nx-2) * (idx + idz)            + &
         c1*(uo(3:nz-2,2:nx-3) + uo(3:nz-2,4:nx-1))*idx + &
         c2*(uo(3:nz-2,1:nx-4) + uo(3:nz-2,5:nx  ))*idx + &
         c1*(uo(2:nz-3,3:nx-2) + uo(4:nz-1,3:nx-2))*idz + &
         c2*(uo(1:nz-4,3:nx-2) + uo(5:nz  ,3:nx-2))*idz

</syntaxhighlight>

1. Inject source wavelet: $\left[\Delta U-f(t)\right]$ <syntaxhighlight lang="fortran">

    ud = ud - ww(it) * rr

</syntaxhighlight>

1. Scale by velocity: $\left[\Delta U-f(t)\right]v^{2}$ <syntaxhighlight lang="fortran">

    ud= ud *vv*vv

</syntaxhighlight>

1. Time step: $U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}$ <syntaxhighlight lang="fortran">

    up = 2*uo - um + ud * dt2
    um =   uo
    uo =   up

</syntaxhighlight>

Python program

<syntaxhighlight lang="python">

1. !/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

1. setup I/O files

Fw=m8r.Input() Fv=m8r.Input ("vel") Fr=m8r.Input ("ref") Fo=m8r.Output()

1. 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)

1. 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)

1. 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')

1. 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>

Wave propagation snapshot.

1. Open files for input/output. <syntaxhighlight lang="python">

Fw=m8r.Input() Fv=m8r.Input ("vel") Fr=m8r.Input ("ref") Fo=m8r.Output() </syntaxhighlight>

1. 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>

1. 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>

1. 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>

1. Loop over time. <syntaxhighlight lang="python">

for it in range(nt): </syntaxhighlight>

1. Compute Laplacian: $\Delta U$. <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>

1. Inject source wavelet: $\left[\Delta U-f(t)\right]$ <syntaxhighlight lang="python">

   ud = ud - ww[it] * rr

</syntaxhighlight>

1. Scale by velocity: $\left[\Delta U-f(t)\right]v^{2}$ <syntaxhighlight lang="python">

   ud= ud *vv*vv

</syntaxhighlight>

1. Time step: $U_{i+1}=\left[\Delta U-f(t)\right]v^{2}\Delta t^{2}+2U_{i}-U_{i-1}$ <syntaxhighlight lang="python">

   up = 2*uo - um + ud * dt2
   um =   uo
   uo =   up

</syntaxhighlight>