void pt2dwrite1(sf_file F, 
		pt2d   *v, 
		size_t  n1, 
		int     k)
/*< output point2d 1-D vector >*/
{
    int i1;
    float **w;
    w=sf_floatalloc2(k,n1);

    for( i1=0; i1<(int)n1; i1++) {
	;        w[i1][0] = v[i1].x; 
	;        w[i1][1] = v[i1].z; 
	if(k==3) w[i1][2] = v[i1].v;
    }
    sf_floatwrite(w[0],k*n1,F);

    free( *w); free(w);
}

int main (int argc, char* argv[])
{
    bool scale;
    int n1,n2, i1,i2, n;
    float d1, *dat, *der;
    sf_file in=NULL, out=NULL;

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

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

    dat = sf_floatalloc(n1);
    der = sf_floatalloc(n1);

    if (!sf_getint("order",&n)) n=6;
    /* filter order */

    if (!sf_getbool("scale",&scale) || 
	!sf_histfloat(in,"d1",&d1)) scale=false;
    /*(scale=n if scale by 1/dx )*/

    sf_deriv_init(n1, n, 0.);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(dat,n1,in);
	sf_deriv(dat,der);

	if (scale) {
	    for (i1=0; i1 < n1; i1++) {
		der[i1] /= d1;
	    }
	}

	sf_floatwrite(der,n1,out);
    }


    exit(0);
}

