//.........这里部分代码省略.........
			/* load source */
#ifdef _OPENMP 
#pragma omp parallel for \
			private(ix,iz) \
			shared(term,rr,padnx,padnz,it)
#endif
			for(ix=0; ix<padnx; ix++){
				for(iz=0; iz<padnz; iz++){
					term[ix][iz] += rr[ix*padnz+iz]*array->ww[it];
				}
			}

			/* update */
#ifdef _OPENMP 
#pragma omp parallel for \
			private(ix,iz) \
			shared(p0,p1,p2,vv,term,padnx,padnz,dt2)
#endif
			for(ix=4; ix<padnx-4; ix++){
				for(iz=4; iz<padnz-4; iz++){
					p2[ix][iz]=2*p1[ix][iz]-p0[ix][iz]+vv[ix][iz]*vv[ix][iz]*dt2*term[ix][iz];
				}
			}
			
			/* swap wavefield pointer of different time steps */
			tmparray=p0; p0=p1; p1=p2; p2=tmparray;

			/* boundary condition */
			apply_sponge(p0, acpar->bc, padnx, padnz, nb);
			apply_sponge(p1, acpar->bc, padnx, padnz, nb);
		} // end of time loop

		/* check */
		if(wit != wnt) sf_error("Incorrect number of wavefield snapshots");
		wit--;

		/* read data */
		sf_seek(Fdat, is*nr*nt*sizeof(float), SEEK_SET);
		sf_floatread(dd[0], nr*nt, Fdat);

		/* initialization */
		memset(p0[0], 0., padnzx*sizeof(float));
		memset(p1[0], 0., padnzx*sizeof(float));
		memset(p2[0], 0., padnzx*sizeof(float));
		memset(term[0], 0., padnzx*sizeof(float));
		
		/* backward propagation */
		for(it=nt-1; it>=0; it--){
			if(verb) sf_warning("Backward propagation is=%d; it=%d;", is+1, it);

			/* laplacian operator */
			laplace(p1, term, padnx, padnz, dx2, dz2);
			
			/* load data */
			for(ir=0; ir<acpar->nr2[is]; ir++){
				rx=acpar->r0_v[is]+ir*acpar->dr_v;
				term[rx][rz] += dd[acpar->r02[is]+ir][it];
			}

			/* update */
#ifdef _OPENMP 
#pragma omp parallel for \
			private(ix,iz) \
			shared(p0,p1,p2,vv,term,padnx,padnz,dt2)
#endif
			for(ix=4; ix<padnx-4; ix++){

int main(int argc, char* argv[])
{
    int n1, n2, n3, i2, i3, iw, n, dim;
    float o1, d1, m1, a1, m2, a2, f;
    sf_complex *data, *data2=NULL;
    char key[5];
    sf_file in, in2, ma, ma2, out, out2=NULL;

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

    if (NULL != sf_getstring("in2")) {
	/* optional second input file */
	in2 = sf_input("in2");
    } else {
	in2 = NULL;
    }
    
    ma = sf_input("ma");
    ma2 = sf_input("ma2");

    out = sf_output("out");

    if (NULL != in2) out2 = sf_output("out2");

    if ((SF_COMPLEX != sf_gettype (in)) ||
	(NULL != in2 && 
	 SF_COMPLEX != sf_gettype (in2))) sf_error("Need complex data");

    if (!sf_histint(in,"n1",&n1)) n1=1;

    if (NULL != in2 && sf_histint(in2,"n1",&n) && n != n1)
	sf_error("Size mismatch in in2: %d != %d",n,n1);

    if (!sf_getint("dim",&dim)) dim=1;
    /* data dimensionality */
    sprintf(key,"n%d",dim+1);

    if (!sf_histint(in,key,&n3)) n3=1;

    if (NULL != in2 && sf_histint(in2,key,&n) && n != n3)
	sf_error("Size mismatch in in2 [%s]: %d != %d",key,n,n3);

    n2 = sf_leftsize(in,1);
    n2 /= n3;

    /* n3 is the number of windows, n2xn1 is the window size */

    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&o1)) o1=0.;
    
    if (!sf_histint(ma,"n1",&n) || n != 2)
	sf_error("Wrong n1= in ma");
    if (!sf_histint(ma2,"n1",&n) || n != 2)
	sf_error("Wrong n1= in ma2");
    if (!sf_histint(ma,"n2",&n) || n != n3)
	sf_error("Wrong n2= in ma");
    if (!sf_histint(ma2,"n2",&n) || n != n3)
	sf_error("Wrong n2= in ma2");
    
    data = sf_complexalloc(n1);
    if (NULL != in2) data2 = sf_complexalloc(n1);

    for (i3=0; i3 < n3; i3++) { /* loop over windows */
	sf_floatread(&m1,1,ma);
	sf_floatread(&a1,1,ma);

	sf_floatread(&m2,1,ma2);
	sf_floatread(&a2,1,ma2);
	
	for (i2=0; i2 < n2; i2++) { /* loop over traces in a window */
	    sf_complexread(data,n1,in);
	    if (NULL != in2) sf_complexread(data2,n1,in2);
	
	    /* Frequency-domain Gaussian scaling and smoothing */
	    for (iw=0; iw < n1; iw++) {
		f = iw*d1;
		f *= f;
		if (m1 > m2) {
		    f = (a2*m1)/(a1*m2)*expf(f*(1./m1-1./m2));
#ifdef SF_HAS_COMPLEX_H
		    data[iw] *= f;
#else
		    data[iw] = sf_crmul(data[iw],f);
#endif
		} else if (NULL != in2) {
		    f = (a1*m2)/(a2*m1)*expf(f*(1./m2-1./m1));
#ifdef SF_HAS_COMPLEX_H
		    data2[iw] *= f;
#else
		    data2[iw] = sf_crmul(data2[iw],f);
#endif
		}
	    }
	    
	    sf_complexwrite(data,n1,out);
	    if (NULL != in2) sf_complexwrite(data2,n1,out2);
	}
    }
    
//.........这里部分代码省略.........

int main(int argc,char *argv[])
{
    float x0,y0,x00,y00,t0,d1,d2,d3,o1,o2,o3,k11,k21,k31,k41,k12,k22,k32,k42;
    float kxi,gamma,omega,phi;/*parameter define*/
    int i,j,k,n1,n2,n3,sn1,pow1,pow2;
    double h;
    float  *inpsg,*xt,*xt1,*re;
    sf_complex *out1;
    sf_file in, out, restor=NULL;
    bool verb,cosine;
    
    sf_init (argc, argv); 
    in = sf_input("in");
    out = sf_output("out");

    if (NULL != sf_getstring ("restor")) {
        restor = sf_input("restor");
    } else {
	restor = NULL;
	sn1=1;
    }
    
   
    sf_setformat(out,"native_complex");
    
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&o1)) sf_error("No o1= in input");
    if (!sf_histint(in,"n2",&n2)) n2=1;
    if (!sf_histfloat(in,"d2",&d2)) d2=1;
    if (!sf_histfloat(in,"o2",&o2)) o2=1;
    if (!sf_histint(in,"n3",&n3)) n3=1;
    if (!sf_histfloat(in,"d3",&d3)) d3=1;
    if (!sf_histfloat(in,"o3",&o3)) o3=1;
    
    if(!sf_getfloat("gamma",&gamma)) gamma=0.75;
    /*strength of external force*/
    if(!sf_getfloat("omega",&omega)) omega=1;
    /*angular frequence of external force*/
    if(!sf_getfloat("kxi",&kxi)) kxi=1;
    /*adjustment for input signal*/
    if(!sf_getfloat("x0",&x0)) x0=0;
    /*initial value of x0*/
    if(!sf_getfloat("y0",&y0)) y0=0;
    /*initial value of y0*/
    if(!sf_getint("pow1",&pow1)) pow1=1;
    /*power of first non-linear restitution term*/
    if(!sf_getint("pow2",&pow2)) pow2=3;
    /*power of second non-linear restitution term*/
    if(!sf_getfloat("phi",&phi)) phi=0.;
    /*phase of cosine signal unit=pi*/
    if(!sf_getbool("verb",&verb)) verb = false;
    /* verbosity flag */
    if(!sf_getbool("cosine",&cosine)) cosine = true;
    /* if n need extenal input for periodic restoring force */
    
    if (NULL != restor){
	
	if (!sf_histint(restor,"n1",&sn1)) sf_error("No n1= in s");
	if (n1>sn1) sf_error("insufficient lenth of external force");
    }
    
    
    h=d1;
    /*step lenth of R-K4*/
    
    re = sf_floatalloc(sn1);
    inpsg = sf_floatalloc(n1*n2);
    xt = sf_floatalloc(n1*n2);
    xt1 = sf_floatalloc(n1*n2);
    out1 = sf_complexalloc(n1*n2);
    
    
    if (NULL != restor) {
	sf_floatread(re,sn1,restor);
    }else {
	for (i=0;i<sn1; i++)
	    re[i] = 0;
    }

    
    if(cosine) {
	for(k = 0; k < n3 ; k++) {
	    sf_floatread(inpsg,n1*n2,in);
	    x00=x0,y00=y0;
	    t0=0;
	    if(verb) sf_warning("step %d of %d;",k,n3);
	    for(j = 0;j < n2;j++) {
		x00=x0,y00=y0;
		t0=0;
		for(i = 0;i < n1;i++){
		    
		    k11 = function1(t0,x00,y00,omega);
		    k12 = function2(t0,x00,y00,inpsg[i+j*n1],kxi,
				    gamma,omega,pow1,pow2,phi);
		    k21 = function1(t0+h/2,x00+k11*h/2,y00+k12*h/2,omega);
		    k22 = function2(t0+h/2,x00+k11*h/2,
				    y00+k12*h/2,inpsg[i+j*n1],kxi,
				    gamma,omega,pow1,pow2,phi);
		    k31 = function1(t0+h/2,x00+k21*h/2,y00+k22*h/2,omega);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool verb;
    int i1, i3, n12, n3, niter, n[2], rect[2];
    float *u1, *u2, *u3, norm;
    sf_file inp, out, den;

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

    if (SF_FLOAT != sf_gettype(inp)) sf_error("Need float input");
    
    if (!sf_histint(inp,"n1",n)) sf_error("No n1= in input");
    if (!sf_histint(inp,"n2",n+1)) sf_error("No n2= in input");
    n3 = sf_leftsize(inp,2);
    n12 = n[0]*n[1];

    if (!sf_getint("rect1",&rect[0])) rect[0]=1;
    if (!sf_getint("rect2",&rect[1])) rect[1]=1;
    /* smoothing radius */

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

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

    u1 = sf_floatalloc(n12);
    u2 = sf_floatalloc(n12);
    u3 = sf_floatalloc(n12);

	runtime_init(n12*sizeof(float));
    for (i3=0; i3 < n3; i3++) 
	{
		sf_floatread(u1, n12, inp);
		sf_floatread(u2, n12, den);

    	sf_divn_init(2, n12, n, rect, niter, verb);
		/* smooth division */
		norm = 0.;
		for (i1=0; i1 < n12; i1++) 
		    norm += u2[i1] * u2[i1];
		norm = sqrtf(n12/norm);

		for (i1=0; i1 < n12; i1++) 
		{
			u1[i1] *= norm;
			u2[i1] *= norm;
		}
		sf_divn (u1, u2, u3);

		sf_floatwrite(u3, n12, out);

		sf_divn_close();
		norm = runtime(1);
		sf_warning("%d of %d, %f MB/sec;", i3, n3, norm);
	}

	free(u1);
	free(u2);
	free(u3);
    exit(0);
}

int main(int argc, char* argv[])
{
    int nt, nt2, nx, i1, i2, ch, n12, n122, fk;
    bool adj, sm, domod;
    float dt, dt2, dx, ot, ot2, ox, epst2;
    float v_1, v_2, v_3, v_4, eps, passthr;
    float * data, * output, * datat2, * outputt2, * smooth, * model;
    sf_file inp, out, pifk;
    /* smoothing variables */
    int nrep, dim, dim1, n[SF_MAX_DIM], rect[SF_MAX_DIM], s[SF_MAX_DIM], i0, i, j, nvar;
    bool diff[SF_MAX_DIM], box[SF_MAX_DIM];
    int irep;
    char key[6];
    sf_triangle tr;
    /* kirchhoff params */
    bool hd;
    int sw;
    float *vrms, v0;
    char *test;
    sf_file vel;

    //MADAGASCAR C API
    /* initialize */
    sf_init(argc,argv);
    inp = sf_input("in");
    out = sf_output("out");

    /* get dimensions from input */
    if (!sf_histint(inp,"n1",&nt)) sf_error("No n1= in inp");
    if (!sf_histint(inp,"n2",&nx)) sf_error("No n2= in inp");
    if (!sf_histfloat(inp,"d1",&dt)) sf_error("No d1= in inp");
    if (!sf_histfloat(inp,"d2",&dx)) sf_error("No d2= in inp");
    if (!sf_histfloat(inp,"o1",&ot)) sf_error("No o1= in inp");
    if (!sf_histfloat(inp,"o2",&ox)) sf_error("No o2= in inp");
    
    /* kirchhoff parameters */
    if (!sf_getbool("hd",&hd)) hd=true;
    if (!sf_getbool("domod",&domod)) domod=true;
    /* if y, apply half-derivative filter */
    if (!sf_getint("sw",&sw)) sw=0;
    /* if > 0, select a branch of the antialiasing operation */
    
    /* smoothing part - determines various params including dimension along
    which smoothing is performed */
    dim = sf_filedims (inp,n);
    dim1 = -1;
    for (i=0; i < dim; i++) {
	snprintf(key,6,"rect%d",i+1);
	if (!sf_getint(key,rect+i)) rect[i]=1;
	/*( rect#=(1,1,...) smoothing radius on #-th axis )*/ 
	if (rect[i] > 1) dim1 = i;
	snprintf(key,6,"diff%d",i+1);
	if (!sf_getbool(key,diff+i)) diff[i]=false;
	/*( diff#=(n,n,...) differentiation on #-th axis )*/
	snprintf(key,5,"box%d",i+1);
	if (!sf_getbool(key,box+i)) box[i]=false;
	/*( box#=(n,n,...) box (rather than triangle) on #-th axis )*/
    }
    
    /* creating parameters for smoothing filter */
    s[0] = 1;
    s[1] = nt;
    //s[2] = ny; for 3D case
    nvar = nt*nx; // 2D 
    // nvar = nt*nx*ny; // 3D
    
    /* to output f-k pi filter */
    if (NULL != sf_getstring("pifk")) {
        pifk = sf_output("pifk");
    } else {
        pifk = NULL;
    }    

    /* get parameters from command line */
    if (!sf_getbool("adj",&adj)) adj=false;

    /* if perform derivative filtering */

    if (!sf_getbool("sm",&sm)) sm=true;

    /* if y, do adjoint integration */
    
    if (!sf_getfloat("v_1",&v_1)) sf_error("No integration range specified"); 
    if (!sf_getfloat("v_2",&v_2)) sf_error("No integration range specified"); 
    if (!sf_getfloat("v_3",&v_3)) sf_error("No integration range specified"); 
    if (!sf_getfloat("v_4",&v_4)) sf_error("No integration range specified");  
    if (!sf_getfloat("passthr",&passthr)) passthr = 0.001; // threshold for tail elimination
    
    if (!sf_getfloat("eps",&eps)) eps = 0.001; // damper for pi
    if (!sf_getfloat("epst2",&epst2)) epst2 = 0.001; // damper for t2warp
    
    /* new axis length */
	if (!sf_getint("pad",&nt2)) nt2=nt; /* output time samples */
	
	if (!sf_getint("repeat",&nrep)) nrep=1;
    /* repeat filtering several times */
	
	n12 = nt2*nx;
    
    data = sf_floatalloc(nt*nx);
//.........这里部分代码省略.........

static void
useek_pipe_rw_double (const char * ext, SF_INFO * psfinfo_write, SF_INFO * psfinfo_read)
{	static double buffer [PIPE_TEST_LEN] ;
	static double data [PIPE_TEST_LEN] ;
	SNDFILE *outfile ;
	SNDFILE *infile_piped ;

	int k, status = 0 ;
	int pipefd [2] ;
	pid_t pida ;

	for (k = 0 ; k < PIPE_TEST_LEN ; k++)
		data [k] = PIPE_INDEX (k) ;

	/*
	** Create the pipe.
	*/
	pipe (pipefd) ;

	/*
	** Attach the write end of the pipe to be written to.
	*/
	if ((outfile = sf_open_fd (pipefd [1], SFM_WRITE, psfinfo_write, SF_TRUE)) == NULL)
	{	printf ("\n\n%s %d : unable to create unseekable pipe for write type \"%s\".\n", __func__, __LINE__, ext) ;
		printf ("\t%s\n\n", sf_strerror (outfile)) ;
		exit (1) ;
		} ;

	if (sf_error (outfile) != SF_ERR_NO_ERROR)
	{	printf ("\n\n%s %d : unable to open unseekable pipe for write type \"%s\".\n\n", __func__, __LINE__, ext) ;
		exit (1) ;
		} ;

	/*
	** Attach the read end of the pipe to be read from.
	*/
	if ((infile_piped = sf_open_fd (pipefd [0], SFM_READ, psfinfo_read, SF_TRUE)) == NULL)
	{	printf ("\n\n%s %d : unable to create unseekable pipe for read type. \"%s\".\n\n", __func__, __LINE__, ext) ;
		exit (1) ;
		} ;

	if (sf_error (infile_piped) != SF_ERR_NO_ERROR)
	{	printf ("\n\n%s %d : unable to open unseekable pipe for read type \"%s\".\n\n", __func__, __LINE__, ext) ;
		exit (1) ;
		} ;

	/* Fork a child process that will write directly into the pipe. */
	if ((pida = fork ()) == 0) /* child process */
	{	test_writef_double_or_die (outfile, 0, data, PIPE_TEST_LEN, __LINE__) ;
		exit (0) ;
		} ;

	/* In the parent process, read from the pipe and compare what is read
	** to what is written, if they match everything went as planned.
	*/
	test_readf_double_or_die (infile_piped, 0, buffer, PIPE_TEST_LEN, __LINE__) ;
	if (memcmp (buffer, data, sizeof (buffer)) != 0)
	{	printf ("\n\n%s %d : unseekable pipe test failed for file type \"%s\".\n\n", __func__, __LINE__, ext) ;
		exit (1) ;
		} ;

	/* Wait for the child process to return. */
	waitpid (pida, &status, 0) ;
	status = WEXITSTATUS (status) ;
	sf_close (outfile) ;
	sf_close (infile_piped) ;

	if (status != 0)
	{	printf ("\n\n%s %d : status of child process is %d for file type %s.\n\n", __func__, __LINE__, status, ext) ;
		exit (1) ;
		} ;

	return ;
} /* useek_pipe_rw_double */

int main (int argc, char* argv[]) 
{
    int ix, curxapp, leftShift;
   
/* Initialize RSF */
    sf_init (argc,argv);
/* Input files */
    inDags_   = sf_input("in");
    inDagsSq_ = sf_input("dataSq");

/* check that the input is float */
    if ( SF_FLOAT != sf_gettype (inDags_) )   sf_error ("Need float input: dip-angle gathers");
    /* dip-angle gathers - stacks in the scattering-angle direction */
    if ( SF_FLOAT != sf_gettype (inDagsSq_) ) sf_error ("Need float input: dip-angle gathers in squares");
    /* stacks of amplitude squares in the scattering-angle direction */

/* Output file */
    sembFile_ = sf_output("out");

/* Depth/time axis */
    if ( !sf_histint   (inDags_, "n1", &zNum_) )   sf_error ("Need n1= in input");
    if ( !sf_histfloat (inDags_, "d1", &zStep_) )  sf_error ("Need d1= in input");
    if ( !sf_histfloat (inDags_, "o1", &zStart_) ) sf_error ("Need o1= in input");
/* Dip angle axis */
    if ( !sf_histint   (inDags_, "n2", &dipNum_) )   sf_error ("Need n2= in input");
    if ( !sf_histfloat (inDags_, "d2", &dipStep_) )  sf_error ("Need d2= in input");
    if ( !sf_histfloat (inDags_, "o2", &dipStart_) ) sf_error ("Need o2= in input");
/* x axis */
    if ( !sf_histint   (inDags_, "n3", &xNum_) )     sf_error ("Need n3= in input");
    if ( !sf_histfloat (inDags_, "d3", &xStep_) )    sf_error ("Need d3= in input");
    if ( !sf_histfloat (inDags_, "o3", &xStart_) )   sf_error ("Need o3= in input");
	
    if ( !sf_getint ("xapp",    &xapp_) )    xapp_ = 1;
    /* number of CIGs in the inline-direction processed simultaneously */
    if (!xapp_) {sf_warning ("xapp value is changed to 1"); xapp_ = 1;}

    if ( !sf_getint ("dipapp",    &xdipapp_) ) xdipapp_ = 11;
    /* number of traces in the x-dip direction processed simultaneously */
    if (!xdipapp_) {sf_warning ("dipapp value is changed to 11"); xdipapp_ = 11;}

    if ( !sf_getint ("coher",   &coher_) )   coher_ = 11;
    /* height of a vertical window for semblance calculation */
    if (!coher_) {sf_warning ("coher value is changed to 1"); coher_ = 1;}

    if ( !sf_getint ("scatnum", &scatnum_) ) scatnum_ = 1;
    /* shows how many traces were stacked in the scattering angle direction; 
       if the stack was normalized use the default value 
    */ 

	makeWeight_ = true;
    if ( !sf_getfloat ("s1", &s1_) ) {s1_ = -1.f; makeWeight_ = false; }
    /* minimum semblance value */ 
    if ( !sf_getfloat ("s2", &s2_) ) {s2_ = -1.f; makeWeight_ = false; }
    /* maximum semblance value */ 
	ds_ = s2_ - s1_;

    dagSize_ = zNum_ * dipNum_;
    halfCoher_ = coher_ / 2;    /* yes - this is the integer division */	
    halfXapp_  = xapp_ / 2;     /* this is the integer division too	*/

    for (ix = 0; ix < xNum_; ++ix) {
		
		sf_warning ("CIG %d of %d;", ix + 1, xNum_);
		
		/* xapp for the currect core CIG; it can be changed by the checkBoundary () */
		curxapp = xapp_; 
		/* distance between the core gather and the left side of the aperture */
		leftShift = 0;	

		ptrToSembPanel_ = sf_floatalloc (dagSize_);
		memset (ptrToSembPanel_, 0, dagSize_ * sizeof (float));

		readBlockAroundPoint (ix, halfXapp_, &curxapp, &leftShift);
		buildFilter (curxapp, leftShift, ptrToSembPanel_);

		free (ptrToDags_);
		free (ptrToDagsSq_);

		free (ptrToData_);
		free (ptrToDataSq_);

		sf_floatwrite (ptrToSembPanel_, dagSize_, sembFile_);
		free (ptrToSembPanel_);
    }

    sf_warning (".");

    exit(0);
}

sf_file sf_output (/*@[email protected]*/ const char* tag)
/*< Create an output file structure.
  ---
  Should do output after the first call to sf_input. >*/
{
    sf_file file;
    char *headname, *dataname, *path, *name, *format;
    size_t namelen;
    extern int mkstemp (char *tmpl);
    extern off_t ftello (FILE *stream);
	
    file = (sf_file) sf_alloc(1,sizeof(*file));
	
    if (NULL == tag || 0 == strcmp(tag,"out")) {
	file->stream = stdout;
	headname = NULL;
    } else {
	headname = sf_getstring (tag);
	if (NULL == headname) {
	    namelen = strlen(tag)+1;
	    headname = sf_charalloc (namelen);
	    memcpy(headname,tag,namelen);
	}
		
	file->stream = fopen(headname,"w+");
	if (NULL == file->stream) 
        {
            free(file);
	    sf_error ("%s: Cannot write to header file %s:",__FILE__,headname);
        }
    }
	
    file->buf = NULL;
    /*    setbuf(file->stream,file->buf); */
	
    file->pars = sf_simtab_init (tabsize);
    file->head = NULL;
    file->headname = NULL;
	
    file->pipe = (bool) (-1 == ftello(file->stream));
    if (file->pipe && ESPIPE != errno) 
    {
        free(file);
	sf_error ("%s: pipe problem:",__FILE__);
    }
 
    dataname = sf_getstring("out");
    if (NULL == dataname)
	dataname = sf_getstring("--out");
	
    if (file->pipe) {
	file->dataname = sf_charalloc (7);
	memcpy(file->dataname,"stdout",7);
    } else if (NULL == dataname) {
	path = getdatapath();
	file->dataname = sf_charalloc (PATH_MAX+NAME_MAX+1);
	strcpy (file->dataname,path);
	name = file->dataname+strlen(path);
	free (path);
	if (getfilename (file->stream,name)) {
	  if(0==strcmp(name,"/dev/null")){
	      file->dataname = sf_charalloc (7);
	      memcpy(file->dataname,"stdout",7);
	    } else {
	      namelen = strlen(file->dataname);
	      file->dataname[namelen]='@';
	      file->dataname[namelen+1]='\0';
	    }
	} else { /* Invent a name */
	  /* stdout is not a pipe, not /dev/null, not a file in this 
	     directory. 
	     One way to get here is to redirect io to a file not in this 
	     directory.  For example >../myfile.rsf.  In this case getfilename
	     cannot find the file from file->stream by looking in the current 
	     directory.  The function mkstemp is used to create a unique name 
	     to contain the binary data. */
	    sprintf(name,"%sXXXXXX",sf_getprog());
	    (void) close(mkstemp(file->dataname));
	    /* (void) unlink(file->dataname); */
	    /* old code for named pipes below */
	    /*
	      if (NULL == headname &&
	      -1L == fseek(file->stream,0L,SEEK_CUR) &&
	      ESPIPE == errno && 
	      0 != mkfifo (file->dataname, S_IRUSR | S_IWUSR))
	      sf_error ("%s: Cannot make a pipe %s:",
	      __FILE__,file->dataname);
	    */
	}  
    } else {
	namelen = strlen(dataname)+1;
	file->dataname = sf_charalloc (namelen);
	memcpy(file->dataname,dataname,namelen);
	free (dataname);
    }
	
    sf_putstring(file,"in",file->dataname);    
	
    file->op = XDR_ENCODE;
	
//.........这里部分代码省略.........

void sf_fileflush (sf_file file, sf_file src)
/*< outputs parameter to a file (initially from source src)
  ---
  Prepares file for writing binary data >*/ 
{
    time_t tm;
    char line[BUFSIZ];
    /* if already flushed, do nothing */
    if (NULL == file->dataname) return;
	
    if (NULL != src && NULL != src->head) {
	rewind(src->head);
		
	while (NULL != fgets(line,BUFSIZ,src->head)) {
	    fputs(line,file->stream);
	}
    }
    
    tm = time (NULL);
    if (0 >= fprintf(file->stream,"%s\t%s\t%s:\t%[email protected]%s\t%s\n",
		     RSF_VERSION,
		     sf_getprog(),
		     sf_getcdir(),
		     sf_getuser(),
		     sf_gethost(),
		     ctime(&tm)))
	sf_error ("%s: cannot flush header:",__FILE__);
	
    switch (file->type) {
	case SF_FLOAT: 
	    switch (file->form) {
		case SF_ASCII:		    
		    sf_putstring(file,"data_format","ascii_float");
		    break;
		case SF_XDR:
		    sf_putstring(file,"data_format","xdr_float");
		    break;
		default:
		    sf_putstring(file,"data_format","native_float");
		    break;
	    }
	    break;
	case SF_COMPLEX:
	    switch (file->form) {
		case SF_ASCII:		    
		    sf_putstring(file,"data_format","ascii_complex");
		    break;
		case SF_XDR:
		    sf_putstring(file,"data_format","xdr_complex");
		    break;
		default:
		    sf_putstring(file,"data_format","native_complex");
		    break;
	    }
	    break;
	case SF_INT:
	    switch (file->form) {
		case SF_ASCII:		    
		    sf_putstring(file,"data_format","ascii_int");
		    break;
		case SF_XDR:
		    sf_putstring(file,"data_format","xdr_int");
		    break;
		default:
		    sf_putstring(file,"data_format","native_int");
		    break;
	    }
	    break;
        case SF_SHORT:
	    switch (file->form) {
		case SF_ASCII:
		    sf_putstring(file,"data_format","ascii_short");
		    break;
		case SF_XDR:
		    sf_putstring(file,"data_format","xdr_short");
		    break;
		default:
		    sf_putstring(file,"data_format","native_short");
		    break;
	    }
	    break;
	case SF_LONG:
	    switch (file->form) {
		case SF_ASCII:
		    sf_putstring(file,"data_format","ascii_long");
		    break;
		case SF_XDR:
		    sf_putstring(file,"data_format","xdr_long");
		    break;
		default:
		    sf_putstring(file,"data_format","native_long");
		    break;
	    }
	    break;
        case SF_DOUBLE:
	    switch (file->form) {
		case SF_ASCII:
		    sf_putstring(file,"data_format","ascii_double");
		    break;
		case SF_XDR:
//.........这里部分代码省略.........

static void sf_input_error(sf_file file, const char* message, const char* name)
{
    if (error) 
	sf_error ("%s: %s %s:",__FILE__,message,name);
    sf_fileclose(file);
}

sf_file sf_input (/*@[email protected]*/ const char* tag)
/*< Create an input file structure >*/
{
    int esize;
    sf_file file;
    char *filename, *format;
    size_t len;
    extern off_t ftello (FILE *stream);
	
    file = (sf_file) sf_alloc(1,sizeof(*file));
    file->dataname = NULL;
    
    if (NULL == tag || 0 == strcmp(tag,"in")) {
	file->stream = stdin;
	filename = NULL;
    } else {
	filename = sf_getstring (tag);
	if (NULL == filename) {
	    /* this option allows you to call function with 
	       sf_input("mydir/myinput.rsf");  Karl  */
	    len = strlen(tag)+1;
	    filename = sf_charalloc(len);
	    /* should change next line to strcpy or strncpy Karl */
	    memcpy(filename,tag,len);
	}
		
	file->stream = fopen(filename,"r+");
	if (NULL == file->stream) {
	    sf_input_error(file,"Cannot read input (header) file",filename);
	    return NULL;
	}
    }
	
    file->buf = NULL;
    /*    setbuf(file->stream,file->buf); */
	
    /* create a parameter table */
    file->pars = sf_simtab_init (tabsize);
    file->head = sf_tempfile(&file->headname,"w+");
	
    /* read the parameter table from the file */
    sf_simtab_input (file->pars,file->stream,file->head);
	
    if (NULL == infiles) {
	infiles = (sf_file *) sf_alloc(1,sizeof(sf_file));
	infiles[0] = NULL;
	nfile=1;
    }
	
    if (NULL == filename) {
	infiles[0] = file;
    } else {
	free (filename);
	ifile++;
	if (ifile >= nfile) {
	    /* grow array */
	    nfile *= 2;
	    infiles = (sf_file *) realloc(infiles, nfile * sizeof(sf_file));
	    if (NULL == infiles) sf_error("%s: reallocation error",__FILE__);
	}
	infiles[ifile] = file;
    }
	
    filename = sf_histstring(file,"in");
    if (NULL == filename) {
    	sf_input_error (file,"No in= in file",tag);
	return NULL;
    }
    len = strlen(filename)+1;
    file->dataname = sf_charalloc(len);
    memcpy(file->dataname,filename,len);
	
    if (0 != strcmp(filename,"stdin")) {
	file->stream = freopen(filename,"r+b",file->stream);
	if (NULL == file->stream) {
	    sf_input_error(file,"Cannot read data file",filename);
	    return NULL;
	}
    }
    free (filename);
	
    file->pipe = (bool) (-1 == ftello(file->stream));
    if (file->pipe && ESPIPE != errno) 
	sf_error ("%s: pipe problem:",__FILE__);
	
    file->op = XDR_DECODE;
	
    format = sf_histstring(file,"data_format");
    if (NULL == format) {
	if (!sf_histint(file,"esize",&esize) || 0 != esize) {
	    sf_input_error (file,"Unknown format in",tag);
	    return NULL;
	}
	sf_setformat(file,"ascii_float");
    } else {    
	sf_setformat(file,format);
	free (format);
    }
	
    return file;
//.........这里部分代码省略.........

void sf_putline (sf_file file, const char* line)
/*< put a string line to a file >*/
{
    if (0 >= fprintf(file->stream,"\t%s\n",line))
	sf_error ("%s: cannot put line '%s':",__FILE__,line);
}

int main(int argc, char* argv[])
{
    int nt, ns, nx, n2, it, nstr, ifix;
    float dt, t0, eps, v0, *v=NULL, *trace=NULL, *out=NULL, *ww=NULL, wsum;
    char buffer[20];
    map4 stolt;
    fint1 istolt;
    bool inv;
    sf_file in=NULL, st=NULL, vel=NULL;

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

    if (!sf_getbool ("inv", &inv)) inv=false;
    /* if y, inverse stretch */

    if (!sf_getint("nstretch",&nstr)) nstr=1;
    /* number of steps */

    if (inv) {
        if (!sf_histint(in,"n1",&ns)) sf_error("No n1= in input");
        if (!sf_histint(vel,"n1",&nt)) sf_error("No n1= in velocity");
        sf_putint(st,"n1",nt);
    } else {
        if (!sf_histint(in,"n1",&nt)) sf_error("No n1= in input");
        if (!sf_getint("pad",&ns)) ns=nt;
        /* time axis padding */
        sf_putint(st,"n1",ns);
        sf_putint(st,"nstretch", nstr);
    }
    nx = sf_leftsize(in,1);
    n2 = sf_leftsize(vel,1);

    if (n2 != 1 && n2 != nx) sf_error("Wrong number of traces in velocity");

    if (!sf_histfloat(in,"d1",&dt)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&t0) || 0. != t0) sf_error("Need o1=0 in input");
    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* stretch regularization */
    if (!sf_getfloat("vel",&v0)) sf_error("Need vel=");
    /* reference velocity */

    trace = sf_floatalloc(nt);
    v = sf_floatalloc(nt);
    str = sf_floatalloc2(nt,nx);
    out = sf_floatalloc(ns);
    ww = sf_floatalloc(nx);

    if (inv) {
        stolt = NULL;
        istolt = fint1_init(4,ns,0);
    } else {
        stolt =  stretch4_init (ns,t0,dt,nt,eps);
        istolt = NULL;
    }

    wsum = 0.;
    for (ix=0; ix < nx; ix++) {
        if (0==ix || n2 > 1) {
            sf_floatread(v,nt,vel);
            ww[ix] = vt2w (nt,v,str[ix]);
        } else if (1 == n2) {
            ww[ix] = ww[0];
            for (it=0; it < nt; it++) {
                str[ix][it] = str[0][it];
            }
        }
        wsum += ww[ix];
    }

    sf_warning("%d traces processed",nx);

    if (1==nstr) {
        sf_putfloat(st,"stretch", wsum/nx);
    } else {
        for (ix=0; ix < nstr; ix++) {
            ifix=SF_MAX(1,(int) (nx*ix/(nstr-1)));
            sf_warning("%d %d",ix,ifix);

            snprintf(buffer,20,"stretch%d",ix);
            sf_putfloat (st,buffer,ww[ifix]);

            snprintf(buffer,20,"node%d",ix);
            sf_putint (st,buffer,ifix);
        }
    }

    for (ix=0; ix < nx; ix++) {
        for (it=0; it < nt; it++) {
            str[ix][it] *= dt/v0;
        }
        if (inv) {
            sf_floatread(out,ns,in);
            fint1_set(istolt,out);
            stretch(istolt, map, ns, dt, t0, nt, dt, t0, trace, 0.);
            sf_floatwrite (trace,nt,st);
        } else {
            sf_floatread (trace,nt,in);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (!sf_histint(infile,parameter,&n_in[iaxis])){
      n_in[iaxis]=1;
    }
    if(verbose>1){
      fprintf(stderr,"on axis=%d got label,o,d,n=%s,%f,%f,%d \n",iaxis,
	      label_in[iaxis],o_in[iaxis],d_in[iaxis],n_in[iaxis]);
    }
  }
  /* compute index_of_keys after calling segy_init */ 






  if(verbose>0)
    fprintf(stderr,"set up output file for tah - should be stdout\n");
  out = sf_output ("out");

  if(verbose>0)fprintf(stderr,"read name of input headers file\n");
  headers_filename=sf_getstring("headers");
  /* \n
     Trace header file name.  Default is the input data file
     name, with the final .rsf changed to _hdr.rsf 
  */

  if(headers_filename==NULL){
    /* compute headers_filename from infile_filename by replacing the final
       .rsf with _hdr.rsf */
    if(!(0==strcmp(infile_filename+strlen(infile_filename)-4,".rsf"))){
	fprintf(stderr,"parameter input, the name of the input file, does\n");
	fprintf(stderr,"not end with .rsf, so header filename cannot be\n");
	fprintf(stderr,"computed by replacing the final .rsf with _hdr.rsf.\n");
	sf_error("default for headers parameter cannot be computed.");
    }
    headers_filename=malloc(strlen(infile_filename)+60);
    strcpy(headers_filename,infile_filename);
    strcpy(headers_filename+strlen(infile_filename)-4,"_hdr.rsf\0");
    if(verbose>2){
      fprintf(stderr,"parameter header defaulted.  Computed to be #%s#\n",
	      headers_filename);
    }
  }
  if(verbose>2)fprintf(stderr,"parameter header input or computed  #%s#\n",
		       headers_filename);

  if(!sf_getbool("makeheader",&makeheader))makeheader=false;
  /* \n
     Option to load headers using the input file axis labels.  If axis 
     label2 through label9 match a header key then that coordinate is
     loaded to the traces header.  This can be used to load the source
     coordinate to the sx header location.  This may require changing
     the axis label because Madagascar axis labels are not the same as
     segy trace headers.  for example axis 2 coordiante can be loaded in
     sx trace header by:
        <spike.rsf sfput label2=sx \\
           | sftahread headers=spike_hdr.rsf makeheader=y \\ 
           | sftahgethw key=sx >/dev/null
 */
  
  inheaders = sf_input(headers_filename);

  if (!sf_histint(infile,"n1",&n1_traces))
    sf_error("input file does not define n1");
  if (!sf_histint(inheaders,"n1",&n1_headers)) 
    sf_error("input headers file does not define n1");

int main(int argc, char* argv[]) 
{
    bool sub;
    int nx, nx2, nk, nt, m, ix, ik, it, im, n2;
    sf_complex *curr, **wave, c;
    float *rcurr, dt;
    sf_complex **lft, **rht, *cwave, *cwavem, **mat;
    sf_file Fw, Fo, right, left, prop;

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

    if (SF_COMPLEX != sf_gettype(Fw)) sf_error("Need complex input");

    sf_settype(Fo,SF_FLOAT);

    /* Read/Write axes */

    if (!sf_getbool("sub",&sub)) sub=false;
    /* if -1 is included in the matrix */

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

    if (!sf_getint("nt",&nt)) sf_error("No nt= in input");
    if (!sf_getfloat("dt",&dt)) sf_error("No dt= in input");

    sf_putint(Fo,"n2",nt);
    sf_putfloat(Fo,"d2",dt);
    sf_putfloat(Fo,"o2",0.);
    sf_putstring(Fo,"label2","Time");
    sf_putstring(Fo,"unit2","s");

    nk = cfft1_init(nx,&nx2);
    sf_warning("nk=%d\n", nk);

    curr = sf_complexalloc(nx2);
    rcurr= sf_floatalloc(nx2);
    cwave = sf_complexalloc(nk);
 
    if (NULL != sf_getstring("right")) {
	left = sf_input("left");   /* Left matrix */
	right = sf_input("right"); /* Right matrix */
	
	if (SF_COMPLEX != sf_gettype(left) ||
	    SF_COMPLEX != sf_gettype(right)) sf_error("Need complex left and right");
	
	if (!sf_histint(left,"n1",&n2) || n2 != nx) sf_error("Need n1=%d in left",nx);
	if (!sf_histint(left,"n2",&m)) sf_error("Need n2= in left");
	
	if (!sf_histint(right,"n1",&n2) || n2 != m) sf_error("Need n1=%d in right",m);
	if (!sf_histint(right,"n2",&n2) || n2 != nk) sf_error("Need n2=%d in right, now n2=%d",nk,n2);
	
	lft = sf_complexalloc2(nx,m);
	rht = sf_complexalloc2(m,nk);
	
	sf_complexread(lft[0],nx*m,left);
	sf_complexread(rht[0],nk*m,right);
	
	sf_fileclose(left);
	sf_fileclose(right);
	
	cwavem = sf_complexalloc(nk);
	wave = sf_complexalloc2(nx2,m);

	mat = NULL;
    } else { /* Use propagator matrix (for testing) */
	 prop = sf_input("prop");
	 if (SF_COMPLEX != sf_gettype(prop)) sf_error("Need complex prop");

	lft = NULL;
	rht = NULL;

	mat = sf_complexalloc2(nk,nx);
	sf_complexread(mat[0],nx*nk,prop);
	sf_fileclose(prop);

	cwavem = NULL;
	wave = NULL;
    }
	
    sf_complexread (curr,nx,Fw);

    for (ix = nx; ix < nx2; ix++) {
	curr[ix] = sf_cmplx(0.,0.);
    }

    for (ix = 0; ix < nx2; ix++) {
	rcurr[ix]=0.;
    }

    /* propagation in time */
    for (it=0; it < nt; it++) {
	sf_warning("it=%d;",it);

        /* FFT: curr -> cwave */
	cfft1(curr,cwave);

	if (NULL == mat) {
	    for (im = 0; im < m; im++) {
//.........这里部分代码省略.........

int main (int argc, char *argv[])
{
    bool left;
    int ir, nr, n1,n2,n3, m1, m2, m3, n12, nw, nj1, nj2, i3;
    float *u1, *u2, *p;
    sf_file in, out, dip;
    allpass ap;

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

    if (SF_FLOAT != sf_gettype(in) ||
	SF_FLOAT != sf_gettype(dip)) sf_error("Need float type");

    if (!sf_histint(in,"n1",&n1)) sf_error("Need n1= in input");
    if (!sf_histint(in,"n2",&n2)) n2=1;
    if (!sf_histint(in,"n3",&n3)) n3=1;
    n12 = n1*n2;
    nr = sf_leftsize(in,3);

    if (!sf_getbool("left",&left)) left=true;
    /* if using left or right side of PWD */

    if (!sf_histint(dip,"n1",&m1) || m1 != n1) 
	sf_error("Need n1=%d in dip",n1);
    if (1 != n2 && (!sf_histint(dip,"n2",&m2) || m2 != n2)) 
	sf_error("Need n2=%d in dip",n2);
    if (1 != n3 && (!sf_histint(dip,"n3",&m3) || m3 != n3)) 
	sf_error("Need n3=%d in dip",n3);

    if (!sf_getint("order",&nw)) nw=1;
    /* accuracy */
    if (!sf_getint("nj1",&nj1)) nj1=1;
    /* in-line aliasing */
    if (!sf_getint("nj2",&nj2)) nj2=1;
    /* cross-line aliasing */

    for (ir=0; ir < nr; ir++) {
	u1 = sf_floatalloc(n12);
	u2 = sf_floatalloc(n12);
	p  = sf_floatalloc(n12);
	    
	for (i3=0; i3 < n3; i3++) {
	    /* read data */
	    sf_floatread(u1,n12,in);
	    
	    /* read t-x dip */
	    sf_floatread(p,n12,dip);
		
	    ap = allpass_init (nw,nj1,n1,n2,1,p);
		
	    /* apply */
	    if (left) {
		left1(false, false, ap, u1, u2);
	    } else {
		right1(false, false, ap, u1, u2);
	    }
		
	    /* write t-x destruction */
	    sf_floatwrite(u2,n12,out);
	}
	
	free(u1);
	free(u2);
	free(p);
	    
    }
    
    exit (0);
}

int main(int argc, char* argv[])
{
    bool inv;
    int nt, nx, nv, it, ix, iv, nw, n3, i3, ntr, ntm, im;
    float *trace=NULL, *modl=NULL, *r=NULL;
    float vmin, vmax, dv, dx, x0, t0, t, dt, tp, xp;
    char *unit=NULL, *space=NULL, *time=NULL;
    size_t len;
    sf_bands spl=NULL;
    sf_file in, out;

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

    if (!sf_getbool("inv",&inv)) inv=false;
    /* if y, do inverse transform */

    if (!sf_getint("nw",&nw)) nw=2;
    /* accuracy level */

    if (!sf_histint(in,"n2",&nt)) sf_error("No n2= in input");
    if (!sf_histfloat(in,"d2",&dt)) sf