C++ sf_putstring函数代码示例

int main(int argc, char* argv[])
{
    int ix, iy, nx, ny;
    float dx, dy, zx, zy, x, y, *trace;
    sf_file angle;
    
    sf_init(argc,argv);
    angle=sf_output("out");
    sf_setformat(angle,"native_float");

    if (!sf_getint("nx",&nx)) nx=451;
    if (!sf_getint("ny",&ny)) ny=451;

    if (!sf_getfloat("dx",&dx)) dx=0.1;
    if (!sf_getfloat("dy",&dy)) dy=0.1;

    if (!sf_getfloat("zx",&zx)) zx=0.;
    if (!sf_getfloat("zy",&zy)) zy=0.;

    zx = tanf(SF_PI*zx/180.);
    zy = tanf(SF_PI*zy/180.);

    sf_putint(angle,"n1",2*nx-1);
    sf_putfloat(angle,"o1",-(nx-1)*dx);
    sf_putfloat(angle,"d1",dx);
    sf_putstring(angle,"label1","In-line Offset Slope");
    sf_putstring(angle,"unit1","degrees");

    sf_putint(angle,"n2",2*ny-1);
    sf_putfloat(angle,"o2",-(ny-1)*dy);
    sf_putfloat(angle,"d2",dy);
    sf_putstring(angle,"label2","Cross-line Offset Slope");
    sf_putstring(angle,"unit2","degrees");

    trace = sf_floatalloc(2*nx-1);
    
    for (iy=-ny+1; iy < ny; iy++) {
	y = tanf(iy*dy*SF_PI/180.);
	for (ix=-nx+1; ix < nx; ix++) {
	    x = tanf(ix*dx*SF_PI/180.);
	    x = x*x*(1.+zx*zx) + 2.*x*y*zx*zy + y*y*(1.+zy*zy);
	    x /= (1.+zx*zx+zy*zy);
	    if (x > 0.) {
		trace[ix+nx-1] = atanf(sqrtf(x)) * 180./SF_PI;
	    } else {
		trace[ix+nx-1] = -1.;
	    }
	}
	sf_floatwrite(trace,2*nx-1,angle);
    }

    exit(0);
}

void puthead4x(sf_file Fo, int n1, int n2, int n3, int n4, float d1, float d2, float d3, float d4, float o1, float o2, float o3, float o4)
/*<  put head for (y,x,z,t) domain float-type 3D data sets>*/
{
        /* Read/Write axes */
        sf_putint(Fo,"n1",n1);
        sf_putint(Fo,"n2",n2);
        sf_putint(Fo,"n3",n3);
        sf_putint(Fo,"n4",n4);
        sf_putfloat(Fo,"d1",d1);
        sf_putfloat(Fo,"d2",d2);
        sf_putfloat(Fo,"d3",d3);
        sf_putfloat(Fo,"d4",d4);
        sf_putfloat(Fo,"o1",o1);
        sf_putfloat(Fo,"o2",o2);
        sf_putfloat(Fo,"o3",o3);
        sf_putfloat(Fo,"o4",o4);
        sf_putstring(Fo,"label1","z");
        sf_putstring(Fo,"label2","x");
        sf_putstring(Fo,"label3","y");
        sf_putstring(Fo,"label4","t");
        sf_putstring(Fo,"unit1","km");
        sf_putstring(Fo,"unit2","km");
        sf_putstring(Fo,"unit3","km");
        sf_putstring(Fo,"unit4","s");
}