int main(int argc, char* argv[])
{
    bool adj, inv;
    int nn,n1,n2,i1,i2;
    float rho, *pp=NULL;
    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");
    n2 = sf_leftsize(in,1);

    if (!sf_getbool("adj",&adj)) adj=false;
    /* If y, apply adjoint */
    if (!sf_getbool("inv",&inv)) inv=false;
    /* If y, do differentiation instead of integration */
    if (!sf_getfloat("rho",&rho)) rho = 1.-1./n1;
    /* Leaky integration constant */

    if (inv) {
	sf_warning("%s half-order differentiation",adj? "anticausal":"causal");
    } else {
	sf_warning("%s half-order integration",adj? "anticausal":"causal");
    }

    nn = 2*kiss_fft_next_fast_size((n1+1)/2);
    pp = sf_floatalloc(nn);

    sf_halfint_init (inv, nn, rho);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread (pp,n1,in);
	for (i1=n1; i1 < nn; i1++) {
	    pp[i1]=0.;
	}
	sf_halfint (adj, pp);
	sf_floatwrite (pp,n1,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int nt, ns, nx, ix, is, it;
    float *trace=NULL, *picks=NULL, **semb=NULL, s0, ds, s;
    sf_file in=NULL, out=NULL, pick=NULL;

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

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

    if (!sf_histint(in,"n2",&ns)) sf_error("No n2= in input");
    if (!sf_histfloat(in,"o2",&s0)) sf_error("No o2= in input");
    if (!sf_histfloat(in,"d2",&ds)) sf_error("No d2= in input");

    nx = sf_unshiftdim(in,out,2);

    semb = sf_floatalloc2(nt,ns);
    picks = sf_floatalloc (nt);
    trace = sf_floatalloc (nt);

    for (ix=0; ix < nx; ix++) {
	sf_floatread (semb[0],nt*ns,in);
	sf_floatread (picks,nt,pick);

	for (it=0; it < nt; it++) {
	    s = (picks[it] - s0)/ds;
	    is = floorf(s); s -= is;
	    if (is >= 0 && is < ns-1) {
		trace[it] = s * semb[is+1][it] + (1.-s) * semb[is][it];
	    } else {
		trace[it] = 0.;
	    }
	}

	sf_floatwrite (trace,nt,out);
    }

    exit (0);
}

int main(int argc, char* argv[])
{
    int  i1, i2, n1,n2,n3,n12;
    float  s=0, *pp1, *pp2;
    sf_file inp1, inp2,  dif;

    sf_init(argc,argv);
    inp1 = sf_input("in");
    inp2 = sf_input("match");
    dif = sf_output("out");
    
    sf_putfloat(dif,"o2",0);
    sf_putfloat(dif,"o1",0);
    sf_putint(dif,"n1",1);  
    sf_putint(dif,"n2",1);    
	
    if (!sf_histint(inp1,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(inp1,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(inp1,"n3",&n3)) n3=1;

    n2 = sf_leftsize(inp1,1); /* left dimensions after the first one */
    n12 = n1*n2;
    
    if (n3!=1) {sf_putint(dif,"n3",1); sf_putint(dif,"d3",1); }
    
    pp1 = sf_floatalloc(n12);
    pp2 = sf_floatalloc(n12);
    sf_floatread(pp1,n12,inp1);
    sf_floatread(pp2,n12,inp2);

    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    s = s + pow((pp1[i2*n1+i1]-pp2[i2*n1+i1]),2);
	}
    }
    
    sf_warning("The difference is %f", s );
    sf_floatwrite(&s,1,dif);
    exit(0);
}

int main(int argc, char* argv[])
{
    int nx, nz, ix, iz;
    float* trace;
    float dx, dz, x[2];
    sf_file mod;

    sf_init (argc,argv);
    mod = sf_output("out");

    if (!sf_getint("nx",&nx)) nx=400; 
    /* horizontal dimension */
    if (!sf_getint("nz",&nz)) nz=800; 
    /* vertical dimension */

    dx = 1./(nx-1);
    dz = 2./(nz-1);

    sf_putint   (mod,"n1",nz); 
    sf_putfloat (mod,"d1",dz); 
    sf_putfloat (mod,"o1",0.);
    sf_putint   (mod,"n2",nx); 
    sf_putfloat (mod,"d2",dx); 
    sf_putfloat (mod,"o2",0.);
    sf_setformat (mod,"native_float");

    trace = sf_floatalloc(nz);

    for (ix = 0; ix < nx; ix++) {
	x[1] = ix*dx;
	for (iz = 0; iz < nz; iz++) {
	    x[0] = iz*dz;
	    trace[iz] = 1./sqrtf(symes_vel(NULL,x));
	}
	sf_floatwrite(trace,nz,mod);
    }

    exit (0);
}

void pt2dwrite2(sf_file F, 
		pt2d  **v, 
		size_t  n1, 
		size_t  n2,
		int     k)
/*< output point2d 2-D vector >*/
{
    int i1,i2;
    float ***w;
    w=sf_floatalloc3(k,n1,n2);

    for( i2=0; i2<(int)n2; i2++) {
	for( i1=0; i1<(int)n1; i1++) {
	    ;        w[i2][i1][0] = v[i2][i1].x; 
	    ;        w[i2][i1][1] = v[i2][i1].z; 
	    if(k==3) w[i2][i1][2] = v[i2][i1].v;
	}
    }
    sf_floatwrite(w[0][0],k*n1*n2,F);
    
    free(**w); free(*w); free(w);
}

int main(int argc, char* argv[])
{
    int n1, n2, i2;
    float *pp, *qq;
    bool square, adj;
    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");
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    n2 = sf_leftsize(in,1);

    if (!sf_getbool("square",&square)) square=false;
    /* if y, use gradient squared */

    if (!sf_getbool("adj",&adj)) adj=false;
    /* adjoint flag */

    pp = sf_floatalloc(n1);
    qq = sf_floatalloc(n1);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(pp,n1,in);
	if (square) {
	    sf_grad2 (n1,pp,qq);
	} else if (adj) {
	    sf_igrad1_lop (true,false,n1,n1,qq,pp);
	} else {
	    sf_igrad1_lop (false,false,n1,n1,pp,qq);
	}
	sf_floatwrite(qq,n1,out);
    }

    exit(0);
}

int main(int argc,char *argv[])
{
    int i,j,n1,n2,n3;
    float o1,d1,o2,d2,max;
    float *orig,*output;
    sf_file in,out;
    
    sf_init (argc, argv); 
    in = sf_input("in");
    out = sf_output("out");
    
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&o1)) sf_error("No o1= in input");
    if (!sf_histint(in,"n2",&n2)) n2=1;
    if (!sf_histfloat(in,"d2",&d2)) d2=1;
    if (!sf_histfloat(in,"o2",&o2)) o2=1;
    n3=sf_leftsize(in,2);
    
    orig = sf_floatalloc(n1*n2);
    output = sf_floatalloc(n1*n2);
    for(j = 0; j < n3; j++){
	
	sf_floatread(orig,n1*n2,in);
	max=0;
	for(i=0;i<n1*n2;i++){
	    if(orig[i]>max) max=orig[i];
	}
	
	for(i=0;i<n1*n2;i++){
	    output[i]=max-orig[i];
	}
	
	sf_floatwrite(output,n1*n2,out);
    }
    exit(0);
    
}

