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

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

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

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

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

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

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

    sf_putint (out, "n1", ESC3_NUM);
    sf_putfloat (out, "o1", 0.0);
    sf_putfloat (out, "d1", 1.0);
    sf_putstring (out, "label1", "Escape variable");
    sf_putstring (out, "unit1", "");
    sf_putint (out, "n2", nb);
    sf_putfloat (out, "d2", db*180.0/SF_PI);
    sf_putfloat (out, "o2", ob*180.0/SF_PI);
    sf_putstring (out, "label2", "Inclination");
    sf_putstring (out, "unit2", "Degrees");
    sf_putint (out, "n3", na);
    sf_putfloat (out, "d3", da*180.0/SF_PI);
    sf_putfloat (out, "o3", oa*180.0/SF_PI);
    sf_putstring (out, "label3", "Azimuth");
    sf_putstring (out, "unit3", "Degrees");

    sf_putint (out, "n4", nz);
    sf_putfloat (out, "o4", oz);
    sf_putfloat (out, "d4", dz);
    if (!spdom) {
        sf_putstring (out, "label4", "Depth");
        sf_putstring (out, "unit4", "");
    }
    sf_putint (out, "n5", nx);
    sf_putfloat (out, "o5", ox);
    sf_putfloat (out, "d5", dx);
    if (!spdom) {
        sf_putstring (out, "label5", "X");
        sf_putstring (out, "unit5", "");
    }
    sf_putint (out, "n6", ny);
    sf_putfloat (out, "o6", oy);
    sf_putfloat (out, "d6", dy);
    if (!spdom) {
        sf_putstring (out, "label6", "Y");
        sf_putstring (out, "unit6", "");

int main(int argc, char* argv[])
{
    bool verb,pas,adj,abc;          /* execution flags */
    int ix, iz, it;                 /* index variables */
    int nt, nx, nz, depth, nzxpad, nb, n2, snap;
    float ox, oz, dx, dz, dt, dt2, idz2, idx2, cb;

    int nxpad, nzpad;
    float **vvpad;
    
    float **dd, **mm, **vv, ***ww;
    float **u0, **u1, **u2, **tmp;  /* temporary arrays */

    sf_file in, out, vel, wave;     /* I/O files */

    /* initialize Madagascar */
    sf_init(argc,argv);
    
    /* initialize OpenMP support */
#ifdef _OPENMP
    omp_init();
#endif

    if(!sf_getbool("verb", &verb)) verb=false;
    /* verbosity flag */
    if(!sf_getbool("adj", &adj)) adj=false;
    /* adjoint flag, 0: modeling, 1: migration */
    if(!sf_getbool("pas", &pas)) pas=false;
    /* passive flag, 0: exploding reflector rtm, 1: passive seismic imaging */
    if(!sf_getbool("abc",&abc)) abc = false;
    /* absorbing boundary condition */
    if(!sf_getint("snap", &snap)) snap=0;
    /* wavefield snapshot flag */
    if(!sf_getint("depth", &depth)) depth=0;
    /* surface */
	
    /* setup I/O files */
    in  = sf_input("in");
    out = sf_output("out");
    vel = sf_input("velocity");
    /* velocity model */
    
    /* Dimensions */
    if(!sf_histint  (vel, "n1", &nz)) sf_error("No n1= in velocity");
    if(!sf_histint  (vel, "n2", &nx)) sf_error("No n2= in velocity");
    if(!sf_histfloat(vel, "o1", &oz)) sf_error("No o1= in velocity");
    if(!sf_histfloat(vel, "o2", &ox)) sf_error("No o2= in velocity");
    if(!sf_histfloat(vel, "d1", &dz)) sf_error("No d1= in velocity");
    if(!sf_histfloat(vel, "d2", &dx)) sf_error("No d2= in velocity");

    if(adj){ /* migration */
        if(!sf_histint(in, "n1", &nt)) sf_error("No n1= in data");
        if(!sf_histfloat(in, "d1", &dt)) sf_error("No d1= in data");
        if(!sf_histint(in, "n2", &n2) || n2!=nx) sf_error("Need n2=%d in data", nx);

        sf_putint   (out, "n1", nz);
        sf_putfloat (out, "o1", oz);
        sf_putfloat (out, "d1", dz);
        sf_putstring(out, "label1", "Depth");
        sf_putstring(out, "unit1" , "km");
        sf_putint   (out, "n2", nx);
        sf_putfloat (out, "o2", ox);
        sf_putfloat (out, "d2", dx);
        sf_putstring(out, "label2", "Distance");
        sf_putstring(out, "unit2" , "km");
        if (pas) {
            sf_putint   (out, "n3", nt);
            sf_putfloat (out, "d3", dt);
            sf_putfloat (out, "o3", 0.0f);
            sf_putstring(out, "label3", "Time");
            sf_putstring(out, "unit3" , "s");
        }
    }else{ /* modeling */
        if(!sf_getint("nt", &nt)) sf_error("Need nt=");
        if(!sf_getfloat("dt", &dt)) sf_error("Need dt=");
        
        sf_putint   (out, "n1", nt);
        sf_putfloat (out, "d1", dt);
        sf_putfloat (out, "o1", 0.0);
        sf_putstring(out, "label1", "Time");
        sf_putstring(out, "unit1" , "s");
        sf_putint   (out, "n2", nx);
        sf_putfloat (out, "o2", ox);
        sf_putfloat (out, "d2", dx);
        sf_putstring(out, "label2", "Distance");
        sf_putstring(out, "unit2" , "km");
        if (pas) {
            sf_putint   (out, "n3", 1);
        }
    }
    
    /* dimension of padded boundary */
    if(!sf_getint("nb", &nb) || nb<NOP) nb = NOP;
    if(!sf_getfloat("cb", &cb)) cb = 0.0f;
    nxpad = nx+2*nb;
    nzpad = nz+2*nb;
    nzxpad = nzpad*nxpad;
    depth = depth+nb;

    /* set Laplacian coefficients */
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
	int   i;
        float *data;
        char *fn;

        sf_init(argc,argv);

        FILE *Fi;
        sf_file Fo;
        
        cjbsegy *tr;
 
        tr = calloc(sizeof(cjbsegy), 1);

        int nx, ny, nazim, nang, ntau;
        float dazim, dang, dtau, fazim=0.0, fang=0.0, ftau=0;
        int dx, dy, fx=0, fy=0;

        if (!sf_getint("nx",&nx)) nx=101;
        if (!sf_getint("ny",&ny)) ny=101;
        if (!sf_getint("nazim",&nazim)) nazim=8;
        if (!sf_getint("nang",&nang)) nang=21;
        if (!sf_getint("ntau",&ntau)) ntau=101;
        if (!sf_getint("dx",&dx)) dx=1;
        if (!sf_getint("dy",&dy)) dy=1;
        if (!sf_getfloat("dazim",&dazim)) dazim=22.5;
        if (!sf_getfloat("dang",&dang)) dang=2.0;
        if (!sf_getfloat("dtau",&dtau)) dtau=0.002;
        if (!sf_getint("fx",&fx)) fx=1;
        if (!sf_getint("fy",&fy)) fy=1;
        if (!sf_getfloat("ftau",&ftau)) ftau=0;
        if (NULL==(fn=sf_getstring("fn"))) fn="kpstm.ladcig.su.agc";

        /* setup I/O files */
        Fo = sf_output("out");

        if((Fi=fopen(fn,"rb"))==NULL)
        {
           printf("File %s open error!\n",fn);
           exit(0);
        }

        fread(tr,sizeof(cjbsegy),1,Fi);
        int iline0=tr->ep;
        sf_warning("ns=%d dt=%f iLineNo=%d ",tr->ns, tr->dt,iline0);


        if(fseek(Fi, 0L, 2) ==-1)
          printf("input file size unknown; Please specify n2\n");
        int nxy=(int) (ftell(Fi)/((60+ntau)*sizeof(float)));
       
        sf_warning("nxy=%d nx=%d ny=%d ",nxy, nx,ny);
        sf_warning("nazim=%d nang=%d ntau=%d",nazim,nang,ntau);
        sf_warning("dx=%d dy=%d dazim=%f dang=%f dtau=%f",dx,dy,dazim,dang,dtau);
        sf_warning("fx=%d fy=%d fazim=%f fang=%f ftau=%f",fx,fy,fazim,fang,ftau);

        if(nxy!=nx*ny*nazim*nang) {
          sf_warning("nx * ny * nazim * nang != nxy ");
          exit(0);  
         };

        sf_putint(Fo,"n1",ntau);
        sf_putint(Fo,"n2",nang);
        sf_putint(Fo,"n3",nazim);
        sf_putint(Fo,"n4",nx);
        sf_putint(Fo,"n5",ny);
        sf_putfloat(Fo,"d1",dtau);
        sf_putfloat(Fo,"o1",ftau);
        sf_putfloat(Fo,"d2",dang);
        sf_putfloat(Fo,"o2",fang);
        sf_putfloat(Fo,"d3",dazim);
        sf_putfloat(Fo,"o3",fazim);
        sf_putfloat(Fo,"d4",dx);
        sf_putfloat(Fo,"o4",fx);
        sf_putfloat(Fo,"d5",dy);
        sf_putfloat(Fo,"o5",fy);
        sf_putstring(Fo,"label1","z");
        sf_putstring(Fo,"label2","angle");
        sf_putstring(Fo,"label3","azimuth");
        sf_putstring(Fo,"label4","x");
        sf_putstring(Fo,"label5","y");
        sf_putstring(Fo,"unit1","ms");
        sf_putstring(Fo,"unit2","degree");
        sf_putstring(Fo,"unit3","degree");
        sf_putstring(Fo,"unit4","m");
        sf_putstring(Fo,"unit5","m");

        data = sf_floatalloc(ntau);

        rewind(Fi);
        for(i=0;;i++)
        {
          fread(tr,sizeof(cjbsegy),1,Fi);
          if(tr->ep != iline0){
            sf_warning("Read iLineNo=%d finished",iline0);
            iline0=tr->ep;
          }
          fread(data,sizeof(float),ntau,Fi);
          if(feof(Fi))break;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool adj,timer,verb,gmres;
    int nt, nx, nz, nx2, nz2, nzx, nzx2, ntx, pad1, snap, gpz, wfnt, i;
    int m2, n2, nk, nth=1;
    int niter, mem;
    float dt, dx, dz, ox;
    sf_complex *img, *imgout, *dat, **lt1, **rt1, **lt2, **rt2, ***wvfld;
    sf_file data, image, leftf, rightf, leftb, rightb, snaps;
    double time=0.,t0=0.,t1=0.;
    geopar geop;

    sf_init(argc,argv);

    /* essentially doing imaging */
    adj = true;

    if (!sf_getbool("gmres",&gmres)) gmres=false;
    if (gmres) {
      if (!sf_getint("niter",&niter)) niter=10;
      if (!sf_getint("mem",&mem)) mem=20;
    }

    if(! sf_getbool("timer",&timer)) timer=false;
    if (!sf_getbool("verb",&verb)) verb=false;
    if (!sf_getint("snap",&snap)) snap=0;
    /* interval for snapshots */
    if (!sf_getint("pad1",&pad1)) pad1=1;
    /* padding factor on the first axis */
    if(!sf_getint("gpz",&gpz)) gpz=0;
    /* geophone surface */

    /* adj */
    if (!sf_getint("nz",&nz)) sf_error("Need nz=");
    /* depth samples */
    if (!sf_getfloat("dz",&dz)) sf_error("Need dz=");
    /* depth sampling */

    /* for */
    if (!sf_getint("nt",&nt)) sf_error("N