int main (int argc, char *argv[]) 
{    
    int m, k, l, seed;
    sf_file  bshuffle,ashuffle;
    sf_axis ax,at,ap,av;
    int nx,nt,np,nv,iteration,*a1, *a2;
    float ***bsh, ***bsh2;
    sf_init(argc,argv);

    bshuffle=sf_input("in");
    ashuffle=sf_output("out");

    if (!sf_getint("iteration",&iteration)) iteration=1;
    if (!sf_getint("seed",&seed)) seed=2012;
    at=sf_iaxa( bshuffle,1); nt=sf_n(at);
    ax=sf_iaxa( bshuffle,3); nx=sf_n(ax);
    ap=sf_iaxa( bshuffle,2); np=sf_n(ap);
    av=sf_iaxa( bshuffle,4); nv=sf_n(av);
    bsh=sf_floatalloc3(nt,np,nx);
    bsh2=sf_floatalloc3(nt,np,nx);
    a1=sf_intalloc(np);
    a2=sf_intalloc(np);
    sf_warning("ntpx=%d",nt*np*nx);

    srand(seed);


    for (m=0; m<np; m++) {
	a1[m]=rand();
    }

    bubble(a1, a2, np);


    for (k=0; k<nv; k++) {
	sf_floatread(bsh[0][0],nt*np*nx,bshuffle);
	for(l=0; l<nx; l++) {
	    for (m=0; m<np; m++) {
		memcpy(bsh2[l][m], bsh[l][a2[m]], nt*sizeof(float));
	    }
	}
	sf_floatwrite(bsh2[0][0],nt*np*nx,ashuffle);
    }

    free(bsh[0][0]);
    free(bsh[0]);
    free(bsh);
    free(bsh2[0][0]);
    free(bsh2[0]);
    free(bsh2);
    free(a1);
    free(a2);
    exit (0);
}

int main(int argc, char* argv[])
{
    sf_init(argc,argv);
    sf_file Fin,Fout;
    long long int n1,n2,n3;
    float d1,o1,d2,o2;
    Fin=sf_input ("in" );
    sf_axis a1,a2,a3;
    a1 = sf_iaxa(Fin,1); n1 = sf_n(a1); d1=sf_d(a1);o1=sf_o(a1); 
    a2=sf_iaxa(Fin,2);n2=sf_n(a2);d2=sf_d(a2);o2=sf_o(a2);
    fprintf(stderr,"%lld %lld %f %f %f %f\n",n1,n2,d1,d2,o1,o2);
    int maxoffset;
    if (!sf_getint("maxoffset",&maxoffset)) maxoffset=100000000;
        
    Fout=sf_output("out");
    sf_putint(Fout,"n2",15);
    float temp[15];
    int i,j;
    int ntrace=0;
    for(i=0;i<n2;i++)
    for(j=0;j<n2;j++)
    {
        if((abs(j-i)%2==0)&&(fabs(j-i)*d2)<=maxoffset)
        {
            ntrace++;
        }
    }
    sf_putint(Fout,"n1",ntrace);
    ntrace=0;
    for(i=0;i<n2;i++)
    for(j=0;j<n2;j++)
    {
        if((abs(j-i)%2==0)&&(fabs(j-i)*d2)<=maxoffset)
        {
            float s=o2+i*d2;
            float r=o2+j*d2;
            float h=(r-s)/2.0;
            float x=(s+r)/2;
            temp[0]=i;
            temp[1]=s;
            temp[2]=r;
            temp[3]=0;
            temp[4]=0;
            temp[5]=x;
            temp[6]=h;
	        temp[14]=ntrace;
	        ntrace++;
            fprintf(stderr,"temp[0]=%f %f %f %f %f ntrace=%f\n",temp[0],s,r,h,x,temp[14]);
            sf_floatwrite(temp,15,Fout);
        }   
    }
    sf_close();
    return 0;
}

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

    sf_file Fu,Fw; /* I/O files */
    sf_axis a1,a2,a3; /* cube axes */

    float ***uu=NULL, ***ww=NULL;
    int *ii,ik;

    int nh1,nh2,nh3;
    int nb1,nb2,nb3;

    int ih1,ih2,ih3;
    int ib1,ib2,ib3;
    int  n3;
    int  i3;
    int  j1, j2, j3;
    int  k1, k2, k3;

    int ompchunk;

    int lo1,hi1;
    int lo2,hi2;
    int lo3,hi3;

    /*------------------------------------------------------------*/
    /* init RSF */
    sf_init(argc,argv);

    if(! sf_getint("ompchunk",&ompchunk)) ompchunk=1;  /* OpenMP data chunk size */
    if(! sf_getbool("verb",&verb)) verb=false;         /* verbosity flag */

    Fu = sf_input ("in" ); /*  input field */
    Fw = sf_output("out"); /* wigner distribution */

    /* read axes */
    a1=sf_iaxa(Fu,1); sf_setlabel(a1,"a1"); if(verb) sf_raxa(a1);
    a2=sf_iaxa(Fu,2); sf_setlabel(a2,"a2"); if(verb) sf_raxa(a2);
    a3=sf_iaxa(Fu,3); sf_setlabel(a3,"a3"); if(verb) sf_raxa(a3);
    n3 = sf_n(a3);

    if(! sf_getint("nh1",&nh1)) nh1=0;
    if(! sf_getint("nh2",&nh2)) nh2=0;
    if(! sf_getint("nh3",&nh3)) nh3=0;
    if(n3<=1) nh3=0;
    if(verb) sf_warning("nh1=%d nh2=%d nh3=%d",2*nh1+1,2*nh2+1,2*nh3+1);

    nb1 = sf_n(a1);
    nb2 = sf_n(a2);
    nb3=2*nh3+1;

    uu=sf_floatalloc3(nb1,nb2,nb3);
    ww=sf_floatalloc3(nb1,nb2,nb3); 

    ii = sf_intalloc(nb3);
    for(ib3=0;ib3<nb3;ib3++) {
	ii[ib3]=ib3;
    }

    /*------------------------------------------------------------*/
    /* low end on axis 3 */
    /*------------------------------------------------------------*/
    for(         ib3=0; ib3<nb3; ib3++) {
	for(     ib2=0; ib2<nb2; ib2++) {
	    for( ib1=0; ib1<nb1; ib1++) {
		ww[ib3][ib2][ib1] = 0.;
	    }
	}
    }

    sf_floatread(uu[0][0],nb1*nb2*nb3,Fu);

    for(        ih3=-nh3; ih3<nh3+1; ih3++) { lo3=SF_ABS(ih3); hi3=nb3-lo3;
	for(    ih2=-nh2; ih2<nh2+1; ih2++) { lo2=SF_ABS(ih2); hi2=nb2-lo2;
	    for(ih1=-nh1; ih1<nh1+1; ih1++) { lo1=SF_ABS(ih1); hi1=nb1-lo1;
		
		for(        ib3=lo3; ib3<nh3+1;ib3++) { j3=ii[ib3-ih3]; k3=ii[ib3+ih3];
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic,ompchunk) \
    private(ib1,ib2,j2,j1,k2,k1) \
    shared(ib3,j3,k3,ih2,ih1,lo2,hi2,lo1,hi1,uu,ww)
#endif
		    for(    ib2=lo2; ib2<hi2;  ib2++) { j2=ib2-ih2; k2=ib2+ih2;
			for(ib1=lo1; ib1<hi1;  ib1++) { j1=ib1-ih1; k1=ib1+ih1;
			    
			    ww[ib3][ib2][ib1] += uu[j3][j2][j1] 
				*                uu[k3][k2][k1];
			} /* nb1  */
		    }     /* nb2  */
		}         /* nb3  */

	    } /* nh1 */
	}     /* nh2 */
    }         /* nh3 */
    
    for(ih3=0;ih3<nh3+1;ih3++) {
	sf_floatwrite(ww[ih3][0],nb1*nb2,Fw);
    }