int main(int argc, char* argv[])
{ 
    int dim, m[SF_MAX_DIM], nd, i, niter, rect[SF_MAX_DIM];
    float **inp, *rat;
    char key[6];
    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,m);

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

    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */
    
    inp = sf_floatalloc2 (m[0],nd/m[0]);
    rat = sf_floatalloc (nd);

    envcorr_init(dim,m,rect,niter);
    
    sf_floatread(inp[0],nd,in);

    envcorr(inp,rat);

    sf_floatwrite(rat,nd,out);

    exit (0);
}

int main(int argc, char** argv)
{
    int n1, n2, i1, i2;
    float* vint=NULL;
    float o1, d1, v0, alpha, t;
    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\n");
    if (!sf_histint (in,"n2",&n2)) sf_error ("No n2= in input\n");

    if (!sf_histfloat (in,"o1",&o1)) o1=0.;
    if (!sf_histfloat (in,"d1",&d1)) sf_error ("No d1= in input\n");

    if (!sf_getfloat ("v0",&v0)) v0 = 1.5;
    /* initial velocity */
    if (!sf_getfloat ("alpha",&alpha)) alpha = 0.5;
    /* velocity gradient */

    vint = sf_floatalloc(n1);
    for (i1=0; i1 < n1; i1++) {
	t = alpha*(o1+i1*d1);
	if (t > 0.) {
	    vint[i1] = v0 * sqrtf((expf(t)-1.)/t);
	} else {
	    vint[i1] = v0;
	}
    }
    for (i2=0; i2 < n2; i2++) {
	sf_floatwrite(vint,n1,out);
    }

    exit (0);
}

int main (int argc, char *argv[])
{
    int n1, n2, n3, i2, i3, order;
    float *trace;
    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");
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);
    
    if (!sf_getint("order",&order)) order=4;
    /* Accuracy order */

    trace = sf_floatalloc (n1);

    pick0_init (n1, n2, order);

    for (i3=0; i3 < n3; i3++) {
	for (i2=0; i2 < n2; i2++) {
	    sf_floatread(trace,n1,in);
	    pick0_set (i2, trace);
	}

	for (i2=0; i2 < n2; i2++) {
	    pick0_delta (i2, trace);
	    sf_floatwrite(trace,n1,out);
	}
    }
    

    exit (0);
}

int main(int argc, char *argv[])
{
    int i, m, n, k, n1a, n1b, n2a, n2b, n3a, n3b; /*n1 is trace length, n2 is the number of traces, n3 is the number of 3th axis*/
    sf_file in, mat, out;
    float *a, *b, *c;
    char *type;
    sf_init(argc,argv);

    in =  sf_input("in");    
    /*  in is the input matrix A, with m rows and n columns */
    mat= sf_input("mat");
    /* mat is the input matrix B, with n rows and k columns */
    out = sf_output("out");
    /* out is the output matrix C, with m rows and k columns */

    if (NULL == (type=sf_getstring("type"))) type="mul";
    /* [mul, add, sub, dotmul] operation type, the default is mul  */


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

    if (!sf_histint(mat,"n1",&n1b)) sf_error("No n1= in input");
    if (!sf_histint(mat,"n2",&n2b)) sf_error("No n2= in input");
    n3b = sf_leftsize(mat,2);
    
    if (type[0]=='m')
	{
	if (n1a!=n2b) sf_error("n1 of A must be equal to n2 of B for multiplication");
	m=n2a;n=n1a;k=n1b;
	}

    if (type[0]!='m')
	{
	if (n1a!=n1b || n2a!=n2b) sf_error(" n1 of A must be equal to n1 of B and n2 of A must be equal to n2 of B");
	m=n2a;n=n1a;k=n;	
	}
    if (n3a != n3b) sf_error("n3 must be equal");
    
    sf_putint(out, "n1", k);   
    sf_putint(out, "n2", m);
    sf_putint(out, "n3", n3a);

    a = sf_floatalloc(m*n);
    b = sf_floatalloc(n*k);
    c = sf_floatalloc(m*k);

    for(i=0;i<n3a;i++)  {
	sf_floatread(a, m*n, in);

	if(type[0]!='m')
	sf_floatread(b, m*n, mat);
	else
	sf_floatread(b,n*k,mat);

	switch(type[0])
	{
		case 'm': mul(a, b, m, n, k, c); break;
		case 'a': add(a, b, m, n, c); break;
		case 's': sub(a, b, m, n, c); break;
		case 'd': dotmul(a, b, m, n, c);
	}
	sf_floatwrite(c, m*k, out);
    }

    exit(0);
}

