Guide to programming with madagascar - Revision history

Fomels: typos

2024-11-26T15:34:35Z

typos

← Older revision		Revision as of 15:34, 26 November 2024
Line 9:		Line 9:
	finite-difference modeling <ref>"Hello world" of seismic imaging.</ref>		finite-difference modeling <ref>"Hello world" of seismic imaging.</ref>
	written with the RSF library.		written with the RSF library.
	The ~~program is demonstrated with the~~ C, C++ and Fortran 90 interfaces.		The C, C++, and Fortran 90 interfaces demonstrate the program.
	The acoustic wave-equation		The acoustic wave-equation
	<center><math>		<center><math>
Line 23:		Line 23:
	<math>v</math> is the velocity, and		<math>v</math> is the velocity, and
	<math>U</math> is a scalar wavefield.		<math>U</math> is a scalar wavefield.
	A discrete time-step involves the following computations:		A discrete-time step involves the following computations:
	<center><math>		<center><math>
	U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1} \;,		U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1} \;,
Line 157:		Line 157:
	sf_axis at,az,ax; /* cube axes */		sf_axis at,az,ax; /* cube axes */
	</syntaxhighlight>		</syntaxhighlight>
	#Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="c">		#Declare multi-dimensional arrays for input, output, and computations. <syntaxhighlight lang="c">

	float ww,vv,rr; / I/O arrays*/		float ww,vv,rr; / I/O arrays*/
Line 176:		Line 176:
	sf_oaxa(Fo,at,3);		sf_oaxa(Fo,at,3);
	</syntaxhighlight>		</syntaxhighlight>
	#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="c">		#Allocate arrays and read wavelet, velocity, and reflectivity. <syntaxhighlight lang="c">
	ww=sf_floatalloc(nt); sf_floatread(ww ,nt ,Fw);		ww=sf_floatalloc(nt); sf_floatread(ww ,nt ,Fw);
	vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);		vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
Line 326:		Line 326:
	</syntaxhighlight>		</syntaxhighlight>

	#Declare input, output and auxiliary file cubes (of type <tt>CUB</tt>). <syntaxhighlight lang="cpp">		#Declare input, output, and auxiliary file cubes (of type <tt>CUB</tt>). <syntaxhighlight lang="cpp">
	CUB Fw("in", "i"); Fw.headin(); //Fw.report();		CUB Fw("in", "i"); Fw.headin(); //Fw.report();
	CUB Fv("vel","i"); Fv.headin(); //Fv.report();		CUB Fv("vel","i"); Fv.headin(); //Fv.report();
Line 404:		Line 404:
	call from_par("verb",verb,.false.)		call from_par("verb",verb,.false.)

	! ~~setup~~ I/O files		! Setup I/O files
	Fw=rsf_input ("in")		Fw=rsf_input ("in")
	Fv=rsf_input ("vel")		Fv=rsf_input ("vel")
Line 464:		Line 464:
	</syntaxhighlight>		</syntaxhighlight>

	#Declare input, output and auxiliary file tags. <syntaxhighlight lang="fortran">		#Declare input, output, and auxiliary file tags. <syntaxhighlight lang="fortran">
	type(file) :: Fw,Fv,Fr,Fo ! I/O files		type(file) :: Fw,Fv,Fr,Fo ! I/O files
	</syntaxhighlight>		</syntaxhighlight>
Line 470:		Line 470:
	type(axa) :: at,az,ax ! cube axes		type(axa) :: at,az,ax ! cube axes
	</syntaxhighlight>		</syntaxhighlight>
	#Declare multi-dimensional arrays for input, output and computations. <syntaxhighlight lang="fortran">		#Declare multi-dimensional arrays for input, output, and computations. <syntaxhighlight lang="fortran">
	real, allocatable :: vv(:,:),rr(:,:),ww(:) ! I/O arrays		real, allocatable :: vv(:,:),rr(:,:),ww(:) ! I/O arrays
	real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays		real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays
Line 489:		Line 489:
	call oaxa(Fo,at,3)		call oaxa(Fo,at,3)
	</syntaxhighlight>		</syntaxhighlight>
	#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="fortran">		#Allocate arrays and read wavelet, velocity, and reflectivity. <syntaxhighlight lang="fortran">
	allocate(ww(nt)); call rsf_read(Fw,ww)		allocate(ww(nt)); call rsf_read(Fw,ww)
	allocate(vv(nz,nx)); call rsf_read(Fv,vv)		allocate(vv(nz,nx)); call rsf_read(Fv,vv)
