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, div;
    int     mm, nvx, nvz, ns;
    int     hnkx, hnkz, nkx, nkz, nxz, nkxz;
    int     hnkx1=1, hnkz1=1, 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 **apvx, **apvz, **apvxx, **apvzz;    /* projection deviation operator of P-wave for a location */

    float ****ex=NULL, ****ez=NULL;                      /* operator for whole model for P-wave*/
    float **exx=NULL, **ezz=NULL;                        /* operator for constant model for P-wave*/

    float **vp0, **vs0, **epsi, **del, **theta;         /* velocity model */
    float **p1, **p2, **p3, **q1, **q2, **q3, **p3c=NULL, **q3c=NULL, **sum=NULL;  /* wavefield array */
    float *kx, *kz, *kkx, *kkz, *kx2, *kz2, **taper;

    clock_t t2, t3, t4, t5;
    float   timespent, fx,fz; 
    char    *tapertype;

    int     isep=1;
    int     ihomo=1;

    double  vp2, vs2, ep2, de2, the;

    sf_file Fo1, Fo2, Fo3, Fo4, Fo5, Fo6, Fo7, Fo8;
    sf_file Fvp0, Fvs0, Feps, Fdel, Fthe;

    sf_axis az, ax;
       
    sf_init(argc,argv);

    /* 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 */
    Fvp0 = sf_input ("in");  /* vp0 using standard input */
    Fvs0 = sf_input ("vs0");  /* vs0 */
    Feps = sf_input ("epsi");  /* epsi */
    Fdel = sf_input ("del");  /* delta */
    Fthe = sf_input ("the");  /* theta */

    /* Read/Write axes */
    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);
    fz=sf_o(az);

    /* 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("dx=%f dz=%f",dx,dz);

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

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

    nxz=nx*nz;
    mm=2*_m+1;
//.........这里部分代码省略.........

int main (int argc, char *argv[])
{
    int ir, nr, n1,n2,n3, m1, m2, m3, n12, nw, nj1, nj2, i3;
    float *u1, *u2, *p, *ani;
    sf_file in, out, dip, aniso;
    bool verb;
    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;
    nr = sf_leftsize(in,2);

    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 (NULL != sf_getstring("aniso")) aniso = sf_input("aniso"); /*aniso field */
	else aniso = NULL ;

    if (!sf_getint("order",&nw)) nw=1;
    /* accuracy */

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

    if (!sf_getint("nj1",&nj1)) nj1=1;
    /* aliasing */

    nj2 = 1;
    n3 = 1;

    for (ir=0; ir < nr; ir++) {
	if (verb) sf_warning("slice %d of %d", ir+1, nr);
	u1  = sf_floatalloc(n12);
	u2  = sf_floatalloc(n12);
	p   = sf_floatalloc(n12);
	ani = sf_floatalloc(n12);
	
	for (i3=0; i3 < n3; i3++) {

            int i12;

	    /* read data */
	    sf_floatread(u1,n12,in);
	    
	    /* read t-x dip */
	    sf_floatread(p,n12,dip);
		if (NULL != aniso)  
			sf_floatread(ani,n12,aniso);		
		else 
			for (i12=0; i12<n12; i12++) {
				ani[i12]=1;
			}
		
  
	    ap = allpass_init (nw,nj1,n1,n2,1,p,ani);
	    
	    /* apply */
	    allpass1(false, ap, u1, u2);
	    
	    /* write t-x destruction */
	    sf_floatwrite(u2,n12,out);
	}
	
	free(u1);
	free(u2);
	free(p);
	
    }
    
    
    exit (0);
}

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

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 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 b1, b2, b3, n1, n2, n3, i, nshot, ndim, is,n123, /* **p, */ j;
    float br1, br2, br3, o1, o2, o3, d1, d2, d3;
    float **s, **t, **v, **vv,**vs,**q,**a,**sg,**bet;
    char *sfile, *file;
    bool plane[3];
    sf_file vvf, bv,  eta, time, shots;

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

    if (SF_FLOAT != sf_gettype(bv)) 
	sf_error("Need float input");
    if (NULL != (file = sf_getstring("vv"))) {
	vvf = sf_input(file);
	free(file);
    } else {
	vvf = NULL;
    }
    if (NULL != (file = sf_getstring("eta"))) {
	eta = sf_input(file);
	free(file);
    } else {
	eta = NULL;
    }
    if(!sf_histint(bv,"n1",&n1)) sf_error("No n1= in input");
    if(!sf_histint(bv,"n2",&n2)) sf_error("No n2= in input");
    if(!sf_histint(bv,"n3",&n3)) n3=1;
    if (n1 <2 || n2<2) sf_error("n1 and n2 muxt be bigger than 1");

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

    if(!sf_histfloat(bv,"o1",&o1)) o1=0.;
    if(!sf_histfloat(bv,"o2",&o2)) o2=0.;
    if(!sf_histfloat(bv,"o3",&o3)) o3=0.;

    /* if y, the input is background time; n, Velocity */

    /* The value of the constant eta */

    if(!sf_getfloat("br1",&br1)) br1=d1;    
    if(!sf_getfloat("br2",&br2)) br2=d2; 
    if(!sf_getfloat("br3",&br3)) br3=d3;
    /* Constant-velocity box around the source (in physical dimensions) */
 
    if(!sf_getbool("plane1",&plane[2])) plane[2]=false;
    if(!sf_getbool("plane2",&plane[1])) plane[1]=false;
    if(!sf_getbool("plane3",&plane[0])) plane[0]=false;
    /* plane-wave source */

    if(!sf_getint("b1",&b1)) b1= plane[2]? n1: (int) (br1/d1+0.5); 
    if(!sf_getint("b2",&b2)) b2= plane[1]? n2: (int) (br2/d2+0.5); 
    if(!sf_getint("b3",&b3)) b3= plane[0]? n3: (int) (br3/d3+0.5); 
    /* Constant-velocity box around the source (in samples) */

    if( b1<1 ) b1=1;  
    if( b2<1 ) b2=1;  
    if( b3<1 ) b3=1;

    sfile = sf_getstring("shotfile");
    /* File with shot locations (n2=number of shots, n1=3) */

    if(NULL != sfile) {
	shots = sf_input("shotfile");

	if (SF_FLOAT != sf_gettype(shots)) 
	    sf_error("Need float shotfile");
	if(!sf_histint(shots,"n2",&nshot)) 
	    sf_error("No n2= in shotfile");
	if(!sf_histint(shots,"n1",&ndim) || ndim != 3) 
	    sf_error("Need n1=3 in shotfile");
  
	s = sf_floatalloc2 (ndim,nshot);
	sf_floatread(s[0],nshot*ndim,shots);
	sf_fileclose(shots);
    
	sf_putint (time,"n4",nshot);
	free (sfile);
    } else {
	nshot = 1;
	ndim = 3;
    
	s = sf_floatalloc2 (ndim,nshot);     

	if(!sf_getfloat("zshot",&s[0][0])  ) s[0][0]=0.; 
	/* Shot location (used if no shotfile) */
	if(!sf_getfloat("yshot",&s[0][1])) s[0][1]=o2 + 0.5*(n2-1)*d2;
	if(!sf_getfloat("xshot",&s[0][2])) s[0][2]=o3 + 0.5*(n3-1)*d3;

	sf_warning("Shooting from zshot=%g yshot=%g xshot=%g",
		   s[0][0],s[0][1],s[0][2]);
    }

    n123 = n1*n2*n3;

    t  = sf_floatalloc2(n1,n2);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
 
    bool hermite_false, hermite_true;
    int n1, n2, npml, pad1, pad2, ns, nw;
    float d1, d2, **v, ds, os, dw, ow;
    sf_complex ****f,  ****obs; 
    sf_file in, out, misfit, source, receiver, record;
    char *order;
    int uts, mts, is, i, j, iw, iter, niter;
    float **recloc;
    float **m_old, **m_new;
    float **d_new, **d_old;
    float **g_old, **g_new; 
    sf_complex ****r_new, ****r_old;
    sf_complex ****Fg; 
    float alpha, beta, gnorm, rnorm; 
    float *datamisfit;

    sf_init(argc, argv);

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

    if (!sf_getint("uts",&uts)) uts=0;