int main (int argc, char *argv[])
{
    bool verb;
    int n1,n2,n3, n12, n23, ref2, ref3, i2,i3,i1, ud, lr, order;
    float eps, ***dat, ***p, ***q, **p2, **q2, *trace;
    sf_file dip, out, seed, cost;

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

    if (!sf_histint(dip,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(dip,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(dip,"n3",&n3)) sf_error("No n3= in input");
    n23 = n2*n3;
    n12 = n1*n23;

    if (!sf_getbool("verb",&verb)) verb=false;
    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* regularization */

    if (!sf_getint("ref2",&ref2)) ref2=0;
    if (!sf_getint("ref3",&ref3)) ref3=0;
    /* reference trace */

    sf_putint(out,"n4",1);

    p = sf_floatalloc3(n1,n2,n3);
    q = sf_floatalloc3(n1,n2,n3);
    dat = sf_floatalloc3(n1,n2,n3);

    sf_floatread(p[0][0],n12,dip);
    sf_floatread(q[0][0],n12,dip);

    p2 = sf_floatalloc2(n2,n3);
    q2 = sf_floatalloc2(n2,n3);

    sf_floatread(p2[0],n23,cost);
    sf_floatread(q2[0],n23,cost);

    dijkstra_init(n2,n3,p2,q2);
    dijkstra_source(ref2,ref3);
    sf_floatread(dat[ref3][ref2],n1,seed);
    if (verb) sf_warning("%d %d",ref2,ref3);

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

    predict_init(n1,n2, eps*eps, order, 1, false);
    
    while (dijskstra_step(&i2,&i3,&ud,&lr)) {
	if (verb) sf_warning("%d %d",i2,i3);

	trace = dat[i3][i2];
	for (i1=0; i1 < n1; i1++) {
	    trace[i1] = dat[i3-lr][i2-ud][i1];
	} 
	if (ud > 0) {
	    predict_step(false,true,trace,p[i3][i2-ud]);
	} else if (ud < 0) {
	    predict_step(false,false,trace,p[i3][i2]);
	}

	if (lr > 0) {
	    predict_step(false,true,trace,q[i3-lr][i2]);
	} else if (lr < 0) {
	    predict_step(false,false,trace,q[i3][i2]);
	}
    }
	
    sf_floatwrite(dat[0][0],n12,out);

    exit (0);
}

//.........这里部分代码省略.........
    sf_putfloat(out,"o2",o2); 
    sf_putfloat(out,"o3",0.0); 

    nxb = nx + nxl + nxr;
    nzb = nz + nbt + nbb;
    
	sf_putint(out,"n1",nzb);
    sf_putint(out,"n2",nxb);
    sf_putint(out,"n3",nsnap);


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

    old    =  sf_floatalloc2(nzb,nxb);
    cur    =  sf_floatalloc2(nzb,nxb);
    aa     =  sf_floatalloc3(nzb,nxb,3);
    

    bd2_init(nx,nz,nxl,nxr,nbt,nbb,cxl,cxr,czt,czb);

    /*input & extend velocity model*/
    v = sf_floatalloc2(nzb,nxb);
    vtmp = sf_floatalloc2(nz,nx);
    sf_floatread(vtmp[0],nx*nz,vel);
	v = extmodel(vtmp, nz, nx, nbt);

    sigma = sf_floatalloc2(nzb,nxb);
    sf_floatread(vtmp[0],nx*nz,fsigma);
	sigma = extmodel(vtmp, nz, nx, nbt);

    delta = sf_floatalloc2(nzb,nxb);
    sf_floatread(vtmp[0],nx*nz,fdelta);
	delta = extmodel(vtmp, nz, nx, nbt);
    
    seta = sf_floatalloc2(nzb,nxb);
    sf_floatread(vtmp[0],nx*nz,fseta);
	seta = extmodel(vtmp, nz, nx, nbt);

    v0 =0.0;
    for (ix=0; ix < nxb; ix++) {
        for (iz=0; iz < nzb; iz++) {
            v0 += v[ix][iz]*v[ix][iz];
         }
    }

    v0 = sqrtf(v0/(nxb*nzb));

	fprintf(stderr, "v0=%f\n\n", v0);


    for (ix=0; ix < nxb; ix++) {
        for (iz=0; iz < nzb; iz++) {
            cur[ix][iz] = 0.0;
            old[ix][iz] = 0.0; 
        }
    }

    /* propagation in time */
    pamstep2_init(nzb,nxb,dz,dx,opt);

    for (it=0; it < nt; it++) {
		fprintf(stderr, "\b\b\b\b\b%d", it);

        pamstep2(old, cur, nzb, nxb, dz, dx, v0, v, sigma, delta, seta, dt); 
        old[isx+nxl][isz+nbt] += wav[it];
     
		bd2_decay(old); 
        bd2_decay(cur); 
        tmp = old;
        old = cur;
        cur = tmp;

		if (it%nsnap==0)
		sf_floatwrite(cur[0], nxb*nzb, out);

    }

    pamstep2_close();
    bd2_close();

    free(**aa);
    free(*aa);
    free(aa);
    free(*v);     
    free(*sigma);     
    free(*delta);     
    free(*seta);     
	free(*vtmp);
    free(*cur);     
    free(*old);     
    free(v);     
    free(sigma);     
    free(delta);     
    free(seta);     
	free(vtmp);
    free(cur);     
    free(old);     
    exit(0); 
}           

//.........这里部分代码省略.........
    aa=sf_maxa(1,0,1); 
    sf_setlabel(aa,""); 
    sf_setunit (aa,""); 

    if(verb) {
	sf_raxa(az);
	sf_raxa(ax);
	sf_raxa(at);
    }

    /* write axes */
    sf_oaxa(Fi,az,1);
    sf_oaxa(Fi,ax,2);
    sf_oaxa(Fi,aa,3);
    
    /*------------------------------------------------------------*/
    /* allocate work arrays */
    ii = sf_floatalloc2(nz,nx); 
    us = sf_floatalloc2(nz,nx);
    ur = sf_floatalloc2(nz,nx);
    for    (ix=0; ix<nx; ix++) {
	for(iz=0; iz<nz; iz++) {
	    ii[ix][iz]=0.;
	}
    }

    /*------------------------------------------------------------*/
    if(isreversed) { /* receiver wavefield is reversed */

	if(verb) fprintf(stderr,"nt\n");
	for (it=0; it<nt; it++) {
	    if(verb) fprintf(stderr,"\b\b\b\b\b\b\b\b\b\b%04d",it);
	    
	    sf_floatread(us[0],nz*nx,Fs);
	    sf_floatread(ur[0],nz*nx,Fr);
	    
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic)				\
    private(ix,iz)							\
    shared (ii,us,ur,nx,nz)
#endif
	    for    (ix=0; ix<nx; ix++) {
		for(iz=0; iz<nz; iz++) {
		    ii[ix][iz] += us[ix][iz]*ur[ix][iz];
		}
	    }
	    
	} /* it */
	if(verb) fprintf(stderr,"\n");
	
    } else { /* receiver wavefield is NOT reversed */

	if(verb) fprintf(stderr,"nt\n");
	for (it=0; it<nt; it++) {
	    if(verb) fprintf(stderr,"\b\b\b\b\b%d",(nt-it-1));
	    
	    sf_floatread(us[0],nz*nx,Fs);
	    iseek=(off_t)(nt-1-it)*nz*nx*sizeof(float);
	    sf_seek(Fr,iseek,SEEK_SET);
	    sf_floatread(ur[0],nz*nx,Fr);
	    
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic)				\
    private(ix,iz)							\
    shared (ii,us,ur,nx,nz)
#endif
	    for    (ix=0; ix<nx; ix++) {
		for(iz=0; iz<nz; iz++) {
		    ii[ix][iz] += us[ix][iz]*ur[ix][iz];
		}
	    }

	} /* it */
	if(verb) fprintf(stderr,"\n");

    } /* end "is reversed" */
    /*------------------------------------------------------------*/
       
    /*------------------------------------------------------------*/
    /* scale image */
    scale = 1./nt;
    for    (ix=0; ix<nx; ix++) {
	for(iz=0; iz<nz; iz++) {
	    ii[ix][iz] *=scale;
	}
    }

    /*------------------------------------------------------------*/
    /* write image */
    sf_floatwrite(ii[0],nx*nz,Fi);
    
    /*------------------------------------------------------------*/
    /* deallocate arrays */
    free(*ii); free(ii);
    free(*us); free(us);
    free(*ur); free(ur);
    /*------------------------------------------------------------*/
    
    exit (0);
}

