int main(int argc, char* argv[])
{
    int i, dim, nm, niter, n[SF_MAX_DIM], rect[SF_MAX_DIM];
    char key[6];
    float *rough, *smooth, *limit, lsum, li;
    sf_file inp, out, lim;

    sf_init(argc,argv);
    inp = sf_input("in");
    out = sf_output("out");
    lim = sf_input("limit");
    /* limiter */

    if (SF_FLOAT != sf_gettype(inp) ||
	SF_FLOAT != sf_gettype(lim)) sf_error("Need float input");
    dim = sf_filedims (inp,n);

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

		if (n[i] > 1) {
			snprintf(key,6,"rect%d",i+1);
			if (!sf_getint(key,rect+i)) rect[i]=1;
			/*( rect#=(1,1,...) smoothing radius on #-th axis )*/ 
		} else {
			rect[i]=1;
		}
    }

    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */
    
    sf_divn_init(dim,nm,n,rect,niter,true);
    
    rough = sf_floatalloc(nm);
    smooth = sf_floatalloc(nm);
    limit = sf_floatalloc(nm);

    sf_floatread(rough,nm,inp);
    sf_floatread(limit,nm,lim);

    lsum = 0.;
    for (i = 0; i < nm; i++) {
		li = limit[i];
		lsum += li*li;
    }
    lsum = sqrtf (lsum/nm);

    for (i=0; i < nm; i++) {
		limit[i] /= lsum;
		rough[i] *= limit[i];
    }
    
    sf_divn(rough,limit,smooth);

    sf_floatwrite(smooth,nm,out);

    exit(0);
}

int main(int argc, char* argv[])
{
    int nx, nf, i2, n2, lag, niter;
    bool single, verb;
    float eps;
    sf_complex *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_COMPLEX != sf_gettype(in) ||
	SF_COMPLEX != 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_complexalloc(nx);
    ff = sf_complexalloc(nf);

    if (!sf_getbool("single",&single)) single=true;
    /* single channel or multichannel */

    if (!sf_getint("lag",&lag)) lag=1;
    /* lag for internal convolution */

    yy = sf_complexalloc(nx);

    if (!sf_getfloat("eps",&eps)) eps=1.0f;
    /* regularizion parameter */

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

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

    if (!single) sf_complexread (ff,nf,filt);

    for (i2=0; i2 < n2; i2++) {
        if (single) sf_complexread (ff,nf,filt);
        sf_complexread (xx,nx,in);

	cicai1_init(nf,ff,lag);

	sf_csolver_reg(sf_ccopy_lop,sf_ccgstep,cicai1_lop,nx,nx,nx,yy,xx,niter,eps,"verb",verb,"end");
	sf_ccgstep_close();
	
	sf_complexwrite (yy,nx,out);
    }
    
    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 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);
}

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[])
{
    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 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 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);
}

size_t sf_esize(sf_file file)
/*< return element size >*/
{
    sf_datatype type;
	
    type = sf_gettype (file);
    switch (type) {
	case SF_FLOAT: 
	    return sizeof(float);
	    break;
	case SF_INT:
	    return sizeof(int);
	    break;
	case SF_COMPLEX:
	    return 2*sizeof(float);
	    break;
        case SF_SHORT:
            return sizeof(short);
            break;
	case SF_LONG:
	    return sizeof(off_t);
	    break;
        case SF_DOUBLE:
            return sizeof(double);
            break;
	default:
	    return sizeof(char);
	    break;
    }
}

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 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, 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,o1,d1,d2,i1,i2,ix,iz; 

    float  ***dat,**sum;
    float sum1; 
    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");
    
    /* parameters from input file*/
    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"o1",&o1)) sf_error("No o1= in input");
    if (!sf_histint(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histint(in,"d2",&d2)) sf_error("No d2= in input");
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(in,"n3",&n3)) sf_error("No n3= in input");
    fprintf(stderr,"values: %d %d\n",n1,n2,n3);
   sum = sf_floatalloc2(n3,n3);
     dat = sf_floatalloc3 (n1,n2,n3);

    fprintf(stderr,"values: %d %d\n",n1,n2);
  dat = sf_floatalloc2 (n1,n2);

   for (ix=0; ix<n2; ++ix) {
      for (iz=0; iz<n1; ++iz) {

           dat[ix][iz]=0.0;
      }
     }
sum1=0;

    fprintf(stderr,"values1: %f\n",sum1);
    sf_floatread(dat[0],n1*n2,in); 

//sum1=0;

    fprintf(stderr,"values1: %f\n",sum1);
      for (i2=o1; i2<n2; ++i2){
            for (i1=o1; i1<n1; ++i1){
      sum1 += dat[i2][i1];
        }
        }
    fprintf(stderr,"values: %f\n",sum1);
   sum[0][0]=sum1; 
    sf_floatwrite(sum[0],1*1,out);
      

  exit(0);
}

