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

static void mig_grab_par(int *op,int *mig_nz,float *mig_dz,float *fd_a,float *fd_b,float *trick)
{
    if (!sf_getint  ("operator",op))    *op=1;
    if (!sf_getfloat("mig_dz",mig_dz))  *mig_dz=4.0;
    if (!sf_getint  ("mig_nz",mig_nz))  *mig_nz=1000;
    if (!sf_getfloat("fd_a",fd_a))      *fd_a=0.5;
    if (!sf_getfloat("fd_b",fd_b))      *fd_b=0.0;
    if (!sf_getfloat("fd_trick",trick)) *trick=1.0/6.0;
}

int main(int argc, char* argv[])
{
   int j, nx, nt, n;
   float dt, dx, lambda; /* dt corresponds to k, dx corresponds to h */
   float *uo, *up, *u_true, *dif, *t; /*uo->u(n), up->u(n+1)*/
   sf_file out;
   sf_init(argc,argv);
   out=sf_output("out");
   
  if(!sf_getint("nt",&nt)) sf_error("No nt in command line!");
  /* number of temporal points */
  if(!sf_getfloat("dt",&dt)) sf_error("No dt in command line!");
  /* temporal sampling */

  if(!sf_getint("nx",&nx)) sf_error("No nx in command line!");
  /* number of spatial points */
  if(!sf_getfloat("dx",&dx)) sf_error("No dx in command line!");
  /* spatial sampling */
   
  sf_putint(out,"n1",nx);
  sf_putfloat(out,"d1",dx);
  sf_putfloat(out,"o1",0);
 
   lambda=dt/(dx*dx); /* correpsonding to a^2k/h^2 */

   uo=sf_floatalloc(nx);
   up=sf_floatalloc(nx);
   u_true=sf_floatalloc(nx);
   dif=sf_floatalloc(nx);

/* initialize it=0 */
   for(j=0;j<nx;j++)
	uo[j]=f(dx*j);

/* loop over it until finish computing up[nt-1] */
   for(n=0;n<nt-1;n++)  /* (nt-1) iterations */
	{
	up[0]=0;up[nx-1]=0;
	for(j=1;j<nx-1;j++)
		up[j]=lambda*uo[j-1]+(1-2*lambda)*uo[j]+lambda*uo[j+1];
	t=uo;
	uo=up;
	up=t;
	}
   for(j=0;j<nx;j++)
	{u_true[j]=f_true(j*dx,(nt-1)*dt); dif[j]=fabs(u_true[j]-uo[j]);}

   for(j=0;j<nx;j++)
   {
	sf_warning("%.1f   %.8f   %.8f   %.3e",j*dx, u_true[j], uo[j], dif[j]);
   }

   /* Here uo refers to the final result. */
   sf_floatwrite(uo, nx, out);
   exit(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 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);
}

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 n1, n2, i1, i2;
	sf_file in, out;
	float *u1, *v, *p, *q, lamda, sigma;

	sf_init(argc, argv);

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

	if(!sf_getfloat("lamda", &lamda)) lamda=0.9;
	if(!sf_getfloat("sigma", &sigma)) sigma=1.0;
	if(!sf_histint(in, "n1", &n1)) sf_error("n1 needed in input");
	n2 = sf_leftsize(in, 1);

	u1 = sf_floatalloc(n1);
	v = sf_floatalloc(n1);
	p = sf_floatalloc(n1);
	q = sf_floatalloc(n1);

	for(i2=0; i2<n2; i2++)
	{
		sf_floatread(u1, n1, in);
		v[0] = u1[0];
		for(i1=1; i1<n1; i1++)
		{
			p[i1] = (u1[i1] - v[i1-1])/sigma;
			p[i1] *= p[i1];
			p[i1] = expf(-p[i1]);
			v[i1] = lamda * p[i1];
			v[i1] = (1-v[i1])*u1[i1] + v[i1]*v[i1-1];
		}
		v[n1-1] = u1[n1-1];
		for(i1=n1-2; i1>=0; i1--)
		{
			q[i1] = (u1[i1] - v[i1+1])/sigma;
			q[i1] *= q[i1];
			q[i1] = expf(-q[i1]);
			v[i1] = lamda * q[i1];
			v[i1] = (1-v[i1])*u1[i1] + v[i1]*v[i1+1];
		}
		for(i1=1; i1<n1-1; i1++) v[i1] = (p[i1] * q[i1]);
		v[0] = v[1]; v[n1-1] = v[n1-2];
		sf_floatwrite(v, n1, out);
	}

	free(u1);
	free(p);
	free(q);
	free(v);
}