//.........这里部分代码省略.........
	*(vel[0][0]+ix) *= *(vel[0][0]+ix)*dt2;
    if (fsrf && !hybrid) {
	for (iy=0; iy<nypad; iy++)
	    for (ix=0; ix<nxpad; ix++)
		memset(vel[iy][ix],0,sizeof(float)*fdm->nb);
    }

    for (it=0; it<nt; it++) {
	if (verb)  sf_warning("it=%d;",it+1);
#if defined _OPENMP && _DEBUG
	tic=omp_get_wtime();
#endif
    
	step_forward(u0,u1,vel,rho,fdcoef_d2,fdcoef_d1,nop,nzpad,nxpad,nypad);
    
	if (adj) { /* backward inject source wavelet */
	    if (expl) {
		sf_seek(file_wav,(off_t)(nt-it-1)*sizeof(float),SEEK_SET);
		sf_floatread(ws,1,file_wav);
		ws[0] *= dt2;
		if (sinc) sinc3d_inject1(u0,ws[0],cssinc);
		else      lint3d_inject1(u0,ws[0],cslint);
	    } else { 
		sf_seek(file_wav,(off_t)(nt-it-1)*ns*sizeof(float),SEEK_SET);
		sf_floatread(ws,ns,file_wav);
		for (is=0; is<ns; is++) ws[is] *= dt2;
		if (sinc) sinc3d_inject(u0,ws,cssinc);
		else      lint3d_inject(u0,ws,cslint);
	    }
	} else { /* forward inject source wavelet */
	    if (expl) {
		sf_floatread(ws,1,file_wav);
		ws[0] *= dt2;
		if (sinc) sinc3d_inject1(u0,ws[0],cssinc);
		else      lint3d_inject1(u0,ws[0],cslint);
	    } else {
		sf_floatread(ws,ns,file_wav);
		for (is=0; is<ns; is++) ws[is] *= dt2;
		if (sinc) sinc3d_inject(u0,ws,cssinc);
		else      lint3d_inject(u0,ws,cslint);
	    }
	}

	/* apply abc */
	if (dabc) {
	    if (hybrid) apply_abc(u0,u1,nz,nx,ny,nbd,abc,nop,damp);
	    else {
		abcone3d_apply(u0,u1,nop,abc,fdm);
		sponge3d_apply(u0,spo,fdm);
		sponge3d_apply(u1,spo,fdm);
	    }
	}

	/* loop over pointers */
	ptr_tmp = u0;  u0 = u1;  u1 = ptr_tmp;

	/* extract snapshot */
	if (snap && it%jsnap==0) {
	    int fy = (floor)((sf_o(acy)-fdm->oypad)/fdm->dy);
	    int jy = floor(sf_d(acy)/fdm->dy);
	    float **ptr_slice;
	    for (iy=0; iy<sf_n(acy); iy++) {
		ptr_slice = u0[fy+iy*jy];
		cut3d_slice(ptr_slice,oslice,fdm,acz,acx);
		sf_floatwrite(oslice[0],sf_n(acz)*sf_n(acx),file_wfl);
	    }
	}

	/* extract receiver data */
	if (sinc) sinc3d_extract(u0,u_dat,crsinc);
	else      lint3d_extract(u0,u_dat,crlint);

	sf_floatwrite(u_dat,nr,file_dat);

#if defined _OPENMP && _DEBUG
	toc=omp_get_wtime(); 
	fprintf(stderr,"%5.2gs",(float)(toc-tic));
#endif
    }
