void sf_close(void)
/*< Remove temporary files >*/
{
    int i;
    sf_file file;
	
    if (NULL == infiles) return;
    
    for (i=0; i <= (int) ifile; i++) {
	file = infiles[i];
		
	if (NULL != file && 
	    NULL != file->dataname && 
	    file->pipe &&
	    0 != strcmp("stdin",file->dataname) && 
	    0 != unlink(file->dataname))
	    sf_warning ("%s: trouble removing %s:",__FILE__,file->dataname);
	    
	sf_fileclose(file);
	free(file);
    }
    free(infiles);
    infiles=NULL;
    ifile=nfile=0;
}

int main(int argc, char* argv[])
{
    int nd, m1, na, nr, niter;
    float *rr, *dd, *coord, o1, d1, eps;
    char *header;
    sf_file in, out, head;

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

    /* create model */
    if (!sf_getint ("nx",&m1)) sf_error("Need n1=");
    /* number of bins */
    sf_putint(out,"n1",m1);

    if (!sf_getfloat("x0",&o1)) sf_error("Need o1=");
    /* grid origin */
    sf_putfloat (out,"o1",o1);

    if (!sf_getfloat("dx",&d1)) sf_error("Need d1=");
    /* grid sampling */
    sf_putfloat (out,"d1",d1);

    if (!sf_getint("niter",&niter)) niter=1+m1*3/2; niter *= 2;
    /* number of conjugate-gradient iterations */
    if (!sf_getfloat("eps",&eps)) eps=0.2;
    /* regularization parameter */

    /* create filter */
    if (!sf_getint("na",&na)) na=3;
    nr = m1 + na;

    rr = sf_floatalloc(nr);

    if (!sf_histint(in,"n1",&nd)) nd=1;
    coord = sf_floatalloc(nd);
    dd = sf_floatalloc(nd);

    header = sf_getstring("head");
    if (NULL == header) { 
	header = sf_histstring(in,"head");
	if (NULL == header) sf_error("Need head=");
    }

    head = sf_input(header);
    if (SF_FLOAT != sf_gettype(head)) sf_error("Need float head");
    sf_floatread (coord,nd,head);
    sf_fileclose (head);

    sf_floatread (dd,nd,in);

    levint1 (niter, m1, nr, nd, coord, dd, o1, d1, rr, eps);
    
    sf_floatwrite (rr,m1,out);


    exit(0);
}

int main(int argc, char* argv[])
{
    int rank, nodes, ndim, job, axis, axis2, jobs;
    off_t n[SF_MAX_DIM];
    char **commands, cmdline[SF_CMDLEN], *iname;
    sf_file inp=NULL, out=NULL;
    MPI_Status stat;

    MPI_Init(&argc,&argv);

    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nodes);
    if (nodes < 2) {
	fprintf(stderr,"Need at least two nodes!\n");
	MPI_Finalize();
	exit(1);
    }

    if (!rank) { /* master node */
	sf_init(argc,argv);

	inp = sf_input("--input");
	out = sf_output("--output");

	ndim = sf_largefiledims (inp,n);

	if (!sf_getint("split",&axis)) axis=ndim;
	/* axis to split */
	
	commands = sf_split(inp,axis,nodes,&jobs,ndim,n,argc,argv);  

	for (job=0; job < jobs; job++) {
	    strncpy(cmdline,commands[job],SF_CMDLEN);
	    MPI_Send(cmdline, SF_CMDLEN, MPI_CHAR, job+1, 0, MPI_COMM_WORLD);
	}

	iname = sf_getstring("--input");

	if (!sf_getint("join",&axis2)) axis2=axis;
	/* axis to join (0 means add) */

	sf_out(out,axis2,iname);
	
	for (job=0; job < jobs; job++) {
	    MPI_Recv(&rank,1, MPI_INT, job+1, 1, MPI_COMM_WORLD,&stat);
	    if (axis2 > 0) sf_join(out,job);
	}
	if (0==axis2) sf_add(out,jobs);

	sf_fileclose(inp);
    } else { /* slave nodes */
	MPI_Recv(cmdline, SF_CMDLEN, MPI_CHAR, 0, 0, MPI_COMM_WORLD,&stat);
	fprintf(stderr,"CPU %d: %s\n",rank,cmdline); 
        sf_system(cmdline);
	MPI_Send(&rank,1,MPI_INT,0,1,MPI_COMM_WORLD);
    }
    
    MPI_Finalize();
}

