C++ sf_input函数代码示例

//.........这里部分代码省略.........
    /*------------------------------------------------------------*/
    /* wavefield cut params                                       */
    /*------------------------------------------------------------*/
    sf_axis   acz=NULL,acx=NULL,acy=NULL;
    int       nqz,nqx,nqy;
    float     oqz,oqx,oqy;
    float     dqz,dqx,dqy;
    sf_complex***uc=NULL; /* tmp array for output wavefield snaps */

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

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

    /*------------------------------------------------------------*/
    /* read execution flags                                       */
    /*------------------------------------------------------------*/
    if(! sf_getbool("verb",&verb)) verb=false; /* verbosity flag */
    if(! sf_getbool("snap",&snap)) snap=false; /* wavefield snapshots flag */
    if(! sf_getbool("free",&fsrf)) fsrf=false; /* free surface flag */
    if(! sf_getbool("back",&back)) back=false; /* backward extrapolation flag (for rtm) */
    if(! sf_getbool("esou",&esou)) esou=false; /* explosive force source */

    /*------------------------------------------------------------*/
    /* I/O files                                                  */
    /*------------------------------------------------------------*/
    Fwav = sf_input ("in" ); /* wavelet   */
    Fdat = sf_output("out"); /* data      */
    Fsou = sf_input ("sou"); /* sources   */
    Frec = sf_input ("rec"); /* receivers */
    Fccc = sf_input ("ccc"); /* stiffness */
    Frkp = sf_input ("rkp"); /* app. rank */
    Fltp = sf_input ("ltp"); /* left mat  */
    Frtp = sf_input ("rtp"); /* right mat */
    Fwfp = sf_output("wfp"); /* wavefield */
    Frks = sf_input ("rks"); /* app. rank */
    Flts = sf_input ("lts"); /* left mat  */
    Frts = sf_input ("rts"); /* right mat */
    Fwfs = sf_output("wfs"); /* wavefield */

    /*------------------------------------------------------------*/
    /* axes                                                       */
    /*------------------------------------------------------------*/
    at = sf_iaxa(Fwav,4); sf_setlabel(at,"t"); if(verb) sf_raxa(at); /* time */
    az = sf_iaxa(Fccc,1); sf_setlabel(az,"z"); if(verb) sf_raxa(az); /* depth */
    ax = sf_iaxa(Fccc,2); sf_setlabel(ax,"x"); if(verb) sf_raxa(ax); /* space x */
    ay = sf_iaxa(Fccc,3); sf_setlabel(ay,"y"); if(verb) sf_raxa(ay); /* space y */

    asx = sf_iaxa(Fsou,2); sf_setlabel(asx,"sx"); if(verb) sf_raxa(asx); /* sources x */
    asy = sf_iaxa(Fsou,3); sf_setlabel(asy,"sy"); if(verb) sf_raxa(asy); /* sources y */
    arx = sf_iaxa(Frec,2); sf_setlabel(arx,"rx"); if(verb) sf_raxa(arx); /* receivers x */
    ary = sf_iaxa(Frec,3); sf_setlabel(ary,"ry"); if(verb) sf_raxa(ary); /* receivers y */

    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(asx)*sf_n(asy);