int main(int argc, char* argv[])
{
    float x;	    /* trace value       */
    int ns=512;	    /* number of samples */
    int ntr=8;	    /* number of traces  */
    int dtime=16;   /* one way time in samples through "ocean" layer */
    float rbot=0.8; /* reflection strength of "ocean" bottom         */
    int h=100;	    /* location in samples of two way reverb train 	 */
    float amp=0.2;  /* strength of reflector */
    int loc=170;    /* location of reflector on trace 1 in time samples */
    int dip=12;	    /* dip of reflector in time samples */
    int i, j, k, sgn;
    float *trace;
    sf_file traces;

    sf_init(argc,argv);
    traces = sf_output("out");
    sf_setformat(traces,"native_float");
    sf_putint(traces,"n1",ns);
    sf_putint(traces,"n2",ntr);

    sf_putfloat(traces,"o1",0);
    sf_putfloat(traces,"d1",0.004);
    sf_putfloat(traces,"o2",1);
    sf_putfloat(traces,"d2",1);

    sf_putstring(traces,"label1","Time");
    sf_putstring(traces,"unit1","s");
    sf_putstring(traces,"label2","Trace");

    trace = sf_floatalloc(ns);

    for (j = 0; j < ntr; j++) {
        for (i = 0; i < ns; i++) {
            if (i >= h && ((i-h) % dtime == 0)) {
                k = (i-h)/dtime;
                sgn = (ISODD(k) ? -1 : 1);
                x = sgn * (k+1) * powf(rbot, k);
            } else {
                x = 0.0f;
            }
            if (i == loc + j*dip) x += amp;
            trace[i] = x;
        }
        sf_floatwrite(trace,ns,traces);
    }

    exit(0);
}

void puthead2kx(sf_file Fo, int n1, int n2, float d1, float d2, float o1, float o2)
/*<  put head for (kx,kz) domain float-type 2D data sets>*/
{
        /* Read/Write axes */
        sf_putint(Fo,"n1",n1);
        sf_putint(Fo,"n2",n2);
        sf_putfloat(Fo,"d1",d1);
        sf_putfloat(Fo,"d2",d2);
        sf_putfloat(Fo,"o1",o1);
        sf_putfloat(Fo,"o2",o2);
        sf_putstring(Fo,"label1","kz");
        sf_putstring(Fo,"label2","kx");
        sf_putstring(Fo,"unit1","2*pi/m");
        sf_putstring(Fo,"unit2","2*pi/m");
}

void FwiParams::putOutputParams() {
  sf_putint(vupdates, "n1", nz);
  sf_putint(vupdates, "n2", nx);
  sf_putfloat(vupdates, "d1", dz);
  sf_putfloat(vupdates, "d2", dx);
  sf_putstring(vupdates, "label1", "Depth");
  sf_putstring(vupdates, "label2", "Distance");
  sf_putstring(vupdates, "label3", "Iteration");
  sf_putint(vupdates, "n3", niter);
  sf_putint(vupdates, "d3", 1);
  sf_putint(vupdates, "o3", 1);
  sf_putint(grads, "n1", nz);
  sf_putint(grads, "n2", nx);
  sf_putint(grads, "n3", niter);
  sf_putfloat(grads, "d1", dz);
  sf_putfloat(grads, "d2", dx);
  sf_putint(grads, "d3", 1);
  sf_putint(grads, "o3", 1);
  sf_putstring(grads, "label1", "Depth");
  sf_putstring(grads, "label2", "Distance");
  sf_putstring(grads, "label3", "Iteration");
  sf_putint(illums, "n1", nz);
  sf_putint(illums, "n2", nx);
  sf_putfloat(illums, "d1", dz);
  sf_putfloat(illums, "d2", dx);
  sf_putint(illums, "n3", niter);
  sf_putint(illums, "d3", 1);
  sf_putint(illums, "o3", 1);
  sf_putint(objs, "n1", niter);
  sf_putint(objs, "n2", 1);
  sf_putfloat(objs, "d1", 1);
  sf_putfloat(objs, "o1", 1);
}

bool sf_fft_label(int axis          /* axis number */,
		  const char* label /* input label */, 
		  sf_file out       /* output file */) 
/*< Choose an output label appropriately. 
  Returns false if the label name is not recognized. >*/
{
    char varname[8];

    snprintf(varname,8,"label%d",axis);
    if (0==strcmp(label,"Time")) {
	sf_putstring(out,varname,"Frequency");
    } else if (0==strcmp(label,"time")) {
	sf_putstring(out,varname,"frequency");
    } else if (0==strcmp(label,"Frequency")) {
	sf_putstring(out,varname,"Time");
    } else if (0==strcmp(label,"frequency")) {
	sf_putstring(out,varname,"time");
    } else {
	return false;
    }
    return true;
}