//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (!sf_getbool("verb",&verb))  verb=false; /* Verbosity flag */
    if (!sf_getbool("snap",&snap))  snap=false; /* Wavefield snapshots flag */
    if (!sf_getbool("expl",&expl))  expl=false; /* Multiple sources, one wvlt*/
    if (!sf_getbool("dabc",&dabc))  dabc=false; /* Absorbing BC */
    if (!sf_getbool("cden",&cden))  cden=false; /* Constant density */
    if (!sf_getbool("adj",&adj))    adj=false; /* adjoint flag */

    if (!sf_getbool("free",&fsrf) && !sf_getbool("fsrf",&fsrf)) fsrf=false; /* Free surface flag */

    if (!sf_getint("nbell",&nbell)) nbell=5; /* gaussian for source injection */

    if (!sf_getbool("optfd",&optfd))  optfd=false; /* optimized FD coefficients flag */
    if (!sf_getint("fdorder",&fdorder))  fdorder=4; /* spatial FD order */
    if (!sf_getbool("hybridbc",&hybrid))  hybrid=false;  /* hybrid Absorbing BC */
    if (!sf_getbool("sinc",&sinc)) sinc=false; /* sinc source injection */
  
    /* Initialize variables */
    file_wav = sf_input("in"); /* wavelet */
    file_vel = sf_input("vel"); /* velocity */ 
    file_src = sf_input("sou"); /* sources */
    file_rec = sf_input("rec"); /* receivers */
    file_dat = sf_output("out"); /* data */

    if (snap)  file_wfl = sf_output("wfl"); /* wavefield */
    if (!cden) {
	if (sf_getstring("cden")) {
	    file_den = sf_input ("den"); /* density */
	} else {
	    cden = true;
	    if (verb) sf_warning("No density file provided, running with constant density");
	}
    }
  
    at = sf_iaxa(file_wav,2); sf_setlabel(at,"t"); if(verb) sf_raxa(at); /* time */
    az = sf_iaxa(file_vel,1); sf_setlabel(az,"z"); if(verb) sf_raxa(az); /* depth */
    ax = sf_iaxa(file_vel,2); sf_setlabel(ax,"x"); if(verb) sf_raxa(ax); /* space */
    ay = sf_iaxa(file_vel,3); sf_setlabel(ay,"y"); if(verb) sf_raxa(ay); /* space */

    as = sf_iaxa(file_src,2); sf_setlabel(as,"s"); if(verb) sf_raxa(as); /* sources */
    ar = sf_iaxa(file_rec,2); sf_setlabel(ar,"r"); if(verb) sf_raxa(ar); /* receivers */

    nt = sf_n(at); dt = sf_d(at); 
    nz = sf_n(az); dz = sf_d(az); 
    nx = sf_n(ax); dx = sf_d(ax); 
    ny = sf_n(ay); dy = sf_d(ay); 
    ns = sf_n(as);
    nr = sf_n(ar);

    /* other execution parameters */
    if (snap) {
	if (!sf_getint("jsnap",&jsnap))  jsnap=nt;
	/* # of t steps at which to save wavefield */
    }
    if (!sf_getint("jdata",&jdata)) jdata=1;
    /* # of t steps at which to save receiver data */

    /* setup output data header */
    sf_oaxa(file_dat,ar,1);
    sf_setn(at,(nt-1)/jdata+1);
    sf_setd(at,dt*jdata);
    sf_oaxa(file_dat,at,2);

    /* wavefield cut params */
    /* setup output wavefield header */
    if (snap) {
	if (!sf_getint  ("nqz",&nqz)) nqz=sf_n(az); /* Saved wfld window nz */

int main(int argc, char* argv[])
{
    bool verb;  /* verbosity flag */
    float clip; /* threshold (clip value) */
    sf_file Fc; /* cube file */
    sf_file Fl; /* list file */
    extern int fseeko(FILE *stream, off_t offset, int whence);

    sf_axis   ax,ay,az;
    sf_axis   aa;
    int   ix,iy,iz;
    int   nx,ny,nz,nj,na;
    int   nk=0,jk;

    float **cube;
    float dx,dy,dz;
    float x0,y0,z0;

    FILE* tfile;
    char* tname;
    float t2[3],t3[4];

    /*------------------------------------------------------------*/

    /* init RSF */
    sf_init(argc,argv);
    if(! sf_getbool("verb",&verb)) verb=false;
    if(! sf_getfloat("clip",&clip)) clip=0;

    Fc = sf_input ( "in"); /* input  cube */
    Fl = sf_output("out"); /* output list */

    /* read axes*/
    az=sf_iaxa(Fc,1);
    sf_setlabel(az,"z");
    ax=sf_iaxa(Fc,2);
    sf_setlabel(ax,"x");
    ay=sf_iaxa(Fc,3);
    sf_setlabel(ay,"y");

    nz=sf_n(az);
    z0=sf_o(az);
    dz=sf_d(az);
    nx=sf_n(ax);
    x0=sf_o(ax);
    dx=sf_d(ax);
    ny=sf_n(ay);
    y0=sf_o(ay);
    dy=sf_d(ay);

    na=0;
    if(ny>1) {
        if(verb) sf_warning("initiating 3D points");
        nj=4;
    } else {
        if(verb) sf_warning("initiating 2D points");
        nj=3;
    }
    /*------------------------------------------------------------*/

    cube = sf_floatalloc2(nz,nx);

    tfile = sf_tempfile(&(tname), "w+b");

    for (iy=0; iy<ny; iy++) {
        /*	if(verb) sf_warning("iy=%d",iy);*/

        sf_floatread(cube[0],nz*nx,Fc);

        nk=0;
        for (ix=0; ix<nx; ix++) {
            for (iz=0; iz<nz; iz++) {
                if( fabs(cube[ix][iz]) > clip) {
                    nk++;
                }
            }
        }

        if(ny>1) {
            jk=0;
            for (ix=0; ix<nx; ix++) {
                for (iz=0; iz<nz; iz++) {
                    if( fabs(cube[ix][iz]) > clip) {
                        t3[0] = x0 + ix * dx;
                        t3[1] = y0 + iy * dy;
                        t3[2] = z0 + iz * dz;
                        t3[3] = cube[ix][iz];

                        fseeko(tfile,jk*4*sizeof(float),SEEK_SET);
                        fwrite(   t3,     sizeof(float),4,tfile);
                        jk++;
                    }
                }
            }
        } else {
            jk=0;
            for (ix=0; ix<nx; ix++) {
                for (iz=0; iz<nz; iz++) {
                    if( fabs(cube[ix][iz]) > clip) {
                        t2[0] = x0 + ix * dx;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool verb,correct;
    int nt,nz,nx,m2,nk,nzx,nz2,nx2,n2,pad1;
    int snap, wfnt;
    float dt, wfdt;
    char *mode;

    sf_complex **ini, **cc, **dd;
    
    sf_file Fi,Fo,Fs,Fa,Fb;    /* I/O files */
    sf_axis az,ax,at;    /* cube axes */

    sf_complex **lt, **rt,***wvfld,*alpha,*beta;
    sf_file left, right;

    sf_init(argc,argv);
    if(!sf_getbool("verb",&verb)) verb=true; /* verbosity */
    if(!sf_getint("nt",&nt)) sf_error("Need nt!"); 
    if(!sf_getfloat("dt",&dt)) sf_error("Need dt!");
    if(!sf_getint("snap",&snap)) snap=0;     /* interval for snapshots */
    if(!sf_getbool("correct",&correct)) correct=false; /*jingwei's correction*/
    mode=sf_getstring("mode"); /* default mode is pspi */
    
    if (mode[0]=='p') {sf_warning(">>>>> Using PSPI+PSPI! <<<<< \n");}
    else if (mode[0]=='n') {sf_warning(">>>>> Using NSPS+NSPS! <<<<< \n");}
    else if (mode[0]=='m') {sf_warning(">>>>> Using NSPS+PSPI! <<<<< \n");}
    else if (mode[0]=='x') {sf_warning(">>>>> Using PSPI+NSPS! <<<<< \n");}
    else sf_error("Specify mode!");

    /* setup I/O files */
    Fi = sf_input ("in" );
    Fo = sf_output("out");
    sf_settype(Fo,SF_COMPLEX);

    /* Read/Write axes */
    az = sf_iaxa(Fi,1); nz = sf_n(az); 
    ax = sf_iaxa(Fi,2); nx = sf_n(ax); 
    at = sf_iaxa(Fi,3);

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

    if (snap > 0) {
      Fs = sf_output("snaps");	/* (optional) snapshot file */
      wfnt = (int)(nt-1)/snap + 1;
      wfdt = dt*snap;
      sf_oaxa(Fs,az,1); 
      sf_oaxa(Fs,ax,2); 
      sf_setn(at,wfnt);
      sf_setd(at,wfdt);
      sf_oaxa(Fs,at,3);
      sf_settype(Fs,SF_COMPLEX);
    } else Fs=NULL;

    if (!sf_getint("pad1",&pad1)) pad1=1; /* padding factor on the first axis */

    nz2 = kiss_fft_next_fast_size(nz*pad1);
    nx2 = kiss_fft_next_fast_size(nx);
    nk = nz2*nx2; /*wavenumber*/
    nzx = nz*nx;

    /* propagator matrices */
    left = sf_input("left");
    right = sf_input("right");

    if (!sf_histint(left,"n1",&n2) || n2 != nzx) sf_error("Need n1=%d in left",nzx);
    if (!sf_histint(left,"n2",&m2))  sf_error("Need n2= in left");
    
    if (!sf_histint(right,"n1",&n2) || n2 != m2) sf_error("Need n1=%d in right",m2);
    if (!sf_histint(right,"n2",&n2) || n2 != nk) sf_error("Need n2=%d in right",nk);
  
    lt = sf_complexalloc2(nzx,m2);
    rt = sf_complexalloc2(m2,nk);

    sf_complexread(lt[0],nzx*m2,left);
    sf_complexread(rt[0],m2*nk,right);

    sf_fileclose(left);
    sf_fileclose(right);

    ini=sf_complexalloc2(nz,nx);
    cc=sf_complexalloc2(nz,nx);
    dd=sf_complexalloc2(nz,nx);

    sf_complexread(ini[0],nzx,Fi);
    sf_fileclose(Fi);

    if (snap>0) wvfld = sf_complexalloc3(nz,nx,wfnt);
    else wvfld = NULL;

    if (correct) {
	Fa = sf_input("alpha");
	Fb = sf_input("beta");
	
	if (!sf_histint(Fa,"n1",&n2) || n2 != nzx) sf_error("Need n1=%d in alpha",nzx);
	if (!sf_histint(Fb,"n1",&n2) || n2 != nk) sf_error("Need n1=%d in beta",nk);
	
	alpha = sf_complexalloc(nzx);
	beta = sf_complexalloc(nk);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
	int	ix, iz, jx, jz,ixx,izz,ixf,izf,i,j,im, jm,nx,nz,nxf,nzf,nxpad,nzpad,it,ii,jj;
	float   kxmax,kzmax;

        float   f0, t, t0, dx, dz, dxf, dzf,dt, dkx, dkz, dt2;
        int     mm, nvx, nvz, ns;
        int     hnkx, hnkz, nkx, nkz, nxz, nkxz;
        int     hnkx1, hnkz1, nkx1, nkz1;
        int     isx, isz, isxm, iszm; /*source location */

        int     itaper; /* tapering or not for spectrum of oprtator*/
        int     nstep;            /* every nstep in spatial grids to calculate filters sparsely*/

        float   *coeff_1dx, *coeff_1dz, *coeff_2dx, *coeff_2dz; /* finite-difference coefficient */

        float **apx, **apz, **apxx, **apzz;        /* polarization operator of P-wave for a location */
        float **apxs, **apzs, **apxxs, **apzzs;    /* polarization operator of SV-wave for a location */

        float ****ex, ****ez;                      /* operator for whole model for P-wave*/
        float ****exs, ****ezs;                    /* operator for whole model for SV-wave*/
        float **exx, **ezz;                        /* operator for constant model for P-wave*/
        float **exxs, **ezzs;                      /* operator for constant model for SV-wave*/

        float **vp0, **vs0, **epsi, **del, **theta;         /* velocity model */
        float **p1, **p2, **p3, **q1, **q2, **q3, **p3c, **q3c, **sum;  /* wavefield array */

        float *kx, *kz, *kkx, *kkz, *kx2, *kz2, **taper;

        clock_t t1, t2, t3, t4;
        float   timespent; 
        float   A, fx, fz;

        int     isep=1;
        int     ihomo=1;

        char    *tapertype;

	double  vp2, vs2, ep2, de2, the;

        sf_init(argc,argv);

        sf_file Fo1, Fo2, Fo3, Fo4, Fo5, Fo6, Fo7, Fo8, Fo9, Fo10, Fo11, Fo12;
       
        t1=clock();
 
        /*  wavelet parameter for source definition */
        f0=30.0;                  
        t0=0.04;                  
        A=1.0;                  

        /* time samping paramter */
        if (!sf_getint("ns",&ns)) ns=301;
        if (!sf_getfloat("dt",&dt)) dt=0.001;
        if (!sf_getint("isep",&isep)) isep=0;             /* if isep=1, separate wave-modes */ 
        if (!sf_getint("ihomo",&ihomo)) ihomo=0;          /* if ihomo=1, homogeneous medium */
        if (NULL== (tapertype=sf_getstring("tapertype"))) tapertype="D"; /* taper type*/
        if (!sf_getint("nstep",&nstep)) nstep=1; /* grid step to calculate operators: 1<=nstep<=5 */

        sf_warning("isep=%d",isep);
        sf_warning("ihomo=%d",ihomo);
        sf_warning("tapertype=%s",tapertype);
        sf_warning("nstep=%d",nstep);

        sf_warning("ns=%d dt=%f",ns,dt);
        sf_warning("read velocity model parameters");

        /* setup I/O files */
        sf_file Fvp0, Fvs0, Feps, Fdel, Fthe;

        Fvp0 = sf_input ("in");  /* vp0 using standard input */
        Fvs0 = sf_input ("vs0");  /* vs0 */
        Feps = sf_input ("epsi");  /* epsi */
        Fdel = sf_input ("del");  /* delta */
        Fthe = sf_input ("the");  /* theta */

        /* Read/Write axes */
        sf_axis az, ax;
        az = sf_iaxa(Fvp0,1); nvz = sf_n(az); dz = sf_d(az)*1000.0;
        ax = sf_iaxa(Fvp0,2); nvx = sf_n(ax); dx = sf_d(ax)*1000.0;
        fx=sf_o(ax)*1000.0;
        fz=sf_o(az)*1000.0;

        /* source definition */
        isx=nvx/2;
        isz=nvz/2;
        //isz=nvz*2/5;

        /* wave modeling space */
	nx=nvx;
	nz=nvz;
        nxpad=nx+2*m;
        nzpad=nz+2*m;

        sf_warning("fx=%f fz=%f dx=%f dz=%f",fx,fz,dx,dz);

        sf_warning("nx=%d nz=%d nxpad=%d nzpad=%d", nx,nz,nxpad,nzpad);

	vp0=sf_floatalloc2(nz,nx);	
	vs0=sf_floatalloc2(nz,nx);	
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
     	int	i,j,im,jm,nx,nz,nxpad,nzpad,it,ii,jj;

        float   f0, t, t0, dx, dz,dt, dt2;
        int     mm, nvx, nvz, ns;
        int     isx, isz, isxm, iszm; /*source location */

        float   *coeff_1dx, *coeff_1dz, *coeff_2dx, *coeff_2dz; /* finite-difference coefficient */

        float **vp0, **vs0, **epsi, **del, **theta;         /* velocity model */
        float **p1, **p2, **p3, **q1, **q2, **q3;  /* wavefield array */

        float   A, fx, fz;
		clock_t t1, t2, t3;

        sf_init(argc,argv);

        sf_file Fo1, Fo2;
		float timespent;
       
		t1 = clock();

        /*  wavelet parameter for source definition */
        f0=30.0;                  
        t0=0.04;                  
        A=1.0;                  

        /* time samping paramter */
        if (!sf_getint("ns",&ns)) ns=301;
        if (!sf_getfloat("dt",&dt)) dt=0.001;

        sf_warning("ns=%d dt=%f",ns,dt);
        sf_warning("read velocity model parameters");

        /* setup I/O files */
        sf_file Fvp0, Fvs0, Feps, Fdel, Fthe;

        Fvp0 = sf_input ("in");  /* vp0 using standard input */
        Fvs0 = sf_input ("vs0");  /* vs0 */
        Feps = sf_input ("epsi");  /* epsi */
        Fdel = sf_input ("del");  /* delta */
        Fthe = sf_input ("the");  /* theta */

        /* Read/Write axes */
        sf_axis az, ax;
        az = sf_iaxa(Fvp0,1); nvz = sf_n(az); dz = sf_d(az)*1000.0;
        ax = sf_iaxa(Fvp0,2); nvx = sf_n(ax); dx = sf_d(ax)*1000.0;
        fx=sf_o(ax)*1000.0;
        fz=sf_o(az)*1000.0;

        /* source definition */
        isx=nvx/2;
        isz=nvz/2;
        //isz=nvz*2/5;

        /* wave modeling space */
	nx=nvx;
	nz=nvz;
        nxpad=nx+2*_m;
        nzpad=nz+2*_m;

        sf_warning("fx=%f fz=%f dx=%f dz=%f",fx,fz,dx,dz);

        sf_warning("nx=%d nz=%d nxpad=%d nzpad=%d", nx,nz,nxpad,nzpad);

	vp0=sf_floatalloc2(nz,nx);	
	vs0=sf_floatalloc2(nz,nx);	
	epsi=sf_floatalloc2(nz,nx);	
	del=sf_floatalloc2(nz,nx);	
	theta=sf_floatalloc2(nz,nx);	

        int nxz=nx*nz;
        mm=2*_m+1;

        dt2=dt*dt;
        isxm=isx+_m;  /* source's x location */
        iszm=isz+_m;  /* source's z-location */

        /* read velocity model */
        sf_floatread(vp0[0],nxz,Fvp0);
        sf_floatread(vs0[0],nxz,Fvs0);
        sf_floatread(epsi[0],nxz,Feps);
        sf_floatread(del[0],nxz,Fdel);
        sf_floatread(theta[0],nxz,Fthe);

        for(i=0;i<nx;i++)
        for(j=0;j<nz;j++)
           theta[i][j] *= SF_PI/180.0;

        Fo1 = sf_output("out"); /* Elastic-wave x-component */
        Fo2 = sf_output("Elasticz"); /* Elastic-wave z-component */
        /* setup I/O files */
        puthead3(Fo1, nz, nx, 1, dz/1000.0, dx/1000.0, dt, fz/1000.0, fx/1000.0, 0.0);
        puthead3(Fo2, nz, nx, 1, dz/1000.0, dx/1000.0, dt, fz/1000.0, fx/1000.0, 0.0);

	/****************begin to calculate wavefield****************/
	/****************begin to calculate wavefield****************/
        sf_warning("==================================================");
        sf_warning("==      Porpagation Using Elastic Wave Eq.      ==");
//.........这里部分代码省略.........

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

    fint1 sft;
    int  ext;

    float a,n,f,da,a0,t0,dt,s;
    int   fint,na,nx,nz,nt;

    sf_axis ax,az,at,aa;
    int ix,iz,it,ia;

    float   **stk=NULL, **ang=NULL, *tmp=NULL, *vel=NULL;
    sf_file  Fstk=NULL,  Fang=NULL, velocity=NULL;

    sf_init (argc,argv);

    /*------------------------------------------------------------*/
    Fstk = sf_input("in");
    Fang = sf_output("out");

    if (SF_FLOAT != sf_gettype(Fstk)) sf_error("Need float input");

    az=sf_iaxa(Fstk,1); nz=sf_n(az);
    at=sf_iaxa(Fstk,2); nt=sf_n(at); t0=sf_o(at); dt=sf_d(at);
    ax=sf_iaxa(Fstk,3); nx=sf_n(ax);

    if (!sf_getint  ("na",&na)) na=nt;          /* number of angles*/
    if (!sf_getfloat("da",&da)) da=90/(nt-1);   /* angle sampling */
    if (!sf_getfloat("a0",&a0)) a0=0.;          /* angle origin */
    aa = sf_maxa(na,a0,da);
    sf_oaxa(Fang,aa,2);

    if (!sf_getint("extend",&ext)) ext=4;       /* tmp extension */
    /*------------------------------------------------------------*/

    if (!sf_getbool("top",&top)) top=false;     /* velocity scaling option */

    if (top) {
	velocity = sf_input("velocity");
	vel = sf_floatalloc(nz);
    } else {
	velocity = NULL;
	vel = NULL;
    }

    stk = sf_floatalloc2(nz,nt);
    ang = sf_floatalloc2(nz,na);
    tmp = sf_floatalloc(nt);

    sft = fint1_init(ext,nt, 0);

    for (ix = 0; ix < nx; ix++) {
	sf_floatread(stk[0],nz*nt,Fstk);
	if (top) sf_floatread(vel,nz,velocity);

	/*------------------------------------------------------------*/
	for (iz = 0; iz < nz; iz++) {
	    for (it = 0; it < nt; it++) {
		tmp[it] = stk[it][iz];
	    }
	    fint1_set(sft,tmp);

	    for (ia=0; ia < na; ia++) {
		a = a0+ia*da;          /* ang or p */

		if (top) {
		    s = a*vel[iz];
		    if (s >= 1.) {
			n = t0 - 10.*dt;
		    } else {
			n = s/sqrtf(1.0-s*s);
		    }
		} else {
		    n = 1.0/(sinf(a/180.0*SF_PI)); /* 1/sin : no angle close to 0 */
		}

		f = (n - t0) / dt;
		fint = f;

		if (fint >= 0 && fint < nt) {
		    ang[ia][iz] = fint1_apply(sft,fint,f-fint,false);
		} else {
		    ang[ia][iz] = 0.;
		}
	    }
	}
	/*------------------------------------------------------------*/

	sf_floatwrite(ang[0],nz*na,Fang);
    }

    exit (0);
}

int main (int argc, char ** argv)
{
   /* RSF variables */
   bool verb;
   float *dat=NULL, **datf=NULL, **image=NULL, **tr=NULL, **ts=NULL, **trs=NULL;
   float **px=NULL, **pz=NULL;
   sf_file Fin=NULL, Fout=NULL, Ftt=NULL;
   sf_axis a1,a2,a1t,a2t,a3t,ax,az;
   int     nf,ntaper,i1,i2,i2t,itt;
   float   fmax,df,theta,dtheta,tg,tgtap,tmin,xs,xr,eps;
   aalias  aa;
   fslice  tabtt=NULL;
   int     n1,n2,n1t,n2t,n3t,nx,nz;
   float   o1,d1,o2,d2,o1t,d1t,o2t,d2t,o3t,d3t,ox,oz,dx,dz;

   sf_init(argc,argv);

   Fin  = sf_input  ("in" );
   Fout = sf_output ("out" );
   Ftt  = sf_input  ("ttfile" );

   /* axes */
   a1  = sf_iaxa(Fin,1);   /* data */
   a2  = sf_iaxa(Fin,2);   /* data */
   a1t = sf_iaxa(Ftt,1);   /* traveltime */
   a2t = sf_iaxa(Ftt,2);   /* traveltime */
   a3t = sf_iaxa(Ftt,3);   /* traveltime */

   o1  = sf_o(a1);  n1   = sf_n(a1);  d1  = sf_d(a1);
   o2  = sf_o(a2);  n2   = sf_n(a2);  d2  = sf_d(a2);
   o1t = sf_o(a1t); n1t  = sf_n(a1t); d1t = sf_d(a1t);
   o2t = sf_o(a2t); n2t  = sf_n(a2t); d2t = sf_d(a2t);
   o3t = sf_o(a3t); n3t  = sf_n(a3t); d3t = sf_d(a3t);

   /* migration parameters */
   if(! sf_getbool("verb",&verb))     verb = false; /* verbosity flag */
   if(! sf_getfloat("theta",&theta))  theta = 30.;  /* maximum dip */
   if(! sf_getfloat("dtheta",&dtheta))  dtheta = theta/3;  /* taper zone */
   if(dtheta>theta) dtheta=theta;
   if(! sf_getfloat("df",&df))        df = 5.;    /* anti-aliasing sampling */
   if(!sf_getfloat("fmax",&fmax)) fmax=.5/d1;
   if(fmax>(.5/d1)) fmax=.5/d1;
   if (!sf_getint("ntaper",&ntaper)) ntaper=11;
   if(!sf_getfloat("tmin",&tmin)) tmin=3*d1;
   if(!sf_getfloat("xs",&xs)) sf_error("missing xs parameter\n");

   /* image parameters */
   if (!sf_getint("nx",&nx))   nx=n2t;
   if(!sf_getfloat("ox",&ox))  ox=o2t;
   if(!sf_getfloat("dx",&dx))  dx=d2t;
   if (!sf_getint("nz",&nz))   nz=n1t;
   if(!sf_getfloat("oz",&oz))  oz=o1t;
   if(!sf_getfloat("dz",&dz))  dz=d1t;

   /* checking dimensions */
   if((dx!=d2t)||(dz!=d1t)) 
     sf_error("sampling interval have to be the same in:\n"
	      " image and traveltime file\n");
   if(ox<o2t) ox=o2t; 
   if(oz<o1t) oz=o1t;
   if((ox+(nx-1)*dx)>(o2t+(n2t-1)*d2t)) nx=floor(((o2t+(n2t-1)*d2t)-ox)/dx)+1;
   if((oz+(nz-1)*dz)>(o1t+(n1t-1)*d1t)) nz=floor(((o1t+(n1t-1)*d1t)-oz)/dz)+1;

   /* output axis */
   ax = sf_maxa(nx,ox,dx); if(verb) sf_raxa(ax);
   az = sf_maxa(nz,oz,dz); if(verb) sf_raxa(az);
   sf_oaxa(Fout,az,1);
   sf_oaxa(Fout,ax,2);

   /* anti-aliasing */
   /* df = fmax; */
   nf = initAalias(-1,verb,fmax,df,n1,d1,&aa);
   /* fprintf(stderr,"forcing nf=%d df=%f\n",nf,df); */

   /* aperture angle */
   tg    = tan(SF_PI*theta/180);
   tgtap = tan(SF_PI*(theta-dtheta)/180);
   if(verb) sf_warning("tgmax=%f tgtap=%f",tg,tgtap);

   /* allocating */
   dat   = sf_floatalloc(n1);
   image = sf_floatalloc2(nz,nx);
   datf  = sf_floatalloc2 (n1,nf);
   ts    = sf_floatalloc2(n1t,n2t);
   tr    = sf_floatalloc2(n1t,n2t);
   trs   = sf_floatalloc2(n1t,n2t);
   px    = sf_floatalloc2(n1t,n2t);
   pz    = sf_floatalloc2(n1t,n2t);

   if(verb) sf_warning("initializing traveltime loading");
   /* initializing traveltime maps */
   tabtt = fslice_init(n1t*n2t,n3t,sizeof(float));
   fslice_load(Ftt,tabtt,SF_FLOAT);
   if(verb) sf_warning("traveltime loading has finished");

   /* reading the source-slice from traveltime table */
   eps = .01*d2;
   itt = floor((xs+eps-o3t)/d3t);
   fslice_get(tabtt,itt,ts[0]);
   if(verb) sf_warning("traveltime table from source was read");
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    /* Laplacian coefficients */
    float c0=-30./12.,c1=+16./12.,c2=- 1./12.;
 
    bool verb,free, ifoneway, ifsponge;           /* 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,nb;        /* 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 */

#ifdef _OPENMP
    /* Testing for OpenMP */
    double start_time, end_time;
#endif
 
    sf_init(argc,argv);

    /* OMP parameters */
#ifdef _OPENMP
    omp_init();
#endif

    if(! sf_getbool("verb",&verb)) verb=0;
    if(! sf_getbool("free",&free)) free=false;
    if(! sf_getbool("ifoneway",&ifoneway)) ifoneway=true;
    if(! sf_getbool("ifsponge",&ifsponge)) ifsponge=true;
    if(! sf_getint("nb",&nb)) nb=5;
 
    /* 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 (iz=0; iz<nz; iz++) {
	for (ix=0; ix<nx; ix++) {
	    um[ix][iz]=0;
	    uo[ix][iz]=0;
	    up[ix][iz]=0;
	    ud[ix][iz]=0;
	}
    }

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

	/* Starting timer */
#ifdef _OPENMP
    start_time = omp_get_wtime();
#endif

    for (it=0; it<nt; it++) 
{
	if(verb) fprintf(stderr,"\b\b\b\b\b%d",it);
 
	/* 4th order laplacian */
	if(ifoneway)
	{
#ifdef _OPENMP
#pragma omp parallel for default(none)						\
    private(ix,iz)						\
    shared(ud,nb,uo,c0,c1,c2,nx,nz,idx,idz)
#endif
	for (iz=nb; iz<nz-nb; iz++) {
	    for (ix=nb; ix<nx-nb; ix++) {
		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;	  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    Ffden  = sf_input("fden"); /*forward density*/
    Fbvel  = sf_input("bvel"); /*backward velocity*/
    Fbden  = sf_input("bden"); /*backward velocity*/
    if (wantrecord) {
	Frcd   = sf_input("rec");   /*record*/
    } else {
	Frcd   = sf_output("rec");   /*record*/
    }
    Fimg1  = sf_output("out");   /*Imaging*/
    Fimg2  = sf_output("img2");  /*Imaging*/

    if (wantwf) {
	Ftmpwf  = sf_output("tmpwf");/*wavefield snap*/
	Ftmpbwf = sf_output("tmpbwf");
    }

    FGx = sf_input("Gx"); 
    FGz = sf_input("Gz");
    Fsxx = sf_input("sxx");
    Fsxz = sf_input("sxz");
    Fszx = sf_input("szx");
    Fszz = sf_input("szz");
    
    if (SF_FLOAT != sf_gettype(Ffvel)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(Ffden)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(Fbvel)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(Fbden)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(Fsrc))  sf_error("Need float input");

   
    /*--- parameters of source ---*/
    srcpar srcp;
    srcp = createsrc();
    at = sf_iaxa(Fsrc, 1); nt = sf_n(at);  dt = sf_d(at);      
    if (!sf_getbool("srcdecay", &srcdecay)) srcdecay=SRCDECAY;
    /*source decay*/
    if (!sf_getint("srcrange", &srcrange)) srcrange=SRCRANGE;
    /*source decay range*/
    if (!sf_getfloat("srctrunc", &srctrunc)) srctrunc=SRCTRUNC;
    /*trunc source after srctrunc time (s)*/
    srcp->nt = nt; srcp->dt = dt; 
    srcp->decay = srcdecay; srcp->range=srcrange; srcp->trunc=srctrunc;
    loadsrc(srcp, Fsrc);

    /*--- parameters of SG LFD Coefficient ---*/
    ax = sf_iaxa(Ffvel, 2); nxb = sf_n(ax); dx = sf_d(ax); ox = sf_o(ax);
    az = sf_iaxa(Ffvel, 1); nzb = sf_n(az); dz = sf_d(az); oz = sf_o(az);

    if (!sf_histint(FGx, "n1", &nxz)) sf_error("No n1= in input");
    if (nxz != nxb*nzb) sf_error (" Need nxz = nxb*nzb");
    if (!sf_histint(FGx,"n2", &lenx)) sf_error("No n2= in input");
    if (!sf_histint(FGz,"n2", &lenz)) sf_error("No n2= in input");
    
    initsglfdcoe(nxb, nzb, lenx, lenz);
    loadcoe(nzb, nxb, FGx, FGz);
    loadschm(Fsxx, Fsxz, Fszx, Fszz);
    marg = getmarg();

    /* pml parameters */
    pmlpar pmlp;
    pmlp = creatpmlpar();
    if (!sf_getint("pmlsize", &pmlout)) pmlout=PMLOUT;
    /* size of PML layer */
    if (!sf_getint("pmld0", &pmld0)) pmld0=PMLD0;
    /* PML parameter */
    if (!sf_getbool("decay",&decay)) decay=DECAY_FLAG;

int main(int argc, char* argv[])
{
    bool verb;        
    int it,iz,im,ik,ix,i,j;     /* index variables */
    int nt,nz,nx, m2, nk, nzx, nz2, nx2, nzx2, n2, pad1;
    sf_complex c;

    float  *rr;      /* I/O arrays*/
    sf_complex *ww, *cwave, *cwavem;

    sf_complex **wave, *curr;
    float *rcurr;
    
    sf_file Fw,Fr,Fo;    /* I/O files */
    sf_axis at,az,ax;    /* cube axes */

    sf_complex **lt, **rt;
    sf_file left, right;


    sf_init(argc,argv);
    if(!sf_getbool("verb",&verb)) verb=false; /* verbosity */

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

    if (SF_COMPLEX != sf_gettype(Fw)) sf_error("Need complex input");
    if (SF_FLOAT != sf_gettype(Fr)) sf_error("Need float ref");

    sf_settype(Fo,SF_FLOAT);

    /* Read/Write axes */
    at = sf_iaxa(Fw,1); nt = sf_n(at); 
    az = sf_iaxa(Fr,1); nz = sf_n(az); 
    ax = sf_iaxa(Fr,2); nx = sf_n(ax); 

    sf_oaxa(Fo,az,1); 
    sf_oaxa(Fo,ax,2); 
    sf_oaxa(Fo,at,3);
    
    if (!sf_getint("pad1",&pad1)) pad1=1; /* padding factor on the first axis */

    nk = cfft2_init(pad1,nz,nx,&nz2,&nx2);

    nzx = nz*nx;
    nzx2 = nz2*nx2;

    /* propagator matrices */
    left = sf_input("left");
    right = sf_input("right");

    if (!sf_histint(left,"n1",&n2) || n2 != nzx) sf_error("Need n1=%d in left",nzx);
    if (!sf_histint(left,"n2",&m2))  sf_error("Need n2= in left");
    
    if (!sf_histint(right,"n1",&n2) || n2 != m2) sf_error("Need n1=%d in right",m2);
    if (!sf_histint(right,"n2",&n2) || n2 != nk) sf_error("Need n2=%d in right",nk);
//    if (!sf_histint(Fw,"n1",&nxx)) sf_error("No n1= in input");
  
    lt = sf_complexalloc2(nzx,m2);
    rt = sf_complexalloc2(m2,nk);

    sf_complexread(lt[0],nzx*m2,left);
    sf_complexread(rt[0],m2*nk,right);

//    sf_fileclose(left);
//    sf_fileclose(right);

    /* read wavelet & reflectivity */
    ww=sf_complexalloc(nt);  
    sf_complexread(ww,nt ,Fw);

    rr=sf_floatalloc(nzx); 
    sf_floatread(rr,nzx,Fr);

    curr   = sf_complexalloc(nzx2);
    rcurr  = sf_floatalloc(nzx2);

    cwave  = sf_complexalloc(nk);
    cwavem = sf_complexalloc(nk);
    wave   = sf_complexalloc2(nzx2,m2);


    for (iz=0; iz < nzx2; iz++) {
	curr[iz] = sf_cmplx(0.,0.);
	rcurr[iz]= 0.;
    }

    /* MAIN LOOP */
    for (it=0; it<nt; it++) {
	if(verb) sf_warning("it=%d;",it);

	/* matrix multiplication */
	cfft2(curr,cwave);

	for (im = 0; im < m2; im++) {
	    for (ik = 0; ik < nk; ik++) {
#ifdef SF_HAS_COMPLEX_H
		cwavem[ik] = cwave[ik]*rt[ik][im];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	sf_histfloat(Fd_v,"d1",&dt);
	sf_histint(Fd_v,"n2",&gpl_v);
      } else Fd_v = NULL;
      src = -1; n_srcs = -1;
      spx = NULL; spz = NULL;
      f0 = NULL; t0 = NULL; A = NULL;
    } else {
      Fd = NULL;
      if (!sf_getint("nt",&nt)) sf_error("Need nt!");
      if (!sf_getfloat("dt",&dt)) sf_error("Need dt!");
      if (!sf_getint("gpl",&gpl)) gpl = -1; /* geophone length */
      if (!sf_getint("gpl_v",&gpl_v)) gpl_v = -1; /* geophone height */
      if (!sf_getint("src",&src)) src=0; /* source type */
      if (!sf_getint("n_srcs",&n_srcs)) n_srcs=1; /* source type */
      spx = sf_intalloc(n_srcs);
      spz = sf_intalloc(n_srcs);
      f0  = sf_floatalloc(n_srcs);
      t0  = sf_floatalloc(n_srcs);
      A   = sf_floatalloc(n_srcs);
      if (!sf_getints("spx",spx,n_srcs)) sf_error("Need spx!"); /* shot position x */
      if (!sf_getints("spz",spz,n_srcs)) sf_error("Need spz!"); /* shot position z */
      if (!sf_getfloats("f0",f0,n_srcs)) sf_error("Need f0! (e.g. 30Hz)");   /*  wavelet peak freq */
      if (!sf_getfloats("t0",t0,n_srcs)) sf_error("Need t0! (e.g. 0.04s)");  /*  wavelet time lag */
      if (!sf_getfloats("A",A,n_srcs)) sf_error("Need A! (e.g. 1)");     /*  wavelet amplitude */
    }
    if (!sf_getint("gpx",&gpx)) gpx = -1; /* geophone position x */
    if (!sf_getint("gpz",&gpz)) gpz = -1; /* geophone position z */
    if (!sf_getint("gpx_v",&gpx_v)) gpx_v = -1; /* geophone position x */
    if (!sf_getint("gpz_v",&gpz_v)) gpz_v = -1; /* geophone position z */

    if (SF_FLOAT != sf_gettype(Fi)) sf_error("Need float input");

    /* Read/Write axes */
    az = sf_iaxa(Fi,1); nz = sf_n(az); dz = sf_d(az);
    ax = sf_iaxa(Fi,2); nx = sf_n(ax); dx = sf_d(ax);
    nz1 = nz-nbt-nbb;
    nx1 = nx-nbl-nbr;
    if (gpx==-1) gpx = nbl;
    if (gpz==-1) gpz = nbt;
    if (gpl==-1) gpl = nx1;
    if (gpx_v==-1) gpx_v = nbl;
    if (gpz_v==-1) gpz_v = nbt;
    if (gpl_v==-1) gpl_v = nz1;
    ntsnap=0;
    if (snap)
        for (it=0;it<nt;it++)
            if (it%snap==0) ntsnap++;

    if (mig) { /*output final wavefield*/
      sf_setn(az,nz1);
      sf_setn(ax,nx1);
      sf_oaxa(Fo,az,1);
      sf_oaxa(Fo,ax,2);
      sf_settype(Fo,SF_FLOAT);
    } else { /*output data*/
      sf_setn(ax,gpl);
      /*output horizontal data is mandatory*/
      sf_putint(Fo,"n1",nt);
      sf_putfloat(Fo,"d1",dt);
      sf_putfloat(Fo,"o1",0.);
      sf_putstring(Fo,"label1","Time");
      sf_putstring(Fo,"unit1","s");
      sf_oaxa(Fo,ax,2);
      sf_settype(Fo,SF_FLOAT);
      /*output vertical data is optional*/
      if (NULL!=sf_getstring("dat_v")) {

int main (int argc, char *argv[])
{
    int  adj;             /* forward or adjoint */
    int  eic;             /* EIC or CIC */
    bool verb;            /* verbosity */
    float eps;            /* dip filter constant */
    int   nrmax;          /* number of reference velocities */
    float dtmax;          /* time error */
    int   pmx,pmy;        /* padding in the k domain */
    int   tmx,tmy;        /* boundary taper size */

    int nhx, nhy, nhz, nht, nc;
    int nhx2,nhy2,nhz2,nht2;
    float dht, oht;

    sf_axis amx,amy,az;
    sf_axis alx,aly;
    sf_axis aw,ae,ac,aa;
    sf_axis ahx,ahy,ahz,aht;         

    /* I/O files */
    sf_file Bws=NULL;   /*  background wavefield file Bws */
    sf_file Bwr=NULL;   /*  background wavefield file Bwr */
    sf_file Bs=NULL;    /*  background slowness file Bs   */
    sf_file Ps=NULL;    /*  slowness perturbation file Ps */
    sf_file Pi=NULL;    /*  image perturbation file Pi    */
    sf_file Fc=NULL;    /*  CIP coordinates               */
    sf_file Pws=NULL;   /*  perturbed wavefield file Pws */
    sf_file Pwr=NULL;   /*  perturbed wavefield file Pwr */
    sf_file Pti=NULL;  

    int ompnth=1;

    wexcub3d cub; /* wavefield hypercube */
    wexcip3d cip; /* CIP gathers */
    wextap3d tap; /* tapering */
    wexssr3d ssr; /* SSR operator */
    wexlsr3d lsr; /* LSR operator */
    wexslo3d slo; /* slowness */

    wexmvaop3d mva;
    float dsmax;

    /*------------------------------------------------------------*/
    sf_init(argc,argv);

    /* OMP parameters */
#ifdef _OPENMP
    ompnth=omp_init();
#endif

    if (!sf_getbool(  "verb",&verb ))  verb =  true; /* verbosity flag */
    if (!sf_getint(   "adj",&adj   ))  sf_error("Specify adjoint!"); /* y=ADJ Back-projection; n=FWD Forward Scattering */
    if (!sf_getint(  "feic",&eic   ))  sf_error("Specify EIC!");     /* extended I.C. flag */
    if (!sf_getfloat(  "eps",&eps  ))   eps =  0.01; /* stability parameter */
    if (!sf_getint(  "nrmax",&nrmax)) nrmax =     1; /* max number of refs */
    if (!sf_getfloat("dtmax",&dtmax)) dtmax = 0.004; /* max time error */
    if (!sf_getint(    "pmx",&pmx  ))   pmx =     0; /* padding on x */
    if (!sf_getint(    "pmy",&pmy  ))   pmy =     0; /* padding on y */
    if (!sf_getint(    "tmx",&tmx  ))   tmx =     0; /* taper on x   */
    if (!sf_getint(    "tmy",&tmy  ))   tmy =     0; /* taper on y   */

    ae  = sf_maxa(1,0,1);
    nhx=nhy=nhz=nht=nc=nhx2=nhy2=nhz2=nht2=0;
    oht = dht = 0.0;
    /*------------------------------------------------------------*/
    /* slowness */

    Bs = sf_input("slo");
    alx = sf_iaxa(Bs,1); sf_setlabel(alx,"lx");
    aly = sf_iaxa(Bs,2); sf_setlabel(aly,"ly");
    az =  sf_iaxa(Bs,3); sf_setlabel(az, "z");

    /*------------------------------------------------------------*/
    /* input file */
    if(adj)
        Pi = sf_input("in");
    else
        Ps = sf_input("in");

    /*------------------------------------------------------------*/
    /* wavefield */
    Bws = sf_input("swfl");
    Bwr = sf_input("rwfl");

    amx = sf_iaxa(Bws,1); sf_