int main(int argc, char* argv[]) 
{
    int nx, nt, nkx, nkz, ix, it, ikx, ikz, nz, iz, nbt, nbb, nbl, nbr, nxb, nzb, isx, isz;
    float dt, dx, dkx, kx, dz, dkz, kz, tmpdt, pi=SF_PI, o1, o2, kx0, kz0;
    float **nxt,  **old,  **cur, **ukr, **dercur, **derold, *wav;
    float **vx, vx2, vx0, vx02, **vz, vz2, vz0, vz02, **yi, yi0, **se, se0;
    float ***aa, dx2, dz2, dx4, dz4, ct, cb, cl, cr; //top, bottom, left, right 
    float w1, w10, w2, w20, w3, w30, h1, h10, h2, h20, h3, h30;
    float cosg, cosg0, cosg2, cosg02, sing, sing0, sing2, sing02;
    float vk, vk2, tmpvk, k2, err, dt2, kx1, kz1;
    kiss_fft_cpx **uk, *ctracex, *ctracez;
    kiss_fft_cfg cfgx, cfgxi, cfgz, cfgzi;
    sf_file out, velx, velz, source, yita, seta;
    bool opt;    /* optimal padding */
   // #ifdef _OPENMP
   // int nth;
   // #endif
     

    sf_init(argc,argv);
    out = sf_output("out");
    velx = sf_input("velx");   /* velocity */
    velz = sf_input("velz");   /* velocity */
    yita = sf_input("yita");   /* anistropic parameter*/
    source = sf_input("in");   /* source wavlet*/
    seta = sf_input("seta");   /* TTI angle*/

//    if (SF_FLOAT != sf_gettype(inp)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(velx)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(velz)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(source)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(seta)) sf_error("Need float input");

    if (!sf_histint(velx,"n1",&nx)) sf_error("No n1= in input");
    if (!sf_histfloat(velx,"d1",&dx)) sf_error("No d1= in input");
    if (!sf_histint(velx,"n2",&nz)) sf_error("No n2= in input");
    if (!sf_histfloat(velx,"d2",&dz)) sf_error("No d2= in input");
    if (!sf_histfloat(velx,"o1",&o1)) o1=0.0;
    if (!sf_histfloat(velx,"o2",&o2)) o2=0.0;
  //  if (!sf_histint(inp,"n2",&nt)) sf_error("No n2= in input");
  //  if (!sf_histfloat(inp,"d2",&dt)) sf_error("No d2= in input");
    if (!sf_getbool("opt",&opt)) opt=true;
    /* if y, determine optimal size for efficiency */
    if (!sf_getfloat("dt",&dt)) sf_error("Need dt input");
    if (!sf_getint("nt",&nt)) sf_error("Need nt input");
    if (!sf_getint("isx",&isx)) sf_error("Need isx input");
    if (!sf_getint("isz",&isz)) sf_error("Need isz input");
    if (!sf_getfloat("err",&err)) err = 0.0001;

    if (!sf_getint("nbt",&nbt)) nbt=44;
    if (!sf_getint("nbb",&nbb)) nbb=44;
    if (!sf_getint("nbl",&nbl)) nbl=44;
    if (!sf_getint("nbr",&nbr)) nbr=44;

    if (!sf_getfloat("ct",&ct)) ct = 0.01; /*decaying parameter*/
    if (!sf_getfloat("cb",&cb)) cb = 0.01; /*decaying parameter*/
    if (!sf_getfloat("cl",&cl)) cl = 0.01; /*decaying parameter*/
    if (!sf_getfloat("cr",&cr)) cr = 0.01; /*decaying parameter*/

    

    sf_putint(out,"n1",nx);
    sf_putfloat(out,"d1",dx);
//    sf_putfloat(out,"o1",x0);
    sf_putint(out,"n2",nz);
    sf_putfloat(out,"d2",dz);
    sf_putint(out,"n3",nt);
    sf_putfloat(out,"d3",dt);
    sf_putfloat(out,"o1",o1); 
    sf_putfloat(out,"o2",o2); 
    sf_putfloat(out,"o3",0.0); 

    nxb = nx + nbl + nbr;
    nzb = nz + nbt + nbb;


    nkx = opt? kiss_fft_next_fast_size(nxb): nxb;
    nkz = opt? kiss_fft_next_fast_size(nzb): nzb;
    if (nkx != nxb) sf_warning("nkx padded to %d",nkx);
    if (nkz != nzb) sf_warning("nkz padded to %d",nkz);
    dkx = 1./(nkx*dx);
    kx0 = -0.5/dx;
    dkz = 1./(nkz*dz);
    kz0 = -0.5/dz;
    cfgx = kiss_fft_alloc(nkx,0,NULL,NULL);
    cfgxi = kiss_fft_alloc(nkx,1,NULL,NULL);
    cfgz = kiss_fft_alloc(nkz,0,NULL,NULL);
    cfgzi = kiss_fft_alloc(nkz,1,NULL,NULL);


    uk = (kiss_fft_cpx **) sf_complexalloc2(nkx,nkz);
    ctracex = (kiss_fft_cpx *) sf_complexalloc(nkx);
    ctracez = (kiss_fft_cpx *) sf_complexalloc(nkz);

    wav    =  sf_floatalloc(nt);
    sf_floatread(wav,nt,source);

    old    =  sf_floatalloc2(nxb,nzb);
    cur    =  sf_floatalloc2(nxb,nzb);
    nxt    =  sf_floatalloc2(nxb,nzb);
//.........这里部分代码省略.........

int main(int argc, char* argv[]) 
{
    int nx, nt, ix, it, isx;
    float dt, dx;
    float *old, *nxt, *cur, *sig, *a, *b1, *b2, *b3, *b4, *b5;
    sf_file in, out, Gmatrix, vel;
    int im,im2,im3,im4,im5,ip,ip2,ip3,ip4,ip5;

    sf_init(argc,argv);
    in  = sf_input("in");
    Gmatrix = sf_input("G");   /* velocity */
    vel = sf_input("vel");   /* velocity */
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(Gmatrix)) sf_error("Need float input");
    if (!sf_histint(vel,"n1",&nx)) sf_error("No n1= in input");
    if (!sf_histfloat(vel,"d1",&dx)) sf_error("No d1= in input");
    if (!sf_getint("isx",&isx)) isx=(int)(nx/2);
    if (!sf_getint("nt",&nt)) sf_error("No nt in input");
    if (!sf_getfloat("dt",&dt)) sf_error("No dt in input");

    sf_putint(out,"n1",nx);
    sf_putfloat(out,"d1",dx);
    sf_putint(out,"n2",nt);
    sf_putfloat(out,"d2",dt);
    sf_putfloat(out,"o2",0.0); 

    sig = sf_floatalloc(nx);
    old = sf_floatalloc(nx);
    nxt = sf_floatalloc(nx);
    cur = sf_floatalloc(nx);
    a = sf_floatalloc(nx);
    b1 = sf_floatalloc(nx);
    b2 = sf_floatalloc(nx);
    b3 = sf_floatalloc(nx);
    b4 = sf_floatalloc(nx);
    b5 = sf_floatalloc(nx);
    
 
    sf_floatread(a,nx,Gmatrix);
    sf_floatread(b1,nx,Gmatrix);
    sf_floatread(b2,nx,Gmatrix);
 /*   sf_floatread(b3,nx,Gmatrix);
    sf_floatread(b4,nx,Gmatrix);
    sf_floatread(b5,nx,Gmatrix); */
    sf_floatread(sig,nx,in);		

	/* initial conditions */
    for (ix=0; ix < nx; ix++){
        cur[ix] =  sig[ix];
        old[ix] =  0.0; 
	nxt[ix] = 0.;
    }

    /* propagation in time */
    for (it=0; it < nt; it++) {
	sf_floatwrite(cur,nx,out);
	/* Stencil */
	for (ix=0; ix < nx; ix++) {
            im = ix-1 < 0? ix-1+nx:ix-1; 
            im2 = ix-2 < 0? ix-2+nx:ix-2; 
            im3 = ix-3 < 0? ix-3+nx:ix-3; 
            im4 = ix-4 < 0? ix-4+nx:ix-4; 
            im5 = ix-5 < 0? ix-5+nx:ix-5; 
            ip = ix+1 > nx-1? ix+1-nx:ix+1;
            ip2 = ix+2 > nx-1? ix+2-nx:ix+2;
            ip3 = ix+3 > nx-1? ix+3-nx:ix+3;
            ip4 = ix+4 > nx-1? ix+4-nx:ix+4;
            ip5 = ix+5 > nx-1? ix+5-nx:ix+5;

	    nxt[ix] = ( 0.5* (cur[im]+cur[ip])*b1[ix] +  0.5*(cur[im2]+cur[ip2])*b2[ix]) 
                       - old[ix] + 2.0*cur[ix];

	    old[ix] = cur[ix];
	    cur[ix] = nxt[ix];
	}
    }
    exit(0);
}

