int main (int argc, char* argv[])
{
    bool verb, fsrf, snap, expl, dabc; 
    int  jsnap, ntsnap, jdata;

    /* I/O files */
    sf_file Fwav = NULL; /* wavelet   */
    sf_file Fsou = NULL; /* sources   */
    sf_file Frec = NULL; /* receivers */
    sf_file Fvel = NULL; /* velocity  */
    sf_file Fden = NULL; /* density   */
    sf_file Fdat = NULL; /* data      */
    sf_file Fwfl = NULL; /* wavefield */

    /* cube axes */
    sf_axis at, az, ax;
    sf_axis as, ar;

    int cpuid;
    int nt, nz, nx, ns, nr;
    int it, ia;
    float dt, dz, dx;

    /* FDM structure */
    /* Wee keep these structures for compatibility,
       as soon as PETSC version is finished,
       this will be removed */
    fdm2d fdm = NULL;

    /* I/O arrays */
    float *ww = NULL;           /* wavelet   */
    pt2d  *ss = NULL;           /* sources   */
    pt2d  *rr = NULL;           /* receivers */
    float *dd = NULL;           /* data      */

    float **u,**v;

    /* linear interpolation weights/indices */
    lint2d cs,cr;

    /* wavefield cut params */
    sf_axis acz = NULL, acx = NULL;
    int   nqz, nqx;
    float oqz, oqx;
    float dqz, dqx;
    float **uc = NULL;

    sf_petsc_aimplfd2 aimplfd;
    PetscErrorCode ierr;
    /* PETSc Initialization */
    ierr = PetscInitialize (&argc, &argv, 0, 0); CHKERRQ(ierr);
    MPI_Comm_rank (MPI_COMM_WORLD, &cpuid);

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

    /*------------------------------------------------------------*/

    if (!sf_getbool ("verb", &verb)) verb = false; /* verbosity flag */
    if (!sf_getbool ("snap", &snap)) snap = false; /* wavefield snapshots flag */
    if (!sf_getbool ("free", &fsrf)) fsrf = false; /* free surface flag */
    if (!sf_getbool ("expl", &expl)) expl = false; /* "exploding reflector" */
    if (!sf_getbool ("dabc", &dabc)) dabc = false; /* absorbing BC */
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    /* I/O files */
    Fwav = sf_input ("in" ); /* wavelet   */
    Fvel = sf_input ("vel"); /* velocity  */
    Fsou = sf_input ("sou"); /* sources   */
    Frec = sf_input ("rec"); /* receivers */
    Fwfl = sf_output("wfl"); /* wavefield */
    Fdat = sf_output("out"); /* data      */
    Fden = sf_input ("den"); /* density   */

    /*------------------------------------------------------------*/
    /* axes */
    at = sf_iaxa (Fwav,2); sf_setlabel (at,"t"); if (verb && 0 == cpuid) sf_raxa (at); /* time */
    az = sf_iaxa (Fvel,1); sf_setlabel (az,"z"); if (verb && 0 == cpuid) sf_raxa (az); /* depth */
    ax = sf_iaxa (Fvel,2); sf_setlabel (ax,"x"); if (verb && 0 == cpuid) sf_raxa (ax); /* space */

    as = sf_iaxa (Fsou, 2); sf_setlabel (as, "s"); if (verb && 0 == cpuid) sf_raxa (as); /* sources */
    ar = sf_iaxa (Frec, 2); sf_setlabel (ar, "r"); if (verb && 0 == cpuid) sf_raxa (ar); /* receivers */

    nt = sf_n (at); dt = sf_d (at);
    nz = sf_n (az); dz = sf_d (az);
    nx = sf_n (ax); dx = sf_d (ax);

    ns = sf_n (as);
    nr = sf_n (ar);
    /*------------------------------------------------------------*/

    /*------------------------------------------------------------*/
    /* other execution parameters */
    if (!sf_getint ("jdata", &jdata)) jdata = 1;
    if (snap) {  /* save wavefield every *jsnap* time steps */
        if (!sf_getint ("jsnap", &jsnap)) jsnap = nt;
    }
    /*------------------------------------------------------------*/
//.........这里部分代码省略.........

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

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

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

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

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

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

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

int main (int argc, char *argv[])
{
    bool verb, up2, up3;
    unsigned char update;
    int n1,n2,n3, i1,i2,i3, ns2, ns3, ip, np2, np3, n23;
    int order, np, i4, n4, k2, k3, j2, j3, i, jp, j;
    float eps, ***u, **p1, **p2, **cost, *trace, *q2=NULL, *q3=NULL;
    sf_file inp, out, dip;

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

    if (!sf_histint(inp,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(inp,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(inp,"n3",&n3)) sf_error("No n3= in input");
    n23 = n2*n3;
    n4 = sf_leftsize(inp,3);

    if (!sf_getbool("verb",&verb)) verb=false;
    /* verbosity */
    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* regularization */
    
    if (!sf_getint("order",&order)) order=1;
    /* accuracy order */

    if (!sf_getint("ns2",&ns2)) sf_error("Need ns1=");
    if (!sf_getint("ns3",&ns3)) sf_error("Need ns2=");
    /* spray radius */
    np2 = 2*ns2+1;
    np3 = 2*ns3+1;
    np = np2*np3;

    sf_putint(out,"n2",np);
    sf_shiftdim(inp, out, 2);

    cost = sf_floatalloc2(np2,np3);
    for (i3=0; i3 < np3; i3++) {
	for (i2=0; i2 < np2; i2++) {
	    cost[i3][i2] = hypotf(i2-ns2,i3-ns3);
	}
    }

    predict_init (n1, n2, eps*eps, order, 1, true);
    update_init(np2,np3,*cost);

    u = sf_floatalloc3(n1,np,n23);
    for (i3=0; i3 < n23; i3++) {
	for (ip=0; ip < np; ip++) {
	    for (i1=0; i1 < n1; i1++) {
		u[i3][ip][i1] = 0.;
	    }
	}
    }

    p1 = sf_floatalloc2(n1,n23);
    p2 = sf_floatalloc2(n1,n23);

    for (i=0; i < n23; i++) { 
	sf_floatread(p1[i],n1,dip);
    }

    for (i=0; i < n23; i++) { 
	sf_floatread(p2[i],n1,dip);
    }

    for (i4=0; i4 < n4; i4++) {
	for (i=0; i < n23; i++) { 
	    sf_floatread(u[i][ns3*np2+ns2],n1,inp);
	    
	    i2 = i%n2;
	    i3 = i/n2;

	    for (ip=0; ip < np; ip++) {
		update = get_update(ip,&up2,&up3,&jp);
		
		/* from jp to j */
		k2 = jp%np2;
		k3 = jp/np2;
		
		j2 = i2+k2-ns2;
		j3 = i3+k3-ns3;

		if (j2 < 0 || j2 >= n2 || 
		    j3 < 0 || j3 >= n3) continue;

		j = j2+j3*n2;
		trace = u[j][jp];

		if (update & 1) {		
		    if (up2) {
			if (j2==0) continue;
			j2 = j-1;
			q2 = p1[j2];
			k2 = jp-1;
		    } else {
			if (j2==n2-1) continue;
			j2 = j+1;
//.........这里部分代码省略.........

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

    exit (0);
}

int main (int argc, char *argv[])
{
    int  adj;             /* forward or adjoint */
    int  eic;             /* EIC or CIC */
    bool verb;            /* verbosity */
    float eps;            /* dip filter constant */
    int   nrmax;          /* number of reference velocities */
    float dtmax;          /* time error */
    int   pmx,pmy;        /* padding in the k domain */
    int   tmx,tmy;        /* boundary taper size */

    int nhx, nhy, nhz, nht, nc;
    int nhx2,nhy2,nhz2,nht2;
    float dht, oht;

    sf_axis amx,amy,az;
    sf_axis alx,aly;
    sf_axis aw,ae,ac,aa;
    sf_axis ahx,ahy,ahz,aht;         

    /* I/O files */
    sf_file Bws=NULL;   /*  background wavefield file Bws */
    sf_file Bwr=NULL;   /*  background wavefield file Bwr */
    sf_file Bs=NULL;    /*  background slowness file Bs   */
    sf_file Ps=NULL;    /*  slowness perturbation file Ps */
    sf_file Pi=NULL;    /*  image perturbation file Pi    */
    sf_file Fc=NULL;    /*  CIP coordinates               */
    sf_file Pws=NULL;   /*  perturbed wavefield file Pws */
    sf_file Pwr=NULL;   /*  perturbed wavefield file Pwr */
    sf_file Pti=NULL;  

    int ompnth=1;

    wexcub3d cub; /* wavefield hypercube */
    wexcip3d cip; /* CIP gathers */
    wextap3d tap; /* tapering */
    wexssr3d ssr; /* SSR operator */
    wexlsr3d lsr; /* LSR operator */
    wexslo3d slo; /* slowness */

    wexmvaop3d mva;
    float dsmax;

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

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

    if (!sf_getbool(  "verb",&verb ))  verb =  true; /* verbosity flag */
    if (!sf_getint(   "adj",&adj   ))  sf_error("Specify adjoint!"); /* y=ADJ Back-projection; n=FWD Forward Scattering */
    if (!sf_getint(  "feic",&eic   ))  sf_error("Specify EIC!");     /* extended I.C. flag */
    if (!sf_getfloat(  "eps",&eps  ))   eps =  0.01; /* stability parameter */
    if (!sf_getint(  "nrmax",&nrmax)) nrmax =     1; /* max number of refs */
    if (!sf_getfloat("dtmax",&dtmax)) dtmax = 0.004; /* max time error */
    if (!sf_getint(    "pmx",&pmx  ))   pmx =     0; /* padding on x */
    if (!sf_getint(    "pmy",&pmy  ))   pmy =     0; /* padding on y */
    if (!sf_getint(    "tmx",&tmx  ))   tmx =     0; /* taper on x   */
    if (!sf_getint(    "tmy",&tmy  ))   tmy =     0; /* taper on y   */

    ae  = sf_maxa(1,0,1);
    nhx=nhy=nhz=nht=nc=nhx2=nhy2=nhz2=nht2=0;
    oht = dht = 0.0;
    /*------------------------------------------------------------*/
    /* slowness */

    Bs = sf_input("slo");
    alx = sf_iaxa(Bs,1); sf_setlabel(alx,"lx");
    aly = sf_iaxa(Bs,2); sf_setlabel(aly,"ly");
    az =  sf_iaxa(Bs,3); sf_setlabel(az, "z");

    /*------------------------------------------------------------*/
    /* input file */
    if(adj)
        Pi = sf_input("in");
    else
        Ps = sf_input("in");

    /*------------------------------------------------------------*/
    /* wavefield */
    Bws = sf_input("swfl");
    Bwr = sf_input("rwfl");

    amx = sf_iaxa(Bws,1); sf_setlabel(amx,"mx");
    amy = sf_iaxa(Bws,2); sf_setlabel(amy,"my");
    aw  = sf_iaxa(Bws,4); sf_setlabel(aw ,"w" );

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

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

//.........这里部分代码省略.........
    Frkp = sf_input ("rkp"); /* app. rank */
    Fltp = sf_input ("ltp"); /* left mat  */
    Frtp = sf_input ("rtp"); /* right mat */
    Fwfp = sf_output("wfp"); /* wavefield */
    Frks = sf_input ("rks"); /* app. rank */
    Flts = sf_input ("lts"); /* left mat  */
    Frts = sf_input ("rts"); /* right mat */
    Fwfs = sf_output("wfs"); /* wavefield */

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

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

    nt = sf_n(at); dt = sf_d(at);
    nz = sf_n(az); dz = sf_d(az);
    nx = sf_n(ax); dx = sf_d(ax);
    ny = sf_n(ay); dy = sf_d(ay);

    ns = sf_n(asx)*sf_n(asy);
    nr = sf_n(arx)*sf_n(ary);

    /*------------------------------------------------------------*/
    /* other execution parameters                                 */
    /*------------------------------------------------------------*/
    if(! sf_getint("nbell",&nbell)) nbell=NOP;  /* bell size */
    if(verb) sf_warning("nbell=%d",nbell);
    if(! sf_getint("jdata",&jdata)) jdata=1;
    if(snap) {  /* save wavefield every *jsnap* time steps */
	if(! sf_getint("jsnap",&jsnap)) jsnap=nt;
    }
    if(back) {
        shft = (nt-1)%jsnap;
        sf_warning("For backward extrapolation, make sure nbell(%d)=0",nbell);
    } else shft = 0;

    /*------------------------------------------------------------*/
    /* expand domain for FD operators and ABC                     */
    /*------------------------------------------------------------*/
    if( !sf_getint("nb",&nb)) nb=NOP;

    fdm=fdutil3d_init(verb,fsrf,az,ax,ay,nb,1);
    if(nbell) fdbell3d_init(nbell);

    sf_setn(az,fdm->nzpad); sf_seto(az,fdm->ozpad); if(verb) sf_raxa(az);
    sf_setn(ax,fdm->nxpad); sf_seto(ax,fdm->oxpad); if(verb) sf_raxa(ax);
    sf_setn(ay,fdm->nypad); sf_seto(ay,fdm->oypad); if(verb) sf_raxa(ay);

    /*------------------------------------------------------------*/
    /* 3D vector components                                       */
    /*------------------------------------------------------------*/
    nc=3;
    ac=sf_maxa(nc,0,1); /* output 3 cartesian components */

    /*------------------------------------------------------------*/
    /* setup output data header                                   */
    /*------------------------------------------------------------*/
    sf_settype(Fdat,SF_COMPLEX);

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

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

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

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

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


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

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

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

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

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

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

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

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

	    r2 += r*r;
	}

        y3 *= x5;

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

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

int main(int argc, char* argv[])
{
    int	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[])
{
    bool verb;
    bool rays;

    sf_axis az,ax;  /* Cartesian coordinates */
    sf_axis at,ag;  /* Ray coordinates */
    int it,ig;
    int nz,nx,nt,ng;
    float dt,dg,ot,og;

    float xsou,zsou; /* source coordinates */

    sf_file Fv=NULL; /* velocity file */
    sf_file Fw=NULL; /* wavefronfs file */

    float **vv=NULL; /* velocity       */
    pt2d   *wm=NULL; /* wavefront it-1 */
    pt2d   *wo=NULL; /* wavefront it   */
    pt2d   *wp=NULL; /* wavefront it+1 */

    pt2d    Ro;    /* point  on wft it-1 */
    pt2d Pm,Po,Pp; /* points on wft it   */
    pt2d    Qo;    /* point  on wft it+1 */

    /*------------------------------------------------------------*/
    sf_init(argc,argv);
    if(! sf_getbool("verb",&verb)) verb=false;
    if(! sf_getbool("rays",&rays)) rays=false;

    /* velocity file */
    Fv = sf_input ("in");
    az = sf_iaxa(Fv,1); sf_setlabel(az,"z"); nz=sf_n(az); if(verb) sf_raxa(az);
    ax = sf_iaxa(Fv,2); sf_setlabel(ax,"x"); nx=sf_n(ax); if(verb) sf_raxa(ax);

    vv=sf_floatalloc2(nz,nx); 
    sf_floatread(vv[0],nz*nx,Fv);

    /* source location */
    if(! sf_getfloat("xsou",&xsou)) xsou=sf_o(ax) + nx*sf_d(ax)/2;
    if(! sf_getfloat("zsou",&zsou)) zsou=sf_o(az) + nz*sf_d(az)/2;
    if(verb) fprintf(stderr,"xsou=%f zsou=%f\n",xsou,zsou);

    /* time axis */
    if(! sf_getint  ("nt",&nt)) nt=100;
    if(! sf_getfloat("ot",&ot)) ot=0;
    if(! sf_getfloat("dt",&dt)) dt=0.001;
    at = sf_maxa(nt,ot,dt); sf_setlabel(at,"t");

    /* shooting angle axis */
    if(! sf_getint  ("ng",&ng)) ng= 360;
    if(! sf_getfloat("og",&og)) og=-180;
    if(! sf_getfloat("dg",&dg)) dg= 1;
    ag = sf_maxa(ng,og,dg); sf_setlabel(ag,"g");

    /*------------------------------------------------------------*/
    /* wavefronts file (g,t) */
    Fw = sf_output("out");
    sf_oaxa(Fw,ag,1); if(verb) sf_raxa(ag);
    sf_oaxa(Fw,at,2); if(verb) sf_raxa(at);

    /* set the output to complex */
    sf_putint(Fw,"esize",8);
    sf_settype(Fw,SF_COMPLEX);

    /*------------------------------------------------------------*/
    /* allocate wavefronts */
    wm = pt2dalloc1(ng);
    wo = pt2dalloc1(ng);
    wp = pt2dalloc1(ng);

    /* initialize wavefronts */
    for( ig=0; ig<ng; ig++) {
	wm[ig].x=wo[ig].x=wp[ig].x=0;
	wm[ig].z=wo[ig].z=wp[ig].z=0;
	wm[ig].v=wo[ig].v=wp[ig].v=0;
    }

    /*------------------------------------------------------------*/
    /* init HWT */
    hwt2d_init(vv,az,ax,at,ag);

    /*------------------------------------------------------------*/
    /* construct it=0 wavefront */
    it=0;
    for( ig=0; ig<ng; ig++) {
	wm[ig].x=xsou;
	wm[ig].z=zsou;
	wm[ig].v=hwt2d_getv(wm[ig]);
    }
    pt2dwrite1(Fw,wm,ng,2); /* write wavefront it=0 */

    /*------------------------------------------------------------*/
    /* construct it=1 wavefront */
    it=1;
    for( ig=0; ig<ng; ig++) {
	double d,g;

	d = dt * hwt2d_getv(wm[ig]);
	g = (og+ig*dg) * SF_PI/180;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool impresp;
    int nt,nx,nz,nw,init,i,padfactor,nfilt,nkol,it,ix,iz,iw;
    float v,dx,dz,lambda,sixth,gamma,epsdamp,pi2,dw,dt, w,wov;
    sf_complex wov2, a, b, c, d, cshift;
    float ***ppp;
    sf_complex *pp, *qq;
    cfilter aa, fac1, fac2;
    sf_file out, imp=NULL;

    sf_init(argc,argv);
    out = sf_output("out");
    sf_setformat(out,"native_float");

    if (!sf_getint("nz",&nz)) nz=96;
    if (!sf_getint("nx",&nx)) nx=48;
    if (!sf_getint("nt",&nt)) nt=12;
    if (!sf_getint("nw",&nw)) nw=2;
    if (!sf_getint("init",&init)) init=1;

    if (!sf_getfloat("v",&v)) v=1.;
    if (!sf_getfloat("dz",&dz)) dz=1.;
    if (!sf_getfloat("dx",&dx)) dx=2.;
    if (!sf_getfloat("lambda",&lambda)) lambda=nz*dz/4.;

    sf_putint(out,"n1",nz);
    sf_putint(out,"n2",nx);
    sf_putint(out,"n3",nt);

    aa = allocatechelix(9);

    if (!sf_getfloat("sixth",&sixth)) sixth=0.0833;
    if (!sf_getfloat("gamma",&gamma)) gamma=0.667;
    if (!sf_getfloat("epsdamp",&epsdamp)) epsdamp=0.01;
    if (!sf_getint("padfactor",&padfactor)) padfactor=1024;
    if (!sf_getint("nfilt",&nfilt)) nfilt=nx+2;
    if (!sf_getbool("impresp",&impresp)) impresp=false;

    if (impresp) {
	imp = sf_output("imp");
	sf_setformat(imp,"native_complex");
	sf_putint(imp,"n1",2*nx);
	sf_putint(imp,"n2",2);
    }

    ppp = sf_floatalloc3(nz,nx,nt);
    pp = sf_complexalloc(nx*nz);
    qq = sf_complexalloc(nx*nz);

    pi2 = 2.*SF_PI;
    dw = v*pi2/lambda;
    dt =   pi2/(nt*dw);

    nkol=pad2(padfactor*nx);
    /* dkol=pi2/nkol; */

    for (it=0; it < nt; it++) {
	for (ix=0; ix < nx; ix++) {
	    for (iz=0; iz < nz; iz++) {
		ppp[it][ix][iz] = 0.;
	    }
	}
    }

    fac1 = allocatechelix(nfilt);
    fac2 = allocatechelix(nfilt);
    helimakelag(fac1,nx,nz);
    helimakelag(fac2,nx,nz);

    xkolmog_init(nkol);

    for (iw=0; iw < nw; iw++) { /* frequency loop */
	w = (iw+1)*dw;

	if (impresp) w=dw*nw/2;

	wov = w/v;
	wov2 = sf_cmplx(epsdamp,wov);
#ifdef SF_HAS_COMPLEX_H
	wov2 = -wov2*wov2;
#else
	wov2 = sf_cneg(sf_cmul(wov2,wov2));
#endif

	sf_warning("%g %g (%d of %d)",crealf(wov2),cimagf(wov2),iw,nw);

	init_wave(init,nx,dx,nz,dz,pp,wov,nw,iw);

	for (iz=0; iz < nx*nz; iz++) {
	    qq[iz]=sf_cmplx(0.,0.);
	}

	/* isotropic laplacian = 5-point laplacian */
	a= sf_cmplx(0.,0.);
#ifdef SF_HAS_COMPLEX_H
	b= gamma*(1+sixth*wov2)* (-1./(dz*dz));
	c= gamma*(1+sixth*wov2)* (-1./(dx*dx));
	d= gamma*(1+sixth*wov2)* (2/(dx*dx) + 2/(dz*dz))  -wov2;
#else
//.........这里部分代码省略.........

int main(int  argc,char **argv)
{
    sf_init(argc,argv);

    sf_file Fix, Fiy, Fiz;
    sf_file Fi1, Fi2, Fi3;
    sf_file Fo1, Fo2;

    clock_t t1, t2;
    float   timespent;

    t1=clock();

   int isep;

   if (!sf_getint("isep",&isep)) isep=1;

    /* setup I/O files */
    Fix = sf_input("in");
    Fiy = sf_input("apvyy");
    Fiz = sf_input("apvzz");

    Fi1 = sf_input("PseudoPurePx"); /* pseudo-pure P-wave x-component */
    Fi2 = sf_input("PseudoPurePy"); /* pseudo-pure P-wave y-component */
    Fi3 = sf_input("PseudoPurePz"); /* pseudo-pure P-wave z-component */
    Fo1 = sf_output("out");               /* pseudo-pure scalar P-wave */
    Fo2 = sf_output("PseudoPureSepP");    /* separated scalar P-wave */

    int     nx, ny, nz, nxf, nyf, nzf, hnx, hny, hnz, i, j, k;
    float   fx, fy, fz, dx, dy, dz, dxf, dyf, dzf;

    /* Read/Write axes */
    sf_axis az, ax, ay;

    az = sf_iaxa(Fi1,1); nz = sf_n(az); dz = sf_d(az)*1000.0;
    ax = sf_iaxa(Fi1,2); nx = sf_n(ax); dx = sf_d(ax)*1000.0;
    ay = sf_iaxa(Fi1,3); ny = sf_n(ay); dy = sf_d(ay)*1000.0;
    fy=sf_o(ay)*1000.0;
    fx=sf_o(ax)*1000.0;
    fz=sf_o(az)*1000.0;

    sf_warning("nx=%d ny=%d nz=%d dx=%f dy=%f dz=%f",nx,ny,nz,dx,dy,dz);

    puthead3x(Fo1, nz, nx, ny, dz/1000.0, dx/1000.0, dy/1000.0, 0.0, 0.0, 0.0);
    puthead3x(Fo2, nz, nx, ny, dz/1000.0, dx/1000.0, dy/1000.0, 0.0, 0.0, 0.0);

    float*** px=sf_floatalloc3(nz,nx,ny);
    float*** py=sf_floatalloc3(nz,nx,ny);
    float*** pz=sf_floatalloc3(nz,nx,ny);
    float*** p=sf_floatalloc3(nz,nx,ny);

    int iy, ix, iz;

    for(iy=0;iy<ny;iy++)
    for(ix=0;ix<nx;ix++)
    {
      sf_floatread(px[iy][ix],nz,Fi1);
      sf_floatread(py[iy][ix],nz,Fi2);
      sf_floatread(pz[iy][ix],nz,Fi3);
    }

    for(iy=0;iy<ny;iy++)
    for(ix=0;ix<nx;ix++){
       for(iz=0;iz<nz;iz++)
          p[iy][ix][iz] = px[iy][ix][iz] + py[iy][ix][iz] + pz[iy][ix][iz];
       sf_floatwrite(p[iy][ix],nz,Fo1);
    }

if(isep==1){
    /* Read axes */
    sf_axis azf, axf, ayf;
    azf = sf_iaxa(Fix,1); nzf = sf_n(azf); dzf = sf_d(azf)*1000.0;
    axf = sf_iaxa(Fix,2); nxf = sf_n(axf); dxf = sf_d(axf)*1000.0;
    ayf = sf_iaxa(Fix,3); nyf = sf_n(ayf); dyf = sf_d(ayf)*1000.0;

    if(dx!=dxf){
       sf_warning("dx= %f dxf=%f",dx,dxf);
       sf_warning("filter and data spatial sampling don't match in x-axis");
       exit(0);
    }
    if(dy!=dyf){
       sf_warning("dy= %f dyf=%f",dy,dyf);
       sf_warning("filter and data spatial sampling don't match in y-axis");
       exit(0);
    }
    if(dz!=dzf){
       sf_warning("dz= %f dzf=%f",dz,dzf);
       sf_warning("filter and data spatial sampling don't match in z-axis");
       exit(0);
    }

    float*** apvxx=sf_floatalloc3(nzf,nxf,nyf);
    float*** apvyy=sf_floatalloc3(nzf,nxf,nyf);
    float*** apvzz=sf_floatalloc3(nzf,nxf,nyf);

    hnx=nxf/2;
    hny=nyf/2;
    hnz=nzf/2;

    sf_warning("nxf=%d nyf=%d nzf=%d dxf=%f dyf=%f dzf=%f",nxf,nyf,nzf,dxf,dyf,dzf);
//.........这里部分代码省略.........

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    sx;
  float    sy;
  float    gx;
  float    gy;
  float    offset;
  float offx, offy, midx, midy;
  float scalco, scale;
  int indx_sx;
  int indx_sy;
  int indx_gx;
  int indx_gy;
  int indx_offset;
  int indx_scalco;

  float v, dx, dy, x0, y0, t0, t0xy, t, off_dotprod_dip;
  int it;
  float d1, o1;
  
 /*****************************/
  /* initialize verbose switch */
  /*****************************/
  sf_init (argc,argv);

  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);
 
  if(!sf_getfloat("v",&v))sf_error("v, a required parameter not input");
  if(!sf_getfloat("dx",&dx))sf_error("dx, a required parameter not input");
  if(!sf_getfloat("dy",&dy))sf_error("dy, a required parameter not input");
  if(!sf_getfloat("x0",&x0))sf_error("x0, a required parameter not input");
  if(!sf_getfloat("y0",&y0))sf_error("y0, a required parameter not input");
  if(!sf_getfloat("t0",&t0))sf_error("t0, a required parameter not input");

  /******************************************/
  /* 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);
 
  if(verbose>0)fprintf(stderr,"allocate intrace.  n1_traces=%d\n",n1_traces);
  intrace= sf_floatalloc(n1_traces);



  /* maybe I should add some validation that n1== n1_traces+n1_headers+2
     and the record length read in the second word is consistent with 
     n1.  */

  /**********************************************************/
  /* end code block for standard tah Trace And Header setup */
  /* continue with any sf_puthist this tah program calls to */
  /* add to the history file                                */
  /**********************************************************/

  /* put the history from the input file to the output */
  sf_fileflush(out,in);

  /********************************************************/
  /* continue initialization specific to this tah program */
  /********************************************************/
  if (!sf_histfloat(in,"d1",&d1))
    sf_error("input data not define d1");
  if (!sf_histfloat(in,"o1",&o1))
    sf_error("input data not define o1");

  /* segy_init gets the list header keys required by segykey function  */
  segy_init(n1_headers,in);

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

/* initialize the length of axis and the number of tics */
void vp_axis_init (const sf_file in)
{
    float ftmp, num;
    bool want;
    char *stmp, *label, *unit;
    size_t len;
    struct Axis axis;

    if (!sf_getint ("axisfat",&axisfat)) axisfat=0;
    if (!sf_getint ("axiscol",&axiscol)) axiscol=7;
    if (!sf_getint ("labelfat",&labelfat)) labelfat=0;
    if (!sf_getfloat ("labelsz",&labelsz)) labelsz=8.;
    
    if ((NULL == (label=sf_getstring("label1"))) &&
	(NULL == (label=sf_histstring(in,"label1")))) {  
	axis1.label = blank;
    } else if ((NULL == (unit=sf_getstring("unit1"))) &&
	       (NULL == (unit=sf_histstring(in,"unit1")))) {
	axis1.label = label;
    } else {
	len = strlen(label)+strlen(unit)+4;
	axis1.label = sf_charalloc(len);
	snprintf(axis1.label,len,"%s (%s)",label,unit);
	free(label);
	free(unit);
    }

    if ((NULL == (label=sf_getstring("label2"))) &&
	(NULL == (label=sf_histstring(in,"label2")))) {
	axis2.label = blank;
    } else if ((NULL == (unit=sf_getstring("unit2"))) &&
	       (NULL == (unit=sf_histstring(in,"unit2")))) {
	axis2.label = label;           
    } else {
	len = strlen(label)+strlen(unit)+4;
	axis2.label = sf_charalloc(len);
	snprintf(axis2.label,len,"%s (%s)",label,unit);
	free(label);
	free(unit);
    }

    where1 = (NULL != (stmp = sf_getstring ("wherexlabel")) &&
	      'b' == *stmp);
    where2 = (NULL != (stmp = sf_getstring ("whereylabel")) &&
	      'r' == *stmp);
    
    /* checking to see if wantaxis is fetched */
    if (!sf_getbool ("wantaxis", &want)) {
	/* setting the default to put the axes on the plot */
        want = true;
	axis1.want = true;
	axis2.want = true;
    } else if (!want) {
	axis1.want = false;
	axis2.want = false;
    } else {
	if (!sf_getbool ("wantaxis1",&(axis1.want)))
	    axis1.want = true;
	if (!sf_getbool ("wantaxis2",&(axis2.want)))
	    axis2.want = true;
    }

    /* frame or axis */
    wheretics = NULL != (stmp = sf_getstring ("wheretics")) &&
	'a' == *stmp;

    if (!sf_getfloat ("axisor1",&(axis1.or)))
	axis1.or = where1? min2: max2;

    if (!sf_getfloat ("axisor2",&(axis2.or)))
	axis2.or = where2? min1: max1;

    if (!sf_getint ("n1tic",&(axis1.ntic))) axis1.ntic = 1;
    if (!sf_getint ("n2tic",&(axis2.ntic))) axis2.ntic = 1;

    if (!sf_getfloat ("d1num", &(axis1.dnum))) getscl (&axis1);
    if (!sf_getfloat ("d2num", &(axis2.dnum))) getscl (&axis2);

    if (0. == axis1.dnum) sf_error("%s: zero d1num",__FILE__);
    if (0. == axis2.dnum) sf_error("%s: zero d2num",__FILE__);

    if (!sf_getfloat("o1num",&(axis1.num0))) {
	ftmp = (min1 < max1)? min1: max1;
	for (num = floorf(ftmp / axis1.dnum) * axis1.dnum - axis1.dnum; 
	     num < ftmp; num += axis1.dnum) ;
	axis1.num0 = num;
    }

    if (!sf_getfloat("o2num",&(axis2.num0))) {
	ftmp = (min2 < max2)? min2: max2;
	for (num = floorf(ftmp / axis2.dnum) * axis2.dnum - axis2.dnum; 
	     num < ftmp; num += axis2.dnum) ;
	axis2.num0 = num;
    }

    axis1.dtic = axis1.dnum / axis1.ntic ;
    ftmp = (min1 < max1)? min1: max1;
    for (num = axis1.num0 - axis1.ntic * axis1.dtic;
	 num < ftmp; num += axis1.dtic);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool verb;
    int n[SF_MAX_DIM], m[SF_MAX_DIM], rect[SF_MAX_DIM];
    int ndim, mdim, nd, ns, n12, i, j, niter, na, ia, i4, n4=1;
    float *d, *f, *g, mean, a0;
    char key[6], *peffile, *lagfile;
    sf_filter aa;
    sf_file dat, flt, mat, pre, pef, lag;

    sf_init(argc,argv);

    dat = sf_input("in");
    mat = sf_input("match");
    flt = sf_output("out");

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

    ndim = 3;
    mdim = 2;


    if (!sf_histint(dat,"n1",&n[0])) sf_error("No n1= in input");
    if (!sf_histint(dat,"n2",&n[1])) sf_error("No n2= in input");
    if (!sf_histint(dat,"n3",&n[2])) sf_error("No n3= in input");
    n[3] = sf_leftsize(dat,3);

    if (!sf_histint(mat,"n1",&m[0])) sf_error("No n1= in match");
    if (!sf_histint(mat,"n2",&m[1])) sf_error("No n2= in match");
    m[2] = sf_leftsize(mat,2);

    if (mdim > ndim) 
	sf_error("Wrong dimensions: %d > %d",mdim,ndim);
 
    if (m[2] != n[3]) {
	sf_error("Size mismatch [n%d] : %d != %d",3,m[2],n[3]);
    } else {
	n4 = n[3];
    }

    nd = 1;
    for (j=0; j < mdim; j++) {
	if (m[j] != n[j]) 
	    sf_error("Size mismatch [n%d]: %d != %d",
		     j+1,m[j],n[j]);
	nd *= m[j];
	snprintf(key,6,"rect%d",j+1);
	if (!sf_getint(key,rect+j)) rect[j]=1;
    }
    for (ns = 1; j < ndim; j++) {
	ns *= n[j];
	if (pre) {
	    snprintf(key,6,"n%d",j+1);
	    sf_putint(pre,key,n4);
	    snprintf(key,6,"n%d",j+2);
	    sf_putint(pre,key,1);
	}
    }
    n12 = nd*ns;

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

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

    if (NULL == (peffile = sf_getstring("pef"))) { 
        /* signal PEF file (optional) */
	aa = NULL;
    } else {
	pef = sf_input(peffile);
	if (!sf_histint(pef,"n1",&na)) 
	    sf_error("No n1= in pef");
	aa = sf_allocatehelix (na);
	
	if (!sf_histfloat(pef,"a0",&a0)) a0=1.;
	sf_floatread (aa->flt,na,pef);
	for( ia=0; ia < na; ia++) {
	    aa->flt[ia] /= a0;
	}
	if (NULL != (lagfile = sf_getstring("lag")) 
	    || 
	    NULL != (lagfile = sf_histstring(pef,"lag"))) {
            /* file with PEF lags (optional) */
	    lag = sf_input(lagfile);
	    sf_intread(aa->lag,na,lag);
	    sf_fileclose(lag);
	} else {
	    for( ia=0; ia < na; ia++) {
		aa->lag[ia] = ia+1;
	    }
	}
    }

    d = sf_floatalloc(n12);
    f = sf_floatalloc(n12);
//.........这里部分代码省略.........

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;