//#ifdef _OPENMP
//    mts = omp_get_max_threads();
//#else
    mts = 1;
//#endif
    uts = (uts < 1)? mts: uts;

    hermite_false=false;
    hermite_true=true;
    /* Hermite operator */
    
    if (!sf_getint("npml",&npml)) npml=20;
    /* PML width */
    
    if (!sf_getint("niter",&niter)) niter=0; 
    /* Number of iterations */

    if (NULL == (order = sf_getstring("order"))) order="c";
    /* discretization scheme (default optimal 9-point) */

    fdprep_order(order);

    /* read input dimension */
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input.");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= 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.");

    v = sf_floatalloc2(n1,n2);
    sf_floatread(v[0],n1*n2,in);
    
	/* PML padding */
	pad1 = n1+2*npml;
	pad2 = n2+2*npml;

    /* read receiver */ 
    if (NULL == sf_getstring("receiver")) sf_error("Need receiver="); 
    receiver = sf_input("receiver"); 
    recloc=sf_floatalloc2(n1,n2);
    sf_floatread(recloc[0],n1*n2,receiver);

    /* read source */
    if (NULL == sf_getstring("source")) sf_error("Need source=");
    source = sf_input("source");

    if (!sf_histint(source,"n3",&ns)) sf_error("No ns=.");
    if (!sf_histfloat(source,"d3",&ds)) ds=d2;
    if (!sf_histfloat(source,"o3",&os)) os=0.;

    /* read observed data */
    if (NULL == sf_getstring("record")) sf_error("Need record=");
    record = sf_input("record");

    if (!sf_histint(record,"n4",&nw)) sf_error("No nw=.");
    if (!sf_histfloat(record,"d4",&dw)) sf_error("No dw=.");
    if (!sf_histfloat(record,"o4",&ow)) sf_error("No ow=."); 

    f = sf_complexalloc4(n1,n2,ns,nw);
    obs = sf_complexalloc4(n1,n2,ns,nw);

    sf_complexread(f[0][0][0],n1*n2*ns*nw,source);
    sf_complexread(obs[0][0][0],n1*n2*ns*nw,record);
   
    /* allocate variables */
    m_old = sf_floatalloc2(n1,n2);
    m_new = sf_floatalloc2(n1,n2);
    d_old = sf_floatalloc2(n1,n2);
    d_new = sf_floatalloc2(n1,n2);
    g_old = sf_floatalloc2(n1,n2);
    g_new = sf_floatalloc2(n1,n2);
    
    r_old = sf_complexalloc4(n1,n2,ns,nw);
    r_new = sf_complexalloc4(n1,n2,ns,nw);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{ 
    int i, j, is, ip, dim, n[SF_MAX_DIM], ii[SF_MAX_DIM];
    int nsp, **k=NULL, **l=NULL, n1, n2, i1, i2, kk, ll;
    char key[7];
    const char *label, *unit;
    float f, *trace, *mag=NULL, **p=NULL, pp;
    sf_file in, spike;

    sf_init (argc,argv);

    if (!sf_stdin()) { /* no input file in stdin */
	in = NULL;
    } else {
	in = sf_input("in");
    }

    spike = sf_output("out");
    
    if (NULL == in) {
	sf_setformat(spike,"native_float");
    } else if (SF_FLOAT != sf_gettype(in)) {
	sf_error("Need float input");
    }
    
    /* dimensions */
    for (i=0; i < SF_MAX_DIM; i++) {
	snprintf(key,3,"n%d",i+1);
	if (!sf_getint(key,n+i) && 
	    (NULL == in || !sf_histint(in,key,n+i))) break;
	/*( n# size of #-th axis )*/  
	sf_putint(spike,key,n[i]);
    }

    if (0==i) sf_error("Need n1=");
    dim=i;
    
    /* basic parameters */
    for (i=0; i < dim; i++) {
	snprintf(key,3,"o%d",i+1);
	if (!sf_getfloat(key,&f) && 
	    (NULL == in || !sf_histfloat(in,key,&f))) f=0.;
	/*( o#=[0,0,...] origin on #-th axis )*/  
	sf_putfloat(spike,key,f);

	snprintf(key,3,"d%d",i+1);
	if (!sf_getfloat(key,&f) &&
	    (NULL == in || !sf_histfloat(in,key,&f))) f = (i==0)? 0.004: 0.1;
	/*( d#=[0.004,0.1,0.1,...] sampling on #-th axis )*/  
	sf_putfloat(spike,key,f);

	snprintf(key,7,"label%d",i+1);
	if (NULL == (label = sf_getstring(key)) &&
	    (NULL == in || NULL == (label = sf_histstring(in,key))))
	    label = (i==0)? "Time":"Distance";
	/*( label#=[Time,Distance,Distance,...] label on #-th axis )*/  
	if (*label != '\0' && (*label != ' ' || *(label+1) != '\0')) 	
	    sf_putstring(spike,key,label);

	snprintf(key,6,"unit%d",i+1);
	if (NULL == (unit = sf_getstring(key)) &&
	    (NULL == in || NULL == (unit = sf_histstring(in,key))))
	    unit = (i==0)? "s":"km";
        /*( unit#=[s,km,km,...] unit on #-th axis )*/  
	if (*unit != '\0' && (*unit != ' ' || *(unit+1) != '\0')) 	
	    sf_putstring(spike,key,unit);
    }
	
    if (NULL != (label = sf_getstring("title")))
	sf_putstring(spike,"title",label);
    /* title for plots */

    if (!sf_getint("nsp",&nsp)) nsp=1;
    /* Number of spikes */

    if (nsp >= 1) { 
	mag = sf_floatalloc (nsp);
	k = sf_intalloc2 (nsp,dim);
	l = sf_intalloc2 (nsp,dim);
	p = sf_floatalloc2 (nsp,dim);
    
	for (i=0; i < dim; i++) {
	    snprintf(key,3,"k%d",i+1);
	    if ( !sf_getints(key,k[i],nsp)) {
		/*( k#=[0,...] spike starting position )*/
		for (is=0; is < nsp; is++) {
		    k[i][is]=-1;
		}
	    } else {
		for (is=0; is < nsp; is++) {
		    if (k[i][is] > n[i]) 
			sf_error("Invalid k%d[%d]=%d > n%d=%d",
				 i+1,is+1,k[i][is],i+1,n[i]);
		    k[i][is]--; /* C notation */
		}
	    }
	    snprintf(key,3,"l%d",i+1);
	    if (!sf_getints(key,l[i],nsp)) {
		/*( l#=[k1,k2,...] spike ending position )*/
		for (is=0; is < nsp; is++) {
//.........这里部分代码省略.........

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[])
{

    /*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[])
{
    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);
	}
    }

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

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

#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 */

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

    e = sf_floatalloc4 (ESC3_NUM, nb, na, nc*mp);

    if (!sf_getstring ("vspl")) sf_error ("Need vspl=");
    /* Spline coefficients for velocity model */
    vspline = sf_input ("vspl");

    if (!sf_getstring ("scgrid")) sf_error ("Need scgrid=");
    /* Grid of supercells of local escape solutions */
    scgrid = sf_input ("scgrid");

    if (sf_getstring ("scdaemon")) {
        /* Daemon for distributed computation */
        scdaemon = sf_input ("scdaemon");
    }

    if (!sf_getint ("morder", &morder)) morder = 1;
    /* Order of interpolation accuracy in the angular domain (1-3) */
#ifdef LINUX
    if (!sf_getint ("inet", &inet)) inet = 1;
    /* Network interface index */
#endif
    if (!sf_getint ("tdel", &tdel)) tdel = 0;
    /* Optional delay time before connecting (seconds) */
   
    /* Slowness components module [(an)isotropic] */
    esc_slow = sf_esc_slowness3_init (vspline, verb);

    /* Make room for escape variables in output */
    if (spdom)
        sf_shiftdimn (spdom, out, 1, 3);

    sf_putint (out, "n1", ESC3_NUM);
    sf_putfloat (out, "o1", 0.0);
    sf_putfloat (out, "d1", 1.0);
    sf_putstring (out, "label1", "Escape variable");

int main(int argc, char* argv[])
{
    bool inv, verb;
    int i1, n1, iw, nt, nw, i2, n2, rect0, niter, n12, n1w;
    int m[SF_MAX_DIM], *rect;
    float t, d1, w, w0, dw, mean=0.0f, alpha;
    float *trace, *kbsc, *mkbsc, *sscc, *mm, *ww;
    sf_complex *outp, *cbsc;
    sf_file in, out, mask, weight, basis;

    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_histfloat(in,"d1",&d1)) d1=1.;

    if (!sf_getbool("inv",&inv)) inv=false;
    /* if y, do inverse transform */

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

    if (NULL != sf_getstring("basis")) {
        basis = sf_output("basis");
        sf_settype(basis,SF_COMPLEX);
    } else {
        basis = NULL;
    }

    if (!inv) {
        if (!sf_getint("nw",&nw)) { /* number of frequencies */
            nt = 2*kiss_fft_next_fast_size((n1+1)/2);
            nw = nt/2+1;
            dw = 1./(nt*d1);
            w0 = 0.;
        } else {
            if (!sf_getfloat("dw",&dw)) {
                /* frequency step */
                nt = 2*kiss_fft_next_fast_size((n1+1)/2);
                dw = 1./(nt*d1);
            }
            if (!sf_getfloat("w0",&w0)) w0=0.;
            /* first frequency */
        }
        n2 = sf_leftsize(in,1);
        sf_shiftdim(in, out, 2);
        sf_putint(out,"n2",nw);
        sf_putfloat(out,"d2",dw);
        sf_putfloat(out,"o2",w0);
        sf_putstring(out,"label2","Frequency");
        sf_putstring(out,"unit2","Hz");
        sf_settype(out,SF_COMPLEX);

        if (!sf_getint("rect",&rect0)) rect0=10;
        /* smoothing radius (in time, samples) */
        if (!sf_getint("niter",&niter)) niter=100;
        /* number of inversion iterations */
        if (!sf_getfloat("alpha",&alpha)) alpha=0.;
        /* frequency adaptivity */

        for(i2=0; i2 < SF_MAX_DIM; i2 ++) {
            m[i2] = 1;
        }
        m[0] = n1;
    } else {
        n2 = sf_leftsize(in,2);
        if (!sf_histint(in,"n2",&nw)) sf_error("No n2= in input");
        if (!sf_histfloat(in,"d2",&dw)) sf_error("No d2= in input");
        if (!sf_histfloat(in,"o2",&w0)) sf_error("No o2= in input");
        sf_unshiftdim(in, out, 2);
        sf_settype(out,SF_FLOAT);
    }

    if (NULL != basis) {
        sf_shiftdim(in, basis, 2);
        sf_putint(basis,"n2",nw);
        sf_putfloat(basis,"d2",dw);
        sf_putfloat(basis,"o2",w0);
        sf_putstring(basis,"label2","Frequency");
        sf_putstring(basis,"unit2","Hz");
    }

    n1w = n1*nw;
    n12 = 2*n1w;
    dw *= 2.*SF_PI;
    w0 *= 2.*SF_PI;

    trace = sf_floatalloc(n1);
    kbsc    = sf_floatalloc(n12);
    outp = sf_complexalloc(n1w);
    cbsc = sf_complexalloc(n1w);

    rect = sf_intalloc(2*nw);
    for (iw=0; iw < nw; iw++) {
        rect[iw+nw] = rect[iw] = SF_MAX(1, (int) rect0/(1.0+alpha*iw/nw));
    }

    if (!inv) {
        sscc = sf_floatalloc(n12);
//.........这里部分代码省略.........

int main (int argc,char* argv[]) 
{
    int n[3], ns, ns0, nt, is0, s0;
    float o[3], os, os0, d[3], ds, ds0, **t, **tds, *tempt, ss;
    char *type;
    sf_file in, out, deriv, pattern;

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

    /* read input dimensions */
    if(!sf_histint(in,"n1",n  )) sf_error("No n1= in input");
    if(!sf_histint(in,"n2",n+1)) sf_error("No n2= in input");
    if(!sf_histint(in,"n3",n+2)) n[2]=1;
    if(!sf_histint(in,"n4",&ns)) sf_error("No ns= in input");

    if (ns <= 1) sf_error("Provide at least two shots");

    if(!sf_histfloat(in,"d1",d  )) sf_error("No d1= in input");
    if(!sf_histfloat(in,"d2",d+1)) sf_error("No d2= in input");
    if(!sf_histfloat(in,"d3",d+2)) d[2]=d[1];
    if(!sf_histfloat(in,"d4",&ds)) ds=d[1];

    if(!sf_histfloat(in,"o1",o  )) o[0]=0.;
    if(!sf_histfloat(in,"o2",o+1)) o[1]=0.;
    if(!sf_histfloat(in,"o3",o+2)) o[2]=o[1];
    if(!sf_histfloat(in,"o4",&os)) os=o[1];

    /* read traveltime file */
    nt = n[0]*n[1]*n[2];
    
    t = sf_floatalloc2(nt,ns);
    sf_floatread(t[0],nt*ns,in);

    if (NULL == (type = sf_getstring("type"))) type="hermit";
    /* type of interpolation (default Hermit) */

    if (type[0] != 'h') {
	/* linear interpolation */
	deriv = NULL;
	tds = NULL;
    } else {
	/* read derivative file */
	if (NULL == sf_getstring("deriv"))
	    sf_error("Need derivative deriv=");
	deriv = sf_input("deriv");

	tds = sf_floatalloc2(nt,ns);
	sf_floatread(tds[0],nt*ns,deriv);
	sf_fileclose(deriv);
    }
    
    if (NULL != sf_getstring("pattern")) {
	pattern = sf_input("pattern");
    } else {
	pattern = NULL;
    }
    
    if (!sf_getint("ns",&ns0) && 
	(NULL==pattern ||
	 !sf_histint(pattern,"n4",&ns0))) sf_error("Need ns=");
    /* Output source size */
    if (!sf_getfloat("ds",&ds0) && 
	(NULL==pattern ||
	 !sf_histfloat(pattern,"d4",&ds0))) sf_error("Need ds=");
    /* Output source sampling */
    if (!sf_getfloat("os",&os0) &&
	(NULL==pattern ||
	 !sf_histfloat(pattern,"o4",&os0))) sf_error("Need os=");
    /* Output source origin */
    
    sf_putint(out,"n4",ns0);
    sf_putfloat(out,"d4",ds0);
    sf_putfloat(out,"o4",os0);

    /* allocate temporaty memory */
    tempt = sf_floatalloc(nt);

    /* initialization */
    tinterp_init(nt,ds);

    /* loop over sources */
    for (is0=0; is0 < ns0; is0++) {
	s0 = (os0+is0*ds0-os)/ds;
	ss = os0+is0*ds0-os-s0*ds;

	if (s0 < 0) {
	    s0 = 0; ss = 0.;
	}
	if (s0 >= ns-1) {
	    s0 = ns-2; ss = ds;
	}

	/* do interpolation */
	switch (type[0]) {
	    case 'l': /* linear */
		tinterp_linear(tempt,ss,t[s0],t[s0+1]);
		break;

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

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