int main (int argc, char **argv)
{
    bool verb;                /* verbosity flag */
    bool adj;
    int nt;                   /* number of time samples */
    int nx;                   /* number of offsets */ 
    int nf1, nf2, nf;         /* number of filter size */
    int n3, i3;
    float *x, *y, *filt;      /* input and output */
    sf_file in, out, fil;

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

    if (!sf_getbool("adj",&adj)) adj=true;
    /* if y, perform adjoint operation */

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

    /* read input file parameters */
    if (!sf_histint(in,"n1",&nt)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&nx)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);

    if (!sf_histint(fil,"n1",&nf1)) sf_error("No nf1= in filt");
    if (!sf_histint(fil,"n2",&nf2)) sf_error("No nf2= in filt");
    if (!sf_histint(fil,"n3",&nf)) sf_error("No nf3= in filt");

    if (nt != nf1 || nx != nf2 ) sf_error("Nee n1==nf1 && n2==nf2");

    x = sf_floatalloc (nt*nx);
    y = sf_floatalloc (nt*nx);
    filt = sf_floatalloc (nt*nx*nf);

    for (i3=0; i3 < n3; i3++) {
	if(verb) sf_warning("i=%d of %d",i3+1,n3);

	sf_floatread(filt,nt*nx*nf,fil);

	matching_init(filt,nt,nx,nf);

	if (adj) {
	    sf_floatread(y,nt*nx,in);
	} else { /* modeling */
	    sf_floatread(x,nt*nx,in);
	}
    
	pmatching_lop (adj, false, nt*nx, nt*nx, x, y);
	
	if (adj) {
	    sf_floatwrite(x,nt*nx,out);
	} else { /* modeling */
	    sf_floatwrite(y,nt*nx,out);
	}
    }
    
    exit (0);
}

int main(int argc, char* argv[])
{
    int nd, m1, na, nr, niter;
    float *rr, *dd, *coord, o1, d1, eps;
    char *header;
    sf_file in, out, head;

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

    /* create model */
    if (!sf_getint ("nx",&m1)) sf_error("Need n1=");
    /* number of bins */
    sf_putint(out,"n1",m1);

    if (!sf_getfloat("x0",&o1)) sf_error("Need o1=");
    /* grid origin */
    sf_putfloat (out,"o1",o1);

    if (!sf_getfloat("dx",&d1)) sf_error("Need d1=");
    /* grid sampling */
    sf_putfloat (out,"d1",d1);

    if (!sf_getint("niter",&niter)) niter=1+m1*3/2; niter *= 2;
    /* number of conjugate-gradient iterations */
    if (!sf_getfloat("eps",&eps)) eps=0.2;
    /* regularization parameter */

    /* create filter */
    if (!sf_getint("na",&na)) na=3;
    nr = m1 + na;

    rr = sf_floatalloc(nr);

    if (!sf_histint(in,"n1",&nd)) nd=1;
    coord = sf_floatalloc(nd);
    dd = sf_floatalloc(nd);

    header = sf_getstring("head");
    if (NULL == header) { 
	header = sf_histstring(in,"head");
	if (NULL == header) sf_error("Need head=");
    }

    head = sf_input(header);
    if (SF_FLOAT != sf_gettype(head)) sf_error("Need float head");
    sf_floatread (coord,nd,head);
    sf_fileclose (head);

    sf_floatread (dd,nd,in);

    levint1 (niter, m1, nr, nd, coord, dd, o1, d1, rr, eps);
    
    sf_floatwrite (rr,m1,out);


    exit(0);
}