int main(int argc, char* argv[])
{
    int n1, ia, na, ib, nb;
    float a0, b0, a1, b1, da, db, a, b, *trace, *filt;
    sf_file inp, out;
    
    sf_init(argc,argv);
    inp = sf_input("in");
    out = sf_output("out");

    if (!sf_histint(inp,"n1",&n1)) sf_error("No n1= in input");
    
    trace = sf_floatalloc(n1);
    filt = sf_floatalloc(n1);

    if (!sf_getint("na",&na)) na=1;
    if (!sf_getint("nb",&nb)) nb=1;
    if (!sf_getfloat("a0",&a0)) sf_error("Need a0=");
    if (!sf_getfloat("b0",&b0)) sf_error("Need b0=");
    if (!sf_getfloat("a1",&a1)) a1=a0;
    if (!sf_getfloat("b1",&b1)) b1=b0;

    da = (1 >= na)?0.: (a1-a0)/(na-1);
    db = (1 >= nb)?0.: (b1-b0)/(nb-1);

    sf_putint(out,"n2",na);
    sf_putint(out,"n3",nb);

    sf_putfloat(out,"o2",a0);
    sf_putfloat(out,"o3",b0);

    sf_putfloat(out,"d2",da);
    sf_putfloat(out,"d3",db);

    sf_floatread(trace,n1,inp);

    for (ib=0; ib < nb; ib++) {
	b = b0 + ib*db;
	for (ia=0; ia < na; ia++) {
	    a = a0 + ia*da;
	    locov_init(a,b);
	    locov_filt(n1,trace,filt);

	    sf_floatwrite(filt,n1,out);
	}
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int n1, n2, n12, nv, inv;
    bool adj, add;
    float *modl, *data, o1,d1,o2,d2, vel,v0;
    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");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n12 = n1*n2;

    if (!sf_histfloat(in,"o1",&o1)) o1=0.;
    if (!sf_histfloat(in,"o2",&o2)) o2=0.;

    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"d2",&d2)) sf_error("No d2= in input");

    if (!sf_getfloat("vel",&vel)) vel=0.75;
    /* final velocity */
    if (!sf_getfloat("v0",&v0)) v0=0.;
    /* starting velocity */
    if (!sf_getint("nv",&nv)) nv=n1;
    /* number of steps */

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

    if (!sf_getint("inv",&inv)) inv=0;
    /* amplitude type */

    modl = sf_floatalloc(n12);
    data = sf_floatalloc(n12);

    sf_floatread (adj? data: modl,n12,in);

    velcon_init (inv,vel,v0,o1,n1,n2,nv,d1,d2);
    velcon_lop (adj, add, n12, n12, modl, data);

    sf_floatwrite (adj? modl: data,n12,out);

    exit(0);
}

int main(int argc, char*argv[])
{
    sf_file in, out ;
    int i1, i2, n1, n2, *v;
    float o1, d1, **u;
    char *l1, *u1;
    sf_axis ax;

    sf_init(argc, argv);

    in  = sf_input("in");  /* delay file (int) */
    out = sf_output("out");

    if(!sf_histint(in, "n1", &n2)) sf_error("n1 needed");
    sf_shiftdim(in, out, 1);


    if(!sf_getint("n1", &n1)) n1=1000; /* samples */
    if(!sf_getfloat("o1", &o1)) o1=0.0; /* sampling interval */
    if(!sf_getfloat("d1", &d1)) d1=0.004; /* original  */
    if((l1=sf_getstring("l1")) == NULL) l1="Time"; /* label "Time" */
    if((u1=sf_getstring("u1")) == NULL) u1="s"; /* unit "s" */

    ax = sf_maxa(n1, o1, d1);
    sf_setlabel(ax, l1);
    sf_setunit(ax, u1);
    sf_oaxa(out, ax, 1);
    sf_putint(out, "n2", n2);
    sf_settype(out, SF_FLOAT);

    v = sf_intalloc(n2);
    u = sf_floatalloc2(n1, n2);
    sf_intread(v, n2, in);

    for(i2=0; i2<n2; i2++)
        for(i1=0; i1<n1; i1++)
            if(i1==v[i2])	u[i2][i1] = 1;
            else u[i2][i1] = 0;

    sf_floatwrite(u[0], n1*n2, out);

    free(v);
    free(u[0]);
    free(u);

    return 0;

}