int main(int argc, char* argv[])
{
    int n1, n2, i2, esize;
    off_t pos;
    struct skey *sorted;
    float *unsorted;
    char *trace, *header;
    sf_file in, head, out;

    sf_init (argc,argv);
    in = sf_input ("in");
    out = sf_output ("out");
 
    header = sf_getstring("head");
    /* header file */
    if (NULL == header) { 
	header = sf_histstring(in,"head");
	if (NULL == header) sf_error("Need head=");
    }

    head = sf_input(header);
    if (SF_FLOAT != sf_gettype(head))
	sf_error("Need float header");
    n2 = sf_filesize(head);
 
    unsorted = sf_floatalloc(n2);
    sorted = (struct skey*) sf_alloc(n2,sizeof(struct skey));
    
    sf_floatread(unsorted,n2,head);
    for (i2 = 0; i2 < n2; i2++) {
	sorted[i2].key = unsorted[i2];
	sorted[i2].pos = i2;
    }
    free (unsorted);
    sf_fileclose(head);

    qsort(sorted,n2,sizeof(struct skey),key_compare);
 
    if (!sf_histint(in,"n1",&n1)) n1=1;
    esize = sf_esize(in);
    n1 *= esize;

    trace = sf_charalloc(n1);

    sf_unpipe(in,((off_t) n1)*((off_t) n2));
    sf_fileflush(out,in);
    sf_setform(in,SF_NATIVE);
    sf_setform(out,SF_NATIVE);

    pos = sf_tell(in);
    for (i2=0; i2<n2; i2++) {
	sf_seek(in,pos+(sorted[i2].pos)*n1,SEEK_SET);
	sf_charread(trace,n1,in);
	sf_charwrite(trace,n1,out);
    }


    exit(0);
}

int main(int argc, char* argv[])
{
    int nd, n1, niter;
    float *d, *a, *b, alfa, beta, perc;
    sf_file inp, reg, out;

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

    if (!sf_histint(inp,"n1",&nd)) sf_error("No n1= in input");
    if (!sf_histint(reg,"n1",&n1) || n1 != nd)
	sf_error("Need n1=%d in reg",nd);
    if (!sf_histint(reg,"n2",&n1) || n1 != 2)
	sf_error("Need n2=2 in reg");

    sf_putint(out,"n1",2);
    
    d = sf_floatalloc(nd);
    a = sf_floatalloc(nd);
    b = sf_floatalloc(nd);

    sf_floatread(d,nd,inp);
    sf_floatread(a,nd,reg);
    sf_floatread(b,nd,reg);
    sf_fileclose(reg);

    if (!sf_getint("niter",&niter)) niter=10;
    /* number of POCS iterations */
    if (!sf_getfloat("perc",&perc)) perc=90.0;
    /* percentage for sharpening */

    l1_init(nd,niter,perc,true);

    bil1(d,a,b,&alfa,&beta);
    
    sf_floatwrite(&alfa,1,out);
    sf_floatwrite(&beta,1,out);

    exit(0);
}

void sf_join(sf_file out /* output file */, 
	     int job     /* job number */)
/*< join outputs >*/
{
    int i, chunk;
    off_t i2, left, nbuf;
    char *iname, *oname;
    sf_file in;

    iname = inames[job];
    oname = onames[job];

    in = sf_input(oname);
    sf_histint(in,nkey,&chunk);
    sf_setform(in,SF_NATIVE);
	    
    for (i2=0; i2 < size2; i2++) {
	nbuf=BUFSIZ;
	for (left=chunk*size1; left > 0; left -= nbuf) {
	    if (nbuf > left) nbuf=left;
		    
	    sf_charread(buffer,nbuf,in);
	    sf_charwrite(buffer,nbuf,out);
	}
    }

    sf_fileclose(in);

    sf_rm(iname,true,false,false);
    sf_rm(oname,true,false,false);

    for (i=0; i < inpargc; i++) {
	iname = inpnames[job][i];
	sf_rm(iname,true,false,false);
    }
}

int main(int argc, char* argv[])
{
    int n12, n1, n2, n3, i3, sw, niter, liter;
    bool hd, ps, verb;
    float *data, *modl, *vrms, *error=NULL;
    float o1,d1,o2,d2,eps;
    char *errfile;
    sf_file in, out, vel, err=NULL;

    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);

    if (!sf_getbool("hd",&hd)) hd=true;
    /* if y, apply half-derivative filter */
    if (!sf_getbool("ps",&ps)) ps=true;
    /* if y, apply pseudo-unitary weighting */
    if (!sf_getbool("verb",&verb)) verb=false;
    /* verbosity flag */
    if (!sf_getint("sw",&sw)) sw=0;
    /* if > 0, select a branch of the antialiasing operation */
    if (!sf_getint("niter",&niter)) niter=5;
    /* number of non-linear iterations, when niter=1, it's linear */
    if (!sf_getint("liter",&liter)) liter=5;
    /* number of linear iterations */
    if (!sf_getfloat("eps",&eps)) eps=0.;
    /* regularization parameters */

    if (!sf_histfloat(in,"o1",&o1)) sf_error("No o1= in input");
    if (!sf_histfloat(in,"o2",&o2)) sf_error("No o2= in input");
    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"d2",&d2)) sf_error("No d2= in input");

    vrms = sf_floatalloc(n1);

    vel = sf_input("velocity");
    sf_floatread(vrms,n1,vel);
    sf_fileclose(vel);
    
    n12 = n1*n2;
    data = sf_floatalloc(n12);
    modl = sf_floatalloc(n12);

    if (niter > 0) {
	errfile = sf_getstring("err");
	/* output file for error */
	if (NULL != errfile) {
	    err = sf_output(errfile);
	    sf_putint(err,"n1",liter);
	    sf_putfloat(err,"d1",1);
	    sf_putfloat(err,"o1",1);
	    sf_putstring(err,"label1","Iteration Number");
	    sf_putstring(err,"label2","Relative Squared Error");
	    sf_putint(err,"n2",1);
	}
	error = sf_floatalloc(liter);
    }

    kirchnew_init (vrms, o1, d1, d2, n1, n2, sw, ps, hd);

    for (i3=0; i3 < n3; i3++) {

	sf_floatread (data,n12,in);

	inverts(kirchnew_lop,niter,niter,liter,n12,n12,modl,data,error,verb,eps);

	sf_floatwrite (modl,n12,out);

	if (NULL != err) sf_floatwrite(error,liter,err);
    }

    sf_warning(".");

    exit(0);
}

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

    int n1, n2, n3;
    int n2top, n2bot, n3top, n3bot;
    float d1, o1;
    float badSample;
    int size;
    float *panel, *trace;
    float top, bottom;
    int pind, i3, i2, ind;
    float t, stack;
    sf_file dataFile, hTopFile, hBotFile, stackFile;