#ifdef _OPENMP
    wall_clock_time_e = omp_get_wtime();
#else
    wall_clock_time_e = (double) clock() / CLOCKS_PER_SEC;
#endif
    if (verb)
	fprintf(stderr,"\nElapsed time: %lf s\n",wall_clock_time_e-wall_clock_time_s);

    free(**u0); free(*u0); free(u0);
    free(**u1); free(*u1); free(u1);
    free(**vel); free(*vel); free(vel);
    free(u_dat);
    free(ws);
    free(fdcoef_d2); free(fdcoef_d1);
    if (snap) { free(*oslice); free(oslice); }
    if(!cden) { free(**rho); free(*rho); free(rho); }
    if (hybrid) free(damp);
    free(src3d); free(rec3d);

    return 0;
}

//.........这里部分代码省略.........
    /* END DECLARATIONS */
    
    sf_init(argc,argv);
    
    /* read inputs from command line */
    getinputs(true,&wi);
    
    /* compute number of shots */
    nsrc = (wi.isxend-wi.isxbeg)/(wi.iskip); nsrc++;
    
    /* compute number of spatial grid points */
    nxz=wi.nx * wi.nz;
    
    /* compute number of traces, samples in each record */
    ntr=wi.igxend-wi.igxbeg+1;
    nsam=ntr*wi.nt;
    
    /* allocate, initialize p0, p1, v, traces */
    p0=sf_floatalloc(nxz);
    p1=sf_floatalloc(nxz);
    v =sf_floatalloc(nxz);
    tr=sf_floatalloc(nsam);
    
    /* read velocity */
    sf_floatread(v,nxz,wi.vfile);
    
    /* CFL, sanity checks */
    vmax=fgetmax(v,nxz);
    vmin=fgetmin(v,nxz);
    if (vmax*wi.dt>CFL*fmaxf(wi.dx,wi.dz)) {
	sf_warning("CFL criterion violated");
	sf_warning("vmax=%e dx=%e dz=%e dt=%e\n",vmax,wi.dx,wi.dz,wi.dt);
	sf_error("max permitted dt=%e\n",CFL*fmaxf(wi.dx,wi.dz)/vmax);
    }
    if (vmin<=0.0) 
	sf_error("min velocity nonpositive");
    
    /* only square of velocity array needed from here on */
    fsquare(v,nxz);
    
    /* precalculate some coefficients */
    rz=wi.dt*wi.dt/(wi.dz*wi.dz);
    rx=wi.dt*wi.dt/(wi.dx*wi.dx);
    s =2.0*(rz+rx);