int main(int argc, char* argv[])
{
    int n1,n2,n12;
    bool adj;
    float d1,d2,o1,o2,v, *model=NULL, *dat=NULL;
    sf_file in=NULL, out=NULL;

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

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

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

    if (adj) {
	if (!sf_histint(in,"n2",&n2)) sf_error("No n2= 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");

	sf_putint(out,"n2",1);
    } else {
	if (!sf_getint("n2",&n2)) n2=20;
	/* number of offsets (if inv=n) */
	if (!sf_getfloat("d2",&d2)) d2=200.;
	/* offset sampling (if inv=n) */
	if (!sf_getfloat("o2",&o2)) o2=0.;
	/* offset origin (if inv=n) */

	sf_putint(out,"n2",n2);
	sf_putfloat(out,"d2",d2);
	sf_putfloat(out,"o2",o2);
    }
    n12 = n1*n2;

    if (!sf_getfloat ("v",&v)) v=1000.;
    /* velocity */

    dat = sf_floatalloc(n12);
    model =  sf_floatalloc(n1);

    imospray_init (1./v, o2,d2, o1,d1, n1,n2);

    if (adj) {
	sf_floatread(dat,n12,in);
	imospray_lop (true,false,n1,n12,model,dat);
	sf_floatwrite(model,n1,out);
    } else {	
	sf_floatread(model,n1,in);
	imospray_lop (false,false,n1,n12,model,dat);
	sf_floatwrite(dat,n12,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    bool adj, verb;
    int n1, n2, n12, n3, i3, order, ns;
    float *input, *smooth, **slope, eps;
    sf_file in, out, dip;

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

    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);
    n12 = n1*n2;

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

    if (!sf_getint("ns",&ns)) ns=0;
    /* smoothing radius */

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

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

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

    input = sf_floatalloc(n12);
    smooth = sf_floatalloc(n12);
    slope = sf_floatalloc2(n1,n2);

    pwsmooth_init(ns, n1, n2, order, eps);

    for (i3=0; i3 < n3; i3++) {
	if (verb) sf_warning("slice %d of %d;",i3+1,n3);

	sf_floatread(input,n12,in);
	sf_floatread(slope[0],n12,dip);

	pwsmooth_set(slope);

	if (adj) {
	    pwsmooth_lop(true,false,n12,n12,smooth,input);
	} else {
	    pwsmooth_lop(false,false,n12,n12,input,smooth);
	}

	sf_floatwrite(smooth,n12,out);
    }
    if (verb) sf_warning(".");

    exit(0);
}

int main (int argc, char *argv[])
{
    int ir, nr, n1,n2, m1, m2, n12, nw, n3;
    float *u1, *u2, *p1, *p2;
    sf_file in, out, ang;
    omni2 ap;

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

    if (SF_FLOAT != sf_gettype(in) ||
	SF_FLOAT != sf_gettype(ang)) 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;
    n12 = n1*n2;
    nr = sf_leftsize(in,2);
    
    if (!sf_histint(ang,"n1",&m1) || m1 != n1) 
	sf_error("Need n1=%d in dip",n1);
    if (1 != n2 && (!sf_histint(ang,"n2",&m2) || m2 != n2)) 
	sf_error("Need n2=%d in dip",n2);
    if (!sf_histint(ang,"n3",&n3) || 2 != n3)
	sf_error("Need n3=2 in dip");
    
    if (!sf_getint("order",&nw)) nw=1;
    /* accuracy */

    u1 = sf_floatalloc(n12);
    u2 = sf_floatalloc(n12);
    p1 = sf_floatalloc(n12);
    p2 = sf_floatalloc(n12);

    for (ir=0; ir < nr; ir++) {
	/* read dip */

	sf_floatread(p1,n12,ang);
	sf_floatread(p2,n12,ang);
	
	/* read data */
	sf_floatread(u1,n12,in);
		
	ap = opwd2_init (nw,n1,n2,p1,p2);
		
	/* apply */
	opwd21(false, false, ap, u1, u2);
		
	/* write out */
	sf_floatwrite(u2,n12,out);
    }
	        
    exit (0);
}

int main(int argc, char*argv[])
{
	int n1, n2, n3,n4, i2, i3, i4, m2, m3;
	sf_file in, hor, out;
	float **u, **v, **h, o1, d1;
	char *interp;
	sinterp intp;

	sf_init(argc, argv);

	in = sf_input("in");
	hor = sf_input("horizon");
	out = sf_output("out");

	if(!sf_histint(in, "n1", &n1)) sf_error("n1 needed in input");
	if(!sf_histfloat(in, "o1", &o1)) o1=0.0;
	if(!sf_histfloat(in, "d1", &d1)) d1=1.0;
	if(!sf_histint(in, "n2", &n2)) sf_error("n2 needed in input");
	if(!sf_histint(in, "n3", &n3)) n3=1; 
	if(!sf_histint(in, "n4", &n4)) n4=1; 
	if(!sf_histint(hor, "n1", &m2)) sf_error("n1 needed in horizon file");
	if(!sf_histint(hor, "n2", &m3)) m3=1; 

	if((interp=sf_getstring("interp")) ==NULL ) interp="linear";
	/*< interpolation method: nearest, linear >*/

	intp = sinterp_c2f(interp);
	if(n2!=m2 || n3!= m3) sf_error("horizon file not match");

	sf_unshiftdim(in,out,1);
	u = sf_floatalloc2(n1, n2);
	v = sf_floatalloc2(n2, n3);
	h = sf_floatalloc2(n2, n3);

	sf_floatread(h[0], n3*n2, hor);
	for(i4=0; i4<n4; i4++)
	{
		for(i3=0; i3<n3; i3++)
		{
			sf_floatread(u[0], n1*n2, in);
			for(i2=0; i2<n2; i2++)
			v[i3][i2] = intp(u[i2], (h[i3][i2]-o1)/d1, n1);
		}
		sf_floatwrite(*v, n2*n3, out);
	}

	free(*u);
	free(u);
	free(*v);
	free(v);
	free(*h);
	free(h);
}

int main(int argc, char* argv[])
{
    int nx, nf, ny, i2, n2, lag;
    bool each, trans;
    float *xx, *yy, *ff;
    sf_file in, out, filt;

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

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

    if (!sf_histint(filt,"n1",&nf)) sf_error("No n1= in filtin");

    xx = sf_floatalloc(nx);
    ff = sf_floatalloc(nf);

    if (!sf_getbool("trans",&trans)) trans=false;
    /* if y, transient convolution; if n, internal */
    if (!sf_getbool("each",&each)) each=false;
    /* if y, new filter for each trace */
    if (!sf_getint("lag",&lag)) lag=1;
    /* lag for internal convolution */

    ny = trans? (nx+nf-1)*(nx+nf-1): nx;
    yy = sf_floatalloc(ny);
    if (trans){
        sf_putint(out,"n1",nx+nf-1);
        sf_putint(out,"n2",nx+nf-1);
    }

    if (!each) sf_floatread (ff,nf,filt);

    for (i2=0; i2 < n2; i2++) {
        if (each) sf_floatread (ff,nf,filt);
        sf_floatread (xx,nx,in);
        if (trans) {
	    tcai1_ns_init(nf,ff);
	    tcai1_ns_lop (false, false,nx,ny,xx,yy);
	} else {
	    icai1_init(nf,ff,lag);
	    icai1_lop (false, false,nx,ny,xx,yy);
	}
	sf_floatwrite (yy,ny,out);
    }

    exit(0);
}

bool get_tah(float* trace, float* header, 
	     int n1_traces, int n1_headers, sf_file file)
/*< get tah.  return 1 if eof encountered, 0 otherwise >*/
{
  int input_record_length;
  char type_input_record[5];

  /* sergey suggests
     if (!sf_try_charread("tah",file)) return 1;
     but I cannot find function and I donot think it will do all that I need
     For now I am changing back to my stub kls */
  
  if(4!=sf_try_charread2(type_input_record,4,file)){
    /* must have encounterred eof. */
    return true;
  }
  type_input_record[4]='\0';
  sf_intread(&input_record_length,1,file);
  /* fprintf(stderr,"sf_get_tah type_input_record=\"%s\"\n",
     type_input_record); */
  if(strcmp(type_input_record,"tah ")!=0){
    /* not my kind of record.  Just write it back out */
    /* Right now tah programs only have tah records.  Other type records 
       used to be be helpful to me.  They allowed me to share the pipe.  
       Program looked for the type records and they wanted and ignored the 
       ones it was not interested in an may did not even know.  If I 
       add other record types this is where I need to add the 
       sf_get_tah_buffer allocations kls */
    sf_error("non tah record found. Better write this code segment\n");
  }     
  
  if(sizeof(float)*(n1_traces+n1_headers)!=input_record_length){
      sf_warning("%s: n1_traces, and n1_headers are",__FILE__);
      sf_warning("inconsistent with input record length");
      sf_warning("n1_traces=%d, n1_headers=%d, sizeof(float)=%d",
		  n1_traces   , n1_headers  ,(int)sizeof(float) );
      sf_warning("input_record_length=%d",input_record_length);
      sf_error("sizeof(float)*(n1_traces+n1_headers)!=input_record_length");
  }
  
  /*fprintf(stderr,"read n1_traces=%d floats from input file\n",n1_traces); */

  sf_floatread(trace,n1_traces,file);

  /*fprintf(stderr,"read n1_headers=%d floats from input file\n",n1_headers);*/

  sf_floatread(header,n1_headers,file);

  return false;
}

int main(int argc, char* argv[])
{
    bool verb;
    int i, dim, n[SF_MAX_DIM], nd, rect[SF_MAX_DIM], niter;
    float *noi, *sig, *rat, eps;
    char key[6];
    sf_file fnoi, fsig, flow;

    sf_init(argc,argv);
    fnoi = sf_input("in");
    fsig = sf_input("sig"); /* input signal */

    if (SF_FLOAT != sf_gettype(fnoi) ||
	SF_FLOAT != sf_gettype(fsig)) sf_error("Need float input");

    flow = sf_output("out");

    dim = sf_filedims (fnoi,n);
    nd = 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 )*/ 
	nd *= n[i];
    }

    noi = sf_floatalloc(nd);
    sig = sf_floatalloc(nd);
    rat = sf_floatalloc(nd);    

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

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

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

    sf_divn_init(dim, nd, n, rect, niter, verb);

    sf_floatread(noi,nd,fnoi);
    sf_floatread(sig,nd,fsig);

    sf_divne (noi, sig, rat, eps);

    sf_floatwrite(rat,nd,flow);

    exit(0);
}