// Initialize RSF
    sf_init (argc,argv);
// Input files
    dataFile = sf_input ("in");
// check that the input is float
    if ( SF_FLOAT != sf_gettype (dataFile) ) sf_error ("Need float input: data volume");
    /* data volume */

    if ( NULL != sf_getstring ("top") ) {
        /* top horizon */
        hTopFile = sf_input ("top");
    } else {
        sf_error ("Need float input: top horizon");
    }

    if ( NULL != sf_getstring ("bot") ) {
        /* bottom horizon */
        hBotFile = sf_input ("bot");
    } else {
        sf_error ("Need float input: bottom horizon");
    }
// Output file
    stackFile = sf_output ("out");

    // data parameters
    if ( !sf_histint   (dataFile, "n1", &n1) )  sf_error ("Need n1= in input");
    if ( !sf_histfloat (dataFile, "d1", &d1) )  sf_error ("Need d1= in input");
    if ( !sf_histfloat (dataFile, "o1", &o1) )  sf_error ("Need o1= in input");
    if ( !sf_histint   (dataFile, "n2", &n2) )  sf_error ("Need n2= in input");
    if ( !sf_histint   (dataFile, "n3", &n3) )  sf_error ("Need n3= in input");

    if ( !sf_histint   (hTopFile, "n2", &n2top) )  sf_error ("Need n2= in top horizon file");
    if ( !sf_histint   (hTopFile, "n3", &n3top) )  sf_error ("Need n3= in top horizon file");

    if ( !sf_histint   (hTopFile, "n2", &n2bot) )  sf_error ("Need n2= in bottom horizon file");
    if ( !sf_histint   (hTopFile, "n3", &n3bot) )  sf_error ("Need n3= in bottom horizon file");

    // files-consistency checking
    if (n2bot != n2 || n3bot != n3 || n2top != n2 || n3top != n3)
        sf_error ("Horizons are not consistent with data");

    sf_putint (stackFile, "n1", 1);

    // run parameters
    if (!sf_getfloat ("badSample", &badSample)) badSample = 99999.f;
    /* non-interpreted sample in the horizons */

    size = n2 * n3;
    panel = sf_floatalloc (size);
    trace = sf_floatalloc (n1);

    pind = 0;
    for (i3 = 0; i3 < n3; ++i3) {
        for (i2 = 0; i2 < n2; ++i2, ++pind) {
            panel [pind] = 0.f;

            sf_floatread (&top, 1, hTopFile);
            sf_floatread (&bottom, 1, hBotFile);
            sf_floatread (trace, n1, dataFile);

            if (fabs (top - badSample) < 1.f) continue;    // non-interpreted sample in the horizons
            if (fabs (bottom - badSample) < 1.f) continue;

            t = o1;
            ind = 0;
            while (t < top && ind < n1) {
                t += d1;
                ++ind;
            }

            stack = 0.f;
            while (t < bottom && ind < n1) {
                stack += pow (trace [ind], 2);
                ++ind;
                t += d1;
            }

            panel [pind] = stack;
        }
    }

    sf_floatwrite (panel, size, stackFile);

    sf_fileclose (dataFile);
    sf_fileclose (hTopFile);
    sf_fileclose (hBotFile);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
#ifdef _OPENMP
#pragma omp parallel for private(is,ix,it)
#endif
        for (is=shtnum0; is<shtnum; is++)
            for (ix=0; ix<gpl; ix++)
                for (it=0; it<nt; it++)
                    record[is][ix][it] = sf_cmplx(0.,0.);
    }
    /*read starting model*/
    if (NULL != Fstart) {
        sf_complexread(imginv[0],nx*nz,Fstart);
    } else {
        /*transform the dimension to 1d*/
#ifdef _OPENMP
#pragma omp parallel for private(iz)
#endif
        for (iz=0; iz<nz; iz++) {
            for (ix=0; ix<nx; ix++) {
                imginv[ix][iz] = sf_cmplx(0.,0.);
            }
        }
    }

    /* output RSF files */
    if (cpuid==0) {
        sf_setn(az, nz);
        sf_setn(ax, nx);
        sf_oaxa(Fimg, az, 1);
        sf_oaxa(Fimg, ax, 2);
        sf_settype(Fimg,SF_COMPLEX);
    }

    /*close RSF files*/
    sf_fileclose(Fvel);
    sf_fileclose(Fsrc);
    sf_fileclose(left);
    sf_fileclose(right);
    sf_fileclose(leftb);
    sf_fileclose(rightb);

    /*load geopar elements*/
    geop->nx  = nx;
    geop->nz  = nz;
    geop->nxb = nxb;
    geop->nzb = nzb;
    geop->dx  = dx;
    geop->dz  = dz;
    geop->ox  = ox;
    geop->oz  = oz;
    geop->snpint = snpint;
    geop->spz = spz;
    geop->gpz = gpz;
    geop->gpl = gpl;
    geop->top = top;
    geop->bot = bot;
    geop->lft = lft;
    geop->rht = rht;
    geop->nt = nt;
    geop->dt = dt;
    geop->trunc = trunc;
    /*geop->adj   = adj; */
    geop->verb  = verb;
    geop->illum = illum;
    geop->m2 = m2;
    geop->m2b = m2b;
    geop->pad1 = pad1;