Line 616:		Line 616:
	Fo.putaxis(at,3)		Fo.putaxis(at,3)
	</syntaxhighlight>		</syntaxhighlight>
	#Allocate arrays and read wavelet, velocity and reflectivity. <syntaxhighlight lang="python">		#Allocate arrays and read wavelet, velocity, and reflectivity. <syntaxhighlight lang="python">
	ww = numpy.zeros(nt,'f'); Fw.read(ww)		ww = numpy.zeros(nt,'f'); Fw.read(ww)
	vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)		vv = numpy.zeros([nz,nx],'f'); Fv.read(vv)

Sergey: update + python

2021-05-06T02:29:52Z

update + python

Show changes

Sfomel: /* C++ program */

2021-02-08T00:20:38Z

C++ program

← Older revision		Revision as of 00:20, 8 February 2021
Line 291:		Line 291:
	if(verb) cerr << endl;		if(verb) cerr << endl;
	for (int it=0; it<nt; it++) {		for (int it=0; it<nt; it++) {
	if(verb) cerr << "~~��~~" << it;		if(verb) cerr << "\b\b\b\b\b" << it;

	// 4th order laplacian		// 4th order laplacian

Sfomel: /* C++ program */

2021-02-08T00:19:13Z

C++ program

← Older revision		Revision as of 00:19, 8 February 2021
Line 240:		Line 240:
	#include <rsf.hh>		#include <rsf.hh>
	#include <cub.hh>		#include <cub.hh>
	#include <vai.hh>
			#include "vai.hh"

	using namespace std;		using namespace std;

Line 256:		Line 258:
	CUB Fv("vel","i"); Fv.headin(); //Fv.report();		CUB Fv("vel","i"); Fv.headin(); //Fv.report();
	CUB Fr("ref","i"); Fr.headin(); //Fr.report();		CUB Fr("ref","i"); Fr.headin(); //Fr.report();
	CUB Fo("out","o"); Fo.setup(3~~,Fv.esize()~~);		CUB Fo("out","o"); Fo.setup(3);

	// Read/Write axes		// Read/Write axes
Line 289:		Line 291:
	if(verb) cerr << endl;		if(verb) cerr << endl;
	for (int it=0; it<nt; it++) {		for (int it=0; it<nt; it++) {
	if(verb) cerr << "~~\b\b\b\b\b~~" << it;		if(verb) cerr << "��" << it;

	// 4th order laplacian		// 4th order laplacian
	for (int iz=2; iz<nz-2; iz++) {		for (int ix=2; ix<nx-2; ix++) {
	for (int ix=2; ix<nx-2; ix++) {		for (int iz=2; iz<nz-2; iz++) {
	ud[k(iz,ix)] =		ud[k(iz,ix)] =
	c0* uo[ k(iz ,ix )] * (idx+idz) +		c0* uo[ k(iz ,ix )] * (idx+idz) +
Line 323:		Line 325:
	</syntaxhighlight>		</syntaxhighlight>

	*Declare input, output and auxiliary file cubes (of type <tt>CUB</tt>).		#Declare input, output and auxiliary file cubes (of type <tt>CUB</tt>).
	<syntaxhighlight lang="cpp">		<syntaxhighlight lang="cpp">
	CUB Fw("in", "i"); Fw.headin(); //Fw.report();		CUB Fw("in", "i"); Fw.headin(); //Fw.report();
	CUB Fv("vel","i"); Fv.headin(); //Fv.report();		CUB Fv("vel","i"); Fv.headin(); //Fv.report();
	CUB Fr("ref","i"); Fr.headin(); //Fr.report();		CUB Fr("ref","i"); Fr.headin(); //Fr.report();
	CUB Fo("out","o"); Fo.setup(3~~,Fv.esize()~~);		CUB Fo("out","o"); Fo.setup(3);
	</syntaxhighlight>		</syntaxhighlight>
	*Declare, read and write RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.		#Declare, read and write RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.
	<syntaxhighlight lang="cpp">		<syntaxhighlight lang="cpp">
	sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);		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 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);		sf_axis ax = Fv.getax(1); int nx = sf_n(ax); float dx = sf_d(ax);
			</syntaxhighlight>
	~~Fo.putax(0,az);~~		#Declare multi-dimensional <tt>valarrays</tt> for input, output and read data.
	~~Fo.putax(1,ax);~~
	~~Fo.putax(2,at);~~
	~~Fo.headou();~~
	</syntaxhighlight>
	*Declare multi-dimensional <tt>valarrays</tt> for input, output and read data.
	<syntaxhighlight lang="cpp">		<syntaxhighlight lang="cpp">
	valarray<float> ww( nt ); ww=0; Fw >> ww;		valarray<float> ww( nt ); ww=0; Fw >> ww;
	valarray<float> vv( nz*nx ); vv=0; Fv >> vv;		valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
	valarray<float> rr( nz*nx ); rr=0; Fr >> rr;		valarray<float> rr( nz*nx ); rr=0; Fr >> rr;
	</syntaxhighlight>		</syntaxhighlight>
	*Declare multi-dimensional <tt>valarrays</tt> for temporary storage.		#Declare multi-dimensional <tt>valarrays</tt> for temporary storage. <syntaxhighlight lang="cpp">
	<syntaxhighlight lang="cpp">
	valarray<float> um(nz*nx); um=0;		valarray<float> um(nz*nx); um=0;
	valarray<float> uo(nz*nx); uo=0;		valarray<float> uo(nz*nx); uo=0;
Line 354:		Line 350:
	valarray<float> ud(nz*nx); ud=0;		valarray<float> ud(nz*nx); ud=0;
	</syntaxhighlight>		</syntaxhighlight>
	*Initialize multidimensional <tt>valarray</tt> index counter (of type <tt>VAI</tt>).		#Initialize multidimensional <tt>valarray</tt> index counter (of type <tt>VAI</tt>).
	<syntaxhighlight lang="cpp">		<syntaxhighlight lang="cpp">
	VAI k(nz,nx);		VAI k(nz,nx);
	</syntaxhighlight>		</syntaxhighlight>
	*Loop over time.		#Loop over time.
	<syntaxhighlight lang="cpp">		<syntaxhighlight lang="cpp">
	for (int it=0; it<nt; it++) {		for (int it=0; it<nt; it++) {
	</syntaxhighlight>		</syntaxhighlight>
	*Compute Laplacian: <math>\Delta U</math>.		#Compute Laplacian: <math>\Delta U</math>.
	<syntaxhighlight lang="cpp">		<syntaxhighlight lang="cpp">
	for (int iz=2; iz<nz-2; iz++) {		for (int ix=2; ix<nx-2; ix++) {
	for (int ix=2; ix<nx-2; ix++) {		for (int iz=2; iz<nz-2; iz++) {
	ud[k(iz,ix)] =		ud[k(iz,ix)] =
	c0* uo[ k(iz ,ix )] * (idx+idz) +		c0* uo[ k(iz ,ix )] * (idx+idz) +
Line 374:		Line 370:
	}		}
	}		}
	</syntaxhighlight>		</syntaxhighlight>
	*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>		#Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math> <syntaxhighlight lang="cpp">
	<syntaxhighlight lang="cpp">
	ud -= ww[it] * rr;		ud -= ww[it] * rr;
	</syntaxhighlight>		</syntaxhighlight>
	*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>		#Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math> <syntaxhighlight lang="cpp">
	<syntaxhighlight lang="cpp">
	ud = vvvv;		ud = vvvv;
	</syntaxhighlight>		</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>		#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="cpp">		<syntaxhighlight lang="cpp">
	up=(float)2 * uo - um + ud * dt2;		up=(float)2 * uo - um + ud * dt2;

Sfomel: /* C program */

2021-02-08T00:13:18Z

C program

Show changes

Sfomel at 21:05, 9 August 2018

2018-08-09T21:05:40Z

← Older revision		Revision as of 21:05, 9 August 2018
Line 1:		Line 1:
	<center><font size="-1">''This page was created from the LaTeX source in [~~http~~://~~sourceforge~~.~~net~~/p/~~rsf~~/~~code~~/~~HEAD/tree/trunk~~/book/rsf/rsf/demo.tex book/rsf/rsf/demo.tex] using [[latex2wiki]]''</font></center>		<center><font size="-1">''This page was created from the LaTeX source in [https://github.com/ahay/src/blob/master/book/rsf/rsf/demo.tex book/rsf/rsf/demo.tex] using [[latex2wiki]]''</font></center>

	This guide demonstrates a simple time-domain		This guide demonstrates a simple time-domain

Sfomel: /* Fortran 90 program */

2014-04-08T22:57:42Z

Fortran 90 program

← Older revision		Revision as of 22:57, 8 April 2014
Line 395:		Line 395:
	==Fortran 90 program==		==Fortran 90 program==

			<syntaxhighlight lang="fortran">
	<fortran>
	! time-domain acoustic FD modeling		! time-domain acoustic FD modeling
	program AFDMf90		program AFDMf90
Line 476:		Line 475:
	call exit(0)		call exit(0)
	end program AFDMf90		end program AFDMf90
	</~~fortran~~>		</syntaxhighlight>

	*Declare input, output and auxiliary file tags.		*Declare input, output and auxiliary file tags.
	<fortran>		<syntaxhighlight lang="fortran">
	type(file) :: Fw,Fv,Fr,Fo ! I/O files		type(file) :: Fw,Fv,Fr,Fo ! I/O files
	</~~fortran~~>		</syntaxhighlight>
	*Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.		*Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.
	<fortran>		<syntaxhighlight lang="fortran">
	type(axa) :: at,az,ax ! cube axes		type(axa) :: at,az,ax ! cube axes
	</~~fortran~~>		</syntaxhighlight>
	*Declare multi-dimensional arrays for input, output and computations.		*Declare multi-dimensional arrays for input, output and computations.
	<fortran>		<syntaxhighlight lang="fortran">
	real, allocatable :: vv(:,:),rr(:,:),ww(:) ! I/O arrays		real, allocatable :: vv(:,:),rr(:,:),ww(:) ! I/O arrays
	real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays		real, allocatable :: um(:,:),uo(:,:),up(:,:),ud(:,:) ! tmp arrays
	</~~fortran~~>		</syntaxhighlight>
	*Open files for input/output.		*Open files for input/output.
	<fortran>		<syntaxhighlight lang="fortran">
	Fw=rsf_input ("in")		Fw=rsf_input ("in")
	Fv=rsf_input ("vel")		Fv=rsf_input ("vel")
	Fr=rsf_input ("ref")		Fr=rsf_input ("ref")
	Fo=rsf_output("out")		Fo=rsf_output("out")
	</~~fortran~~>		</syntaxhighlight>
	*Read axes from input files; write axes to output file.		*Read axes from input files; write axes to output file.
	<fortran>		<syntaxhighlight lang="fortran">
	call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2)		call iaxa(Fw,at,1); call iaxa(Fv,az,1); call iaxa(Fv,ax,2)
	call oaxa(Fo,az,1); call oaxa(Fo,ax,2); call oaxa(Fo,at,3)		call oaxa(Fo,az,1); call oaxa(Fo,ax,2); call oaxa(Fo,at,3)
	</~~fortran~~>		</syntaxhighlight>
	*Allocate arrays and read wavelet, velocity and reflectivity.		*Allocate arrays and read wavelet, velocity and reflectivity.
	<fortran>		<syntaxhighlight lang="fortran">
	allocate(ww(at%n )); ww=0.; call rsf_read(Fw,ww)		allocate(ww(at%n )); ww=0.; call rsf_read(Fw,ww)
	allocate(vv(az%n,ax%n)); vv=0.; call rsf_read(Fv,vv)		allocate(vv(az%n,ax%n)); vv=0.; call rsf_read(Fv,vv)
	allocate(rr(az%n,ax%n)); rr=0.; call rsf_read(Fr,rr)		allocate(rr(az%n,ax%n)); rr=0.; call rsf_read(Fr,rr)
	</~~fortran~~>		</syntaxhighlight>
	*Allocate temporary arrays.		*Allocate temporary arrays.
	<fortran>		<syntaxhighlight lang="fortran">
	allocate(um(az%n,ax%n)); um=0.		allocate(um(az%n,ax%n)); um=0.
	allocate(uo(az%n,ax%n)); uo=0.		allocate(uo(az%n,ax%n)); uo=0.
	allocate(up(az%n,ax%n)); up=0.		allocate(up(az%n,ax%n)); up=0.
	allocate(ud(az%n,ax%n)); ud=0.		allocate(ud(az%n,ax%n)); ud=0.
	</~~fortran~~>		</syntaxhighlight>
	*Loop over time.		*Loop over time.
	<fortran>		<syntaxhighlight lang="fortran">
	do it=1,at%n		do it=1,at%n
	</~~fortran~~>		</syntaxhighlight>
	*Compute Laplacian: <math>\Delta U</math>.		*Compute Laplacian: <math>\Delta U</math>.
	<fortran>		<syntaxhighlight lang="fortran">
	do iz=2,az%n-2		do iz=2,az%n-2
	do ix=2,ax%n-2		do ix=2,ax%n-2
Line 532:		Line 531:
	end do		end do
	end do		end do
	</~~fortran~~>		</syntaxhighlight>
	*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>		*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>
	<fortran>		<syntaxhighlight lang="fortran">
	ud = ud - ww(it) * rr		ud = ud - ww(it) * rr
	</~~fortran~~>		</syntaxhighlight>
	*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>		*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>
	<fortran>		<syntaxhighlight lang="fortran">
	ud= ud vvvv		ud= ud vvvv
	</~~fortran~~>		</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>		*Time step: <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>
	<fortran>		<syntaxhighlight lang="fortran">
	up = 2uo - um + ud dt2		up = 2uo - um + ud dt2
	um = uo		um = uo
	uo = up		uo = up
	</~~fortran~~>		</syntaxhighlight>

	==References==		==References==
	<references/>		<references/>

Sfomel: /* C++ program */

2014-04-08T22:55:03Z

C++ program

← Older revision		Revision as of 22:55, 8 April 2014
Line 237:		Line 237:
	==C++ program==		==C++ program==

			<syntaxhighlight lang="cpp">
	<cpp>
	// time-domain acoustic FD modeling		// time-domain acoustic FD modeling
	#include <valarray>		#include <valarray>
Line 325:		Line 324:
	exit(0);		exit(0);
	}		}
			</syntaxhighlight>
	</~~cpp~~>

	*Declare input, output and auxiliary file cubes (of type <tt>CUB</tt>).		*Declare input, output and auxiliary file cubes (of type <tt>CUB</tt>).
	<cpp>		<syntaxhighlight lang="cpp">
	CUB Fw("in", "i"); Fw.headin(); //Fw.report();		CUB Fw("in", "i"); Fw.headin(); //Fw.report();
	CUB Fv("vel","i"); Fv.headin(); //Fv.report();		CUB Fv("vel","i"); Fv.headin(); //Fv.report();
	CUB Fr("ref","i"); Fr.headin(); //Fr.report();		CUB Fr("ref","i"); Fr.headin(); //Fr.report();
	CUB Fo("out","o"); Fo.setup(3,Fv.esize());		CUB Fo("out","o"); Fo.setup(3,Fv.esize());
	</~~cpp~~>		</syntaxhighlight>
	*Declare, read and write RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.		*Declare, read and write RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.
	<cpp>		<syntaxhighlight lang="cpp">
	sf_axis at = Fw.getax(0); int nt = sf_n(at); float dt = sf_d(at);		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 az = Fv.getax(0); int nz = sf_n(az); float dz = sf_d(az);
Line 345:		Line 343:
	Fo.putax(2,at);		Fo.putax(2,at);
	Fo.headou();		Fo.headou();
	</~~cpp~~>		</syntaxhighlight>
	*Declare multi-dimensional <tt>valarrays</tt> for input, output and read data.		*Declare multi-dimensional <tt>valarrays</tt> for input, output and read data.
	<cpp>		<syntaxhighlight lang="cpp">
	valarray<float> ww( nt ); ww=0; Fw >> ww;		valarray<float> ww( nt ); ww=0; Fw >> ww;
	valarray<float> vv( nz*nx ); vv=0; Fv >> vv;		valarray<float> vv( nz*nx ); vv=0; Fv >> vv;
	valarray<float> rr( nz*nx ); rr=0; Fr >> rr;		valarray<float> rr( nz*nx ); rr=0; Fr >> rr;
	</~~cpp~~>		</syntaxhighlight>
	*Declare multi-dimensional <tt>valarrays</tt> for temporary storage.		*Declare multi-dimensional <tt>valarrays</tt> for temporary storage.
	<cpp>		<syntaxhighlight lang="cpp">
	valarray<float> um(nz*nx); um=0;		valarray<float> um(nz*nx); um=0;
	valarray<float> uo(nz*nx); uo=0;		valarray<float> uo(nz*nx); uo=0;
	valarray<float> up(nz*nx); up=0;		valarray<float> up(nz*nx); up=0;
	valarray<float> ud(nz*nx); ud=0;		valarray<float> ud(nz*nx); ud=0;
	</~~cpp~~>		</syntaxhighlight>
	*Initialize multidimensional <tt>valarray</tt> index counter (of type <tt>VAI</tt>).		*Initialize multidimensional <tt>valarray</tt> index counter (of type <tt>VAI</tt>).
	<cpp>		<syntaxhighlight lang="cpp">
	VAI k(nz,nx);		VAI k(nz,nx);
	</~~cpp~~>		</syntaxhighlight>
	*Loop over time.		*Loop over time.
	<cpp>		<syntaxhighlight lang="cpp">
	for (int it=0; it<nt; it++) {		for (int it=0; it<nt; it++) {
	</~~cpp~~>		</syntaxhighlight>
	*Compute Laplacian: <math>\Delta U</math>.		*Compute Laplacian: <math>\Delta U</math>.
	<cpp>		<syntaxhighlight lang="cpp">
	for (int iz=2; iz<nz-2; iz++) {		for (int iz=2; iz<nz-2; iz++) {
	for (int ix=2; ix<nx-2; ix++) {		for (int ix=2; ix<nx-2; ix++) {
Line 379:		Line 377:
	}		}
	}		}
	</~~cpp~~>		</syntaxhighlight>
	*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>		*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>
	<cpp>		<syntaxhighlight lang="cpp">
	ud -= ww[it] * rr;		ud -= ww[it] * rr;
	</~~cpp~~>		</syntaxhighlight>
	*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>		*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>
	<cpp>		<syntaxhighlight lang="cpp">
	ud = vvvv;		ud = vvvv;
	</~~cpp~~>		</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>		*Time step: <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>
	<cpp>		<syntaxhighlight lang="cpp">
	up=(float)2 * uo - um + ud * dt2;		up=(float)2 * uo - um + ud * dt2;
	um = uo;		um = uo;
	uo = up;		uo = up;
	</~~cpp~~>		</syntaxhighlight>

	==Fortran 90 program==		==Fortran 90 program==

Sfomel: /* C program */

2014-04-08T22:52:57Z

C program

← Older revision		Revision as of 22:52, 8 April 2014
Line 32:		Line 32:
	==C program==		==C program==

			<syntaxhighlight lang="c">
	<c>
	/* time-domain acoustic FD modeling */		/* time-domain acoustic FD modeling */
	#include <rsf.h>		#include <rsf.h>
Line 144:		Line 143:
	exit (0);		exit (0);
	}		}
	</c>		</syntaxhighlight>



	*Declare input, output and auxiliary file tags: <tt>Fw</tt> for input wavelet, <tt>Fv</tt> for velocity, <tt>Fr</tt> for reflectivity, and <tt>Fo</tt> for output wavefield.		*Declare input, output and auxiliary file tags: <tt>Fw</tt> for input wavelet, <tt>Fv</tt> for velocity, <tt>Fr</tt> for reflectivity, and <tt>Fo</tt> for output wavefield.
	<c>		<syntaxhighlight lang="c">
	sf_file Fw,Fv,Fr,Fo; /* I/O files */		sf_file Fw,Fv,Fr,Fo; /* I/O files */
	</c>		</syntaxhighlight>
	*Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.		*Declare RSF cube axes: <tt>at</tt> time axis, <tt>ax</tt> space axis, <tt>az</tt> depth axis.
	<c>		<syntaxhighlight lang="c">
	sf_axis at,az,ax; /* cube axes */		sf_axis at,az,ax; /* cube axes */
	</c>		</syntaxhighlight>
	*Declare multi-dimensional arrays for input, output and computations.		*Declare multi-dimensional arrays for input, output and computations.
	<c>		<syntaxhighlight lang="c">
	float ww,vv,rr; / I/O arrays*/		float ww,vv,rr; / I/O arrays*/
	</c>		</syntaxhighlight>
	*Open files for input/output.		*Open files for input/output.
	<c>		<syntaxhighlight lang="c">
	Fw = sf_input ("in" );		Fw = sf_input ("in" );
	Fo = sf_output("out");		Fo = sf_output("out");
	Fv = sf_input ("vel");		Fv = sf_input ("vel");
	Fr = sf_input ("ref");		Fr = sf_input ("ref");
	</c>		</syntaxhighlight>
	*Read axes from input files; write axes to output file.		*Read axes from input files; write axes to output file.
	<c>		<syntaxhighlight lang="c">
	at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);		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);		az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
Line 176:		Line 173:
	sf_oaxa(Fo,ax,2);		sf_oaxa(Fo,ax,2);
	sf_oaxa(Fo,at,3);		sf_oaxa(Fo,at,3);
	</c>		</syntaxhighlight>
	*Allocate arrays and read wavelet, velocity and reflectivity.		*Allocate arrays and read wavelet, velocity and reflectivity.
	<c>		<syntaxhighlight lang="c">
	ww=sf_floatalloc(nt); sf_floatread(ww ,nt ,Fw);		ww=sf_floatalloc(nt); sf_floatread(ww ,nt ,Fw);
	vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);		vv=sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
	rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);		rr=sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);
	</c>		</syntaxhighlight>
	*Allocate temporary arrays.		*Allocate temporary arrays.
	<c>		<syntaxhighlight lang="c">
	um=sf_floatalloc2(nz,nx);		um=sf_floatalloc2(nz,nx);
	uo=sf_floatalloc2(nz,nx);		uo=sf_floatalloc2(nz,nx);
	up=sf_floatalloc2(nz,nx);		up=sf_floatalloc2(nz,nx);
	ud=sf_floatalloc2(nz,nx);		ud=sf_floatalloc2(nz,nx);
	</c>		</syntaxhighlight>
	*Loop over time.		*Loop over time.
	<c>		<syntaxhighlight lang="c">
	for (it=0; it<nt; it++) {		for (it=0; it<nt; it++) {
	</c>		</syntaxhighlight>
	*Compute Laplacian: <math>\Delta U</math>.		*Compute Laplacian: <math>\Delta U</math>.
	<c>		<syntaxhighlight lang="c">
	for (iz=2; iz<nz-2; iz++) {		for (iz=2; iz<nz-2; iz++) {
	for (ix=2; ix<nx-2; ix++) {		for (ix=2; ix<nx-2; ix++) {
Line 206:		Line 203:
	}		}
	}		}
	</c>		</syntaxhighlight>
	*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>		*Inject source wavelet: <math>\left[ \Delta U - f(t) \right]</math>
	<c>		<syntaxhighlight lang="c">
	for (iz=0; iz<nz; iz++) {		for (iz=0; iz<nz; iz++) {
	for (ix=0; ix<nx; ix++) {		for (ix=0; ix<nx; ix++) {
Line 214:		Line 211:
	}		}
	}		}
	</c>		</syntaxhighlight>
	*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>		*Scale by velocity: <math>\left[ \Delta U - f(t) \right] v^2</math>
	<c>		<syntaxhighlight lang="c">
	for (iz=0; iz<nz; iz++) {		for (iz=0; iz<nz; iz++) {
	for (ix=0; ix<nx; ix++) {		for (ix=0; ix<nx; ix++) {
Line 222:		Line 219:
	}		}
	}		}
	</c>		</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>		*Time step: <math>U_{i+1} = \left[ \Delta U -f(t) \right] v^2 \Delta t^2 + 2 U_{i} - U_{i-1}</math>
	<c>		<syntaxhighlight lang="c">
	for (iz=0; iz<nz; iz++) {		for (iz=0; iz<nz; iz++) {
	for (ix=0; ix<nx; ix++) {		for (ix=0; ix<nx; ix++) {
Line 236:		Line 233:
	}		}
	}		}
	</c>		</syntaxhighlight>

	==C++ program==		==C++ program==

Sfomel at 20:41, 10 August 2013

2013-08-10T20:41:22Z

← Older revision		Revision as of 20:41, 10 August 2013
Line 1:		Line 1:
	<center><font size="-1">''This page was created from the LaTeX source in [http://~~rsf.svn.~~sourceforge.net/~~viewvc~~/rsf/trunk/book/rsf/rsf/demo.tex~~?view=markup~~ book/rsf/rsf/demo.tex] using [[latex2wiki]]''</font></center>		<center><font size="-1">''This page was created from the LaTeX source in [http://sourceforge.net/p/rsf/code/HEAD/tree/trunk/book/rsf/rsf/demo.tex book/rsf/rsf/demo.tex] using [[latex2wiki]]''</font></center>

	This guide demonstrates a simple time-domain		This guide demonstrates a simple time-domain