void sf_setpars (sf_file file)
/*< change parameters to those from the command line >*/
{
    char *in;
	
    in = sf_histstring(file,"in");
    sf_simtab_close(file->pars);
    file->pars = sf_getpars();
	
    if (NULL != in) {
	sf_putstring(file,"in",in);
	free(in);
    }
}

void sf_fft_unit(int axis         /* axis number */,
		 const char* unit /* input unit */,
		 sf_file out      /* output file */)
/*< Change unit to 1/unit or vice versa >*/
{
    size_t len;
    char *unit2, varname[7];

    snprintf(varname,7,"unit%d",axis);
    if (NULL != unit) {
	if (0==strcmp(unit,"s")) {
	    sf_putstring(out,varname,"Hz");
	} else if (0==strcmp(unit,"Hz")) {
	    sf_putstring(out,varname,"s");
	} else if (0==strncmp(unit,"1/",2)) {
	    sf_putstring(out,varname,unit+2);
	} else {
	    len=strlen(unit)+3;
	    unit2 = sf_charalloc(len);
	    snprintf(unit2,len,"1/%s",unit);
	    sf_putstring(out,varname,unit2);
	}
    }
}

void puthead3kx(sf_file Fo, int n1, int n2, int n3, float d1, float d2, float d3, float o1, float o2, float o3)
/*<  put head for (ky,kx,kz) domain float-type 3D data sets>*/
{
        /* Read/Write axes */
        sf_putint(Fo,"n1",n1);
        sf_putint(Fo,"n2",n2);
        sf_putint(Fo,"n3",n3);
        sf_putfloat(Fo,"d1",d1);
        sf_putfloat(Fo,"d2",d2);
        sf_putfloat(Fo,"d3",d3);
        sf_putfloat(Fo,"o1",o1);
        sf_putfloat(Fo,"o2",o2);
        sf_putfloat(Fo,"o3",o3);
        sf_putstring(Fo,"label1","kz");
        sf_putstring(Fo,"label2","kx");
        sf_putstring(Fo,"label3","ky");
        sf_putstring(Fo,"unit1","2*pi/m");
        sf_putstring(Fo,"unit2","2*pi/m");
        sf_putstring(Fo,"unit3","2*pi/m");
}

void puthead2dcommonshot(sf_file Fo, int n1, int n2, int n3, float d1, float d2, float d3, float o1, float o2, float o3)
/*<  put head for 2D (x,t) domain common-shot data sets>*/
{
        /* Read/Write axes */
        sf_putint(Fo,"n1",n1);
        sf_putint(Fo,"n2",n2);
        sf_putint(Fo,"n3",n3);
        sf_putfloat(Fo,"d1",d1);
        sf_putfloat(Fo,"d2",d2);
        sf_putfloat(Fo,"d3",d3);
        sf_putfloat(Fo,"o1",o1);
        sf_putfloat(Fo,"o2",o2);
        sf_putfloat(Fo,"o3",o3);
        sf_putstring(Fo,"label1","z");
        sf_putstring(Fo,"label2","x");
        sf_putstring(Fo,"label3","y");
        sf_putstring(Fo,"unit1","km");
        sf_putstring(Fo,"unit2","km");
        sf_putstring(Fo,"unit3","km");
}