int main(int argc, char* argv[])
{
    bool nan;
    int i, n, nbuf;
    float clip, value, *trace;
    sf_file in=NULL, out=NULL; /* 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");

    if (!sf_getfloat("clip",&clip)) sf_error("Need clip=");
    /* clip value */

    if (!sf_getfloat("value",&value)) value=clip;
    /* replacement value */

    /* allocate floating point buffer */
    nbuf = sf_bufsiz(in)/sizeof(float);
    trace = sf_floatalloc (nbuf);

    /* process buffers */
    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++) {
	    nan = (bool) !isfinite(trace[i]);

	    if (nan) trace[i] = SF_SIG(trace[i])*value; 
	    else if (trace[i] >  clip) trace[i]= value;
	    else if (trace[i] < -clip) trace[i]=-value;
	}

	sf_floatwrite(trace,nbuf,out);
    }


    exit(0);
}

int main(int argc, char* argv[]) {
    int n1, n2, n3, i1, i2, i3;
    float t0,v0, dt,dv, t1, v1, slope,tmax,band, t,taper, *data;
    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");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);

    if (!sf_histfloat(in,"o1",&t0)) sf_error("No o1= in input");
    if (!sf_histfloat(in,"d1",&dt)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o2",&v0)) sf_error("No o2= in input");
    if (!sf_histfloat(in,"d2",&dv)) sf_error("No d2= in input");

    if (!sf_getfloat("t1",&t1)) t1=t0; // start time
    if (!sf_getfloat("v1",&v1)) v1=v0+(n2-1)*dv; // end velocity
    if (!sf_getfloat("band",&band)) band=20*dt; // start time
    slope = (t0+(n1-1)*dt-t1)/(v1-v0);

    data = sf_floatalloc(n1);

    for (i3=0; i3 < n3; i3++) {
        for (i2=0; i2 < n2; i2++) {
            tmax = t1 + i2*dv*slope;
            sf_floatread (data,n1,in);

            for (i1=0; i1 < n1; i1++) {
                t = t0+i1*dt;

                if (t > tmax) {
                    data[i1] = 0.0f;
                } else if (t > tmax-band) {
                    taper = sinf(0.5 * SF_PI * (t-tmax)/band);
                    data[i1] *= taper*taper;
                }
            }
            sf_floatwrite (data,n1,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);
}

void vp_minmax (float emin1, float emin2, float emax1, float emax2)
{
    bool btmp;
    float ftmp;

    fmin1 = (!sf_getfloat ("min1",&min1)); if(fmin1) min1=emin1;
    fmin2 = (!sf_getfloat ("min2",&min2)); if(fmin2) min2=emin2;
    fmax1 = (!sf_getfloat ("max1",&max1)); if(fmax1) max1=emax1;
    fmax2 = (!sf_getfloat ("max2",&max2)); if(fmax2) max2=emax2;

    if (transp) {
	ftmp=min1; min1=min2; min2=ftmp;
	btmp=fmin1; fmin1=fmin2; fmin2=btmp;

	ftmp=max1; max1=max2; max2=ftmp;
	btmp=fmax1; fmax1=fmax2; fmax2=btmp;
    }
}

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

void shot_image_init(float dx,float dy,float dz,int nz,struct shot_image_par_type *image_par)
/*< initialize image >*/
{
    sf_file imagex, imagez;
    float *tmpimg;
    float ohx,ohy;
    int i; 

    tmpimg=sf_floatalloc(nz);
 
    vector_value_f(tmpimg,0.0,nz);
    image_par->dx=dx; image_par->dy=dy; image_par->dz=dz;
    image_par->nx=(image_par->max_x-image_par->min_x)/(image_par->dx)+1;
    image_par->ny=(image_par->max_y-image_par->min_y)/(image_par->dy)+1;
    image_par->nz=nz;
    if (!sf_getfloat("imageohx",&ohx)) ohx=-image_par->dx*image_par->nhx/2;
    if (!sf_getfloat("imageohy",&ohy)) ohy=-image_par->dy*image_par->nhy/2;
    image_par->ohx=ohx;   image_par->ohy=ohy;
    
    imagex = sf_output("Image_hx");
    imagez = sf_output("Image_hz");

    sf_putfloat(imagex,"o1",0.0);              sf_putfloat(imagex,"d1",image_par->dz); sf_putint(imagex,"n1",image_par->nz);
    sf_putfloat(imagex,"o2",image_par->ohx);   sf_putfloat(imagex,"d2",image_par->dx); sf_putint(imagex,"n2",image_par->nhx);
    sf_putfloat(imagex,"o3",image_par->ohy);   sf_putfloat(imagex,"d3",image_par->dy); sf_putint(imagex,"n3",image_par->nhy);
    sf_putfloat(imagex,"o4",image_par->min_x); sf_putfloat(imagex,"d4",image_par->dx); sf_putint(imagex,"n4",image_par->nx);
    sf_putfloat(imagex,"o5",image_par->min_y); sf_putfloat(imagex,"d5",image_par->dy); sf_putint(imagex,"n5",image_par->ny);

    sf_putfloat(imagez,"o1",0.0);              sf_putfloat(imagez,"d1",image_par->dz); sf_putint(imagez,"n1",image_par->nz);
    sf_putfloat(imagez,"o2",image_par->ohx);   sf_putfloat(imagez,"d2",image_par->dx); sf_putint(imagez,"n2",image_par->nhx);
    sf_putfloat(imagez,"o3",image_par->ohy);   sf_putfloat(imagez,"d3",image_par->dy); sf_putint(imagez,"n3",image_par->nhy);
    sf_putfloat(imagez,"o4",image_par->min_x); sf_putfloat(imagez,"d4",image_par->dx); sf_putint(imagez,"n4",image_par->nx);
    sf_putfloat(imagez,"o5",image_par->min_y); sf_putfloat(imagez,"d5",image_par->dy); sf_putint(imagez,"n5",image_par->ny);

    for(i=0;i<image_par->ny*image_par->nx*image_par->nhy*image_par->nhx;i++) {
	sf_floatwrite(tmpimg,image_par->nz,imagex);
	sf_floatwrite(tmpimg,image_par->nz,imagez);
    }
    free(tmpimg);

    image_par->imagex = imagex;
    image_par->imagez = imagez;
}

struct shot_image_par_type shot_image_grab_par(void)
/*< initialize >*/
{
    struct shot_image_par_type image_par;
    float min_x,max_x,min_y,max_y;
    int nhx,nhy,ihx;
    if (!sf_getint("nhx",&nhx) ) nhx=1;
    if (!sf_getint("nhy",&nhy) ) nhy=1;
    if (!sf_getfloat("image_min_x",&min_x)) min_x=1000.0;
    if (!sf_getfloat("image_max_x",&max_x)) max_x=10000.0;
    if (!sf_getfloat("image_min_y",&min_y)) min_y=0.0;
    if (!sf_getfloat("image_max_y",&max_y)) max_y=0.0;
    image_par.nhx=nhx; image_par.nhy=nhy;
    image_par.min_x=min_x; image_par.max_x=max_x;
    image_par.min_y=min_y; image_par.max_y=max_y;
    image_par.sincx_table=sf_floatalloc(8*2001);
    for(ihx=0;ihx<2001;ihx++) mksinc((float)(ihx)*0.0005,8,image_par.sincx_table+ihx*8); 
    sf_warning("Image: minx=%f  maxx=%f\n",image_par.min_x,image_par.max_x);
    return image_par; 
}