int main (int argc, char* argv[]) 
{
    bool der;
    int i1, n1, i2, n2, irep, nrep;
    float *data, *data2, rect;
    sf_file in, out;

    sf_init (argc, argv);
    in = sf_input ("in");
    out = 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");
    n2 = sf_leftsize(in,1);

    if (!sf_getbool("der",&der)) der=false;
    /* compute derivative */

    if (!sf_getfloat("rect",&rect)) rect=1;
    /* smoothing radius */ 
    recgauss_init (n1,der,rect);

    data = sf_floatalloc (n1);
    data2 = der? sf_floatalloc (n1): NULL;

    if (!sf_getint("repeat",&nrep)) nrep=1;
    /* repeat filtering several times */

    for (i2=0; i2 < n2; i2++) {
	if (der) {
	    sf_floatread(data2,n1,in);
	    data[0] = data2[1]-data2[0];
	    for (i1=1; i1 < n1-1; i1++) {
		data[i1] = 0.5*(data2[i1+1]-data2[i1-1]);
	    }
	    data[n1-1] = data2[n1-1]-data2[n1-2];
	} else {
	    sf_floatread(data,n1,in);
	}

	for (irep=0; irep < nrep; irep++) {
	    recgauss (data);
	}
	sf_floatwrite(data,n1,out);
    }    

    exit (0);
}