//.........这里部分代码省略.........
     file name, with the final .rsf changed to _hdr.rsf. 
  */  

  if(outheaders_filename==NULL){
    /* compute headers_filename from output_filename by replacing the final
       .rsf with _hdr.rsf */
    if(!(0==strcmp(output_filename+strlen(output_filename)-4,".rsf"))){
	fprintf(stderr,"parameter output, the name of the output file,\n");
	fprintf(stderr,"does not end with .rsf, so header filename cannot\n");
	fprintf(stderr,"be computed by replacing the final .rsf with\n");
	fprintf(stderr,"_hdr.rsf.\n");
	sf_error("default for outheaders parameter cannot be computed.");
    }
    outheaders_filename=malloc(strlen(output_filename)+60);
    strcpy(outheaders_filename,output_filename);
    strcpy(outheaders_filename+strlen(output_filename)-4,"_hdr.rsf\0");
    if(verbose>1)
      fprintf(stderr,"parameter outheader defaulted.  Computed to be #%s#\n",
			 outheaders_filename);
  }
  if(verbose>1)fprintf(stderr,"parameter outheader input or computed  #%s#\n",
		       outheaders_filename);
  outheaders=sf_output(outheaders_filename);
 
  /* get each of the axis information: 
     label2, n2, o2, d2,  
     label3, n3, o3, d3,
     etc
     label1, n1, o1, d1  is always defaulted from input */
  sf_putint   (output    ,"n1"    ,n1_traces );
  sf_putint   (outheaders,"n1"    ,n1_headers);
  sf_putfloat (outheaders,"d1"    ,1         );
  sf_putfloat (outheaders,"o1"    ,0         );
  sf_putstring(outheaders,"label1","none"    );
  sf_putstring(outheaders,"unit1" ,"none"    );
  
  dim_output=1;
  for (iaxis=1; iaxis<SF_MAX_DIM; iaxis++){
    label_argparmread=n_argparmread=o_argparmread=d_argparmread=false;
    sprintf(parameter,"label%d",iaxis+1);
    fprintf(stderr,"try to read %s\n",parameter);
    if ((label[iaxis]=sf_getstring(parameter))) {
      /*(label#=(2,...)  name of each of the axes. 
	   label1 is not changed from input. Each label must be a 
	   header key like cdp, cdpt, or ep.  The trace header 
	   values are used to define the output trace location in
	   the output file. )*/
      fprintf(stderr,"got %s=%s\n",parameter,label[iaxis]);
      sf_putstring(output    ,parameter,label[iaxis]);
      sf_putstring(outheaders,parameter,label[iaxis]);
      label_argparmread=true;
    }
    sprintf(parameter,"n%d",iaxis+1);
    fprintf(stderr,"try to read %s\n",parameter);
    if (sf_getlargeint  (parameter,&n_output[iaxis])) {
      /*( n#=(2,...) number of locations in the #-th dimension )*/ 
      fprintf(stderr,"got %s=%lld\n",parameter,(long long) n_output[iaxis]);
      sf_putint(output    ,parameter,n_output[iaxis]);
      sf_putint(outheaders,parameter,n_output[iaxis]);
      n_argparmread=true;
    }
    sprintf(parameter,"o%d",iaxis+1);
    if (sf_getfloat(parameter,&o_output[iaxis])) {
      /*( o#=(2,...) origin of the #-th dimension )*/ 
      sf_putfloat(output    ,parameter,o_output[iaxis]);
      sf_putfloat(outheaders,parameter,o_output[iaxis]);

//.........这里部分代码省略.........
            sf_warning ("Processing V_x");
        buf2 = sf_floatalloc (n);
        sf_floatread (buf2, n, velx);
        sf_fileclose (velx);
        /* Convert to V_x^2 */
        for (i = 0; i < n; i++)
            buf2[i] *= buf2[i];
        set_multi_UBspline_3d_s (velspline, ic, buf2);
        ic++;
        /* Convert to (V_z*V_x)^2 */
        for (i = 0; i < n; i++)
            buf[i] *= buf2[i];
    }
    if (eta) {
        if (verb)
            sf_warning ("Processing Eta");
        sf_floatread (buf2, n, eta);
        sf_fileclose (eta);
        /* Convert to -8*eta/(1 + 2*eta)*(V_z*V_x)^2 */
        for (i = 0; i < n; i++) {
            buf2[i] = -8.0*buf2[i]/(1.0 + 2.0*buf2[i]);
            buf2[i] *= buf[i];
        }
        set_multi_UBspline_3d_s (velspline, ic, buf2);
        ic++;
    }
    if (theta) {
        if (verb)
            sf_warning ("Processing Theta");
        sf_floatread (buf, n, theta);
        sf_fileclose (theta);
        /* Convert to radians */
        for (i = 0; i < n; i++)
            buf[i] = buf[i]*SF_PI/180.0;
        set_multi_UBspline_3d_s (velspline, ic, buf);
        ic++;
    }
    if (phi) {
        if (verb)
            sf_warning ("Processing Phi");
        sf_floatread (buf, n, phi);
        sf_fileclose (phi);
        /* Convert to radians */
        for (i = 0; i < n; i++)
            buf[i] = buf[i]*SF_PI/180.0;
        set_multi_UBspline_3d_s (velspline, ic, buf);
        ic++;
    }

    if (buf2)
        free (buf2);
    free (buf);

    sz = (size_t)sizeof(multi_UBspline_3d_s) +
         (size_t)velspline->nc;
#ifdef HAVE_SSE
    if (sizeof(multi_UBspline_3d_s) % 64)
        sz += 64 - (sizeof(multi_UBspline_3d_s) % 64);
#endif
    /* Make output a 1-D file of coefficients */
    sf_unshiftdim2 (velz, out, 1);
    /* Set up output */
    sf_settype (out, SF_UCHAR);
    sf_putlargeint (out, "n1", sz);
    sf_putfloat (out, "o1", 0.0);
    sf_putfloat (out, "d1", 1.0);
    sf_putstring (out, "label1", "Spline coefficients");
    sf_putstring (out, "unit1", "");
    sf_putstring (out, "label2", "");
    sf_putstring (out, "unit2", "");

    sf_putint (out, "Nz", nz);
    sf_putfloat (out, "Oz", oz);
    sf_putfloat (out, "Dz", dz);
    sf_putint (out, "Nx", nx);
    sf_putfloat (out, "Ox", ox);
    sf_putfloat (out, "Dx", dx);
    sf_putint (out, "Ny", ny);
    sf_putfloat (out, "Oy", oy);
    sf_putfloat (out, "Dy", dy);

    sf_putint (out, "Nc", nc);
    sf_putstring (out, "splines", "y");

    if (verb) {
        sf_warning ("Number of spline coefficients: %lu",
                    velspline->nc/(size_t)sizeof(float));
        sf_warning ("Writing spline coefficients");
    }

    sf_ucharwrite ((unsigned char*)velspline, (size_t)sizeof(multi_UBspline_3d_s), out);
#ifdef HAVE_SSE
    if (sizeof(multi_UBspline_3d_s) % 64)
        sf_ucharwrite (pad, (size_t)(64 - (sizeof(multi_UBspline_3d_s) % 64)), out);
#endif
    sf_ucharwrite ((unsigned char*)velspline->coefs, (size_t)velspline->nc, out);
    destroy_Bspline (velspline);

    return 0;
}