int main (int argc, char *argv[])
{
    bool both;
    int ir, nr, n1,n2,n3,n4, m1, m2, m3, n12, n123, nw, nj1, nj2, i3;
    float *u1, *u2, *p;
    sf_file in, out, dip;
    off_t pos=0;
    allpass ap;

    sf_init(argc,argv);
    in = sf_input ("in");
    dip = sf_input ("dip");
    out = sf_output ("out");

    if (SF_FLOAT != sf_gettype(in) ||
            SF_FLOAT != sf_gettype(dip)) sf_error("Need float type");

    if (!sf_histint(in,"n1",&n1)) sf_error("Need n1= in input");
    if (!sf_histint(in,"n2",&n2)) n2=1;
    if (!sf_histint(in,"n3",&n3)) n3=1;
    n12 = n1*n2;
    n123 = n12*n3;
    nr = sf_leftsize(in,3);

    if (!sf_histint(dip,"n1",&m1) || m1 != n1)
        sf_error("Need n1=%d in dip",n1);
    if (1 != n2 && (!sf_histint(dip,"n2",&m2) || m2 != n2))
        sf_error("Need n2=%d in dip",n2);
    if (1 != n3 && (!sf_histint(dip,"n3",&m3) || m3 != n3))
        sf_error("Need n3=%d in dip",n3);

    if (!sf_getbool("both",&both)) both=false;
    /* if y, compute both left and right predictions */

    if (1 == n3) {
        n4=0;
        if (both) sf_putint(out,"n3",2);
    } else {
        if (!sf_getint("n4",&n4)) n4=2;
        /* what to compute in 3-D. 0: in-line, 1: cross-line, 2: both */

        if (n4 > 2) n4=2;
        if (2==n4) sf_putint(out,"n4",both? 4:2);

        if (0 != n4 || both) {
            sf_unpipe(in,(off_t) n123*sizeof(float));
            pos = sf_tell(in);
        }
    }

    if (!sf_getint("order",&nw)) nw=1;
    /* accuracy */
    if (!sf_getint("nj1",&nj1)) nj1=1;
    /* in-line aliasing */
    if (!sf_getint("nj2",&nj2)) nj2=1;
    /* cross-line aliasing */

    for (ir=0; ir < nr; ir++) {

        if (1 != n4) { /* in-line */
            u1 = sf_floatalloc(n12);
            u2 = sf_floatalloc(n12);
            p  = sf_floatalloc(n12);

            for (i3=0; i3 < n3; i3++) {
                /* read data */
                sf_floatread(u1,n12,in);

                /* read t-x dip */
                sf_floatread(p,n12,dip);

                ap = allpass_init (nw,nj1,n1,n2,1,p);

                /* apply */
                allpass1(false, false, ap, u1, u2);

                /* write t-x destruction */
                sf_floatwrite(u2,n12,out);
            }

            free(u1);
            free(u2);
            free(p);

        }

        if (0 != n4) { /* cross-line */
            u1 = sf_floatalloc(n123);
            u2 = sf_floatalloc(n123);
            p  = sf_floatalloc(n123);

            /* read data */
            sf_seek(in,pos,SEEK_SET);
            sf_floatread(u1,n123,in);

            /* read t-y dip */
            sf_floatread(p,n123,dip);

            ap = allpass_init(nw,nj2,n1,n2,n3,p);

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

int main(int argc, char* argv[])
{
    bool sym;
    float o1,d1, o2,d2, tol;
    int nd, n1, n2, n12, niter, rect1, rect2, nw;
    float **xy, *z, *m, *m2, *d;
    sf_file in, out, coord, pattern;

    sf_init(argc,argv);
    in = sf_input("in");
    out = sf_output("out");
    coord = sf_input("coord");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"n1",&nd)) sf_error("No n1= in input");

    if (NULL != sf_getstring("pattern")) {
	/* pattern file for output dimensions */
	pattern = sf_input("pattern");
	
	if (!sf_histint(pattern,"n1",&n1)) sf_error("No n1= in pattern");
	if (!sf_histint(pattern,"n2",&n2)) sf_error("No n2= in pattern");
	if (!sf_histfloat(pattern,"d1",&d1)) d1=1.;
	if (!sf_histfloat(pattern,"d2",&d2)) d2=1.;
	if (!sf_histfloat(pattern,"o1",&o1)) o1=0.;
	if (!sf_histfloat(pattern,"o2",&o2)) o2=0.;
	
	sf_fileclose(pattern);
    } else {
	if (!sf_getint("n1",&n1)) sf_error("Need n1=");
	if (!sf_getint("n2",&n2)) sf_error("Need n2=");
	if (!sf_getfloat("d1",&d1)) d1=1.;
	if (!sf_getfloat("d2",&d2)) d2=1.;
	if (!sf_getfloat("o1",&o1)) o1=0.;
	if (!sf_getfloat("o2",&o2)) o2=0.;
    }

    n12 = n1*n2;

    sf_putint(out,"n1",n1);
    sf_putint(out,"n2",n2);
    sf_putfloat(out,"d1",d1);
    sf_putfloat(out,"d2",d2);
    sf_putfloat(out,"o1",o1);
    sf_putfloat(out,"o2",o2);

    if (!sf_getint("niter",&niter)) niter=10;
    /* number of iterations */

    xy = sf_floatalloc2(2,nd);
    sf_floatread(xy[0],nd*2,coord);

    if (!sf_getint("rect1",&rect1)) rect1=1;
    if (!sf_getint("rect2",&rect2)) rect2=1;
    /* smoothing regularization */
    
    if (!sf_getint("nw",&nw)) nw=2;
    /* interpolator size */
    
    if (!sf_getbool("sym",&sym)) sym=false;
    /* if y, use symmetric shaping */
    if (!sf_getfloat("tol",&tol)) tol=1e-3;
    /* tolerance for stopping iteration */
    
    nnshape_init(sym,nd, n1,n2, o1,o2, d1,d2, 
		 rect1,rect2, nw, 2, xy);
    sf_gmres_init(n12,niter); 
 
    z = sf_floatalloc (n12);
    m = sf_floatalloc (n12);
    d = sf_floatalloc(nd);

    sf_floatread(d,nd,in);
 
    /* make right-hand side */
    nnshape_back(d,z);
    nnshape_smooth(z);
    /* invert */
    sf_gmres(z,m,nnshape,NULL,niter,tol,true);
    if (sym) nnshape_smooth(m);
    
    sf_floatwrite (m, n12, out);

    exit(0);
}

int main(int argc, char* argv[])
{
    int n1,n2,n3, jump, i1,i2,i3, j;
    float d, *pp=NULL, *zero=NULL, *one=NULL;
    sf_file in=NULL, out=NULL, mask=NULL;

    sf_init(argc,argv);
    in = sf_input("in");
    out = sf_output("out");
    mask = sf_output("mask");

    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);

    if (!sf_getint("jump",&jump)) jump=2;
    /* aliasing */

    if (n2 > 1) {
	sf_putint(out,"n2",n2*jump);
	sf_putint(mask,"n2",n2*jump);
	if (sf_histfloat(in,"d2",&d)) {
	    sf_putfloat(out,"d2",d/jump);
	    sf_putfloat(mask,"d2",d/jump);
	}
    }
    if (n3 > 1) {
	sf_putint(out,"n3",n3*jump);
	sf_putint(mask,"n3",n3*jump);
	if (sf_histfloat(in,"d3",&d)) {
	    sf_putfloat(out,"d3",d/jump);
	    sf_putfloat(mask,"d3",d/jump);
	}
    }

    pp = sf_floatalloc(n1);
    zero = sf_floatalloc(n1);
    one = sf_floatalloc(n1);

    for (i1=0; i1 < n1; i1++) {
	zero[i1] = 0.;
	one[i1] = 1.;
    }

    for (i3=0; i3 < n3; i3++) {
	for (i2=0; i2 < n2; i2++) {
	    sf_floatread(pp,n1,in);
	    sf_floatwrite(pp,n1,out);
	    sf_floatwrite(one,n1,mask);
	    if (n2 >1) {
		for (j=0; j < jump-1; j++) {
		    sf_floatwrite(zero,n1,out);
		    sf_floatwrite(zero,n1,mask);
		}
	    }
	}
	if (n3 > 1) {
	    for (i2=0; i2 < n2*jump; i2++) {
		sf_floatwrite(zero,n1,out);
		sf_floatwrite(zero,n1,mask);
	    }
	}
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int i, iter, niter, p[6][2], status, *mask;
    float *buf, *buf2, *wht;
    double rn, rnp, alpha, beta;
    pid_t pid[6]={1,1,1,1,1,1};
    off_t nm, nd, msiz, dsiz, pos;
    size_t nbuf, mbuf, dbuf;
    FILE *xfile, *Rfile, *gfile, *sfile, *Sfile;
    char *x, *R, *g, *s, *S, *prog;
    sf_file mod, dat, from, mwt, x0, known;  /* input */
    sf_file to, out; /* output */

    extern int fseeko(FILE *stream, off_t offset, int whence);
    extern off_t ftello (FILE *stream);

    sf_init(argc,argv);

    dat = sf_input("in");
    mod = sf_input("mod");

    if (SF_FLOAT != sf_gettype(mod) ||
	SF_FLOAT != sf_gettype(dat)) 
	sf_error("Need float type in mod and dat");

    for (i=0; i < argc-1; i++) {
	argv[i]=argv[i+1];
    }
    for (i=0; i < argc-1; i++) {	
	/* find the program to run */
	if (NULL == strchr(argv[i],'=')) {
	    /* first one without the '=' */
	    prog = argv[0];
	    argv[0] = argv[i];
	    argv[i] = prog;
	    break;
	}
    }

    argv[argc-1] = sf_charalloc(6);
    snprintf(argv[argc-1],6,"adj=X");

    if (!sf_getint("niter",&niter)) niter=1;
    /* number of iterations */

    Rfile = sf_tempfile(&R,"w+b"); 
    xfile = sf_tempfile(&x,"w+b"); 
    gfile = sf_tempfile(&g,"w+b");
    sfile = sf_tempfile(&s,"w+b");
    Sfile = sf_tempfile(&S,"w+b");

    fclose(Rfile);
    fclose(xfile);
    fclose(gfile);
    fclose(sfile);
    fclose(Sfile);

    nm = sf_filesize(mod);
    nd = sf_filesize(dat);

    /* I/O buffers */
    nbuf = BUFSIZ/sizeof(float);
    buf  = sf_floatalloc(nbuf);
    buf2 = sf_floatalloc(nbuf);

    if (NULL != sf_getstring("mwt")) {
	mwt = sf_input("mwt"); /* model weight */
	wht = sf_floatalloc(nbuf);
    } else {
	mwt = NULL;
	wht = NULL;
    }

    if (NULL != sf_getstring("known")) {
	known = sf_input("known"); /* known model mask */
	if (SF_INT != sf_gettype(known)) sf_error("Need int type in known");
	mask = sf_intalloc(nbuf);
    } else {
	known = NULL;
	mask = NULL;
    }

    if (NULL != sf_getstring("x0")) {
	x0 = sf_input("x0"); /* initial model */
    } else {
	x0 = NULL;
    }

    for (i=0; i < 6; i++) { /* make six pipes */
	if (pipe(p[i]) < 0) sf_error("pipe error:");
    }

    for (iter=0; iter < niter; iter++) {
	for (i=0; i < 6; i++) { /* fork six children */
	    if ((pid[i] = fork()) < 0) sf_error("fork error:");
	    if (0 == pid[i]) break;
	}
	
	if (0 == pid[0]) {	
	    /* feeds rr to p[0] */
//.........这里部分代码省略.........

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

    /*survey parameters*/
    int   nx, nz;
    float dx, dz;
    int   n_srcs;
    int   *spx, *spz;
    int   gpz, gpx, gpl;
    int   gpz_v, gpx_v, gpl_v;
    int   snap;
    /*fft related*/
    bool  cmplx;
    int   pad1;
    /*absorbing boundary*/
    bool abc;
    int nbt, nbb, nbl, nbr;
    float ct,cb,cl,cr;
    /*source parameters*/
    int src; /*source type*/
    int nt,ntsnap;
    float dt,*f0,*t0,*A;
    /*misc*/
    bool verb, ps, mig;
    float vref;

    pspar par;
    int nx1, nz1; /*domain of interest*/
    int it;
    float *vel,**dat,**dat_v,**wvfld,*img; /*velocity profile*/
    sf_file Fi,Fo,Fd,Fd_v,snaps; /* I/O files */
    sf_axis az,ax; /* cube axes */

    sf_init(argc,argv);

    if (!sf_getint("snap",&snap)) snap=0; /* interval for snapshots */
    if (!sf_getbool("cmplx",&cmplx)) cmplx=true; /* use complex fft */
    if (!sf_getint("pad1",&pad1)) pad1=1; /* padding factor on the first axis */
    if(!sf_getbool("abc",&abc)) abc=false; /* absorbing flag */
    if (abc) {
      if(!sf_getint("nbt",&nbt)) sf_error("Need nbt!");
      if(!sf_getint("nbb",&nbb)) nbb = nbt;
      if(!sf_getint("nbl",&nbl)) nbl = nbt;
      if(!sf_getint("nbr",&nbr)) nbr = nbt;
      if(!sf_getfloat("ct",&ct)) sf_error("Need ct!");
      if(!sf_getfloat("cb",&cb)) cb = ct;
      if(!sf_getfloat("cl",&cl)) cl = ct;
      if(!sf_getfloat("cr",&cr)) cr = ct;
    } else {
      nbt = 0; nbb = 0; nbl = 0; nbr = 0;
      ct = 0; cb = 0; cl = 0; cr = 0;
    }
    if (!sf_getbool("verb",&verb)) verb=false; /* verbosity */
    if (!sf_getbool("ps",&ps)) ps=false; /* use pseudo-spectral */
    if (ps) sf_warning("Using pseudo-spectral...");
    else sf_warning("Using pseudo-analytical...");
    if (!sf_getbool("mig",&mig)) mig=false; /* use pseudo-spectral */
    if (mig) sf_warning("Time-reversal propagation");
    else sf_warning("Forward modeling");
    if (!sf_getfloat("vref",&vref)) vref=1500; /* reference velocity (default using water) */

    /* setup I/O files */
    Fi = sf_input ("in");
    Fo = sf_output("out");
    if (mig) {
      gpl = -1;
      gpl_v = -1;
      if (NULL==sf_getstring("dat") && NULL==sf_getstring("dat_v"))
	sf_error("Need Data!");
      if (NULL!=sf_getstring("dat")) {
	Fd = sf_input("dat");
	sf_histint(Fd,"n1",&nt);
	sf_histfloat(Fd,"d1",&dt);
	sf_histint(Fd,"n2",&gpl);
      } else Fd = NULL;
      if (NULL!=sf_getstring("dat_v")) {
	Fd_v = sf_input("dat_v");
	sf_histint(Fd_v,"n1",&nt);
	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 */
//.........这里部分代码省略.........

int main(int argc, char* argv[]) 
{
    int nx, nt, ix, it, isx;
    float dt, dx;
    float *old, *nxt, *cur, *sig, *v;
    sf_file in, out, vel;
    int im,im2,im3,im4,im5,ip,ip2,ip3,ip4,ip5;

    sf_init(argc,argv);
    in  = sf_input("in");
    vel = sf_input("vel");   /* velocity */
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(vel,"n1",&nx)) sf_error("No n1= in input");
    if (!sf_histfloat(vel,"d1",&dx)) sf_error("No d1= in input");
    if (!sf_getint("isx",&isx)) isx=(int)(nx/2);
    if (!sf_getint("nt",&nt)) sf_error("No nt in input");
    if (!sf_getfloat("dt",&dt)) sf_error("No dt in input");

    sf_putint(out,"n1",nx);
    sf_putfloat(out,"d1",dx);
    sf_putint(out,"n2",nt);
    sf_putfloat(out,"d2",dt);
    sf_putfloat(out,"o2",0.0); 

    sig = sf_floatalloc(nx);
    old = sf_floatalloc(nx);
    nxt = sf_floatalloc(nx);
    cur = sf_floatalloc(nx);
    v   = sf_floatalloc(nx);
    
 
    sf_floatread(v,nx,vel);
    sf_floatread(sig,nx,in);		

	/* initial conditions */
    for (ix=0; ix < nx; ix++){
        cur[ix] =  sig[ix];
        old[ix] =  0.0; 
	nxt[ix] = 0.;
    }

    /* propagation in time */
    for (it=0; it < nt; it++) {
	sf_floatwrite(cur,nx,out);
	/* Stencil */
	for (ix=0; ix < nx; ix++) {
            im  =(ix-1+nx)%nx; 
            im2 =(ix-2+nx)%nx; 
            im3 =(ix-3+nx)%nx; 
            im4 =(ix-4+nx)%nx; 
            im5 =(ix-5+nx)%nx; 
            ip = (ix+1+nx)%nx;
            ip2 =(ix+2+nx)%nx;
            ip3 =(ix+3+nx)%nx;
            ip4 =(ix+4+nx)%nx;
            ip5 =(ix+5+nx)%nx;
	    nxt[ix] = dt*dt/(12.0*dx*dx)*(-30.0*cur[ix] +16.0* (cur[im]+cur[ip]) - (cur[im2]+cur[ip2]))*v[ix]*v[ix] 
                      +2.0*cur[ix]- old[ix];
	}
       for (ix=0; ix < nx; ix++) {
	    old[ix] = cur[ix];
	    cur[ix] = nxt[ix];
	}
    }
    exit(0);
}

int main(int argc, char* argv[]) 
{
    int nx, nt, ix, it, what;
    float dt, dx, w, g, a0, a1, a2, a3, a4, a5;
    float *old, *nxt, *cur, *sig, *v, *vx, **a; 
    sf_file inp, out, vel, grad;

    sf_init(argc,argv);
    inp = sf_input("in");
    out = sf_output("out");
    vel = sf_input("vel");   /* velocity */
    grad = sf_input("grad"); /* velocity gradient */

    if (SF_FLOAT != sf_gettype(inp)) sf_error("Need float input");
    if (!sf_histint(vel,"n1",&nx)) sf_error("No n1= in input");
    if (!sf_histfloat(vel,"d1",&dx)) sf_error("No d1= in input");
    if (!sf_histint(inp,"n1",&nt)) sf_error("No n2= in input");
    if (!sf_histfloat(inp,"d1",&dt)) sf_error("No d2= in input");

    sf_putint(out,"n1",nx);
    sf_putfloat(out,"d1",dx);
    sf_putint(out,"n2",nt);
    sf_putfloat(out,"d2",dt);

    if (!sf_getint("what",&what)) what=2; /* 2nd or 4th order for FD*/
    sig = sf_floatalloc(nt);
    old = sf_floatalloc(nx);
    nxt = sf_floatalloc(nx);
    cur = sf_floatalloc(nx);
    v = sf_floatalloc(nx);
    vx = sf_floatalloc(nx);

    sf_floatread(v,nx,vel);
    sf_floatread(vx,nx,grad);

    sf_fileclose(vel);
    sf_fileclose(grad);

    sf_floatread(sig,nt,inp);		
    
    switch(what) {
	case 2:
	    /* 2nd order FD*/
	    a = sf_floatalloc2(3,nx);
	    for (ix=0; ix < nx; ix++) {
		/* dimensionless velocity */
		w = v[ix] * dt/dx;
		/* dimensionless gradient */
		g = 0.5 * vx[ix] * dt;

		a1 = w*w * (1.0 + g*g);
		a2= g*w;
	
		a[ix][0] = a1+a2;
		a[ix][1] = -2*a1;
		a[ix][2] = a1-a2;

		/* initial conditions */
		cur[ix] = 0.;
		nxt[ix] = 0.;
	    }

	    free(v);
	    free(vx);

	    /* propagation in time */
	    for (it=0; it < nt; it++) {


		for (ix=0; ix < nx; ix++) {
		    old[ix] = cur[ix];
		    cur[ix] = nxt[ix];
		}   

		/* Stencil */
		nxt[0] = sig[it]; 
		for (ix=1; ix < nx-1; ix++) {
		    nxt[ix] = cur[ix]*a[ix][1] + cur[ix+1]*a[ix][0] + cur[ix-1]*a[ix][2];
		}
		nxt[nx-1] = cur[nx-1]*a[nx-1][1] + cur[nx-1]*a[nx-1][0] + cur[nx-2]*a[nx-1][2];
	
		for (ix=1; ix < nx; ix++) {
		    nxt[ix] += 2*cur[ix] - old[ix];
		} 

		sf_floatwrite(nxt,nx,out);
	    }
	    break;
        
        case 4:
	    /* 4th order FD*/
	    a = sf_floatalloc2(5,nx);
     
	    for (ix=0; ix < nx; ix++) {
		/* dimensionless velocity */
		w = v[ix] * dt/dx;
		/* dimensionless gradient */
		g = vx[ix] * dt;

		a1 = w*w * (4.0 + g*g);
//.........这里部分代码省略.........

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_setlabel(amx,"mx");
    amy = sf_iaxa(Bws,2); sf_setlabel(amy,"my");
    aw  = sf_iaxa(Bws,4); sf_setlabel(aw ,"w" );

    Pws = sf_tmpfile(NULL); sf_settype(Pws,SF_COMPLEX);
    Pwr = sf_tmpfile(NULL); sf_settype(Pwr,SF_COMPLEX);

    /*------------------------------------------------------------*/
    cub = wex_cube(verb,
                   amx,amy,az,
                   alx,aly,
                   aw,
                   ae,
                   eps,
                   ompnth);
//.........这里部分代码省略.........

int main (int argc, char* argv[]) 
{
    int n1,n2,n3; /*n1 is trace length, n2 is the number of traces, n3 is the number of 3th axis*/
    int nfw;    /*nfw is the filter-window length*/
    int m;
    int i,j,k,ii;
    bool boundary;
    
    float *trace;
    float *tempt;
    float *temp1;
    float *extendt;
    sf_file in, out;
    
    sf_init (argc, argv); 
    in = sf_input("in");
    out = sf_output("out");
    
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);
    /* get the trace length (n1) and the number of traces (n2) and n3*/
    
    if (!sf_getint("nfw",&nfw)) sf_error("Need integer input");
    /* filter-window length (positive and odd integer)*/
    
    if (!sf_getbool("boundary",&boundary)) boundary=false;
    /* if y, boundary is data, whereas zero*/
    
    if (nfw < 1)  sf_error("Need positive integer input"); 
    if (nfw%2 == 0) nfw++;
    m=(nfw-1)/2;
    
    trace = sf_floatalloc(n1*n2);
    tempt = sf_floatalloc(n1*n2);
    temp1 = sf_floatalloc(nfw);
    extendt = sf_floatalloc((n1+2*m)*n2);
    
    for(ii=0;ii<n3;ii++){
	sf_floatread(trace,n1*n2,in);
	
	for(i=0;i<n1*n2;i++){
	    tempt[i]=trace[i];
	}
	
	bound1(tempt,extendt,nfw,n1,n2,boundary);
	
	/************1D median filter****************/
	
	for(i=0;i<n2;i++){
	    for(j=0;j<n1;j++){
		for(k=0;k<nfw;k++){
		    temp1[k]=extendt[(n1+2*m)*i+j+k];
		}
		trace[n1*i+j]=sf_quantile(m,nfw,temp1); 
	    }
	}
	
	sf_floatwrite(trace,n1*n2,out);
    }

    exit (0);
}

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

    float *data /*input [nk] */; 
    int nk      /* data length */;
    float dk    /* data sampling */; 
    float k0    /* initial location */;
    bool verb   /* verbosity flag */;
    int niter   /* number of iterations */; 

    int ik, iter;
    float x1, x2, x3, dx1, dx2, dx3, x10, x20, x30, x4, x5;
    float k, f, s, l, r, p, q, s2, y1, y2, y3, sk, s1, s3;
    float a11, a12, a13, a22, a23, a33, a2, a3, a4, a5, a6, d;
    float b11, b12, b13, b22, b23, b33, eps, r2; 
    
    sf_file in, out, prm;

    /* Estimate shape (Caution: data gets corrupted) */ 

    sf_init (argc,argv);
    in = sf_input("in");
    out = sf_output("out");
    prm = sf_output("prm");


    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"n1",&nk)) sf_error("No n1= in input");
    if (!sf_histfloat(in,"d1",&dk)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&k0)) sf_error("No o1= in input");
 
    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */
    if (!sf_getfloat("x10",&x10)) x10=6.;
    /* initial nonlinear parameter x1 value */
    if (!sf_getfloat("x20",&x20)) x20=-0.5;
    /* initial nonlinear parameter x2 value */
    if (!sf_getfloat("x30",&x30)) x30=200.;
    /* initial nonlinear parameter x3 value */
    if (!sf_getbool("verb",&verb)) verb=false;
    /* verbosity flag */

    sf_putint(prm,"n1",3);
    sf_putint(prm,"nk",nk);
    sf_putfloat(prm,"dk",dk);
    sf_putfloat(prm,"k0",k0);
    sf_fileflush(prm,in);

    data = sf_floatalloc(nk);
    
    sf_floatread(data,nk,in);

    x1 = x10; /* initial x1 */
    x2 = x20; /* initial x2 */
    x3 = x30; /* initial x3 */
    x4 = x3*x3;

    eps = 10.*FLT_EPSILON;
    eps *= eps;
    a11 = (float)nk;

    if (verb) sf_warning("got x10=%g x20=%g x30=%g k0=%g niter=%d nk=%d dk=%g",
			 x10,x20,x30,k0,niter,nk,dk);
	
    /* Gauss iterations */
    for (iter = 0; iter < niter; iter++) {
        x5 = 2.*x2*x3;
        s2 = y1 = y2 = y3 = sk = s1 = s3 = a12 = r2 = 0.;
	for (ik = 0; ik < nk; ik++) {
	    k = ik*dk + k0;
	    k *= k;
            s = 1. + x4*k;
	    l = log(s);
            q = 1./s;
	    p = k*q;
            
            s2 += l*l;
            sk += p;
            a12 += l;

	    f = log(data[ik]);       
	    r = x1 + x2*l - f;

	    y1 += r;
	    y2 += r*l;
	    y3 += r*p;
	    s1 += p*(r + x2*l);
            s3 += p*q*(2.*x2*x4*k + r*(1.-x4*k));

	    r2 += r*r;
	}

        y3 *= x5;

	a13 = x5*sk;
	a22 = s2;
        a23 = 2.*x3*s1;
	a33 = 2.*x2*s3;

	a2 = a12*a12;
//.........这里部分代码省略.........

int main (int argc, char* argv[]) {
    int nz, nx, ny, nb, na, iz, ix, iy, ia, ib, fz, lz;
    int icpu = 0, ncpu = 1, morder = 2, ic, nc = 1, mp = 1, ith = 0, inet = 0, tdel = 0;
    float dz, oz, dx, ox, dy, oy, db, ob, da, oa, aper;
    float z, x, y;
    float ****e;
    sf_file spdom, vspline = NULL, scgrid = NULL, scdaemon = NULL, 
            out;
    char *ext = NULL;
    bool verb, parab, mmaped, rfail;
    sf_esc_slowness3 esc_slow;
    sf_esc_scglstor3 esc_scgrid_lstor;
    sf_esc_tracer3 *esc_tracers;
    sf_esc_scgrid3 *esc_scgrids;
    sf_timer timer;

    sf_init (argc, argv);

    if (!sf_stdin ()) {
        spdom = NULL;
    } else {
        spdom = sf_input ("in");
        /* Spatial (z,x,y) domain */
    }

    out = sf_output ("out");
    /* Escape values */

    /* Spatial dimensions */
    if (spdom) {
        if (!sf_histint (spdom, "n1", &nz)) sf_error ("No n1= in input");
        if (!sf_histint (spdom, "n2", &nx)) sf_error ("No n2= in input");
        if (!sf_histint (spdom, "n3", &ny)) sf_error ("No n3= in input");
        if (!sf_histfloat (spdom, "d1", &dz)) sf_error ("No d1= in input");
        if (!sf_histfloat (spdom, "o1", &oz)) sf_error ("No o1= in input");
        if (!sf_histfloat (spdom, "d2", &dx)) sf_error ("No d2= in input");
        if (!sf_histfloat (spdom, "o2", &ox)) sf_error ("No o2= in input");
        if (!sf_histfloat (spdom, "d3", &dy)) sf_error ("No d3= in input");
        if (!sf_histfloat (spdom, "o3", &oy)) sf_error ("No o3= in input");
        if (!sf_histint (spdom, "icpu", &icpu)) icpu = 0;
        /* Current CPU number */
        if (!sf_histint (spdom, "ncpu", &ncpu)) ncpu = 1;
        /* Total number of CPUs */
    }
    ext = sf_escst3_warnext (spdom);
    if (!sf_getint ("nz", &nz) && !spdom) sf_error ("Need nz=");
    /* Number of samples in z axis */
    if (!sf_getfloat ("oz", &oz) && !spdom) sf_error ("Need oz=");
    /* Beginning of z axis */
    if (!sf_getfloat ("dz", &dz) && !spdom) sf_error ("Need dz=");
    /* Sampling of z axis */
    if (!sf_getint ("nx", &nx) && !spdom) sf_error ("Need nx=");
    /* Number of samples in x axis */
    if (!sf_getfloat ("ox", &ox) && !spdom) sf_error ("Need ox=");
    /* Beginning of x axis */
    if (!sf_getfloat ("dx", &dx) && !spdom) sf_error ("Need dx=");
    /* Sampling of x axis */
    if (!sf_getint ("ny", &ny) && !spdom) sf_error ("Need ny=");
    /* Number of samples in y axis */
    if (!sf_getfloat ("oy", &oy) && !spdom) sf_error ("Need oy=");
    /* Beginning of y axis */
    if (!sf_getfloat ("dy", &dy) && !spdom) sf_error ("Need dy=");
    /* Sampling of y axis */

    if (!sf_getint ("na", &na)) na = 360;
    /* Number of azimuth phase angles */
    da = 2.0*SF_PI/(float)na;
    oa = 0.5*da;

    if (!sf_getint ("nb", &nb)) nb = 180;
    /* Number of inclination phase angles */
    db = SF_PI/(float)nb;
    ob = 0.5*db;

#ifdef _OPENMP
    if (!sf_getint ("mp", &mp)) mp = 1;
    /* Bufferization factor for multicore processing (number of points in buffer = mp*nc) */
    if (!sf_getint ("nc", &nc)) nc = 0;
    /* Number of threads to use for ray tracing (OMP_NUM_THREADS by default) */
    if (nc)
        omp_set_num_threads (nc); /* User override */
    else
        nc = omp_get_max_threads (); /* Current default */
    sf_warning ("%s Using %d threads", ext, omp_get_max_threads ());
    sf_warning ("%s Buffering %d points", ext, nc*mp);
#endif

    if (!sf_getfloat ("aper", &aper)) aper = SF_HUGE;
    /* Maximum aperture in x and y directions from current point (default - up to grid boundaries) */
    if (aper != SF_HUGE)
        aper = fabsf (aper);

    if (!sf_getbool ("parab", &parab)) parab = true;
    /* y - use parabolic approximation of trajectories, n - straight line */

    if (!sf_getbool ("mmaped", &mmaped)) mmaped = true;
    /* n - do not use memory mapping for local data access */

    if (!sf_getbool ("rfail", &rfail)) rfail = true;
    /* n - do not quit if remote processing fails, try local processing */
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    int i1, n1, i2, n2, i3, n3, ip, np, it;
    float **dat=NULL, *trace=NULL, dp, p0, p, t;
    vint1 str;
    sf_file inp=NULL, out=NULL;

    sf_init(argc,argv);
    inp = sf_input("in");
    out = sf_output("out");

    if (!sf_histint(inp,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(inp,"n2",&n2) || n2 < 2) sf_error("No n2= in input");
    n3 = sf_leftsize(inp,2);

    if (!sf_getint("np",&np)) sf_error("Need np=");
    /* number of slopes */
    if (!sf_getfloat("dp",&dp)) sf_error("Need dp=");
    /* slope sampling */
    if (!sf_getfloat("p0",&p0)) sf_error("Need p0=");
    /* first slope */

    sf_putint(out,"n3",np);
    sf_putfloat(out,"d3",dp);
    sf_putfloat(out,"o3",p0);

    sf_shiftdim(inp,out,3);

    dat = sf_floatalloc2(n1,n2);
    trace = sf_floatalloc(n2);
    trace[n2-1] = 0.;

    str = vint1_init (4,n1,n2-1);

    for (i3=0; i3 < n3; i3++) {
	sf_floatread(dat[0],n1*n2,inp);
	vint1_set(str,dat+1);

	for (ip=0; ip < np; ip++) {
	    p = p0 + ip*dp;

	    for (i1=0; i1 < n1; i1++) {
		t = i1+p;
		it = floorf(t);

		if (it < -1 || it > n1) {
		    for (i2=0; i2 < n2-1; i2++) {
			trace[i2]=0.;
		    }
		} else {
		    vint1_apply(str,it,t-it,false,trace);
		}

		for (i2=0; i2 < n2; i2++) {
		    dat[i2][i1] = trace[i2];
		}
	    }
	    sf_floatwrite(dat[0],n1*n2,out);
	}
    }

    exit(0);
}