int main(int argc, char* argv[])
{
    int n2, i2, i1, ic, id, nc, n;
    double xp;
    float *roots, *sol, *dat, x, o, d;
    sf_file inp, coef, out;

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

    if (!sf_histint(inp,"n1",&nc))   sf_error("No n1= in input");
    nc++;
    n2 = sf_leftsize(inp,1);

    if (NULL != sf_getstring("coef")) {
        /* (optional) coefficients */
	coef = sf_output("coef"); 
	sf_putint(coef,"n1",nc);
    } else {
	coef = NULL;
    }

    roots = sf_floatalloc(nc-1);

    out = sf_output("out");
    if (!sf_getint("n1",&n)) sf_error("Need n1=");
    /* number of samples */
    if (!sf_getfloat("d1",&d)) sf_error("Need d1=");
    /* sampling */
    if (!sf_getfloat("o1",&o)) sf_error("Need o1=");
    /* origin */
    sf_putint(out,"n1",n);
    sf_putfloat(out,"d1",d);
    sf_putfloat(out,"o1",o);

    dat = sf_floatalloc(n);
    sol = sf_floatalloc(nc);
    sol[0]=1.0;

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(roots,nc-1,inp);
	
	for (ic=1; ic < nc; ic++) {
	    sol[ic]=0.0;
	}

	for (ic=0; ic < nc-1; ic++) {
	    for (id=0; id <= ic; id++) {
		sol[id+1] -= roots[ic]*sol[id];
	    }
	}

	if (NULL != coef) sf_floatwrite(sol,nc,coef);

	for (i1=0; i1 < n; i1++) {
	    dat[i1] = 0.;
	    x = o+i1*d;
	    xp = 1.0;
	    for (ic=0; ic < nc; ic++) {
		dat[i1] += xp*sol[ic];
		xp *= x;
	    }
	}

	sf_floatwrite(dat,n,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int niter, nd, n1, n2, i1, i2;
    float **vr, **vi, **wt, **v0, **bk, *tmp, wti, perc;
    sf_file vrms, vint, weight, vout, block;

    sf_init(argc,argv);
    vrms = sf_input("in");
    vint = sf_output("out");
    weight = sf_input("weight");
    block = sf_input("block");

    if (!sf_histint(vrms,"n1",&n1)) sf_error("No n1= in input");
    n2 = sf_leftsize(vrms,1);
    nd = n1*n2;

    vr = sf_floatalloc2(n1,n2);
    vi = sf_floatalloc2(n1,n2);
    wt = sf_floatalloc2(n1,n2);
    v0 = sf_floatalloc2(n1,n2);
    bk = sf_floatalloc2(n1,n2);
    tmp = sf_floatalloc(n1);

    sf_floatread(vr[0],nd,vrms);
    sf_floatread(wt[0],nd,weight);
    sf_floatread(bk[0],nd,block);

    if (!sf_getfloat("perc",&perc)) perc=50.0;
    /* percentage for sharpening */

    blockder_init(n1, n2, perc, bk[0], wt[0]);

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

    wti = 0.;
    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    wti += wt[i2][i1]*wt[i2][i1];
	}
    }
    if (wti > 0.) wti = sqrtf(n1*n2/wti);

    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    vr[i2][i1] *= vr[i2][i1]*(i1+1.0f); /* vrms^2*t - data */
	    wt[i2][i1] *= wti/(i1+1.0f); /* decrease weight with time */	 
	    v0[i2][i1] = -vr[i2][0];
	}
	sf_causint_lop(false,true,n1,n1,v0[i2],vr[i2]);
    }
    
    blockder(niter, vr[0], vi[0]);
 
    for (i2=0; i2 < n2; i2++) {
	sf_causint_lop(false,false,n1,n1,vi[i2],tmp);
	for (i1=0; i1 < n1; i1++) {
	    vi[i2][i1] = tmp[i1] - v0[i2][i1];
	}
	sf_causint_lop(false,false,n1,n1,vi[i2],vr[i2]);
    }

    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    vr[i2][i1] = sqrtf(fabsf(vr[i2][i1]/(i1+1.0f)));
	    vi[i2][i1] = sqrtf(fabsf(vi[i2][i1]));
	}
    }

    sf_floatwrite(vi[0],nd,vint);

    if (NULL != sf_getstring("vrmsout")) {
	/* optionally, output predicted vrms */
	vout = sf_output("vrmsout");

	sf_floatwrite(vr[0],nd,vout);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, i1, i2, niter;
    float **b, *h, ***rt, d, maxd;
    sf_complex **data;
    sf_file inp, out, bad;

    sf_init (argc,argv);
    inp = sf_input("in");
    out = sf_output("out");
    
    if (SF_COMPLEX != sf_gettype(inp)) sf_error("Need complex input");
    if (!sf_histint(inp,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(inp,"n2",&n2)) sf_error("No n2= in input");

    sf_settype(out,SF_FLOAT);

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

    if (NULL != sf_getstring("badness")) {
	bad = sf_output("badness");
	/* (optional) badness attribute file */
	sf_settype(bad,SF_FLOAT);
    } else {
	bad = NULL;
    }

    data = sf_complexalloc2(n1,n2);
    sf_complexread(data[0],n1*n2,inp);

    /* normalize */
    maxd = 0;
    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    d = cabsf(data[i2][i1]);
	    if (maxd < d) maxd=d;
	}
    }
    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
#ifdef SF_HAS_COMPLEX_H
	    data[i2][i1] = data[i2][i1]/maxd;