/*    two=2.0;
      nz=wi.nz; */
    
    /* shot loop */
    isrc=0;
    isx=wi.isxbeg;
    while (isx <= wi.isxend) {
	
	/* initialize pressure fields, traces */
	fzeros(p0,nxz);
	fzeros(p1,nxz);
	fzeros(tr,nsam);
	
	/* initialize movie frame counter */
	imf=0;
	
	/* time loop */
	for (it=0;it<wi.nt;it++) {
	    
	    /* construct next time step, overwrite on p0 */
	    
	    step_forward(p0,p1,v,wi.nz,wi.nx,rz,rx,s);
	    
	    /* tack on source */
	    p0[wi.isz+isx*wi.nz]+=fgetrick(it*wi.dt,wi.freq);
	    
	    /* swap pointers */
	    tmp=p0;
	    p0=p1;
	    p1=tmp;
	    
	    /* store trace samples if necessary */
	    if (NULL != wi.tfile) 
		for (ix=0;ix<ntr;ix++) 
		    tr[ix*wi.nt+it]=p1[(wi.igxbeg+ix)*wi.nz+wi.igz];
	    
	    /* write movie snap to file if necessary */
	    if (NULL != wi.mfile && wi.nm && !(it%wi.nm)) {
		sf_floatwrite(p1,nxz,wi.mfile);
		imf++;
	    }
	    
	    /* next t */
	}

	/* write traces to file if necessary */
	if (NULL != wi.tfile) 
	    sf_floatwrite(tr,nsam,wi.tfile);
	
	isx += wi.iskip;
	isrc++;
    } 


    exit(0);
}