int main(int argc, char* argv[])
{
    int j, k, n, n2, i3, n3, iter, niter;
    float **a, *e, **v, s2;
    sf_file mat, val, eig;

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

    if (SF_FLOAT != sf_gettype(mat)) sf_error("Need float input");
    if (!sf_histint(mat,"n1",&n)) sf_error("No n1= in input");
    if (!sf_histint(mat,"n2",&n2) || n2 != n) sf_error("Need n1=n2 in input");
    n3 = sf_leftsize(mat,2);

    sf_putint(val,"n2",1);

    if (!sf_getint("niter",&niter)) niter=10;

    a = sf_floatalloc2(n,n);
    e = sf_floatalloc(n);

    if (NULL != sf_getstring("eig")) {
	eig = sf_output("eig"); /* eigenvectors */
	v = sf_floatalloc2(n,n);
	for (j=0; j < n; j++) {
	    for (k=0; k < n; k++) {
		v[j][k] = (j==k)? 1.0:0.0;
	    }
	}
    } else {
	eig = NULL;
	v = NULL;
    }

    jacobi_init(n);

    for (i3=0; i3 < n3; i3++) {
	sf_floatread(a[0],n*n,mat);
	
	for (iter=0; iter < niter; iter++) {
	    s2 = 0.;
	    for (j=0; j < n-1; j++) {
		for (k=j+1; k < n; k++) {
		    s2 += jacobi(a,j,k,v);
		}
	    }
	    sf_warning("iter=%d s2=%g",iter+1,s2);
	}

	for (j=0; j < n; j++) {
	    e[j]=a[j][j];
	}

	sf_floatwrite(e,n, val);
	if (NULL != v) 	sf_floatwrite(v[0],n*n, eig);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, i2, nw;
    float *trace=NULL;
    sf_bands spl;
    sf_file in=NULL, out=NULL;

    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);

    if (!sf_getint("nw",&nw)) sf_error("Need nw=");
    /* filter size */

    trace = sf_floatalloc(n1);
    spl = sf_spline_init(nw,n1);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(trace,n1,in);
	sf_banded_solve(spl,trace);
	sf_floatwrite(trace,n1,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int i, id, dim, n[SF_MAX_DIM], nd;
    kiss_fft_cpx *din;
    float *dout;
    sf_file fin, fout;

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

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

    dim = sf_filedims (fin,n);
    nd = 1;

    for (i=0; i < dim; i++) nd *= n[i];

    din = (kiss_fft_cpx*) sf_complexalloc(nd);
    dout = sf_floatalloc(nd);

    sf_floatread((float*)din,nd*2,fin);

    for (id=0; id < nd; id++) {
	dout[id] = hypotf(din[id].r,din[id].i);
    }

    sf_settype(fout, SF_FLOAT);
    sf_setform(fout, SF_NATIVE);

    sf_floatwrite(dout,nd,fout);

    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, i2;
    bool adj;
    float *pp, *qq;
    sf_file in=NULL, out=NULL;

    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);

    pp = sf_floatalloc(n1);
    qq = sf_floatalloc(n1);

    if (!sf_getbool("adj",&adj)) adj=false;
    /* if y, do adjoint integration */

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(pp,n1,in);
	if (adj) {
	    sf_causint_lop (true,false,n1,n1,qq,pp);
	} else {
	    sf_causint_lop (false,false,n1,n1,pp,qq);
	}
	sf_floatwrite(qq,n1,out);
    }


    exit(0);
}