int main (int argc, char* argv[]) 
{
    int dim, dim1, i, ider, nrep, i2, n1, n2, nderiv;
    int n[SF_MAX_DIM], rect[SF_MAX_DIM];
    char key[6];
    float* data;
    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");

    dim = sf_filedims (in,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;
    }

    n1 = n2 = 1;
    for (i=0; i < dim; i++) {
	if (i <= dim1) {
	    n1 *= n[i];
	} else {
	    n2 *= n[i];
	}
    }

    data = sf_floatalloc (n1);

    if (!sf_getint("ider",&ider)) ider=0;
    /* direction of the derivative (0 means no derivative) */
    
    if (!sf_getint("repeat",&nrep)) nrep=1;
    /* repeat smoothing several times */

    if (!sf_getint("nderiv",&nderiv)) nderiv=6;
    /* derivative filter accuracy */

    sf_dtrianglen_init(dim1+1,rect,n);
    
    for (i2=0; i2 < n2; i2++) {
	sf_floatread(data,n1,in);

	sf_dtrianglen(ider, nrep, nderiv, data);
	
	sf_floatwrite(data,n1,out);
    }    

    exit (0);
}

int main(int argc, char*argv[])
{
	sf_file in, out;
	int nf, n1, n2, n3, m, n;
	int i3;
	float **wav, ****fb;
	char *interp;

	sf_init(argc, argv);

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

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

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

	sf_shiftdim2(in,out,2);

	if(!sf_getint("m", &m)) m=1;
	/* b[-m, ... ,n] */
	if(!sf_getint("n", &n)) n=1;
	/* b[-m, ... ,n] */
	if ((interp=sf_getstring("interp"))==NULL) interp="maxflat";
	/* interpolation method: maxflat lagrange bspline */

	nf = m+n+1;

	wav = sf_floatalloc2(n1,n2);
	fb  = sf_floatalloc4(n1, n2, nf, nf);

	sf_putint(out, "n3", nf);
	sf_putint(out, "n4", nf);
	sf_putfloat(out, "o3", 0);
	sf_putfloat(out, "d3", 1);
	sf_putfloat(out, "o4", 0);
	sf_putfloat(out, "d4", 1);


	fbank_init(m, n, interp);


	for(i3=0; i3<n3; i3++)
	{
		sf_floatread(wav[0], n1*n2, in);
		fbank2(n1, n2, wav, fb);
		sf_floatwrite(fb[0][0][0], n1*n2*nf*nf, out);
	}

	fbank_close();
	return 0;
}

int main (int argc, char* argv[]) {
    int i, n, nbuf;
    float upper, lower, upperval, lowerval, *trace=NULL;
    sf_file in, out;

    /* 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");

    n = sf_filesize (in);

    if (!sf_getfloat ("upper", &upper)) upper = +FLT_MAX;
    /* upper range limit */
    if (!sf_getfloat ("lower", &lower)) lower = -FLT_MAX;
    /* lower range limit */
    if (!sf_getfloat ("upperval", &upperval)) upperval = +FLT_MAX;
    /* upper range value */
    if (!sf_getfloat ("lowerval", &lowerval)) lowerval = -FLT_MAX;
    /* lower range value */

    nbuf = BUFSIZ/sizeof(float);
    trace = sf_floatalloc (nbuf);

    /* loop over traces */
    for (n = sf_filesize (in); n > 0; n -= nbuf) {
        if (nbuf > n) nbuf=n;

        sf_floatread (trace, nbuf, in);

        for (i = 0; i < nbuf; i++) {
            if (upper > lower) {
                if (trace[i] > upper)
                    trace[i] = upperval;
                if (trace[i] < lower)
                    trace[i] = lowerval;
            } else {
                if (trace[i] > upper &&
                        trace[i] < lower)
                    trace[i] = lowerval;
            }
        }

        sf_floatwrite (trace, nbuf, out);
    }

    return 0;
}

int main(int argc, char* argv[]) 
{
    int nt, nx, ntx, n2, n3, next;
    float v0, v1, dt, dx, t0, kur;
    float **data;
    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",&nt)) sf_error("No n1= in input");
    if (!sf_histfloat(inp,"d1",&dt)) sf_error("No d1= in input");
    if (!sf_histfloat(inp,"o1",&t0)) t0=0.;  

    if (!sf_histint(inp,"n2",&nx)) sf_error("No n2= in input");
    if (!sf_histfloat(inp,"d2",&dx)) sf_error("No d2= in input");

    ntx = nt*nx;

    if (!sf_getfloat("v0",&v0)) v0=SF_EPS;
    /* initial velocity */

    if (!sf_getfloat("v",&v1)) sf_error("Need v=");
    /* final velocity */

    if (!sf_getint("pad",&n2)) n2=nt; /* padding for stretch */
    if (!sf_getint("pad2",&n3)) n3=2*kiss_fft_next_fast_size((n2+1)/2);
    /* padding for FFT */

    if (!sf_getint("extend",&next)) next=4;
    /* trace extension */

    velcon_init(nt,nx,dt,dx,t0,n2,n3,next);

    data = sf_floatalloc2(nt,nx);

    sf_floatread(data[0],ntx,inp);

    kur = kurtosis(ntx,data[0]);
    sf_warning("kurtosis before: %g",kur);

    velcon(data,v0,v1);

    kur = kurtosis(ntx,data[0]);
    sf_warning("kurtosis after: %g",kur);

    sf_floatwrite(data[0],ntx,out);

    exit(0);
}

