|
|
|
| Madagascar Programming Reference Manual | |
|
|
|
Next: Explanation of the code
Up: An example: Finite-Difference modeling
Previous: Introduction
1 /* time-domain acoustic FD modeling */
2 #include <rsf.h>
3
4 int main(int argc, char* argv[])
5 {
6 /* Laplacian coefficients */
7 float c0=-30./12.,c1=+16./12.,c2=- 1./12.;
8
9 bool verb; /* verbose flag */
10 sf_file Fw=NULL,Fv=NULL,Fr=NULL,Fo=NULL; /* I/O files */
11 sf_axis at,az,ax; /* cube axes */
12 int it,iz,ix; /* index variables */
13 int nt,nz,nx;
14 float dt,dz,dx,idx,idz,dt2;
15
16 float *ww,**vv,**rr; /* I/O arrays*/
17 float **um,**uo,**up,**ud;/* tmp arrays */
18
19 sf_init(argc,argv);
20 if(! sf_getbool("verb",&verb)) verb=0;
21
22 /* setup I/O files */
23 Fw = sf_input ("in" );
24 Fo = sf_output("out");
25 Fv = sf_input ("vel");
26 Fr = sf_input ("ref");
27
28 /* Read/Write axes */
29 at = sf_iaxa(Fw,1); nt = sf_n(at); dt = sf_d(at);
30 az = sf_iaxa(Fv,1); nz = sf_n(az); dz = sf_d(az);
31 ax = sf_iaxa(Fv,2); nx = sf_n(ax); dx = sf_d(ax);
32
33 sf_oaxa(Fo,az,1);
34 sf_oaxa(Fo,ax,2);
35 sf_oaxa(Fo,at,3);
36
37 dt2 = dt*dt;
38 idz = 1/(dz*dz);
39 idx = 1/(dx*dx);
40
41 /* read wavelet, velocity & reflectivity */
42 ww = sf_floatalloc(nt); sf_floatread(ww ,nt ,Fw);
43 vv = sf_floatalloc2(nz,nx); sf_floatread(vv[0],nz*nx,Fv);
44 rr = sf_floatalloc2(nz,nx); sf_floatread(rr[0],nz*nx,Fr);
45
46 /* allocate temporary arrays */
47 um = sf_floatalloc2(nz,nx);
48 uo = sf_floatalloc2(nz,nx);
49 up = sf_floatalloc2(nz,nx);
50 ud = sf_floatalloc2(nz,nx);
51
52 for (iz=0; iz<nz; iz++) {
53 for (ix=0; ix<nx; ix++) {
54 um[ix][iz]=0;
55 uo[ix][iz]=0;
56 up[ix][iz]=0;
57 ud[ix][iz]=0;
58 }
59 }
60
61 /* MAIN LOOP */
62 if(verb) fprintf(stderr,"\n");
63 for (it=0; it<nt; it++) {
64 if(verb) fprintf(stderr,"\b\b\b\b\b%d",it);
65
66 /* 4th order laplacian */
67 for (iz=2; iz<nz-2; iz++) {
68 for (ix=2; ix<nx-2; ix++) {
69 ud[ix][iz] =
70 c0* uo[ix ][iz ] * (idx+idz) +
71 c1*(uo[ix-1][iz ] + uo[ix+1][iz ])*idx +
72 c2*(uo[ix-2][iz ] + uo[ix+2][iz ])*idx +
73 c1*(uo[ix ][iz-1] + uo[ix ][iz+1])*idz +
74 c2*(uo[ix ][iz-2] + uo[ix ][iz+2])*idz;
75 }
76 }
77
78 /* inject wavelet */
79 for (iz=0; iz<nz; iz++) {
80 for (ix=0; ix<nx; ix++) {
81 ud[ix][iz] -= ww[it] * rr[ix][iz];
82 }
83 }
84
85 /* scale by velocity */
86 for (iz=0; iz<nz; iz++) {
87 for (ix=0; ix<nx; ix++) {
88 ud[ix][iz] *= vv[ix][iz]*vv[ix][iz];
89 }
90 }
91
92 /* time step */
93 for (iz=0; iz<nz; iz++) {
94 for (ix=0; ix<nx; ix++) {
95 up[ix][iz] = 2*uo[ix][iz]
96 - um[ix][iz]
97 + ud[ix][iz] * dt2;
98
99 um[ix][iz] = uo[ix][iz];
100 uo[ix][iz] = up[ix][iz];
101 }
102 }
103
104 /* write wavefield to output */
105 sf_floatwrite(uo[0],nz*nx,Fo);
106 }
107 if(verb) fprintf(stderr,"\n");
108 sf_close()
109 exit(0);
110 }
2011-07-02