int main (int argc, char *argv[])
{
    int n123, niter, order, nj1,nj2, i, j, liter, dim;
    int n[SF_MAX_DIM], rect[3], nr, ir; 
    float p0, *u, *p, pmin, pmax, eps;
    float **allu, **allp, d1, d2, d3, o1, o2, o3, *sendbuf, *recvbuf;
    bool verb, **mm;
    sf_file in, out, mask, dip0;

    int cpuid, numprocs, nrpad, iturn;
    MPI_Comm comm=MPI_COMM_WORLD;

    sf_init(argc,argv);

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(comm, &cpuid);
    MPI_Comm_size(comm, &numprocs);

    in = sf_input ("--input");
    out = sf_output ("--output");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float type");
    if (cpuid==0) sf_warning("numprocs=%d", numprocs);

    dim = sf_filedims(in,n);
    if (dim < 2) n[1]=1;
    n123 = n[0]*n[1];
    nr = 1;
    for (j=2; j < dim; j++) {
        nr *= n[j];
    }

    n[2]= 1;
    rect[2]=1;
    nj2=1;

    if(nr%numprocs==0)
        nrpad=nr;
    else
        nrpad=(nr/numprocs+1)*numprocs;

    /* set up output dimension */
    if(cpuid==0){
        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_histfloat(in, "d2", &d2)) sf_error("No d2= in input");
        if (!sf_histfloat(in, "o2", &o2)) sf_error("No o2= in input");
        if (!sf_histfloat(in, "d3", &d3)) sf_error("No d3= in input");
        if (!sf_histfloat(in, "o3", &o3)) sf_error("No o3= in input");

        sf_putint(out, "n1", n[0]);
        sf_putfloat(out, "d1", d1);
        sf_putfloat(out, "o1", o1);
        sf_putstring(out, "label1", "Depth");
        sf_putstring(out, "unit1", "m");
        sf_putint(out, "n2", n[1]);
        sf_putfloat(out, "d2", d2);
        sf_putfloat(out, "o2", o2);
        sf_putstring(out, "label2", "Offset");
        sf_putstring(out, "unit2", "m");
        sf_putint(out, "n3", nr);
        sf_putfloat(out, "d3", d3);
        sf_putfloat(out, "o3", o3);
        sf_putstring(out, "label3", "CIGs");
        sf_putstring(out, "unit3", "m");
    }

    if (!sf_getint("niter",&niter)) niter=5;
    /* number of iterations */
    if (!sf_getint("liter",&liter)) liter=20;
    /* number of linear iterations */

    if (!sf_getint("rect1",&rect[0])) rect[0]=1;
    /* dip smoothness on 1st axis */
    if (!sf_getint("rect2",&rect[1])) rect[1]=1;
    /* dip smoothness on 2nd axis */

    if (!sf_getfloat("p0",&p0)) p0=0.;
    /* initial dip */

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

    if (!sf_getbool("verb",&verb)) verb = false;
    /* verbosity flag */
    if (!sf_getfloat("pmin",&pmin)) pmin = -FLT_MAX;
    /* minimum dip */
    if (!sf_getfloat("pmax",&pmax)) pmax = +FLT_MAX;
    /* maximum dip */

    if (!sf_getfloat("eps",&eps)) eps=0.0f;
    /* regularization */

    /* initialize dip estimation */
    dip3_init(n[0], n[1], n[2], rect, liter, eps, false);

    /* initial dip file */
    if(NULL != sf_getstring("dip0")){
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
    int n1, n2, ninf, i1, i2, i, *inter2;
    char *label, *unit;
    float o1, d1, o2, d2, x, z;
    float *v0, *dvdx, *dvdz, *x0, *z0, *trace, **inter;
    sf_file model, surface;

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

    if (SF_FLOAT != sf_gettype(surface)) sf_error("Need float input");

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

    sf_shiftdim(surface, model, 1);
    sf_putint(model,"n3",1);

    if (!sf_histint(surface,"n2",&ninf)) ninf=1; 

    if (!sf_getint("n1",&n1)) sf_error("Need n1=");
    /* Number of samples on the depth axis */
    if (!sf_getfloat("d1",&d1)) sf_error("Need d1=");
    /* Sampling of the depth axis */
    if (!sf_getfloat("o1",&o1)) o1=0.;
    /* Origin of the depth axis */

    sf_putint(model,"n1",n1);
    sf_putfloat(model,"d1",d1);
    sf_putfloat(model,"o1",o1);

    if (NULL == (label = sf_getstring("label1"))) label="Depth";
    /* depth axis label */
    sf_putstring(model,"label1",label);

    if (NULL != (unit = sf_getstring("unit1"))) /* depth axis unit */
	sf_putstring(model,"unit1",unit);

    inter = sf_floatalloc2(n2,ninf);
    inter2 = sf_intalloc(ninf);
    sf_floatread(inter[0],n2*ninf,surface);

    ninf++; /* more layers than interfaces */
    v0 = sf_floatalloc(ninf);
    x0 = sf_floatalloc(ninf);
    z0 = sf_floatalloc(ninf);
    dvdx = sf_floatalloc(ninf);
    dvdz = sf_floatalloc(ninf);

    /* Input layer velocities and velocity derivatives */
    if (!sf_getfloats("x0",x0,ninf)) 
	for(i=0;i< ninf;i++) x0[i] = 0.;
    if (!sf_getfloats("z0",z0,ninf))
	for(i=0;i< ninf;i++) z0[i] = 0.;
    if (!sf_getfloats("v00",v0,ninf))
	for(i=0;i< ninf;i++) v0[i] = 1500.+ 500*i;
    if (!sf_getfloats("dvdx",dvdx,ninf)) 
	for(i=0;i< ninf;i++) dvdx[i] = 0.;
    if (!sf_getfloats("dvdz",dvdz,ninf)) 
	for(i=0;i< ninf;i++) dvdz[i] = 0.;

    trace = sf_floatalloc(n1);

    /* compute linear velocity */
    for(i2=0; i2 < n2; i2++) { 
	x = o2+i2*d2;
	for (i=0; i < ninf-1; i++) {
	    inter2[i] = floorf(0.5+(inter[i][i2]-o1)/d1);
	}
	for(i1=0; i1 < n1; i1++) {
	    z = o1+i1*d1;
	    for (i=0; i < ninf-1; i++) {
		if (i1 < inter2[i]) break;
	    }
	    trace[i1] = v0[i] + (x-x0[i])*dvdx[i] + (z-z0[i])*dvdz[i];
	}
	sf_floatwrite(trace,n1,model);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int n[SF_MAX_DIM], a[SF_MAX_DIM], center[SF_MAX_DIM], gap[SF_MAX_DIM];
    int *pch, *nh, dim, n123, nf, i, niter, nbf, nbp, id, ip, ig, np;
    int *kk, *pp;
    float *dd, eps, dabs, di;
    nfilter aa, bb;
    char varname[6], *lagfile;
    sf_file in, flt, lag, mask, patch, reg;

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

    dim = sf_filedims(in,n);

    if (NULL == (lagfile = sf_getstring("lag"))) sf_error("Need lag=");
    /* output file for filter lags */

    lag = sf_output(lagfile);
    sf_settype(lag,SF_INT);

    sf_putstring(flt,"lag",lagfile);

    sf_putints(lag,"n",n,dim);

    if (!sf_getints("a",a,dim)) sf_error("Need a=");

    if (!sf_getints("center",center,dim)) {
        for (i=0; i < dim; i++) {
            center[i] = (i+1 < dim && a[i+1] > 1)? a[i]/2: 0;
        }
    }

    if (!sf_getints("gap",gap,dim)) {
        for (i=0; i < dim; i++) {
            gap[i] = 0;
        }
    }

    n123 = 1;
    for (i=0; i < dim; i++) {
        n123 *= n[i];
    }

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

    if (NULL != sf_getstring("maskin")) {
        /* optional input mask file */
        mask = sf_input("maskin");

        switch (sf_gettype(mask)) {
        case SF_INT:
            sf_intread (kk,n123,mask);
            break;
        case SF_FLOAT:
            sf_floatread (dd,n123,mask);
            for (i=0; i < n123; i++) {
                kk[i] = (dd[i] != 0.);
            }
            break;
        default:
            sf_error ("Wrong data type in maskin");
            break;
        }

        sf_fileclose (mask);
    } else {
        for (i=0; i < n123; i++) {
            kk[i] = 1;
        }
    }

    sf_floatread(dd,n123,in);

    dabs = fabsf(dd[0]);
    for (i=1; i < n123; i++) {
        di = fabsf(dd[i]);
        if (di > dabs) dabs=di;
    }

    random_init(2004);
    for (i=0; i < n123; i++) {
        dd[i] = dd[i]/dabs+ 100.*FLT_EPSILON*(random0()-0.5);;
    }

    pp = sf_intalloc(n123);
    if (NULL != sf_getstring("pch")) {
        patch = sf_input("pch");
        if (SF_INT != sf_gettype(patch)) sf_error("Need int pch");

        sf_intread(pp,n123,patch);

        np = pp[0];
        for (i=1; i < n123; i++) {
            if (pp[i] > np) np = pp[i];
        }

        sf_fileclose(patch);
//.........这里部分代码省略.........