int main (int argc, char *argv[])
{
    bool  verb;         /* verbosity */
    int ompnth=1;       /* number of threads */
    bool  down;         /* adjoint flag */
    bool causal;

    sf_file Fslo=NULL;  /*  slowness file  S(nlx,nly,nz   ) */
    sf_file Fsou=NULL;  /*      data file Ds(nmx,nmy,   nw) */
    sf_file Fwfl=NULL;  /*      data file Dr(nmx,nmy,   nw) */

    weicub3d cub;       /* hypercube */
    weitap3d tap;       /* tapering */
    weissr3d ssr;       /* SSR operator */
    weislo3d slo;       /* slowness */
    weiop3d weop;       /* WEI operator */

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

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

    if (!sf_getbool(  "verb",&verb  ))  verb = false; /* verbosity flag */
    if (!sf_getbool(  "down",&down  ))  down = true;  /* up/down   flag */
    if (!sf_getbool("causal",&causal))  sf_error("Specify causal!"); /* causality flag */

    /*------------------------------------------------------------*/
    if(verb) fprintf(stderr,"init cube...");
    cub = wei_cube(verb,ompnth);
    if(verb) fprintf(stderr,"OK\n");

    /*------------------------------------------------------------*/
    Fslo = sf_input ("slo");
    if (SF_FLOAT !=sf_gettype(Fslo)) sf_error("Need float slowness");
    weislo_inp(cub,Fslo); /* input slowness */

    /*------------------------------------------------------------*/
    Fsou = sf_input ( "in");
    if (SF_COMPLEX !=sf_gettype(Fsou)) sf_error("Need complex sdat");
    weiwfl_inp(cub,Fsou); /* input adjoint source dimensions */

    /*------------------------------------------------------------*/
    Fwfl = sf_output("out"); sf_settype(Fwfl,SF_COMPLEX);

    /*------------------------------------------------------------*/
    if(verb) fprintf(stderr,"init slo...");
    slo = weislo_init(cub,Fslo);
    if(verb) fprintf(stderr,"OK\n");
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    if(verb) fprintf(stderr,"init tap...");
    tap = weitap_init(cub);
    if(verb) fprintf(stderr,"OK\n");
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    if(verb) fprintf(stderr,"init ssr...");
    ssr = weissr_init(cub,slo);
    if(verb) fprintf(stderr,"OK\n");
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    if(verb) fprintf(stderr,"init weop...");
    weop = weiwfl_init(cub);
    if(verb) fprintf(stderr,"OK\n");
    /*------------------------------------------------------------*/

    adjwfl(weop,cub,ssr,tap,slo,Fsou,Fwfl,down,causal);
    weiwfl_close(weop);

    /*------------------------------------------------------------*/
    if(verb) fprintf(stderr,"close structures...");
    weislo_close(slo);
    weissr_close(cub,ssr);
    weitap_close(tap);
    if(verb) fprintf(stderr,"OK\n");
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    /* close files */
    if (Fsou!=NULL) sf_fileclose(Fsou);
    if (Fwfl!=NULL) sf_fileclose(Fwfl);
    if (Fslo!=NULL) sf_fileclose(Fslo);
    /*------------------------------------------------------------*/

    exit (0);
}