int main (int argc, char* argv[])
{
    int n, n1, n2, i1, i2, iv, nv;
    float dv, vel, v2, k, dk, dt;
    float *t, *d;
    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");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n = n1-1;

    if (!sf_getint("nv",&nv)) sf_error("Need nv=");
    if (!sf_getfloat("vel",&vel)) sf_error("Need vel=");

    v2 = 0.5*vel*vel*SF_SIG(vel);
    dv = v2/(nv-1);

    if (!sf_histfloat(in,"dt",&dt)) sf_error("No dt= in input");
    if (!sf_histfloat(in,"d2",&dk)) sf_error("No dk= in input");
    dk *= 2*SF_PI;

    t = sf_floatalloc(n1);
    d = sf_floatalloc(n);
    tricheb_init(n);

    for (i2=0; i2 < n2; i2++) {
	k = i2*dk;
	k *= k*dt*dv*0.25;

	sf_floatread(t,n1,in);
	tricheb_step(-dv,k);

	for (iv=0; iv < nv; iv++) {
	    tricheb_apply (d,t);
	    tricheb_solve (d,t+1);
	    t[0] = 0.;
	    for (i1=1; i1 < n; i1 +=2) {
		t[0] += t[i1]*SF_SIG(dv) - t[i1+1];
	    }
	}

	sf_floatwrite(t,n1,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    off_t nsiz;
    int nbuf, k, order;
    char *type;
    float *fbuf;
    sf_file inp, out;
    double (*func)(double);

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

    if (SF_FLOAT != sf_gettype(inp)) sf_error("Need float input");
    
    if (NULL == (type = sf_getstring("type"))) type="J";
    if (!sf_getint("order",&order)) order=0;

    switch(type[0]) {
	case 'i':
	case 'I':
	    switch (order) {
		case 0:
		    func = bessel_I0;
		    break;
		case 1:
		    func = bessel_I1; 
		    break;
		default:
		    sf_error("Order %d for type %s is not implemented",order,type);
	    }
	    break;
	default:
	    sf_error("Type %s is not implemented",type);
	    break;
    }

    nbuf = BUFSIZ/sizeof(float);
    fbuf = sf_floatalloc(nbuf);

    for (nsiz = sf_filesize(inp); nsiz > 0; nsiz -= nbuf) {
	if (nbuf > nsiz) nbuf = nsiz;
	sf_floatread(fbuf,nbuf,inp);
	for (k=0; k < nbuf; k++) {
	    fbuf[k] = func(fbuf[k]);
	}
	sf_floatwrite(fbuf,nbuf,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    char        *mode=NULL;             /* Standard types           */
    int         i1,i2,n1,n2;
    float       *list=NULL;

    sf_file     in=NULL,out=NULL;       /* RSF types                */

    sf_init(argc,argv);                 /* Initialize RSF           */

    in  = sf_input("in");               /* Open files               */
    out = sf_output("out");

    sf_fileflush(out,in);               /* Copy header to output    */

                                        /* Check input data type    */
    if (sf_gettype(in) != SF_FLOAT) sf_error("Input must be float.");

                                        /* Get the file sizes       */
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input.");

    n2 = sf_leftsize(in,1);

                                        /* Check mode               */
    if ((mode = sf_getstring("mode")) == NULL) mode = "min";
    /* 'min' (default) or 'max' */

    if (strcmp(mode,"min")!=0 && strcmp(mode,"max")!=0)
        sf_error ("Unknown mode %s.",mode);

    list = sf_floatalloc(n1);           /* Allocate list            */

    for (i2=0; i2<n2; i2++)             /* Process each list        */
    {
        sf_floatread(list,n1,in);       /* Read a list              */

                                        /* Find the min or max      */
        if (strcmp(mode,"min")==0)
            for (i1=1; i1<n1; i1++) list[i1] = SF_MIN(list[i1],list[i1-1]);

        if (strcmp(mode,"max")==0)
            for (i1=1; i1<n1; i1++) list[i1] = SF_MAX(list[i1],list[i1-1]);

        sf_floatwrite(list,n1,out);     /* Write the output         */
    }



    exit (0);
}