int main(int argc, char* argv[])
{
    int p[4][2], i, im, id, status;
    unsigned long mseed, dseed;
    off_t nm, nd, msiz, dsiz;
    size_t nbuf, mbuf, dbuf;
    float *buf;
    double dp;
    pid_t pid[6]={1,1,1,1,1,1};
    sf_file mod=NULL;
    sf_file dat=NULL;
    sf_file pip=NULL;

    sf_init(argc,argv);

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

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

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

    nbuf = BUFSIZ/sizeof(float);
    buf = sf_floatalloc(nbuf);

    mseed = (unsigned long) time(NULL);
    init_genrand(mseed);
    mseed = genrand_int32();
    dseed = genrand_int32();

    for (i=0; i < argc-1; i++) {
	argv[i]=argv[i+1];
    }
    argv[argc-1] = sf_charalloc(6);
    snprintf(argv[argc-1],6,"adj=X");

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

    for (i=0; i < 6; i++) { /* fork six children */
	if ((pid[i] = fork()) < 0) sf_error("fork error:");
	if (0 == pid[i]) break;
    }

    if (0 == pid[0]) {	
	/* makes random model and writes it to p[0] */

	close(p[0][0]);
	close(STDOUT_FILENO);
	DUP(p[0][1]);

	pip = sf_output("out");
	sf_fileflush(pip,mod);

	init_genrand(mseed);
	for (msiz=nm, mbuf=nbuf; msiz > 0; msiz -= mbuf) {
	    if (msiz < mbuf) mbuf=msiz;

	    sf_random(mbuf,buf);

	    sf_floatwrite(buf,mbuf,pip);
	}
    } 

    if (0 == pid[1]) {
	/* reads from p[0], runs the program, and writes to p[1] */

	close(p[0][1]);
	close(STDIN_FILENO);
	DUP(p[0][0]);

	close(p[1][0]);
	close(STDOUT_FILENO);
	DUP(p[1][1]);

	argv[argc-1][4]='0';
	execvp(argv[0],argv);

	_exit(1);
    }

    if (0 == pid[2]) {
	/* reads from p[1] and multiplies it with random data */
	
	close(p[1][1]);
	close(STDIN_FILENO);
	DUP(p[1][0]);

	pip = sf_input("in");

	init_genrand(dseed);
	dp = 0.;
	for (dsiz=nd, dbuf=nbuf; dsiz > 0; dsiz -= dbuf) {
	    if (dsiz < dbuf) dbuf=dsiz;

	    sf_floatread(buf,dbuf,pip);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (!sf_histfloat(in,"o3",&cmp0)) cmp0=0.;

    /* get other parameters */
    if (!sf_getint("w",&w)) w=50;  
    /* sliding window size*/    
    if (!sf_getfloat("sft",&sft)) sft=1;  
    /*weight shift*/
    if (!sf_getfloat("ee",&ee)) ee=1.0;      
    if (!sf_getfloat("esp",&esp)) esp=1.0;       
    
    sf_putint(out,"n1",nt);
    sf_putint(out,"n2",ncmp);

    sf_putfloat(out,"d1",dt);
    sf_putfloat(out,"o1",t0);
    sf_putfloat(out,"d2",dcmp);
    sf_putfloat(out,"o2",cmp0);
    sf_putint(out,"n3",1);
    
    
                  
    indata = sf_floatalloc(nt*nh);
    outdata = sf_floatalloc(nt);
    outweight = sf_floatalloc(nt*nh);
    stack = sf_floatalloc(nt);
    win1 = sf_floatalloc(w);
    win2 = sf_floatalloc(w);
 
    for (icmp=0; icmp < ncmp; icmp++){ 
        
         sf_floatread(indata,nt*nh,in);

         for (it=0; it < nt; it++){
              stack[it] = 0;
              outdata[it] = 0;
	      for (ih=0;ih<nh;ih++){
	           outweight[nt*ih+it]=0;
	      }
         }
	 
	 
         /* computer the directly stack trace */
         for (it=0; it < nt; it++){
              zero=0;
              for(ih=0; ih < nh; ih++){
                 if (indata[ih*nt+it]!=0)
                     zero++;
                 
                 stack[it]+= indata[ih*nt+it];
              }
              if (zero==0)
              stack[it]=stack[it];
              else 
              stack[it]=stack[it]/zero;
         }         
         /* estimate the noise variances */
         for (it=0; it < nt; it++){ 
              zero = 0;
	       sumweight=0;
              for (ih=0; ih < nh; ih++){
		  /* sumwin1 = 0;
                    sumwin2 = 0;
                    sumab = 0; */
                   shift = SF_MAX(0,SF_MIN(nt-w, it-w/2-1));
		   /* weight=0; */
		   var=0;
		   maxwin1=0;
		   maxwin2=0;
                   for (iw=0; iw <w; iw++){
                        win1[iw] = indata[ih*nt+iw+shift];
                        win2[iw] = stack[iw+shift];
			if (fabs(win1[iw])>fabs(maxwin1))
			     maxwin1 = fabs(win1[iw]);
			if (fabs(win2[iw])>fabs(maxwin2))
			     maxwin2 = fabs(win2[iw]);
	           }
		   for (iw=0; iw <w; iw++){		     
			var += fabs(win1[iw]/maxwin1-win2[iw]/maxwin2)*fabs(win1[iw]/maxwin1-win2[iw]/maxwin2);    
                   }
                   outweight[nt*ih+it] =1./(var+ee); 
                   
                   if (indata[ih*nt+it]!=0){
                     zero++;
                     outdata[it] += pow(1./(var+ee),1)*indata[ih*nt+it];
                     sumweight +=  pow(1./(var+ee),1);
		   }
              }
              if (zero==0)
                  outdata[it]=0;              
              else     
                  outdata[it] = outdata[it]/(zero*sumweight+esp);
              
         }
	 
         sf_floatwrite(outdata,nt,out);
         sf_warning("running cmp is = %d of %d",icmp, ncmp); 
     }

    exit(0);
}