int main(int argc, char* argv[])
{
    int nt, nd;
    float dist;
    float **points, *point;
    kd_node tree, near;
    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",&nd)) sf_error("No n1= in input");
    if (!sf_histint(inp,"n2",&nt)) sf_error("No n2= in input");
    sf_putint(out,"n2",1);

    points = sf_floatalloc2(nd,nt);
    sf_floatread(points[0],nd*nt,inp);

    tree = kd_tree(points,nt,nd);

    point = sf_floatalloc(nd);
    if (!sf_getfloats("point",point,nd)) sf_error("Need point=");
    
    dist = SF_HUGE;
    kd_nearest(tree, point, 0, nd, &near, &dist);

    sf_floatwrite(kd_coord(near),nd,out);

    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, n3, i3;
    float **pp, **qq;
    sf_file in=NULL, out=NULL;

    sf_init(argc,argv);

    in = sf_input("in");
    out = sf_output("out");

    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);

    pp = sf_floatalloc2(n1,n2);
    qq = sf_floatalloc2(n1,n2);

    for (i3=0; i3 < n3; i3++) {
	sf_floatread(pp[0],n1*n2,in);
	sf_sobel2(n1,n2,pp,qq);
	sf_floatwrite(qq[0],n1*n2,out);
    }


    exit(0);
}

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

    int i,n1,n2,n22,n;
    float *A, *x, *d;
    bool verb;
    sf_file in, rhs, out;

    sf_init (argc,argv);
    in = sf_input("in");
    rhs= sf_input("rhs");
    out = sf_output("out");
    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float");

    if(!sf_histint(in,"n1",&n1)) sf_error("No n1 in input matrix!");
    if(!sf_histint(in,"n2",&n2)) sf_error("No n2 in input matrix!");
    if(!sf_histint(rhs,"n2",&n22))sf_error("No n2 in input vector!");
    if(!sf_getbool("verb",&verb)) verb=false;
    if((n1!=n2) || (n1!=n22)) sf_error("Dimension mistake!");
    n=n1;

    x = sf_floatalloc(n);
    A=sf_floatalloc(n*n);
    d=sf_floatalloc(n);

    sf_floatread(A,n*n,in);
    sf_floatread(d,n,rhs);

    trid_init(n, A);
    trid_solve(d, x);
    
    if(verb)
    {
	for(i=0;i<n;i++)
	sf_warning("x[%d]=%f",i+1,x[i]);
    }
    
   sf_putint(out,"n2",n);
   sf_putint(out,"n1",1);
   sf_floatwrite(x,n,out);

   trid_close();
   free(A);
   free(d);
   free(x);
   exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, box1,box2,klo1,khi1,klo2,khi2,kmid1,kmid2;
    float  **dat,**nsum,**neg,sumA;
    sf_file in, out; /* Input and output files */

    /* Initialize RSF */
    sf_init(argc,argv);
    /* standard input */
    in  = sf_input("in");
    /* standard output */
    out = sf_output("out");

    /* check that the input is float */
    if (SF_FLOAT != sf_gettype(in)) 
	sf_error("Need float input");

    /* n1 is the fastest dimension (trace length) */
    if (!sf_histint(in,"n1",&n1)) 
	sf_error("No n1= in input");
    /* leftsize gets n2*n3*n4*... (the number of traces) */
    n2 = sf_leftsize(in,1);

    /* parameter from the command line (i.e. box1=50 box2=50 ) */
    if (!sf_getfloat("box1",&box1)) sf_error("Need box1=");
    if (!sf_getfloat("box2",&box2)) sf_error("Nedd box2=");

    /* allocate floating point array */
    dat= sf_floatalloc2 (n1,n2);
	
    /* initialise the size of the searching box*/
    klo1=0;
    khi1=box1;
    klo2=0;
    khi2=box2;

    

    for (int klo1=0, int khi1=box1; khi1<=xmax; ++klo1,++khi1)
	
    /* loop over traces */
    for (i2=0; i2 < n2; i2++) {

	/* read a trace */
	sf_floatread(trace,n1,in);

	/* loop over samples */
	for (i1=0; i1 < n1; i1++) {
	    if      (trace[i1] >  clip) trace[i1]= clip;
	    else if (trace[i1] < -clip) trace[i1]=-clip;
	}

	/* write a trace */
	sf_floatwrite(trace,n1,out);
    }


    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, i2, esize;
    off_t pos;
    struct skey *sorted;
    float *unsorted;
    char *trace, *header;
    sf_file in, head, out;

    sf_init (argc,argv);
    in = sf_input ("in");
    out = sf_output ("out");
 
    header = sf_getstring("head");
    /* header file */
    if (NULL == header) { 
	header = sf_histstring(in,"head");
	if (NULL == header) sf_error("Need head=");
    }

    head = sf_input(header);
    if (SF_FLOAT != sf_gettype(head))
	sf_error("Need float header");
    n2 = sf_filesize(head);
 
    unsorted = sf_floatalloc(n2);
    sorted = (struct skey*) sf_alloc(n2,sizeof(struct skey));
    
    sf_floatread(unsorted,n2,head);
    for (i2 = 0; i2 < n2; i2++) {
	sorted[i2].key = unsorted[i2];
	sorted[i2].pos = i2;
    }
    free (unsorted);
    sf_fileclose(head);

    qsort(sorted,n2,sizeof(struct skey),key_compare);
 
    if (!sf_histint(in,"n1",&n1)) n1=1;
    esize = sf_esize(in);
    n1 *= esize;

    trace = sf_charalloc(n1);

    sf_unpipe(in,((off_t) n1)*((off_t) n2));
    sf_fileflush(out,in);
    sf_setform(in,SF_NATIVE);
    sf_setform(out,SF_NATIVE);

    pos = sf_tell(in);
    for (i2=0; i2<n2; i2++) {
	sf_seek(in,pos+(sorted[i2].pos)*n1,SEEK_SET);
	sf_charread(trace,n1,in);
	sf_charwrite(trace,n1,out);
    }


    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, i2, i, j;
    float *data, sigma, rand, thresh, thresh2;
    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);

    if (!sf_getfloat("sigma",&sigma)) sigma=1.;
    /* noise magnitude */

    if (!sf_getfloat("thresh",&thresh)) thresh=0.93;
    /* noise threshold */

    if (!sf_getfloat("thresh2",&thresh2)) thresh2=0.4;
    /* noise threshold */

    data = sf_floatalloc(n1);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(data,n1,in);

	random_init (121794);
	
	for (i=0; i < n1; i++) {
	    data[i] += (random0()-0.5)*0.001;
	}
	
	random_init (121094);

	i=0;
	while (i < n1) {
	    rand = random0();
	    if (rand < thresh) { /* no noise */
		i++;
	    } else { 
		for (j=i; j < n1; j++) {
		    rand = random0();
		    data[j] += sigma*rand;
	
		    rand = random0();
		    if (rand < thresh2) break;
		}
		i = j + 1;
	    }
	}
	
	sf_floatwrite(data,n1,out);
    }


    exit(0);
}