#else
	    data[i2][i1] = sf_crmul(data[i2][i1],1.0f/maxd);
#endif
	}
    }

    if (NULL != bad) {
	b = sf_floatalloc2(n1,n2);
	rt = sf_floatalloc3(n1,n2,2);

	grad2init (n1, n2, data, rt);

	for (i2=0; i2 < n2; i2++) {
	    b[i2][n1-1] = 0.;
	    b[i2][n1-2] = 0.;
	}

	for (i1=0; i1 < n1; i1++) {
	    b[n2-1][i1] = 0.;
	    b[n2-2][i1] = 0.;
	}
	
	for (i2=0; i2 < n2-2; i2++) {
	    for (i1=0; i1 < n1-2; i1++) {
		b[i2][i1] = 
		    (rt[0][i2+1][i1] - rt[0][i2][i1]) - 
		    (rt[1][i2][i1+1] - rt[1][i2][i1]);
	    }
	}

	sf_floatwrite(b[0],n1*n2,bad);
    }

    if (niter > 0) {
	h = sf_floatalloc(n1*n2);
	unwraper (n1,n2, data, h, niter);
	sf_floatwrite(h,n1*n2,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    bool verb;
    int n[SF_MAX_DIM], n0[SF_MAX_DIM], rect[SF_MAX_DIM];
    int a[SF_MAX_DIM], center[SF_MAX_DIM], gap[SF_MAX_DIM];
    int ndim, dim, n123, n123s, i, ia, ns, i1, niter, na, i4, n4, *kk;
    float *d, *f, *dd;
    double mean;
    char *lagfile, key[6];
    sf_filter aa;
    sf_file in, filt, lag;
 
    sf_init(argc,argv);

    in = sf_input("in");
    filt = sf_output("out");

    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(filt,"lag",lagfile);

    ndim = sf_filedims(in,n);

    if (!sf_getint("dim",&dim)) dim=ndim; /* number of dimensions */

    n4 = sf_leftsize(in,dim);
    
    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_getint("na",&na)) na=0;
    /* filter size */

    if (0 == na) {
	if (!sf_getints("gap",gap,dim)) {
	    for (i=0; i < dim; i++) {
		gap[i] = 0;
	    }
	}
	
	aa = createhelix(dim, n, center, gap, a); /* allocate PEF */
	
	for (i=0; i < dim; i++) {	    
	    n0[i] = n[i];
	}
    } else {
	aa =  sf_allocatehelix (na);
	if (!sf_getints ("lags", aa->lag, na)) sf_error("Need lags=");
	if (!sf_getints ("n", n0, dim)) {
	    for (i=0; i < dim; i++) {	    
		n0[i] = n[i];
	    }
	}
    }

    n123 = 1;
    for (i=0; i < dim; i++) {
	n123 *= n[i];
	
	snprintf(key,6,"rect%d",i+1);
	if (!sf_getint(key,rect+i)) rect[i]=1;
    }

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

    for (i1=0; i1 < n123; i1++) {
	kk[i1] = 1;
    }

    bound (dim, n0, n, a, aa);
    find_mask(n123, kk, aa);

    na = aa->nh;

    snprintf(key,3,"n%d",dim+1);
    sf_putint(filt,key,na);
    sf_shiftdim(in, filt, dim+1);    
    
    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */

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

    n123s = n123*na;
    
    d = sf_floatalloc(n123s);
    f = sf_floatalloc(n123s);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool verb,pas,adj,abc;          /* execution flags */
    int ix, iz, it;                 /* index variables */
    int nt, nx, nz, depth, nzxpad, nb, n2, snap;
    float ox, oz, dx, dz, dt, dt2, idz2, idx2, cb;

    int nxpad, nzpad;
    float **vvpad;
    
    float **dd, **mm, **vv, ***ww;
    float **u0, **u1, **u2, **tmp;  /* temporary arrays */

    sf_file in, out, vel, wave;     /* I/O files */

    /* initialize Madagascar */
    sf_init(argc,argv);
    
    /* initialize OpenMP support */
#ifdef _OPENMP
    omp_init();
#endif

    if(!sf_getbool("verb", &verb)) verb=false;
    /* verbosity flag */
    if(!sf_getbool("adj", &adj)) adj=false;
    /* adjoint flag, 0: modeling, 1: migration */
    if(!sf_getbool("pas", &pas)) pas=false;
    /* passive flag, 0: exploding reflector rtm, 1: passive seismic imaging */
    if(!sf_getbool("abc",&abc)) abc = false;
    /* absorbing boundary condition */
    if(!sf_getint("snap", &snap)) snap=0;
    /* wavefield snapshot flag */
    if(!sf_getint("depth", &depth)) depth=0;
    /* surface */
	
    /* setup I/O files */
    in  = sf_input("in");
    out = sf_output("out");
    vel = sf_input("velocity");
    /* velocity model */
    
    /* Dimensions */
    if(!sf_histint  (vel, "n1", &nz)) sf_error("No n1= in velocity");
    if(!sf_histint  (vel, "n2", &nx)) sf_error("No n2= in velocity");
    if(!sf_histfloat(vel, "o1", &oz)) sf_error("No o1= in velocity");
    if(!sf_histfloat(vel, "o2", &ox)) sf_error("No o2= in velocity");
    if(!sf_histfloat(vel, "d1", &dz)) sf_error("No d1= in velocity");
    if(!sf_histfloat(vel, "d2", &dx)) sf_error("No d2= in velocity");

    if(adj){ /* migration */
        if(!sf_histint(in, "n1", &nt)) sf_error("No n1= in data");
        if(!sf_histfloat(in, "d1", &dt)) sf_error("No d1= in data");
        if(!sf_histint(in, "n2", &n2) || n2!=nx) sf_error("Need n2=%d in data", nx);

        sf_putint   (out, "n1", nz);
        sf_putfloat (out, "o1", oz);
        sf_putfloat (out, "d1", dz);
        sf_putstring(out, "label1", "Depth");
        sf_putstring(out, "unit1" , "km");
        sf_putint   (out, "n2", nx);
        sf_putfloat (out, "o2", ox);
        sf_putfloat (out, "d2", dx);
        sf_putstring(out, "label2", "Distance");
        sf_putstring(out, "unit2" , "km");
        if (pas) {
            sf_putint   (out, "n3", nt);
            sf_putfloat (out, "d3", dt);
            sf_putfloat (out, "o3", 0.0f);
            sf_putstring(out, "label3", "Time");
            sf_putstring(out, "unit3" , "s");
        }
    }else{ /* modeling */
        if(!sf_getint("nt", &nt)) sf_error("Need nt=");
        if(!sf_getfloat("dt", &dt)) sf_error("Need dt=");
        
        sf_putint   (out, "n1", nt);
        sf_putfloat (out, "d1", dt);
        sf_putfloat (out, "o1", 0.0);
        sf_putstring(out, "label1", "Time");
        sf_putstring(out, "unit1" , "s");
        sf_putint   (out, "n2", nx);
        sf_putfloat (out, "o2", ox);
        sf_putfloat (out, "d2", dx);
        sf_putstring(out, "label2", "Distance");
        sf_putstring(out, "unit2" , "km");
        if (pas) {
            sf_putint   (out, "n3", 1);
        }
    }
    
    /* dimension of padded boundary */
    if(!sf_getint("nb", &nb) || nb<NOP) nb = NOP;
    if(!sf_getfloat("cb", &cb)) cb = 0.0f;
    nxpad = nx+2*nb;
    nzpad = nz+2*nb;
    nzxpad = nzpad*nxpad;
    depth = depth+nb;

    /* set Laplacian coefficients */
//.........这里部分代码省略.........

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


  int nt, nx, it, ix, niter, iter, ntfft, nxfft,np, ip, ikt, ikx, iktn, ikxn, ifsnr; /* iktn, ikxn, iNyquist*/
  float dt, dx, pmin, pmax, dp, p, cmax, scalar, sembpmax, num, den;
  float *sembp, *mask, *gy, *fden, *fshift, *SNR;
  float **fdata, **taup, **odata, **tdata, **odatat, **semb; /* tdata is the true data */
  kiss_fft_cpx **cdata, **cdatat;
  char *type;
  sf_file inp, outp, m, spec1, spec2, trued, snr; 


  sf_init(argc,argv);

  inp=sf_input("in");
  m=sf_input("mask");
  outp=sf_output("out");

  if(!sf_histint(inp,"n1",&nt)) sf_warning("No n1 in input");
  if(!sf_histint(inp,"n2",&nx)) sf_warning("No n2 in input");
  if(!sf_histfloat(inp,"d1",&dt)) sf_warning("No n1 in input");
  if(!sf_histfloat(inp,"d2",&dx)) sf_warning("No n2 in input");

  ntfft = 2*kiss_fft_next_fast_size((nt+1)/2);
  nxfft = 2*kiss_fft_next_fast_size((nx+1)/2);
  scalar = 1./(ntfft*nxfft);
  iktn=ntfft/2; ikxn=nxfft/2;
  float dkt = 1.0/(ntfft*dt), fkt = 0.0,kt;
  float dkx = 1.0/(nxfft*dx), fkx = 0.0,kx;

    if (NULL == (type=sf_getstring("type"))) type="amplitude";
    /* [amplitude, semblance] thresholding type, the default is amplitude thresholding  */ 

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

  	if(!sf_getint("ifsnr",&ifsnr)) ifsnr = 0;
	/* If compute SNR during iteration */


	if(type[0]=='s')
	{

  		if(!sf_getfloat("pmin",&pmin)) pmin=-2;
        /* minimum p */		
  		if(!sf_getfloat("pmax",&pmin)) pmax=2;
        /* maximum p */			
  		if(!sf_getint("np",&np)) np=nx;
        /* number of p */
		dp=(pmax-pmin)/(np-1);		
  		sembp =sf_floatalloc(np);
  		semb  =sf_floatalloc2(nt,np);
 		taup  =sf_floatalloc2(nt,np);	
	}

	/* output files */
  	if (NULL!=sf_getstring("spec2")) 
  	{
		spec2=sf_output("spec2");
		sf_putint(spec2, "n1", ntfft);
		sf_putint(spec2, "n2", nxfft);
	}	       		
  	if (NULL!=sf_getstring("spec1")) 
  	{
		spec1=sf_output("spec1");
		sf_putint(spec1, "n1", ntfft);
		sf_putint(spec1, "n2", nxfft);
	}	  
  	if (ifsnr==1 && (NULL!=sf_getstring("true"))) 
  	{
		snr=sf_output("snr");
		trued=sf_input("true");
		tdata=sf_floatalloc2(nt,nx);
		SNR=sf_floatalloc(niter);
		sf_floatread(tdata[0],nt*nx,trued);

		sf_putint(snr,"n1",niter);
		sf_putint(snr,"d1",1);
		sf_putint(snr,"n2",1);
	}	

  /* Allocate memory */
  cdata =(kiss_fft_cpx**) sf_complexalloc2(ntfft,nxfft);
  cdatat =(kiss_fft_cpx**) sf_complexalloc2(ntfft,nxfft); /* temporary file */
  fshift= sf_floatalloc(ntfft);
  fden 	= sf_floatalloc(ntfft);
  gy 	= sf_floatalloc(nxfft);
  mask  =sf_floatalloc(nx);

  odata =sf_floatalloc2(nt,nx); 
  odatat =sf_floatalloc2(ntfft,nxfft); 
  fdata =sf_floatalloc2(ntfft,nxfft);
  memset(&odata[0][0],0,ntfft*nxfft*sizeof(float)); 

  /* Read data */
  sf_floatread(odata[0],nt*nx,inp);
  sf_floatread(mask,nx,m);

	if(type[0]=='s')
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    int n1, n2, n3, i1, i2, i3, j1, j2, k1, k2, nedge, nold, nmin, nmax, n12;
    int **edge;
    float **pp, **ww, **w1, **w2, g1, g2, w, min, max;
    sf_file in=NULL, out=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_getfloat("min",&min)) min=5.0;
    /* minimum threshold */
    if (!sf_getfloat("max",&max)) max=95.0;
    /* maximum threshold */

    n12 = n1*n2;
    nmin = min*0.01*n12;
    if (nmin < 0) nmin=0;
    if (nmin >= n12) nmin=n12-1;
    nmax = max*0.01*n12;
    if (nmax < 0) nmax=0;
    if (nmax >= n12) nmax=n12-1;

    pp = sf_floatalloc2(n1,n2);
    w1 = sf_floatalloc2(n1,n2);
    w2 = sf_floatalloc2(n1,n2);
    ww = sf_floatalloc2(n1,n2);
    edge = sf_intalloc2(n1,n2);

    sf_settype(out,SF_INT);

    for (i3=0; i3 < n3; i3++) {
	sf_floatread(pp[0],n12,in);
	/* gradient computation */
	sf_sobel(n1,n2,pp,w1,w2);

	for (i2=0; i2 < n2; i2++) {
	    for (i1=0; i1 < n1; i1++) {
		/* gradient norm */
		g1 = w1[i2][i1];
		g2 = w2[i2][i1];
		ww[i2][i1] = g1*g1+g2*g2;
	    }
	}
	/* edge thinning */
	for (i2=0; i2 < n2; i2++) {
	    for (i1=0; i1 < n1; i1++) {
		g1 = w1[i2][i1];
		g2 = w2[i2][i1];
		if (fabsf(g1) > fabsf(g2)) {
		    j1=1;
		    if (g2/g1 > 0.5) { 
			j2=1;
		    } else if (g2/g1 < - 0.5) {
			j2=-1; 
		    } else {
			j2=0;
		    }
		} else if (fabsf(g2) > fabsf(g1)) {
		    j2=1;
		    if (g1/g2 > 0.5) { 
			j1=1;
		    } else if (g1/g2 < - 0.5) {
			j1=-1; 
		    } else {
			j1=0;
		    }
		} else {
		    j1=0;
		    j2=0;
		}
		k1 = i1+j1; if (j1 && (k1 < 0 || k1 >= n1)) k1=i1;
		k2 = i2+j2; if (j2 && (k2 < 0 || k2 >= n2)) k2=i2;
		if (ww[i2][i1] <= ww[k2][k1]) {
		    pp[i2][i1] = 0.;
		    continue;
		} 
		k1 = i1-j1; if (k1 < 0 || k1 >= n1) k1=i1;
		k2 = i2-j2; if (k2 < 0 || k2 >= n2) k2=i2;
		if (ww[i2][i1] <= ww[k2][k1]) {
		    pp[i2][i1] = 0.;
		    continue;
		} 
		pp[i2][i1] = ww[i2][i1];
	    }
	}
	/* edge selection */
	max = sf_quantile(nmax,n12,ww[0]);
	min = sf_quantile(nmin,n12,ww[0]);

	nedge=0;
	for (i2=0; i2 < n2; i2++) {
	    for (i1=0; i1 < n1; i1++) {
		w = pp[i2][i1];
//.........这里部分代码省略.........

int main (int argc, char* argv[])
{
    map4 nmo; /* using cubic spline interpolation */
    bool half;
    int it,ix,iy, nt, nx, ny, nw, i4, n4, n;
    float dt, t0, x, y, x0, y0, f, dx, dy, eps;
    float *trace=NULL, *vx=NULL, *vy=NULL, *vxy=NULL, *str=NULL, *out=NULL;
    sf_file cmp=NULL, nmod=NULL, vel=NULL;

    sf_init (argc,argv);
    cmp = sf_input("in");
    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",&nx)) sf_error("No n2= in input");
    if (!sf_histint(cmp,"n3",&ny)) sf_error("No n3= in input");
    n4 = sf_leftsize(cmp,3);

    vel = sf_input("velocity");
    if (SF_FLOAT != sf_gettype(vel)) sf_error("Need float velocity");
    if (!sf_histint(vel,"n1",&n) || nt != n) sf_error("Need n1=%d in velocity",nt);
    if (!sf_histint(vel,"n2",&n) || 3 != n) sf_error("Need n2=3 in velocity",nt);
    if (n4 != sf_leftsize(vel,2)) sf_error("Wrong dimensions in velocity");

    if (!sf_histfloat(cmp,"d2",&dx)) sf_error("No d2= in input");
    if (!sf_histfloat(cmp,"o2",&x0)) sf_error("No o2= in input");
    if (!sf_histfloat(cmp,"d3",&dy)) sf_error("No d3= in input");
    if (!sf_histfloat(cmp,"o3",&y0)) sf_error("No o3= in input");

    if (!sf_getbool("half",&half)) half=true;
    /* if y, the second and third axes are half-offset instead of full offset */

    if (half) { /* get full offset - multiply by 2 */
      dx *= 2.;
      dy *= 2.;
      x0 *= 2.;
      y0 *= 2.;
      sf_warning("Since half=y, offsets were doubled.");
    }else{
      sf_warning("Since half=n, offsets not doubled.");
    }

    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* stretch regularization */

    trace = sf_floatalloc(nt);
    vx = sf_floatalloc(nt);
    vy = sf_floatalloc(nt);
    vxy = sf_floatalloc(nt);

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

    if (!sf_getint("extend",&nw)) nw=8;
    /* trace extension */

    nmo = stretch4_init (nt, t0, dt, nt, eps);

    for (i4=0; i4 < n4; i4++) { /* loop over cmps */
      sf_floatread (vx,nt,vel);
      sf_floatread (vy,nt,vel);
      sf_floatread (vxy,nt,vel);

      for (iy = 0; iy < ny; iy++) {
	y = y0+iy*dy;

	for (ix = 0; ix < nx; ix++) {
	  x = x0 + ix*dx;

	  sf_floatread (trace,nt,cmp);

	  for (it=0; it < nt; it++) {
	    f = t0 + it*dt;
	    f = f*f + x*x*vx[it]+y*y*vy[it]+2*x*y*vxy[it];

	    if (f < 0.) {
	      str[it]=t0-10.*dt;
	    } else {
	      str[it] = sqrtf(f);
	    }
	  }

	  stretch4_define (nmo,str);
	  stretch4_apply (nmo,trace,out);

	  sf_floatwrite (out,nt,nmod);
	} /* ix */
      } /* iy */
    } /* i4 */

    exit (0);
}

int main(int argc, char* argv[])
{
    bool adj, half, verb, normalize;
    int nt, nx, nh, nh2, ix, ih, iy, i, nn, it, **fold, apt;
    float *trace, **image, **v, rho, **stack, *pp, *off;
    float h, x, t, h0, dh, dx, ti, tx, t0, t1, t2, dt, vi, aal, angle;
    sf_file inp, out, vel, gather, offset;

    sf_init (argc,argv);
    inp = sf_input("in");
    vel = sf_input("vel");  
    out = sf_output("out");
       
    if (!sf_getbool("adj",&adj)) adj=true;
    /* adjoint flag (y for migration, n for modeling) */

    if (!sf_getbool("normalize",&normalize)) normalize=true;
    /* normalize for the fold */
 

    if (!sf_histint(inp,"n1",&nt)) sf_error("No n1=");
    if (!sf_histint(inp,"n2",&nx)) sf_error("No n2=");

    if (!sf_histfloat(inp,"o1",&t0)) sf_error("No o1=");
    if (!sf_histfloat(inp,"d1",&dt)) sf_error("No d1=");
    if (!sf_histfloat(inp,"d2",&dx)) sf_error("No d2=");

    if (adj) {
	if (!sf_histint(inp,"n3",&nh)) sf_error("No n3=");
       
	sf_putint(out,"n3",1);
    } else {
	if (!sf_getint("nh",&nh)) sf_error("Need nh=");
	/* number of offsets (for modeling) */
	
	sf_putint(out,"n3",nh);
    }	

    if (NULL != sf_getstring("gather")) {
	gather = sf_output("gather");
    } else {
	gather = NULL;
    }

    if (!sf_getfloat("antialias",&aal)) aal = 1.0;
    /* antialiasing */

    if (!sf_getint("apt",&apt)) apt = nx;
    /* integral aperture */

    if (!sf_getfloat("angle",&angle)) angle = 90.0;
    /* angle aperture */

    angle = fabsf(tanf(angle*SF_PI/180.0));

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

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

    if (!sf_getfloat("rho",&rho)) rho = 1.-1./nt;
    /* Leaky integration constant */

    if (NULL != sf_getstring("offset")) {
	offset = sf_input("offset");
	nh2 = sf_filesize(offset);

	if (nh2 != nh*nx) sf_error("Wrong dimensions in offset");

	off = sf_floatalloc(nh2);	
	sf_floatread (off,nh2,offset);
	sf_fileclose(offset);
    } else {
	if (adj) {
	    if (!sf_histfloat(inp,"o3",&h0)) sf_error("No o3=");
	    if (!sf_histfloat(inp,"d3",&dh)) sf_error("No d3=");
	    sf_putfloat(out,"d3",1.);
	    sf_putfloat(out,"o3",0.);
	} else {
	    if (!sf_getfloat("dh",&dh)) sf_error("Need dh=");
	    /* offset sampling (for modeling) */
	    if (!sf_getfloat("h0",&h0)) sf_error("Need h0=");
	    /* first offset (for modeling) */
	    sf_putfloat(out,"d3",dh);
	    sf_putfloat(out,"o3",h0);
	}
	
	if (!half) dh *= 0.5;

	off = sf_floatalloc(nh*nx);
	for (ix = 0; ix < nx; ix++) {
	    for (ih = 0; ih < nh; ih++) {
		off[ih*nx+ix] = h0 + ih*dh; 
	    }
	}
	offset = NULL;
    }

    v = sf_floatalloc2(nt,nx);
//.........这里部分代码省略.........

    int main(int argc, char* argv[])
{
    int i, n1, n2, n12, nj1, nj2, niter, nw, n3, i3;
    float eps, *d, *s, *nd, **nn, **ss;
    bool verb;
    sf_file in, out, ndip, sdip;

    sf_init (argc,argv);
    in = sf_input("in");
    ndip = sf_input("ndip");
    sdip = sf_input("sdip");
    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");
    n12 = n1*n2;

    n3 = sf_leftsize(in,2);
    if (1==n3) sf_putint (out,"n3",2);

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

    if (!sf_getfloat ("eps",&eps)) eps=1.;
    /* regularization parameter */

    if (!sf_getint("order",&nw)) nw=1;
    /* [1,2,3] accuracy order */
    if (nw < 1 || nw > 3) 
	sf_error ("Unsupported nw=%d, choose between 1 and 3",nw);
    if (!sf_getint("nj1",&nj1)) nj1=1;
    /* antialiasing for noise dip */
    if (!sf_getint("nj2",&nj2)) nj2=1;
    /* antialiasing for signal dip */

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

    nd = sf_floatalloc(2*n12);
    s = sf_floatalloc(n12);
    d = sf_floatalloc(n12);
    nn = sf_floatalloc2(n1,n2);
    ss = sf_floatalloc2(n1,n2);

    for (i=0; i < n12; i++) {
	nd[n12+i] = 0.;
    }

    planesignoi_init (nw, nj1,nj2, n1,n2, nn, ss, eps);

    for (i3=0; i3 < n3; i3++) {
	sf_floatread (d,n12,in);
	sf_floatread (nn[0],n12,ndip);
	sf_floatread (ss[0],n12,sdip);

	allpass22_init(allpass2_init(nw,nj1,n1,n2,nn));
	allpass21_lop (false,false,n12,n12,d,nd);
	

	sf_solver (planesignoi_lop, sf_cgstep, n12, n12*2, s, nd, niter, 
		   "verb", verb, "end");

	sf_cgstep_close();

	sf_floatwrite(s,n12,out);

	if (1==n3) {
	    for (i=0; i < n12; i++) {
		d[i] -= s[i];
	    }
	    sf_floatwrite(d,n12,out);
	}
    }

    exit(0);
}

int main(int argc, char* argv[])
{
  int verbose;
  sf_file in=NULL, out=NULL;
  int n1_traces;
  int n1_headers;

  char* header_format=NULL;
  sf_datatype typehead;
  /* kls do I need to add this?  sf_datatype typein; */
  float* fheader=NULL;
  float* intrace=NULL;
  float* fprevheader=NULL;
  int numkeys;
  int ikey;
  char** list_of_keys;
  int *indx_of_keys;
  char* skey;
  int indx_of_skey;
  int skeyvalue;
  bool pkeychanged;
  int itrace=0;

  /*****************************/
  /* initialize verbose switch */
  /*****************************/
  sf_init (argc,argv);

  /* verbose flag controls ammount of print */
  /*( verbose=1 0 terse, 1 informative, 2 chatty, 3 debug ) */
  /* fprintf(stderr,"read verbose switch.  getint reads command line.\n"); */
  if(!sf_getint("verbose",&verbose))verbose=1;
    /* \n
     flag to control amount of print
     0 terse, 1 informative, 2 chatty, 3 debug
  */
  sf_warning("verbose=%d",verbose);
 
  /******************************************/
  /* input and output data are stdin/stdout */
  /******************************************/

  if(verbose>0)fprintf(stderr,"read in file name\n");  
  in = sf_input ("in");

  if(verbose>0)fprintf(stderr,"read out file name\n");
  out = sf_output ("out");

  if (!sf_histint(in,"n1_traces",&n1_traces))
    sf_error("input data not define n1_traces");
  if (!sf_histint(in,"n1_headers",&n1_headers)) 
    sf_error("input data does not define n1_headers");

  header_format=sf_histstring(in,"header_format");
  if(strcmp (header_format,"native_int")==0) typehead=SF_INT;
  else                                       typehead=SF_FLOAT;

  if(verbose>0)fprintf(stderr,"allocate headers.  n1_headers=%d\n",n1_headers);
  fheader = sf_floatalloc(n1_headers);
  fprevheader = sf_floatalloc(n1_headers);
 
  if(verbose>0)fprintf(stderr,"allocate intrace.  n1_traces=%d\n",n1_traces);
  intrace= sf_floatalloc(n1_traces);

  if(verbose>0)fprintf(stderr,"call list of keys\n");
 
  /* this sf_getstring will create parameter descrpiton in the self doc */
  sf_getstring("pkey"); 
  /* \n
     A comma seperated list of primary header keys to monitor to determine 
     gathers.  The trace number in the gather is counted and put in the
     skey header location.
     \n
  */
  list_of_keys=sf_getnstring("pkey",&numkeys);
  /* List of the primary keys monitored to determine gathers. */

  if(list_of_keys==NULL)
    sf_error("The required parameter \"pkey\" was not found.");
  /* I wanted to use sf_getstrings, but it seems to want a colon seperated
     list of keys (eg key=offset:ep:fldr:cdp) and I wanted a comma seperated
     list of keys (eg key=offset:ep:fldr:cdp).
  numkeys=sf_getnumpars("pkey");
  if(numkeys==0)
    sf_error("The required parameter \"pkey\" was not found.");
  fprintf(stderr,"alloc list_of_keys numkeys=%d\n",numkeys);
  list_of_keys=(char**)sf_alloc(numkeys,sizeof(char*)); 
  sf_getstrings("pkey",list_of_keys,numkeys);
  */
  /* print the list of keys */
  if(verbose>1){
    fprintf(stderr,"numkeys=%d\n",numkeys);
    for(ikey=0; ikey<numkeys; ikey++){
      fprintf(stderr,"list_of_keys[%d]=%s\n",ikey,list_of_keys[ikey]);
    }
  }
  if(NULL==(skey=sf_getstring("skey")))
    sf_error("the required parameter \"skey\" was not found");
  /* The name of the secondary key created by the program. */
  
//.........这里部分代码省略.........

int main (int argc, char* argv[])
{
    int nh, n1,n2, i1,i2, i, n12, niter, dim, n[SF_MAX_DIM], rect[SF_MAX_DIM];
    int shrt, lng;
    float *trace, *hilb, *num, *den, *phase, mean, c;
    char key[6];
    sf_triangle ts, tl;
    sf_file in, out;

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

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    
    dim = sf_filedims (in,n);
    n1 = n[0];
    n12 = 1;
    for (i=0; i < dim; i++) {
	snprintf(key,6,"rect%d",i+1);
	if (!sf_getint(key,rect+i)) rect[i]=1;
	n12 *= n[i];
    }
    n2 = n12/n1;

    num = sf_floatalloc(n12);
    den = sf_floatalloc(n12);
    phase = sf_floatalloc(n12);

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

    if (!sf_getint("order",&nh)) nh=10;
    /* Hilbert transformer order */
    if (!sf_getfloat("ref",&c)) c=1.;
    /* Hilbert transformer reference (0.5 < ref <= 1) */

    sf_hilbert_init(n1, nh, c);

    if (!sf_getint("short",&shrt)) shrt=1;
    /* short smoothing radius */
    if (!sf_getint("long",&lng)) lng=10;
    /* long smoothing radius */

    ts = sf_triangle_init(shrt,n1,false);
    tl = sf_triangle_init(lng,n1,false);

    mean=0.;
    for (i2=0; i2 < n2; i2++) {
	trace = num+n1*i2;
	hilb = den+n1*i2;

	sf_floatread(trace,n1,in);
	sf_hilbert(trace,hilb);

	for (i1=0; i1 < n1; i1++) {
	    trace[i1] = hypotf(trace[i1],hilb[i1]);
	    hilb[i1] = trace[i1];
	}

	sf_smooth2 (ts, 0, 1, false, trace);
	sf_smooth2 (tl, 0, 1, false, hilb);
	
	for (i1=0; i1 < nh; i1++) {
	    trace[i1] = 0.;
	    hilb[i1] = 0.;
	}
	for (i1=nh; i1 < n1-nh; i1++) {
	    mean += hilb[i1]*hilb[i1];
	}
	for (i1=n1-nh; i1 < n1; i1++) {
	    trace[i1] = 0.;
	    hilb[i1] = 0.;
	}
    }
    mean = sqrtf(n12/mean);

    for (i=0; i < n12; i++) {
	num[i] *= mean;
	den[i] *= mean;
    }

    sf_divn_init(dim, n12, n, rect, niter, true);
    sf_divn (num, den, phase);

    sf_floatwrite(phase,n12,out);

    exit(0);
}

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

    float   A, 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;         /* 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, t5;
    float   timespent;
    float   fx, fz;

    int     isep=1;
    int     ihomo=1;

    char    *tapertype;

    double  vp2, vs2, ep2, de2;

    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;

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

    /* 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);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    sf_map4 mo;
    bool inv, each=true;
    int i, nt, n1, i2, n2, nw;
    float o1, d1, t0, dt, eps;
    sf_complex *ctrace, *ctrace2;
    float *trace, *str, *trace2;
    sf_file in, out, warp;

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

    if (!sf_getbool("inv",&inv)) inv=true;
    /* inversion flag */

    if (inv) {
	if (!sf_histint(in,"n1",&nt)) sf_error("No n1= in input");
	
	if (!sf_getint("n1",&n1)) n1=nt; /* output samples - for inv=y */
	if (!sf_getfloat("d1",&d1) && !sf_histfloat(in,"d1",&d1)) d1=1.;
	/*( d1=1 output sampling - for inv=y )*/
	if (!sf_getfloat("o1",&o1) && !sf_histfloat(in,"o1",&o1)) o1=0.;
	/*( o1=0 output origin - for inv=y )*/ 

	sf_putint(out,"n1",n1);
	sf_putfloat(out,"d1",d1);
	sf_putfloat(out,"o1",o1);
    } else {
	if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
	if (!sf_histfloat(in,"d1",&d1)) d1=1.;
	if (!sf_histfloat(in,"o1",&o1)) o1=0.;

	if (!sf_histint(warp,"n1",&nt)) sf_error("No n1= in warp");
	if (!sf_histfloat(warp,"d1",&dt)) dt=d1;
	if (!sf_histfloat(warp,"o1",&t0)) t0=o1;
	
	sf_putint(out,"n1",nt);
	sf_putfloat(out,"d1",dt);
	sf_putfloat(out,"o1",t0);
    }

    n2 = sf_leftsize(in,1);
    nw = sf_leftsize(warp,1);
    if (1 == nw) {
	each = false;
    } else if (n2 != nw) {
	sf_error("Need %d traces in warp, got %d",n2,nw);
    } 
    
    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* stretch regularization */

    trace = sf_floatalloc(nt);
    str = sf_floatalloc(nt);
    trace2 = sf_floatalloc(n1);

    mo = sf_stretch4_init (n1, o1, d1, nt, eps);

    if (SF_COMPLEX == sf_gettype(in)) {
	ctrace = sf_complexalloc(nt);
	ctrace2 = sf_complexalloc(n1);
    } else {
	ctrace = ctrace2 = NULL;
    }

    for (i2=0; i2 < n2; i2++) {
	if (each || 0==i2) {
	    sf_floatread(str,nt,warp);
	    sf_stretch4_define (mo,str);
	}

	if (inv) {
	    if (SF_COMPLEX == sf_gettype(in)) {
		sf_complexread(ctrace,nt,in);
		for (i=0; i < nt; i++) {
		    trace[i] = crealf(ctrace[i]);
		}
		sf_stretch4_apply (false,mo,trace,trace2);
		for (i=0; i < n1; i++) {
		    ctrace2[i] = sf_cmplx(trace2[i],0.0f);
		}
		for (i=0; i < nt; i++) {
		    trace[i] = cimagf(ctrace[i]);
		}
		sf_stretch4_apply (false,mo,trace,trace2);
		for (i=0; i < n1; i++) {
#ifdef SF_HAS_COMPLEX_H
		    ctrace2[i] += sf_cmplx(0.0f,trace2[i]);
#else
		    ctrace2[i] = sf_cadd(ctrace2[i],sf_cmplx(0.0f,trace2[i]));
#endif
		}
		sf_complexwrite (ctrace2,n1,out);
	    } else {
		sf_floatread(trace,nt,in);
		sf_stretch4_apply (false,mo,trace,trace2);
		sf_floatwrite (trace2,n1,out);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool adj,timer,verb,gmres;
    int nt, nx, nz, nx2, nz2, nzx, nzx2, ntx, pad1, snap, gpz, wfnt, i;
    int m2, n2, nk, nth=1;
    int niter, mem;
    float dt, dx, dz, ox;
    sf_complex *img, *imgout, *dat, **lt1, **rt1, **lt2, **rt2, ***wvfld;
    sf_file data, image, leftf, rightf, leftb, rightb, snaps;
    double time=0.,t0=0.,t1=0.;
    geopar geop;

    sf_init(argc,argv);

    /* essentially doing imaging */
    adj = true;

    if (!sf_getbool("gmres",&gmres)) gmres=false;
    if (gmres) {
      if (!sf_getint("niter",&niter)) niter=10;
      if (!sf_getint("mem",&mem)) mem=20;
    }

    if(! sf_getbool("timer",&timer)) timer=false;
    if (!sf_getbool("verb",&verb)) verb=false;
    if (!sf_getint("snap",&snap)) snap=0;
    /* interval for snapshots */
    if (!sf_getint("pad1",&pad1)) pad1=1;
    /* padding factor on the first axis */
    if(!sf_getint("gpz",&gpz)) gpz=0;
    /* geophone surface */

    /* adj */
    if (!sf_getint("nz",&nz)) sf_error("Need nz=");
    /* depth samples */
    if (!sf_getfloat("dz",&dz)) sf_error("Need dz=");
    /* depth sampling */

    /* for */
    if (!sf_getint("nt",&nt)) sf_error("Need nt=");
    /* time samples */
    if (!sf_getfloat("dt",&dt)) sf_error("Need dt=");
    /* time sampling */

    if (adj) { /* migration */
	data = sf_input("in");
	image = sf_output("out");
	sf_settype(image,SF_COMPLEX);

	if (!sf_histint(data,"n1",&nt)) sf_error("No n1= in input");
	if (!sf_histfloat(data,"d1",&dt)) sf_error("No d1= in input");

	if (!sf_histint(data,"n2",&nx)) sf_error("No n2= in input");
	if (!sf_histfloat(data,"d2",&dx)) sf_error("No d2= in input");
	if (!sf_histfloat(data,"o2",&ox)) ox=0.; 

	sf_putint(image,"n1",nz);
	sf_putfloat(image,"d1",dz);
	sf_putfloat(image,"o1",0.);
	sf_putstring(image,"label1","Depth");

	sf_putint(image,"n2",nx);
	sf_putfloat(image,"d2",dx);
	sf_putfloat(image,"o2",ox);
	sf_putstring(image,"label2","Distance");
    } else { /* modeling */
	image = sf_input("in");
	data = sf_output("out");
	sf_settype(data,SF_COMPLEX);

	if (!sf_histint(image,"n1",&nz)) sf_error("No n1= in input");
	if (!sf_histfloat(image,"d1",&dz)) sf_error("No d1= in input");

	if (!sf_histint(image,"n2",&nx))  sf_error("No n2= in input");
	if (!sf_histfloat(image,"d2",&dx)) sf_error("No d2= in input");
	if (!sf_histfloat(image,"o2",&ox)) ox=0.; 	

	if (!sf_getint("nt",&nt)) sf_error("Need nt=");
	/* time samples */
	if (!sf_getfloat("dt",&dt)) sf_error("Need dt=");
	/* time sampling */

	sf_putint(data,"n1",nt);
	sf_putfloat(data,"d1",dt);
	sf_putfloat(data,"o1",0.);
	sf_putstring(data,"label1","Time");
	sf_putstring(data,"unit1","s");

	sf_putint(data,"n2",nx);
	sf_putfloat(data,"d2",dx);
	sf_putfloat(data,"o2",ox);
	sf_putstring(data,"label2","Distance");
    }

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

    nzx = nz*nx;
    nzx2 = nz2*nx2;
//.........这里部分代码省略.........