int main(int argc, char* argv[])
{
    bool velocity, lsint;
    int is, nz, nx, im, nm, order, nshot, ndim, nt, nr, ir, i1, i2, ia;
    float a, b, fb, r0, dr, rmax, xmax, a0, amax, da, r1, r2, tol;
    float dz, dx, z0, x0, *angle, *slow, **shot, *time;
    sf_file shots, vel, out;

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

    /* get 2-D grid parameters */
    if (!sf_histint(vel,"n1",&nz))   sf_error("No n1= in input");
    if (!sf_histint(vel,"n2",&nx))   sf_error("No n2= in input");
    if (!sf_histfloat(vel,"d1",&dz)) sf_error("No d1= in input");
    if (!sf_histfloat(vel,"d2",&dx)) sf_error("No d2= in input");
    if (!sf_histfloat(vel,"o1",&z0)) z0=0.;
    if (!sf_histfloat(vel,"o2",&x0)) x0=0.;
    xmax = x0 + (nx-1)*dx;

    /* additional parameters */
    if(!sf_getbool("vel",&velocity)) velocity=true;
    /* If y, the input is velocity; if n, slowness */
    if(!sf_getint("order",&order)) order=4;
    /* Interpolation order */
    if (!sf_getbool("lsint",&lsint)) lsint=false;
    /* if use least-squares interpolation */

    if (!sf_getint("nt",&nt)) nt=nx*nz;
    /* Maximum number of time steps */

    if (!sf_getint("nr",&nr)) nr=1;
    /* number of recievers */
    if (!sf_getfloat("r0",&r0)) r0=x0;
    /* first receiver */
    if (!sf_getfloat("dr",&dr)) dr=dx;
    /* receiver increment */
    rmax = r0 + (nr-1)*dr;

    /* sanity check */
    if (dr <= 0.) sf_error ("Need positive dr");
    if (r0 < x0 || r0 > xmax || rmax < x0 || rmax > xmax)
	sf_error ("Range of receivers [%f,%f] "
		  "is outside of the model range [%f,%f]",
		  r0,rmax,x0,xmax);
 
    /* get shot locations */
    if (NULL != sf_getstring("shotfile")) {
	/* file with shot locations */
	shots = sf_input("shotfile");
	if (!sf_histint(shots,"n1",&ndim) || 2 != ndim) 
	    sf_error("Must have n1=2 in shotfile");
	if (!sf_histint(shots,"n2",&nshot)) 
	    sf_error("No n2= in shotfile");
  
	shot = sf_floatalloc2 (ndim,nshot);
	sf_floatread(shot[0],ndim*nshot,shots);
	sf_fileclose (shots);
    } else {
	nshot = 1;
	ndim = 2;

	shot = sf_floatalloc2 (ndim,nshot);

	if (!sf_getfloat("zshot",shot[0]))   shot[0][0]=0.;
	/* shot coordinates (if no shotfile) */
	if (!sf_getfloat("yshot",shot[0]+1)) shot[0][1]=x0 + 0.5*(nx-1)*dx;
	
	sf_warning("Shooting from z=%f, x=%f",shot[0][0],shot[0][1]);
    }

    /* specify output dimensions */
    sf_putint (out,"n1",nr);
    sf_putfloat (out,"d1",dr);
    sf_putfloat (out,"o1",r0);
    sf_putint (out,"n2",nshot);
	    
    /* get velocity */
    nm = nz*nx;
    slow = sf_floatalloc(nm);

    sf_floatread(slow,nm,vel);

    if (vel) { /* convert to slowness */
	for(im = 0; im < nm; im++){
	    slow[im] = 1./slow[im];
	}
    }

    /* initialize ray tracing object */
    ct = sf_celltrace_init (lsint, order, nt, nz, nx, dz, dx, z0, x0, slow);
    
    free (slow);

    angle = sf_floatalloc (nr);
    time  = sf_floatalloc (nr);

    if (!sf_getfloat("tol",&tol)) tol=0.01;
    /* Shooting tolerance (in degrees) */
//.........这里部分代码省略.........

int main (int argc, char* argv[]) {
    int nz, nx, ny, ic, nc;
    float dz, oz, dx, ox, dy, oy;
    size_t i, n, sz;
    float *buf = NULL, *buf2 = NULL;
#ifdef HAVE_SSE
    unsigned char pad[64];
#endif
    sf_file velz, velx = NULL, theta = NULL, phi = NULL, eta = NULL, out;
    bool verb;
    Ugrid z_grid, x_grid, y_grid;
    BCtype_s zBC, xBC, yBC;
    multi_UBspline_3d_s *velspline = NULL;

    sf_init (argc, argv);

    velz = sf_input ("in");
    /* Vertical velocity */
    out = sf_output ("out");
    /* Spline coefficients */

    /* Spatial dimensions */
    if (!sf_histint (velz, "n1", &nz)) sf_error ("No n1= in input");
    if (!sf_histint (velz, "n2", &nx)) sf_error ("No n2= in input");
    if (!sf_histint (velz, "n3", &ny)) sf_error ("No n3= in input");
    if (!sf_histfloat (velz, "d1", &dz)) sf_error ("No d1= in input");
    if (!sf_histfloat (velz, "o1", &oz)) oz = 0.0;
    if (!sf_histfloat (velz, "d2", &dx)) sf_error ("No d2= in input");
    if (!sf_histfloat (velz, "o2", &ox)) ox = 0.0;
    if (!sf_histfloat (velz, "d3", &dy)) sf_error ("No d3= in input");
    if (!sf_histfloat (velz, "o3", &oy)) oy = 0.0;

    if (!sf_getbool ("verb", &verb)) verb = false;
    /* verbosity flag */

    n = (size_t)nz*(size_t)nx*(size_t)ny;
    buf = sf_floatalloc (n);
    nc = 1;

    if (sf_getstring ("vx")) {
        /* Horizontal velocity */
        velx = sf_input ("vx");
        nc++;
    }
    if (sf_getstring ("eta")) {
        /* Anellipticity */
        if (NULL == velx)
            sf_error ("Need vx=, if eta= is given");
        eta = sf_input ("eta");
        nc++;
    } else if (velx)
        sf_error ("Need eta=, if vx= is given");
    if (sf_getstring ("theta")) {
        /* Tilt angle elevation */
        if (NULL == velx)
            sf_error ("Need vx=, if theta= is given");
        theta = sf_input ("theta");
        nc++;
    }
    if (sf_getstring ("phi")) {
        /* Tilt angle azimuth */
        if (NULL == theta)
            sf_error ("Need theta=, if phi= is given");
        phi = sf_input ("phi");
        nc++;
    }

    z_grid.start = oz; z_grid.end = oz + (nz - 1)*dz; z_grid.num = nz;
    x_grid.start = ox; x_grid.end = ox + (nx - 1)*dx; x_grid.num = nx;
    y_grid.start = oy; y_grid.end = oy + (ny - 1)*dy; y_grid.num = ny;
    zBC.lCode = zBC.rCode = NATURAL;
    xBC.lCode = xBC.rCode = NATURAL;
    yBC.lCode = yBC.rCode = NATURAL;
    velspline = create_multi_UBspline_3d_s (y_grid, x_grid, z_grid, yBC, xBC, zBC, nc); 

    /* Read data and compute spline coefficients */
    if (verb)
        sf_warning ("Processing V_z");
    sf_floatread (buf, n, velz);
    if (1 == nc) {
        /* Isotropic case - convert velocity to slowness */
        if (verb)
            sf_warning ("Converting to slowness for isotropic case");
        for (i = 0; i < n; i++)
            buf[i] = 1.0/buf[i];
    } else {
        /* Convert to V_z^2 */
        for (i = 0; i < n; i++)
            buf[i] *= buf[i];
    }
    ic = 0;
    set_multi_UBspline_3d_s (velspline, ic, buf);
    ic++;
    if (velx) {
        if (verb)
            sf_warning ("Processing V_x");
        buf2 = sf_floatalloc (n);
        sf_floatread (buf2, n, velx);
        sf_fileclose (velx);
        /* Convert to V_x^2 */
//.........这里部分代码省略.........

int main (int argc,char* argv[]) 
{
    int i, j, ip, n1, n2, np, ndim, order, n123, *pp, n[2], box[2], **shift, b;
    float o1, o2, d1, d2, slow, *dd, **pts, *vv, *h, *bin, *vor, d[2], **rect;
    bool isvel, dist, voro;
    sf_upgrad upg;
    sf_file coord, ord, grid, vel;

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

    if (NULL != sf_getstring("velocity")) {
	vel = sf_input("velocity");

	if(!sf_histint(vel,"n1",&n1)) sf_error("No n1= in vel");
	if(!sf_histint(vel,"n2",&n2)) sf_error("No n2= in vel");
	/* dimensions */
	
	if(!sf_histfloat(vel,"d1",&d1)) sf_error("No d1= in vel");
	if(!sf_histfloat(vel,"d2",&d2)) sf_error("No d2= in vel");
	/* sampling */
	
	if(!sf_histfloat(vel,"o1",&o1)) o1=0.;
	if(!sf_histfloat(vel,"o2",&o2)) o2=0.;
	/* origin */
    } else {
	vel = NULL;

	if(!sf_getint("n1",&n1)) sf_error("Need n1=");
	if(!sf_getint("n2",&n2)) sf_error("Need n2=");
	/* dimensions */
	
	if(!sf_getfloat("d1",&d1)) sf_error("Need d1=");
	if(!sf_getfloat("d2",&d2)) sf_error("Need d2=");
	/* sampling */
	
	if(!sf_getfloat("o1",&o1)) o1=0.;
	if(!sf_getfloat("o2",&o2)) o2=0.;
	/* origin */
    }

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

    n[0]=n1; d[0]=d1;
    n[1]=n2; d[1]=d2;

    if(!sf_getint("order",&order)) order=2;
    /* [1,2] Accuracy order for distance calculation */

    if(!sf_getbool("vel",&isvel)) isvel=true;
    /* if y, the input is velocity; n, slowness squared */

    if (SF_FLOAT != sf_gettype(coord)) sf_error("Need float input");
    if(!sf_histint(coord,"n2",&np)) sf_error("No n2= in input");
    if(!sf_histint(coord,"n1",&ndim) || ndim > 3)  sf_error("Need n1 <= 3 in input");

    pts = sf_floatalloc2 (3,np);
    for (ip=0; ip < np; ip++) {
	sf_floatread(pts[ip],ndim,coord);
	pts[ip][2] = 0.0f;
    }
    
    n123 = n1*n2;

    dd = sf_floatalloc (n123);
    vv = sf_floatalloc (n123);
    pp = sf_intalloc (n123);

    if (NULL != vel) {
	sf_floatread(vv,n123,vel);
	sf_fileclose(vel);

	/* transform velocity to slowness squared */
	if (isvel) {
	    for(i = 0; i < n123; i++) {
		slow = vv[i];
		vv[i] = 1./(slow*slow);
	    }
	} 
    } else {
	for(i = 0; i < n123; i++) {
	    vv[i] = 1.;
	}
    }
    
    /* 1. find distance */
    distance_init (1,n2,n1,np);  
    distance(np,pts,dd,vv,pp,
	     1,n2,n1,
	     0.,o2,o1,
	     1.,d2,d1,
	     order);

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

int main (int argc, char* argv[])
{
    map4 nmo; /* using cubic spline interpolation */
    bool half, slow;
    int it,ix,ih, nt,nx, nh, CDPtype;
    float dt, t0, h, h0, f, dh, eps, dy;
    float *trace, *vel, *off, *str, *out;
    sf_file cmp, nmod, velocity, offset;

    sf_init (argc,argv);
    cmp = sf_input("in");
    velocity = sf_input("velocity");
    nmod = sf_output("out");

    if (SF_FLOAT != sf_gettype(cmp)) sf_error("Need float input");
    if (!sf_histint(cmp,"n1",&nt)) sf_error("No n1= in input");
    if (!sf_histfloat(cmp,"d1",&dt)) sf_error("No d1= in input");
    if (!sf_histfloat(cmp,"o1",&t0)) sf_error("No o1= in input");

    if (!sf_histint(cmp,"n2",&nh)) sf_error("No n2= in input");

    off = sf_floatalloc(nh);

    if (!sf_getbool("half",&half)) half=true;
    /* if y, the second axis is half-offset instead of full offset */


    CDPtype=1;
    if (NULL != sf_getstring("offset")) {
	offset = sf_input("offset");
	sf_floatread (off,nh,offset);
	sf_fileclose(offset);
    } else {
	if (!sf_histfloat(cmp,"d2",&dh)) sf_error("No d2= in input");
	if (!sf_histfloat(cmp,"o2",&h0)) sf_error("No o2= in input");
	
	if (sf_histfloat(cmp,"d3",&dy)) {
	    CDPtype=half? 0.5+dh/dy : 0.5+0.5*dh/dy;
	    if (CDPtype < 1) {
		CDPtype=1;
	    } else if (1 != CDPtype) {
		sf_histint(cmp,"CDPtype",&CDPtype);
	    	sf_warning("CDPtype=%d",CDPtype);
	    }
	} 

	for (ih = 0; ih < nh; ih++) {
	    off[ih] = h0 + ih*dh; 
	}
    }

    if (!sf_getbool("slowness",&slow)) slow=false;
    /* if y, use slowness instead of velocity */

    nx = sf_leftsize(cmp,2);

    if (!sf_getfloat ("h0",&h0)) h0=0.;
    /* reference offset */
    if (half) h0 *= 2.;
    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* stretch regularization */

    trace = sf_floatalloc(nt);
    vel = sf_floatalloc(nt);
    str = sf_floatalloc(nt);
    out = sf_floatalloc(nt);

    nmo = stretch4_init (nt, t0, dt, nt, eps);
    
    for (ix = 0; ix < nx; ix++) {
	sf_floatread (vel,nt,velocity);	

	for (ih = 0; ih < nh; ih++) {
	    sf_floatread (trace,nt,cmp);
	    
	    h = off[ih] + (dh/CDPtype)*(ix%CDPtype); 
	    if (half) h *= 2;
	    h = h*h - h0*h0;
	    
	    for (it=0; it < nt; it++) {
		f = t0 + it*dt;
		if (slow) {
		    f = f*f + h*vel[it]*vel[it];
		} else {
		    f = f*f + h/(vel[it]*vel[it]);
		}
		if (f < 0.) {
		    str[it]=t0-10.*dt;
		} else {
		    str[it] = sqrtf(f);
		}
	    }

	    stretch4_define (nmo,str);
	    stretch4_apply (false,nmo,trace,out);
	    
	    sf_floatwrite (out,nt,nmod);
	}
    }

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

//.........这里部分代码省略.........
        /* Adjoint I.C. operator, dI -> dW */
        wexcip_adj(cub,cip,Bwr,Pws,Pi,eic,1,1); /* Ws dR */
        wexcip_adj(cub,cip,Bws,Pwr,Pi,eic,0,0); /* Wr dR */

        sf_filefresh(Pws);
        sf_filefresh(Pwr);

        /* Adjoint WEMVA operator, dW -> dS */
        wexmva(mva,adj,cub,ssr,lsr,tap,slo,Bws,Bwr,Pws,Pwr,Ps);

    } else {
        /* set up the I/O of output CIP gathers */
        Pi = sf_output("out"); sf_settype(Pi,SF_COMPLEX);

        if(eic){
            /* CIP coordinates */
            Fc = sf_input ("cc" );
            ac = sf_iaxa(Fc,2); sf_setlabel(ac,"cc"); sf_setunit(ac,"");
            nc = sf_n(ac);

            if(! sf_getint("nhx",&nhx)) nhx=0; /* number of lags on the x axis */
            if(! sf_getint("nhy",&nhy)) nhy=0; /* number of lags on the y axis */
            if(! sf_getint("nhz",&nhz)) nhz=0; /* number of lags on the z axis */
            if(! sf_getint("nht",&nht)) nht=0; /* number of lags on the t axis */
            if(! sf_getfloat("dht",&dht)) sf_error("need dht");
            oht = -nht*dht;

            nhx2=2*nhx+1; nhy2=2*nhy+1;
            nhz2=2*nhz+1; nht2=2*nht+1;

            aa=sf_maxa(nhx2,-nhx*cub->amx.d,cub->amx.d);
            sf_setlabel(aa,"hx"); sf_setunit(aa,"");
            if(verb) sf_raxa(aa);
            sf_oaxa(Pi,aa,1);

            aa=sf_maxa(nhy2,-nhy*cub->amy.d,cub->amy.d);
            sf_setlabel(aa,"hy"); sf_setunit(aa,"");
            if(verb) sf_raxa(aa);
            sf_oaxa(Pi,aa,2);

            aa=sf_maxa(nhz2,-nhz*cub->az.d,cub->az.d);
            sf_setlabel(aa,"hz"); sf_setunit(aa,"");
            if(verb) sf_raxa(aa);
            sf_oaxa(Pi,aa,3);

            aa=sf_maxa(nht2,-nht*dht,dht);
            sf_setlabel(aa,"ht"); sf_setunit(aa,"s");
            if(verb) sf_raxa(aa);
            sf_oaxa(Pi,aa,4);
 
            sf_oaxa(Pi,ac,5);

        }
        else{
            sf_oaxa(Pi,amx,1);
            sf_oaxa(Pi,amy,2);
            sf_oaxa(Pi,az, 3);
        }

        cip = wexcip_init(cub,nhx,nhy,nhz,nht,nhx2,nhy2,nhz2,nht2,nc,dht,oht,Fc,eic);
        mva = wexmva_init(cub,cip);
 
        /* WEMVA operator, dS -> dW */
        wexmva(mva,adj,cub,ssr,lsr,tap,slo,Bws,Bwr,Pws,Pwr,Ps);
         
        sf_filefresh(Pws);
        sf_filefresh(Pwr);

        /* I.C. operator, dW -> dI */
        wexcip_for(cub,cip,Bws,Pwr,Pti,eic,0,0); /* CONJ( Ws) dWr */
        sf_seek(Pti,(off_t)0,SEEK_SET);
        wexcip_for(cub,cip,Pws,Bwr,Pti,eic,0,1); /* CONJ(dWs)  Wr */

        sf_filefresh(Pti);
        sf_filecopy(Pi,Pti,SF_COMPLEX);
    }

    /*------------------------------------------------------------*/
    /* close structures */
    wexslo_close(slo);
    wexssr_close(cub,ssr);
    wextap2D_close(tap);
    wexmva_close(mva);
    wexcip_close(cip,eic);
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    /* close files */
    if (Ps!=NULL) sf_fileclose(Ps);
    if (Fc!=NULL) sf_fileclose(Fc);
    if (Pi!=NULL) sf_fileclose(Pi);
    if (Bws!=NULL) sf_fileclose(Bws);
    if (Bwr!=NULL) sf_fileclose(Bwr);
    if (Pws!=NULL) sf_tmpfileclose(Pws);
    if (Pwr!=NULL) sf_tmpfileclose(Pwr);
    if (Pti!=NULL) sf_tmpfileclose(Pti);
    /*------------------------------------------------------------*/

    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 n[SF_MAX_DIM], a[SF_MAX_DIM], center[SF_MAX_DIM], gap[SF_MAX_DIM];
    int *pch, *nh, dim, n123, nf, i, niter, nbf, nbp, id, ip, ig, np;
    int *kk, *pp;
    float *dd, eps, dabs, di;
    nfilter aa, bb;
    char varname[6], *lagfile;
    sf_file in, flt, lag, mask, patch, reg;

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

    dim = sf_filedims(in,n);

    if (NULL == (lagfile = sf_getstring("lag"))) sf_error("Need lag=");
    /* output file for filter lags */

    lag = sf_output(lagfile);
    sf_settype(lag,SF_INT);

    sf_putstring(flt,"lag",lagfile);

    sf_putints(lag,"n",n,dim);

    if (!sf_getints("a",a,dim)) sf_error("Need a=");

    if (!sf_getints("center",center,dim)) {
        for (i=0; i < dim; i++) {
            center[i] = (i+1 < dim && a[i+1] > 1)? a[i]/2: 0;
        }
    }

    if (!sf_getints("gap",gap,dim)) {
        for (i=0; i < dim; i++) {
            gap[i] = 0;
        }
    }

    n123 = 1;
    for (i=0; i < dim; i++) {
        n123 *= n[i];
    }

    dd = sf_floatalloc(n123);
    kk = sf_intalloc(n123);

    if (NULL != sf_getstring("maskin")) {
        /* optional input mask file */
        mask = sf_input("maskin");

        switch (sf_gettype(mask)) {
        case SF_INT:
            sf_intread (kk,n123,mask);
            break;
        case SF_FLOAT:
            sf_floatread (dd,n123,mask);
            for (i=0; i < n123; i++) {
                kk[i] = (dd[i] != 0.);
            }
            break;
        default:
            sf_error ("Wrong data type in maskin");
            break;
        }

        sf_fileclose (mask);
    } else {
        for (i=0; i < n123; i++) {
            kk[i] = 1;
        }
    }

    sf_floatread(dd,n123,in);

    dabs = fabsf(dd[0]);
    for (i=1; i < n123; i++) {
        di = fabsf(dd[i]);
        if (di > dabs) dabs=di;
    }

    random_init(2004);
    for (i=0; i < n123; i++) {
        dd[i] = dd[i]/dabs+ 100.*FLT_EPSILON*(random0()-0.5);;
    }

    pp = sf_intalloc(n123);
    if (NULL != sf_getstring("pch")) {
        patch = sf_input("pch");
        if (SF_INT != sf_gettype(patch)) sf_error("Need int pch");

        sf_intread(pp,n123,patch);

        np = pp[0];
        for (i=1; i < n123; i++) {
            if (pp[i] > np) np = pp[i];
        }

        sf_fileclose(patch);
//.........这里部分代码省略.........