Madagascar Programming Reference Manual

C program

  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     }