int main(int argc, char* argv[])
{
    bool velocity, sym, escvar;
    int is, n[3], im, nm, order, nshot, ndim, two, na, nb;
    int nt, nt1, nr, ir, it, i, ia, ib;
    float t, dt, da=0., a0, amax, db=0., b0, bmax, v0;
    float x[3], p[3], d[3], o[3], **traj, *slow, **s, *a, *b;
    raytrace rt;
    sf_file shots, vel, rays, angles;

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

    /* get 2-D grid parameters */
    if (!sf_histint(vel,"n1",n))     sf_error("No n1= in input");
    if (!sf_histint(vel,"n2",n+1))   sf_error("No n2= in input");
    if (!sf_histint(vel,"n3",n+2))   sf_error("No n3= in input");
    if (!sf_histfloat(vel,"d1",d))   sf_error("No d1= in input");
    if (!sf_histfloat(vel,"d2",d+1)) sf_error("No d2= in input");
    if (!sf_histfloat(vel,"d3",d+2)) sf_error("No d3= in input");
    if (!sf_histfloat(vel,"o1",o))   o[0]=0.;
    if (!sf_histfloat(vel,"o2",o+1)) o[1]=0.;
    if (!sf_histfloat(vel,"o3",o+2)) o[2]=0.;

    /* additional parameters */
    if(!sf_getbool("vel",&velocity)) velocity=true;
    /* If y, input is velocity; if n, slowness */
    if(!sf_getint("order",&order)) order=4;
    /* Interpolation order */

    if (!sf_getint("nt",&nt)) sf_error("Need nt=");
    /* Number of time steps */
    if (!sf_getfloat("dt",&dt)) sf_error("Need dt=");
    /* Sampling in time */

    if (!sf_getbool("sym",&sym)) sym=true;
    /* if y, use symplectic integrator */

    if(!sf_getbool("escvar",&escvar)) escvar=false;
    /* If y - output escape values, n - trajectories */

    /* get shot locations */
    if (NULL != sf_getstring("shotfile")) {
        /* file with shot locations */
        shots = sf_input("shotfile");
        if (!sf_histint(shots,"n1",&ndim) || 3 != ndim) 
            sf_error("Must have n1=2 in shotfile");
        if (!sf_histint(shots,"n2",&nshot)) 
            sf_error("No n2= in shotfile");
  
        s = sf_floatalloc2 (ndim,nshot);
        sf_floatread(s[0],ndim*nshot,shots);
        sf_fileclose (shots);
    } else {
        nshot = 1;
        ndim = 3;

        s = sf_floatalloc2 (ndim,nshot);

        if (!sf_getfloat("zshot",&s[0][0])) s[0][0]=o[0];
        /* shot location in depth (if shotfile is not specified) */
        if (!sf_getfloat("yshot",&s[0][1])) s[0][1]=o[1] + 0.5*(n[1]-1)*d[1];
        /* shot location inline (if shotfile is not specified) */
        if (!sf_getfloat("xshot",&s[0][2])) s[0][2]=o[2] + 0.5*(n[2]-1)*d[2];
        /* shot location crossline (if shotfile is not specified) */

        sf_warning("Shooting from z=%g, y=%g, x=%g",s[0][0],s[0][1],s[0][2]);
    }


    if (NULL != sf_getstring("anglefile")) {
        /* file with initial angles */
        angles = sf_input("anglefile");

        if (!sf_histint(angles,"n1",&nr)) sf_error("No n1= in anglefile");
        if (!sf_histint(angles,"n2",&two) || 2 != two) sf_error("Need n2=2 in anglefile");

        a = sf_floatalloc(nr);
        b = sf_floatalloc(nr);
    } else {
        angles = NULL;

        if (!sf_getint("na",&na)) sf_error("Need na=");
        /* Number of azimuths (if anglefile is not specified) */
        if (!sf_getint("nb",&nb)) sf_error("Need nb=");
        /* Number of inclinations (if anglefile is not specified) */

        if (!sf_getfloat("a0",&a0)) a0 = 0.; 
        /* First azimuth angle in degrees (if anglefile is not specified) */
        if (!sf_getfloat("amax",&amax)) amax=360.;
        /* Maximum azimuth angle in degrees (if anglefile is not specified) */

        if (!sf_getfloat("b0",&b0)) b0 = 0.; 
        /* First inclination angle in degrees (if anglefile is not specified) */
        if (!sf_getfloat("bmax",&bmax)) bmax=180.;
        /* Maximum inclination angle in degrees (if anglefile is not specified) */
        
        /* convert degrees to radians */
        a0 *= SF_PI/180.;
//.........这里部分代码省略.........

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

    off_t pos;
    bool sign, half;
    int   nsx, nsy,   nmx, nmy,   nhx, nhy,   nmhx, nmhy,   nt;
    int   isx, isy,   imx, imy,   ihx, ihy;
    float osx, osy, dsx, dsy, omx, omy, dmx, dmy, ohx, ohy, dhx, dhy, dmhx, dmhy, binx, biny;
    float s_xmin, s_ymin, r_xmin, r_ymin, s_xmax, s_ymax, r_xmax, r_ymax, survey_xmin, survey_ymin, survey_xmax, survey_ymax;
    char *trace, *zero;
    sf_file in, out;


    sf_init(argc,argv);
    in  = sf_input ( "in");
    out = sf_output("out");
    sf_warning("WARNING: This 3-D CMP sorting code is not yet ready for use--this message will be removed when it is. ");

    if (!sf_histint  (in,"n1",&nt)) sf_error("No n1= in input");
    /* Number of samples per trace */

    if (!sf_histint  (in,"n2",&nhx)) sf_error("No n2= in input");
    /* Number of offsets per shot gather in x-direction*/
    if (!sf_histfloat(in,"o2",&ohx)) sf_error("No o2= in input");
    /* First offset x-component */
    if (!sf_histfloat(in,"d2",&dhx)) sf_error("No d2= in input");
    /* Offset increment in x-direction */
    if (!sf_histint  (in,"n3",&nhy)) sf_error("No n3= in input");
    /* Number of offsets per shot gather in y-direction*/
    if (!sf_histfloat(in,"o3",&ohy)) sf_error("No o3= in input");
    /* First offset y-component */
    if (!sf_histfloat(in,"d3",&dhy)) sf_error("No d3= in input");
    /* Offset increment in y-direction */
    if (!sf_histint  (in,"n4",&nsx)) sf_error("No n4= in input");
    /* Number of sources along x-direction*/
    if (!sf_histfloat(in,"d4",&dsx)) sf_error("No d4= in input");
    /* Source spacing in x-direction*/
    if (!sf_histfloat(in,"o4",&osx)) sf_error("No o4= in input");
    /* First source x-coordinate*/
    if (!sf_histint  (in,"n5",&nsy)) sf_error("No n5= in input");
    /* Number of sources along y-direction*/
    if (!sf_histfloat(in,"d5",&dsy)) sf_error("No d5= in input");
    /* Source spacing in y-direction*/
    if (!sf_histfloat(in,"o5",&osy)) sf_error("No o5= in input");
    /* First source y-coordinate*/
    
    if (!sf_getfloat("binx",&binx)) sf_error("No x bin size specified. Please set binx=");
    /*Number of bins along x-direction*/
    if (!sf_getfloat("biny",&biny)) sf_error("No y bin size specified. Please set biny=");
    /*Number of bins along y-direction*/

    if (!sf_getbool("positive",&sign)) sign=true;
    /* initial offset orientation:
       yes is generally for off-end surveys, where the first offsets are positive.  
       no is generally for split-spread surveys with first negative then positive offsets. */

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

    if (!half) {
	dhx /= 2;
	ohx /= 2;
	dhy /= 2;
	ohy /= 2;
    }

    s_xmin = osx;
    s_ymin = osy;
    s_xmax = osx+(nsx*dsx);
    s_ymax = osy+(nsy*dsy);
    
    r_xmin = s_xmin+2*ohx;
    r_ymin = s_ymin+2*ohy;
    r_xmax = s_xmax+2*(ohx+dhx*nhx);
    r_ymax = s_ymax+2*(ohy+dhy*nhy);

    if (s_xmin <= r_xmin){
      survey_xmin = s_xmin;
    }else{
      survey_xmin = r_xmin;
    }

    if (s_ymin <= r_ymin){
      survey_ymin = s_ymin;
    }else{
      survey_ymin = r_ymin;
    }

    if (s_xmax <= r_xmax){
      survey_xmax = r_xmax;
    }else{
      survey_xmax = s_xmax;
    }

    if (s_ymax <= r_ymax){
      survey_ymax = r_ymax;
    }else{
      survey_ymax = s_ymax;
    }
    
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    char *unit, *type;
    bool adj, cig, cmp;
    off_t nzx;
    int nt, nx, ny, ns, nh, nz, i, ix, iz, ih, is, ist, iht, ng;
    float *trace, **out, **table, *stable, *rtable, **tablex, *stablex, *rtablex;
    float ds, s0, x0, y0, dy, s, h, h0, dh, dx, ti, t0, t1, t2, dt, z0, dz, tau;
    float aal, tx, aper;
    sf_file dat, mig, tbl, der;

    sf_init (argc,argv);

    if (!sf_getbool("adj",&adj)) adj=true;
    /* y for migration, n for modeling */

    if (!sf_getbool("cmp",&cmp)) cmp=true;
    /* y for CMP gather, n for shot gather */    

    if (adj) {
	dat = sf_input("in");
	mig = sf_output("out");
    } else {
	mig = sf_input("in");
	dat = sf_output("out");
    }

    tbl = sf_input("table"); /* traveltime table */
    der = sf_input("deriv"); /* source derivative table */

    if (adj) {
	if (!sf_histint(dat,"n1",&nt)) sf_error("No n1= in input");
	if (!sf_histint(dat,"n2",&nh)) sf_error("No n2= in input");
	if (!sf_histint(dat,"n3",&ns)) sf_error("No n3= in input");
	
	if (!sf_histfloat(dat,"o1",&t0)) sf_error("No o1= in input");
	if (!sf_histfloat(dat,"d1",&dt)) sf_error("No d1= in input");
	if (!sf_histfloat(dat,"o2",&h0)) sf_error("No o2= in input");
	if (!sf_histfloat(dat,"d2",&dh)) sf_error("No d2= in input");
	if (!sf_histfloat(dat,"o3",&s0)) sf_error("No o3= in input");
	if (!sf_histfloat(dat,"d3",&ds)) sf_error("No d3= in input");
    } else {	
	if (!sf_getint("nt",&nt)) sf_error("Need nt="); /* time samples */
	if (!sf_getint("nh",&nh)) nh=1; /* offset/receiver samples */
	if (!sf_getint("ns",&ns)) ns=1; /* shot samples */
	
	if (!sf_getfloat("t0",&t0)) t0=0.0; /* time origin */
	if (!sf_getfloat("dt",&dt)) sf_error("Need dt="); /* time sampling */
	if (!sf_getfloat("h0",&h0)) h0=0.0; /* offset/receiver origin */
	if (!sf_getfloat("dh",&dh)) sf_error("Need dh="); /* offset/receiver sampling */
	if (!sf_getfloat("s0",&s0)) s0=0.0; /* shot origin */
	if (!sf_getfloat("ds",&ds)) sf_error("Need ds="); /* shot sampling */
    }

    if (1==nh) dh=0.0;
    if (1==ns) ds=0.0;

    if (!sf_histint(tbl,"n1",&nz)) sf_error("No n1= in table");
    if (!sf_histint(tbl,"n2",&nx)) sf_error("No n2= in table");
    if (!sf_histint(tbl,"n3",&ny)) sf_error("No n3= in table");

    if (!sf_histfloat(tbl,"o1",&z0)) sf_error("No o1= in table");
    if (!sf_histfloat(tbl,"d1",&dz)) sf_error("No d1= in table");

    if (!sf_histfloat(tbl,"o2",&x0)) sf_error("No o2= in table");
    if (!sf_histfloat(tbl,"d2",&dx)) sf_error("No d2= in table");

    if (!sf_histfloat(tbl,"o3",&y0)) sf_error("No o3= in table");
    if (!sf_histfloat(tbl,"d3",&dy)) sf_error("No d3= in table");

    if (!sf_getfloat("tau",&tau)) tau=0.;
    /* static time-shift (in second) */

    if (!sf_getfloat("aperture",&aper)) aper=90.;
    /* migration aperture (in degree) */

    if (!sf_getfloat("antialias",&aal)) aal=1.0;
    /* antialiasing */
    
    if (!sf_getbool("cig",&cig)) cig=false;
    /* y - output common offset/receiver gathers */

    ng = cig? nh: 1;

    if (adj) {
	sf_putint(mig,"n1",nz);
	sf_putint(mig,"n2",nx);
	
	sf_putfloat(mig,"o1",z0);
	sf_putfloat(mig,"d1",dz);
	sf_putfloat(mig,"o2",x0);
	sf_putfloat(mig,"d2",dx);
	sf_putstring(mig,"label1","Depth");
	sf_putstring(mig,"label2","Lateral");
	unit = sf_histstring(dat,"unit2");
	if (NULL != unit) sf_putstring(mig,"unit1",unit);
	if (cig) {
	    sf_putint(mig,"n3",nh);
	    sf_putfloat(mig,"o3",h0);
	    sf_putfloat(mig,"d3",dh);
//.........这里部分代码省略.........

int main (int argc,char* argv[]) 
{
    int i, ip, n1, n2, np, ndim, order, n123, *pp, n[2], box[2], *shift[2], b;
    float o1, o2, d1, d2, slow, *dd, **pts, *vv, *h, *bin, *vor, d[2], *rect[2];
    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);  
    distance(np,pts,dd,vv,pp,
	     1,n2,n1,
	     0.,o2,o1,
	     1.,d2,d1,
	     order);

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

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

{
    float *tabout                            ;
    float  oh, os, oy, or, dh, ds, dy, dr    ;
    float  sloc, hloc, yloc, rloc            ;
    int    nh, ns, ny, nr                    ;
    int    mo, ms, my, mr                    ;
    int    i, j, k, l                        ;
    sf_file inp, out;

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

    if (SF_FLOAT != sf_gettype(inp)) sf_error("Need float input");
    if (!sf_histint(inp, "n1"   , &nh)) sf_error("No n1= in input");
    if (!sf_histint(inp, "n2"   , &ns)) sf_error("No n2= in input");
    if (!sf_histfloat(inp, "o1" , &oh)) sf_error("No o1= in input");
    if (!sf_histfloat(inp, "o2" , &os)) sf_error("No o2= in input");
    if (!sf_histfloat(inp, "d1" , &dh)) sf_error("No d1= in input");
    if (!sf_histfloat(inp, "d2" , &ds)) sf_error("No d2= in input");
		    
    if (!sf_getint("mo", &mo)) sf_error("Need mo=");
    /* offset parameter, a constant offset line will appear in the output every o offset */
    if (!sf_getint("ms", &ms)) sf_error("Need ms=");
    /* source parameter, a constant source line will appear in the output every s source */
    if (!sf_getint("my", &my)) sf_error("Need my=");
    /* midpoint parameter, a constant midpoint line will appear in the output every y midpoint */
    if (!sf_getint("mr", &mr)) sf_error("Need mr=");
    /* receiver parameter, a constant receiver line will appear in the output every r receiver */

    ny = 2 * (ns - 1) + nh ;
    oy = oh/2              ;
    dy = dh/2              ;
    
    nr = ns - 1 + nh ;
    or = os + oh     ;
    dr = dh          ;

    tabout = sf_floatalloc( nh * ns);

    for ( k = 0 ; k < nh * ns ; k++) 
	tabout[k] = 1.0f;

    for ( i = 0 ; i < nh ; i ++ ) {
	for ( j = 0 ; j < ns ; j ++ ) {
	    hloc = oh + i * dh   ;
	    sloc = os + j * ds   ;
	    yloc = sloc + hloc/2 ;
	    rloc = sloc + hloc   ;

	    k = (int) floorf( (yloc - oy)/dy + 0.5 ) ;
	    l = (int) floorf( (rloc - or)/dr + 0.5 ) ;

	    assert( k >= 0 && k < ny ) ;
	    assert( l >= 0 && l < nr ) ;
	    
	    if ( !(i % mo) || !(j % ms) || !(k % my) || !(l % mr) ) {
		tabout[nh * j + i] = 0.0 ;
	    }
	}
    }

    sf_floatwrite(tabout, nh * ns, out);
    
    exit(0);
}

int main(int argc, char* argv[])
{ 
    int i1, n1, i2, m2, m3, n2, n, order, ng, ig, i0, w, nw, iw, jw, nw1, nw2;
    int wi;
    float *coord, **inp, *out, **oth, o1, d1, o2, d2, g0, dg, g;
    float *corr, *win1, *win2, a, b, a2, b2, ab, h, dw;
    bool taper, diff, verb, shift;
    sf_file in, warped, other;

    sf_init (argc, argv);
    in = sf_input("in");
    warped = 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_histfloat(in,"d1",&d1)) sf_error ("No d1= in input");
    if(!sf_histfloat(in,"o1",&o1)) o1 = 0.;

    if(!sf_histint(in,"n2",&m2)) m2 = 1;
    if(!sf_histint(in,"n3",&m3)) m3 = 1;

    if (!sf_getint("ng",&ng)) ng=1;
    /* number of gamma values */
    if (!sf_getfloat("g0",&g0)) sf_error("Need g0=");
    /* gamma origin */
    if (!sf_getfloat("dg",&dg)) dg=g0;
    /* gamma sampling */

    if (!sf_getbool("shift",&shift)) shift=false;
    /* use shift instead of stretch */

    other = sf_input("other");

    if(!sf_histint(other,"n1",&n2)) sf_error ("No n1= in other");
    if(!sf_histfloat(other,"d1",&d2)) sf_error ("No d1= in other");
    if(!sf_histfloat(other,"o1",&o2)) o2 = 0.;

    if (m3 > 1) {
	sf_putint  (warped,"n4",ng);
	sf_putfloat(warped,"d4",dg);
	sf_putfloat(warped,"o4",g0);
    } else if (m2 > 1) {
	sf_putint  (warped,"n3",ng);
	sf_putfloat(warped,"d3",dg);
	sf_putfloat(warped,"o3",g0);
    } else {
	sf_putint  (warped,"n2",ng);
	sf_putfloat(warped,"d2",dg);
	sf_putfloat(warped,"o2",g0);
    }

    m2 *= m3;
    n = n2*m2;

    if(!sf_getint("accuracy",&order)) {
	/* [1-4] interpolation accuracy */
	order = 2;
    } else if (order < 1 || order > 4) {
	sf_error ("accuracy must be between 1 and 4");
    }
    order *= 2;

    if (!sf_getint("nw",&nw)) sf_error ("Need nw=");
    /* number of windows */
    if (!sf_getint("w",&w)) sf_error ("Need w=");
    /* window size */
    if (!sf_getfloat("h",&h)) h=0.5*(w-1);
    /* window overlap */
    if (!sf_getbool("taper",&taper)) taper=true;
    /* window tapering */
    if (!sf_getbool("verb",&verb)) verb=true;
    /* w */

    if (!sf_getbool("diff",&diff)) diff=false;
    /* if y, compute difference power insted of correlation */

    dw = (n2-w)/(nw-1.);
    nw1 = (0.5*w+1.)/dw;
    nw2 = (0.5*w-2.)/dw;
    nw += nw1 + nw2;

    sf_putint(warped,"n1",nw);
    sf_putfloat(warped,"o1",o2+(0.5*w+1.)*d2-nw2*dw*d2);
    sf_putfloat(warped,"d1",dw*d2);

    window1_init (h,dw);

    coord = sf_floatalloc (n2); 
    inp =   sf_floatalloc2 (n1,m2);
    out =   sf_floatalloc (n2);
    oth =   sf_floatalloc2 (n2,m2);

    corr =  sf_floatalloc (nw);

    win1 = sf_floatalloc (w);
    win2 = sf_floatalloc (w);
    win2 = sf_floatalloc (w);

    sf_prefilter_init (order, n1, order*10);     
    for (i2=0; i2 < m2; i2++) {
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    int dim, n[SF_MAX_DIM], rect[SF_MAX_DIM], nd;
    int ns, nt, nm, nx, i, ix, is, it, niter, imax;
    float **slice, *pick0, *pick, *ampl;
    float s0, ds, t, asum, ai, amax, ap, am, num, den, dx;
    char key[6];
    sf_file in, out;

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

    dim = sf_filedims (in,n);
    if (dim < 2) sf_error("Need at least two dimensions");

    nd = 1;
    for (i=0; i < dim; i++) {
	nd *= n[i];
	if (n[i] > 1) {
	    snprintf(key,6,"rect%d",i+1);
	    if (!sf_getint(key,rect+i)) rect[i]=1;
	} else {
	    rect[i]=1;
	}
    }
    nt = n[0];
    ns = n[1];
    nm = nd/ns;
    nx = nm/nt;

    if (!sf_histfloat(in,"o2",&s0)) sf_error("No o2= in input");
    if (!sf_histfloat(in,"d2",&ds)) sf_error("No d2= in input");
    sf_putint(out,"n2",1);

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

    for (i=1; i < dim-1; i++) {
	n[i] = n[i+1];
    }
    dim--;

    divn_init(dim,nm,n,rect,niter);

    slice = sf_floatalloc2(nt,ns);

    pick0 = sf_floatalloc(nm);
    pick = sf_floatalloc(nm);
    ampl = sf_floatalloc(nm);

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

	/* pick blind maximum */
	for (it = 0; it < nt; it++) {
	    i = ix*nt+it;

	    imax = 0;
	    amax = 0.;
	    for (is = 0; is < ns; is++) {
		t = slice[is][it];
		if (t > amax) {
		    amax = t;
		    imax = is;
		}
	    }

	    /* quadratic interpolation for sub-pixel accuracy */
	    if (imax > 0 && imax < ns-1) {
		am = slice[imax-1][it];
		ap = slice[imax+1][it];
		num = 0.5*(am-ap);
		den = am+ap-2.*amax;
		dx = num*den/(den*den + FLT_EPSILON);
		if (fabsf(dx) >= 1.) dx=0.;

		ampl[i] = amax - 0.5*dx*dx*den;
		pick0[i] = s0+(imax+dx)*ds;
	    } else {
		ampl[i] = amax;
		pick0[i] = s0+imax*ds;
	    }
	}
    }
    
    /* normalize amplitudes */
    asum = 0.;
    for (i = 0; i < nm; i++) {
	ai = ampl[i];
	asum += ai*ai;
    }
    asum = sqrtf (asum/nm);

    for(i=0; i < nm; i++) {
	ampl[i] /= asum;
	pick0[i] *= ampl[i];
    }

    divn(pick0,ampl,pick);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    int n1, i1, n2, i2, i, hilbn, it;
    bool norm;
    int fft_size;
    float d1, o1, mt, ee, max, freq, hilbc, *trace, *freqc, *phase, *ric, *hilb, *outtrace;
  
    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");
    n2 = sf_leftsize(in,1);

    /* Frequency and phase*/
    if (NULL == sf_getstring("tphase") || NULL == sf_getstring("tfreq")) {
        if (!sf_getfloat("frequency",&freq)) {
			/* peak frequency for Ricker wavelet (in Hz) */
			if (!sf_getfloat("freq",&freq)) freq=0.2;
				/* peak frequency for Ricker wavelet (as fraction of Nyquist) */
		} else {
			if (!sf_histfloat(in,"d1",&d1)) d1=1.;
			freq *= 2.*d1;
        }

        trace = sf_floatalloc(n1);
        fft_size = 2*kiss_fft_next_fast_size((n1+1)/2);
        ricker_init(fft_size, 0.5*freq, 0);

        for (i2=0; i2 < n2; i2++) {
			sf_floatread(trace,n1,in);
			sf_freqfilt(n1,trace);
			sf_floatwrite(trace,n1,out);
        }
    } else { /*frequency and phase are time-varying*/
        
        sf_file tpha, tfre;
        if (!sf_histfloat(in,"o1",&o1)) o1=0;
        if (!sf_histfloat(in,"d1",&d1)) d1=0.001;
        if (!sf_getfloat("esp",&ee))  ee=0.; 
        /* if norm=y, stable parameter*/
        if (!sf_getbool("norm",&norm)) norm=false;
        tfre = sf_input ("tfreq");
        tpha = sf_input ("tphase");
        if (!sf_histint(tpha,"n1",&i) || i != n1) sf_error("the phase file has not same length");
        if (!sf_histint(tfre,"n1",&i) || i != n1) sf_error("the frequency file has not same length");
        if (!sf_getint("hiborder",&hilbn)) hilbn=6;
	/* Hilbert transformer order */
        if (!sf_getfloat("hibref",&hilbc)) hilbc=1.;
        trace = sf_floatalloc(n1);
        freqc = sf_floatalloc(n1);
        phase = sf_floatalloc(n1);
        outtrace= sf_floatalloc(n1);
        ric   = sf_floatalloc(n1);
        sf_floatread(freqc,n1,tfre);
        sf_floatread(phase,n1,tpha);

         
        sf_hilbert_init(n1, hilbn, hilbc);
        hilb = sf_floatalloc(n1);
        for (i2=0;i2 < n2; i2++) {
			sf_floatread(trace,n1,in);
			for (it=0; it < n1; it++) {
				outtrace[it] = 0.;
			}
			for (i1=0; i1 < n1; i1++) {
	
				max = 0.; 
				for (it=0;it < n1; it++) {

					mt=d1*(i1-it);
					ric[it] = (1-2*(SF_PI*freqc[i1]*mt*SF_PI*freqc[i1]*mt))*exp(-SF_PI*freqc[i1]*mt*SF_PI*freqc[i1]*mt);
	
				}
				sf_hilbert(ric,hilb);
				for (it=0; it < n1; it++) {
					hilb[it] = ric[it]*cosf(phase[i1]*SF_PI/180.) + hilb[it]*sinf(phase[i1]*SF_PI/180.);
					if (hilb[it] > max) max=hilb[it];
				}
				for (it=0;it < n1; it++) {
					if (!norm) {
						max = 1;
						ee  = 0.;
					}
					hilb[it] = hilb[it]*trace[i1]/(max+ee);
					outtrace[it] += hilb[it];
				}
			}
			sf_floatwrite(outtrace,n1,out);
        }
    }
 
    exit(0);
}

int main(int argc, char* argv[])
{
    int n2, i2, nx, ny, niter, nliter;
    float x0, dx, y0, dy, a[3], **data=NULL, **pred=NULL;

    bool verb;
    sf_file in, out, ma;

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

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

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

    if (!sf_getfloat("a0",&a[0])) a[0]=1.;
    /* starting sharpness in xx */
    if (!sf_getfloat("b0",&a[1])) a[1]=0.;
    /* starting sharpness in xy */
    if (!sf_getfloat("c0",&a[2])) a[2]=1.;
    /* starting sharpness in yy */

    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */
    if (!sf_getint("nliter",&nliter)) nliter=1;
    /* number of reweighting iterations */
    if (!sf_getbool("verb",&verb)) verb=false;
    /* verbosity flag */

    sf_putint(ma,"n1",3);
    sf_putint(ma,"n2",1);
    sf_putint(ma,"nx",nx);
    sf_putfloat(ma,"dx",dx);
    sf_putfloat(ma,"x0",x0);
    sf_putint(ma,"ny",ny);
    sf_putfloat(ma,"dy",dy);
    sf_putfloat(ma,"y0",y0);
    sf_fileflush(ma,in);

    data = sf_floatalloc2(nx,ny);
    pred = sf_floatalloc2(nx,ny);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(data[0],nx*ny,in);

	monof2(data,pred,nliter,niter,a,nx,dx,x0,ny,dy,y0,verb);

	sf_floatwrite(a,3,ma);

	sf_floatwrite (pred[0],nx*ny,out);
    }

    exit (0);
}

int main(int argc, char* argv[])
{
    bool adj, velocity, l1norm, plane[3], verb;
    int dim, i, count, n[SF_MAX_DIM], it, nt, **m, nrhs, is, nshot=1, *flag, order, iter, niter, stiter, *k, nfreq, nmem;
    float o[SF_MAX_DIM], d[SF_MAX_DIM], **t, *t0, *s, *temps, **source, *rhs, *ds;
    float rhsnorm, rhsnorm0, rhsnorm1, rate, eps, gama;
    char key[4], *what;
    sf_file sinp, sout, shot, time, reco, rece, topo, grad, norm;
    sf_weight weight=NULL;

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

    if (NULL == (what = sf_getstring("what"))) what="tomo";
    /* what to compute (default tomography) */

    switch (what[0]) {
	case 'l': /* linear operator */

	    if (NULL == sf_getstring("time"))
		sf_error("Need time=");
	    time = sf_input("time");

	    /* read operator dimension from time table */
	    dim = sf_filedims(time,n);
	    
	    nt = 1;
	    for (i=0; i < 3; i++) {
		sprintf(key,"d%d",i+1);
		if (!sf_histfloat(time,key,d+i)) sf_error("No %s= in input",key);
		sprintf(key,"o%d",i+1);
		if (!sf_histfloat(time,key,o+i)) o[i]=0.;
		nt *= n[i]; plane[i] = false;
	    }
	    if (dim < 3) {
		n[2] = 1; o[2] = o[1]; d[2] = d[1]; plane[2] = false;
	    }

	    dim = 2;

	    /* read in shot file */
	    if (NULL == sf_getstring("shot"))
		sf_error("Need source shot=");
	    shot = sf_input("shot");
	    
	    if (!sf_histint(shot,"n2",&nshot)) nshot=1;
	    sf_fileclose(shot);

	    /* read in receiver file */
	    m = sf_intalloc2(nt,nshot);
	    if (NULL == sf_getstring("receiver")) {
		for (is=0; is < nshot; is++) {
		    for (it=0; it < nt; it++) {
			m[is][it] = 1;
		    }
		}
	    } else {
		rece = sf_input("receiver");
		sf_intread(m[0],nt*nshot,rece);
		sf_fileclose(rece);
	    }

	    /* number of right-hand side */
	    nrhs = 0;
	    for (is=0; is < nshot; is++) {
		for (it=0; it < nt; it++) {
		    if (m[is][it] == 1) nrhs++;
		}
	    }
	    rhs = sf_floatalloc(nrhs);

	    t = sf_floatalloc2(nt,nshot);
	    sf_floatread(t[0],nt*nshot,time);

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

	    /* initialize fatomo */
	    fatomo_init(dim,n,d,nshot);

	    /* set operators */
	    fatomo_set(t,m);

	    t0 = sf_floatalloc(nt);

	    if (adj) {
		sf_floatread(rhs,nrhs,sinp);

		fatomo_lop(true,false,nt,nrhs,t0,rhs);

		sf_putint(sout,"n1",nt);
		sf_putint(sout,"n2",1);
		sf_putint(sout,"n3",1);
		sf_floatwrite(t0,nt,sout);
	    } else {
		sf_floatread(t0,nt,sinp);

		fatomo_lop(false,false,nt,nrhs,t0,rhs);
		
//.........这里部分代码省略.........

void velocity_read(sf_file velocity,float *vel,sf_file wave,float theta,int ix_center)
/*< read velocity or anisotropy parameter in tilted coordinates >*/
{
  float *tmpvel;
  int nvel[2],ntmpvel[2];
  int mig_nx, mig_nz, vel_nx, vel_nz;
  float mig_dx,mig_dz,mig_ox,mig_oz;
  float vel_dx,vel_dz,vel_ox,vel_oz;
  int iy,ix,iz,n_block; 
  float tmp_x,tmp_z,tmp_wx,tmp_wz,dtmp_x,dtmp_z;
  int tmp_ix,tmp_iz;
 
  sf_histfloat(wave,"d2",&mig_dx); sf_histfloat(wave,"o2",&mig_ox); sf_histint(wave,"n2",&mig_nx); 
  sf_histfloat(wave,"d1",&mig_dz); sf_histfloat(wave,"o1",&mig_oz); sf_histint(wave,"n1",&mig_nz); 

  sf_histfloat(velocity,"d2",&vel_dx); sf_histfloat(velocity,"o2",&vel_ox); sf_histint(velocity,"n2",&vel_nx); 
  sf_histfloat(velocity,"d1",&vel_dz); sf_histfloat(velocity,"o1",&vel_oz); sf_histint(velocity,"n1",&vel_nz); 

  d3(mig_nx,mig_nz,nvel);
  d3(vel_nx,vel_nz,ntmpvel);

  tmpvel=sf_floatalloc(vel_nz*vel_nx);
  if (theta >=0)
  {
    for(iy=0;iy<1;iy++){    //do iy=1,ny
      n_block=iy*vel_nx;
      sf_seek(velocity,n_block*vel_nz*sizeof(float),SEEK_SET);
      sf_floatread(tmpvel,vel_nz*vel_nx,velocity);

      for(ix=0;ix<mig_nx;ix++){
        dtmp_x=ix*mig_dx-(ix_center-1.0)*mig_dx;
        for(iz=0;iz<mig_nz;iz++){
          dtmp_z=iz*mig_dz;

          tmp_x=mig_ox+dtmp_x*cos(theta)+dtmp_z*sin(theta);
          tmp_x=(tmp_x-vel_ox)/vel_dx;
          tmp_ix=(int)(tmp_x);  tmp_wx=tmp_x-tmp_ix;
          if (tmp_x <0){
            tmp_ix=0;  tmp_wx=0.5;
          }
          if (tmp_ix >= vel_nx-1 ){
            tmp_ix=vel_nx-2;  tmp_wx=0.5;
          }

          tmp_z=mig_oz-dtmp_x*sin(theta)+dtmp_z*cos(theta);
          tmp_z=(tmp_z- vel_oz)/vel_dz;
          tmp_iz=(int)(tmp_z); tmp_wz=tmp_z-tmp_iz;
          if (tmp_z < 0){
            tmp_iz=0;   tmp_wz=0.5;
          }
          if (tmp_iz >= vel_nz-1 ){ 
            tmp_iz=vel_nz-2; tmp_wz=0.5;
          }
         
          vel[i3(iy,ix,iz,nvel)]=tmpvel[i3(0,tmp_ix,tmp_iz,ntmpvel)]*(1-tmp_wz)*(1-tmp_wx)+ 
                                 tmpvel[i3(0,tmp_ix,tmp_iz+1,ntmpvel)]*tmp_wz*(1-tmp_wx)  + 
                                 tmpvel[i3(0,tmp_ix+1,tmp_iz,ntmpvel)]*(1-tmp_wz)*tmp_wx  + 
                                 tmpvel[i3(0,tmp_ix+1,tmp_iz+1,ntmpvel)]*tmp_wz*tmp_wx;    
        }
      }
    }
  }
  else
  
  {
    theta=-theta;
    for(iy=0;iy<1;iy++){
      n_block=iy*vel_nx;
      sf_seek(velocity,n_block*vel_nz*sizeof(float),SEEK_SET);
      sf_floatread(tmpvel,vel_nz*vel_nx,velocity);
        
      for(ix=mig_nx-1;ix>=0;ix--){
        dtmp_x=(mig_nx-1-ix)*mig_dx-(ix_center-1)*mig_dx;
        for(iz=0;iz<mig_nz;iz++){
          dtmp_z=iz*mig_dz;

          tmp_x=mig_ox-dtmp_x*cos(theta)-dtmp_z*sin(theta);
          tmp_x=(tmp_x-vel_ox)/vel_dx;
          tmp_ix=(int)(tmp_x);  tmp_wx=tmp_x-tmp_ix;
          if (tmp_x <0){ 
            tmp_ix=0;  tmp_wx=0.5;
          }
          if (tmp_ix >= vel_nx-1 ){
            tmp_ix=vel_nx-2;  tmp_wx=0.5;
          }
          tmp_z=mig_oz-dtmp_x*sin(theta)+dtmp_z*cos(theta);
          tmp_z=(tmp_z-vel_oz)/vel_dz;
          tmp_iz=(int)(tmp_z); tmp_wz=tmp_z-tmp_iz;
          if (tmp_z < 0){
            tmp_iz=0;   tmp_wz=0.5;
          }
          if (tmp_iz >= vel_nz-1 ){
            tmp_iz=vel_nz-2; tmp_wz=0.5;
          }
          vel[i3(iy,(mig_nx-1)-ix,iz,nvel)]=tmpvel[i3(0,tmp_ix,tmp_iz,ntmpvel)]*(1-tmp_wz)*(1-tmp_wx)+ 
                                      tmpvel[i3(0,tmp_ix,tmp_iz+1,ntmpvel)]*tmp_wz*(1-tmp_wx)  + 
                                      tmpvel[i3(0,tmp_ix+1,tmp_iz,ntmpvel)]*(1-tmp_wz)*tmp_wx  + 
                                      tmpvel[i3(0,tmp_ix+1,tmp_iz+1,ntmpvel)]*tmp_wz*tmp_wx;
        }
      }        
//.........这里部分代码省略.........

int main (int argc, char* argv[])
{
    fint1 nmo;
    bool half, squared, slow;
    int ix,ih,it,iv, nt,nx,nw, nh, nh2, m, CDPtype, mute, *mask, nv, *fold;
    float dt, t0, h0,dh,h, dy, str, v0,dv,v;
    float **traces, *trace, *off, *stack;
    double *dstack;
    mapfunc nmofunc;
    sf_file cmp, stk, offset, msk;

    sf_init (argc,argv);
    cmp = sf_input("in");
    stk = sf_output("out");
    nmofunc =  nmo_map;
    
    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");
    nx = sf_leftsize(cmp,2);

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

    if (!sf_getint("mute",&mute)) mute=12;
    /* mute zone */

    if (!sf_getint("nv",&nv)) sf_error("Need nv=");
    /* number of velocities */
    
    if (!sf_getfloat("v0",&v0)) sf_error("Need v0=");
    /* first velocity */
    
    if (!sf_getfloat("dv",&dv)) sf_error("Need dv=");
    /* step in velocity */

    sf_putint(stk,"n2",nv);
    sf_putfloat(stk,"o2",v0);
    sf_putfloat(stk,"d2",dv);

    CDPtype=1;
    if (NULL != sf_getstring("offset")) {
	offset = sf_input("offset");
	if (SF_FLOAT != sf_gettype(offset)) sf_error("Need float offset");
	nh2 = sf_filesize(offset);
	if (nh2 != nh && nh2 != nh*nx) sf_error("Wrong dimensions in offset");

	off = sf_floatalloc(nh2);	
	sf_floatread (off,nh2,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) && !sf_getint("CDPtype",&CDPtype)) {
	    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);
	    }
	} 	    

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

	offset = NULL;
    }
    
    if (NULL != sf_getstring("mask")) {
	msk = sf_input("mask");
	if (SF_INT != sf_gettype(msk)) sf_error("Need integer mask");
	nh2 = sf_filesize(msk);
	if (nh2 != nh && nh2 != nh*nx) sf_error("Wrong dimensions in mask");
	mask = sf_intalloc(nh2);
	sf_intread (mask,nh2,msk);
	sf_fileclose(msk);
    } else {
	msk = NULL;
	mask = NULL;
    }

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

    if (!sf_getbool("squared",&squared)) squared=false;
    /* if y, the slowness or velocity is squared */

    if (!sf_getfloat ("h0",&h0)) h0=0.;
    /* reference offset */
    if (half) h0 *= 2.;
    if (!sf_getint("extend",&nw)) nw=4;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool verb, adj, hessian;
    int n1, n2, nf, i;
    float df, of;
    sf_complex **img, *coe0, *coe, *dx, *dr;
    int iter, niter, cgiter, count;
    float rhsnorm, rhsnorm0, rhsnorm1, rate, gamma;
    char *what;
    sf_file in, out, image, coef, grad;

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

    if (NULL == (what = sf_getstring("what"))) what="tomo";
    /* what to compute (default tomography) */

    switch (what[0]) {
	case 'l': /* linear operator */

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

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

	    if (!sf_histint(image,"n1",&n1)) sf_error("No n1 in image.");
	    if (!sf_histint(image,"n2",&n2)) sf_error("No n2 in image.");

	    if (!sf_histint(image,"n3",&nf)) sf_error("No n3 in image.");
	    if (!sf_histfloat(image,"o3",&of)) sf_error("No o3 in image.");
	    if (!sf_histfloat(image,"d3",&df)) sf_error("No d3 in image.");

	    img = NULL;

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

	    coe0 = sf_complexalloc(n1*n2*2);
	    sf_complexread(coe0,n1*n2*2,coef);
	    sf_fileclose(coef);
	    
	    /* allocate temporary memory */
	    dx  = sf_complexalloc(n1*n2*2);
	    dr  = sf_complexalloc(n1*n2*nf);

	    /* read input and write output header */
	    if (adj) {
		sf_complexread(dr,n1*n2*nf,in);

		sf_putint(out,"n3",2);
	    } else {
		sf_complexread(dx,n1*n2*2,in);

		sf_putint(out,"n3",nf);
	    }

	    /* initialize operator */
	    freqmig0_init(n1,n2, of,df,nf, img);

	    /* set operator */
	    freqmig0_set(coe0);

	    /* linear operator */
	    if (adj) {
		freqmig0_oper(true,false,n1*n2*2,n1*n2*nf,dx,dr);
	    } else {
		freqmig0_oper(false,false,n1*n2*2,n1*n2*nf,dx,dr);
	    }
	    	    
	    /* write output */
	    if (adj) {
		sf_complexwrite(dx,n1*n2*2,out);
	    } else {
		sf_complexwrite(dr,n1*n2*nf,out);
	    }

	    break;

	case 't': /* tomography */

	    /* read model 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.");
	    
	    /* read initial guess */
	    coe0 = sf_complexalloc(n1*n2*2);
	    sf_complexread(coe0,n1*n2*2,in);
	    
	    /* read image */
	    if (NULL == sf_getstring("image"))
		sf_error("Need image image=");
	    image = sf_input("image");

	    if (!sf_histint(image,"n3",&nf)) sf_error("No n3 in image.");
//.........这里部分代码省略.........

int main(int argc, char* argv[]) 
{
    int nx, nt, ix, it, nz, iz, isx, isz, nxz, na;
    float dt, dx, dz;
    float **nxt,  **old,  **cur, *wav;
    float ***coef, plap; 
    int *stmp, *s1, *s2, ik, SIZE=3;
    sf_file out, vel, source, G;

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

/*    if (SF_FLOAT != sf_gettype(inp)) sf_error("Need float input"); */
    if (SF_FLOAT != sf_gettype(vel)) sf_error("Need float input");
    if (SF_FLOAT != sf_gettype(source)) sf_error("Need float input");
    if (!sf_histint(vel,"n1",&nz)) sf_error("No n1= in input");
    if (!sf_histfloat(vel,"d1",&dz)) sf_error("No d1= in input");
    if (!sf_histint(vel,"n2",&nx)) sf_error("No n2= in input");
    if (!sf_histfloat(vel,"d2",&dx)) sf_error("No d2= in input");
    if (!sf_getfloat("dt",&dt)) sf_error("Need dt input");
    if (!sf_getint("nt",&nt)) sf_error("Need nt input");
    if (!sf_getint("isx",&isx)) sf_error("Need isx input");
    if (!sf_getint("isz",&isz)) sf_error("Need isz input");
    if (!sf_histint(G,"n1",&nxz)) sf_error("No n1= in input");
    if (!sf_histint(G,"n2",&na)) sf_error("No n2= in input");
    if (nx*nz != nxz) sf_error("nx*nz != nxz");

    sf_putint(out,"n1",nz);
    sf_putfloat(out,"d1",dz);
    sf_putint(out,"n2",nx);
    sf_putfloat(out,"d2",dx);
    sf_putint(out,"n3",nt);
    sf_putfloat(out,"d3",dt);
    sf_putfloat(out,"o3",0.0); 

    sf_warning("nt=%d",nt);
    wav    =  sf_floatalloc(nt);
    sf_floatread(wav,nt,source);
    sf_warning("dt=%g",dt);

    old    =  sf_floatalloc2(nz,nx);
    cur    =  sf_floatalloc2(nz,nx);
    nxt    =  sf_floatalloc2(nz,nx);

    coef   =  sf_floatalloc3(nz,nx,na);