int main (int argc, char *argv[]) 
{    
    int m, k, l, seed;
    sf_file  bshuffle,ashuffle;
    sf_axis ax,at,ap,av;
    int nx,nt,np,nv,iteration,*a1, *a2;
    float ***bsh, ***bsh2;
    sf_init(argc,argv);

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

    if (!sf_getint("iteration",&iteration)) iteration=1;
    if (!sf_getint("seed",&seed)) seed=2012;
    at=sf_iaxa( bshuffle,1); nt=sf_n(at);
    ax=sf_iaxa( bshuffle,3); nx=sf_n(ax);
    ap=sf_iaxa( bshuffle,2); np=sf_n(ap);
    av=sf_iaxa( bshuffle,4); nv=sf_n(av);
    bsh=sf_floatalloc3(nt,np,nx);
    bsh2=sf_floatalloc3(nt,np,nx);
    a1=sf_intalloc(np);
    a2=sf_intalloc(np);
    sf_warning("ntpx=%d",nt*np*nx);

    srand(seed);


    for (m=0; m<np; m++) {
	a1[m]=rand();
    }

    bubble(a1, a2, np);


    for (k=0; k<nv; k++) {
	sf_floatread(bsh[0][0],nt*np*nx,bshuffle);
	for(l=0; l<nx; l++) {
	    for (m=0; m<np; m++) {
		memcpy(bsh2[l][m], bsh[l][a2[m]], nt*sizeof(float));
	    }
	}
	sf_floatwrite(bsh2[0][0],nt*np*nx,ashuffle);
    }

    free(bsh[0][0]);
    free(bsh[0]);
    free(bsh);
    free(bsh2[0][0]);
    free(bsh2[0]);
    free(bsh2);
    free(a1);
    free(a2);
    exit (0);
}

/*------------------------------------------------------------*/
lint2d lint2d_make(	int na, 
					pt2d *aa, 
					fdm2d fdm)
		   
/*< init 2D linear interpolation >*/
{
    lint2d ca;
    int    ia;
    float f1,f2;
    
    ca = (lint2d) sf_alloc(1,sizeof(*ca));

    ca->n = na;

    ca->w00 = sf_floatalloc(na);
    ca->w01 = sf_floatalloc(na);
    ca->w10 = sf_floatalloc(na);
    ca->w11 = sf_floatalloc(na);

    ca->jz  = sf_intalloc(na);
    ca->jx  = sf_intalloc(na);

    for (ia=0; ia<na; ia++) {
	
		if(aa[ia].z >= fdm->ozpad && 
			aa[ia].z <  fdm->ozpad + (fdm->nzpad-1)*fdm->dz &&
			aa[ia].x >= fdm->oxpad && 
			aa[ia].x <  fdm->oxpad + (fdm->nxpad-1)*fdm->dx   ) {
				
			ca->jz[ia] = NINT((aa[ia].z-fdm->ozpad)/fdm->dz);
			ca->jx[ia] = NINT((aa[ia].x-fdm->oxpad)/fdm->dx);
				
			f1 = (aa[ia].z-fdm->ozpad)/fdm->dz - ca->jz[ia];
			f2 = (aa[ia].x-fdm->oxpad)/fdm->dx - ca->jx[ia];
		}
		else {
			ca->jz[ia] = 0; 
			ca->jx[ia] = 0;
			
			f1 = 1; 
			f2 = 0;
		}
		
		ca->w00[ia] = (1-f1)*(1-f2);
		ca->w01[ia] = (  f1)*(1-f2);
		ca->w10[ia] = (1-f1)*(  f2);
		ca->w11[ia] = (  f1)*(  f2);
		
    }
	
    return ca;
}

/*------------------------------------------------------------*/
lint2d lint2d_make(int    na, 
		   pt2d*  aa, 
		   fdm2d fdm)
/*< init 2D linear interpolation >*/
{
    lint2d ca;
    int    ia;
    float f1,f2;
    
    ca = (lint2d) sf_alloc(1,sizeof(*ca));

    ca->n = na;

    ca->w00 = sf_floatalloc(na);
    ca->w01 = sf_floatalloc(na);
    ca->w10 = sf_floatalloc(na);
    ca->w11 = sf_floatalloc(na);

    ca->j1  = sf_intalloc(na);
    ca->j2  = sf_intalloc(na);

    for (ia=0;ia<na;ia++) {
	
	if(aa[ia].z >= fdm->o1pad && 
	   aa[ia].z <  fdm->o1pad + (fdm->n1pad-1)*fdm->d1 &&
	   aa[ia].x >= fdm->o2pad && 
	   aa[ia].x <  fdm->o2pad + (fdm->n2pad-1)*fdm->d2   ) {
	    
	    ca->j1[ia] = (int)( (aa[ia].z-fdm->o1pad)/fdm->d1);
	    ca->j2[ia] = (int)( (aa[ia].x-fdm->o2pad)/fdm->d2);
	    
	    f1 = (aa[ia].z-fdm->o1pad)/fdm->d1 - ca->j1[ia];
	    f2 = (aa[ia].x-fdm->o2pad)/fdm->d2 - ca->j2[ia];
	} else {
	    ca->j1[ia] = 0; 
	    ca->j2[ia] = 0;
	    
	    f1 = 1; 
	    f2 = 0;
	}

	ca->w00[ia] = (1-f1)*(1-f2);
	ca->w01[ia] = (  f1)*(1-f2);
	ca->w10[ia] = (1-f1)*(  f2);
	ca->w11[ia] = (  f1)*(  f2);
    }

    return ca;
}

float wmedian(float *temp,float *weight,int nfw)
/*< get a weighted median value >*/
{   
    int i, j, pass, *index, b;
    float wm, a, *data;
    
    data = sf_floatalloc(nfw);
    index = sf_intalloc(nfw);
    for (j=0; j < nfw; j++) {
	data[j] = temp[j]*weight[j];
	index[j] = j;
    }
    for(pass=1; pass < nfw; pass++) {
	for(i=0; i < nfw-pass; i++) {
	    if(data[i] > data[i+1]) {
		a = data[i];
		b = index[i];
		data[i] = data[i+1];
		index[i] = index[i+1];
		data[i+1] = a;
		index[i+1] = b;
	    }
	}
    }
    wm = temp[index[nfw/2]];
    free(index);
    free(data);
    return wm;
}

map4 stretch4_init (int n1, float o1, float d1 /* regular axis */, 
		    int nd                     /* data samples */, 
		    float eps                  /* regularization */)
/*< initialize >*/
{
    int i;
    map4 str;
    
    str = (map4) sf_alloc (1, sizeof(*str));

    str->nt = n1; 
    str->t0 = o1; 
    str->dt = d1; 
    str->nd = nd; 
    str->eps = eps;
    
    str->x = sf_intalloc (nd);
    str->m = sf_boolalloc (nd);
    str->w = sf_floatalloc2 (4,nd);
    str->diag = sf_floatalloc (str->nt);
    
    for (i = 0; i < 3; i++) {
	str->offd[i] = sf_floatalloc (str->nt-1-i);
    }
  
    str->slv = sf_banded_init (str->nt,3);
    str->tslv = sf_spline4_init(str->nt);

    return str;
}

int mutingf(int nt, int nx, float dt, float dx, float dz, int isx, int isz, int gpz, float vel, int wd, float **dat)
/*< muting function >*/
{
  float z2; /* source location */
  float tt;
  int *t;
  int ix,it,cut;

  t = sf_intalloc(nx);

  z2= powf((isz-gpz)*dz,2);

#ifdef _OPENMP
#pragma omp parallel for default(shared) private(ix,tt,cut)
#endif
  for (ix = 0; ix < nx; ix++) {
    tt = sqrtf( powf((isx-ix)*dx,2) + z2 ) / vel;
    cut = (int) (tt/dt) + wd; /* why cannot use int to convert */
    if (cut > nt) cut = nt;
    t[ix] = cut;
  }

#ifdef _OPENMP
#pragma omp parallel for default(shared) private(ix,it)
#endif
  for (ix = 0; ix < nx; ix++) {
    for (it = 0; it < t[ix]; it ++) {
      dat[ix][it] = 0.0f;
    }
  }

  return 0;
}

sf_map sf_stretch_init (int n1, float o1, float d1 /* regular axis */, 
			int nd                     /* data length */, 
			float eps                  /* regularization */, 
			bool narrow                /* if zero boundary */)
/*< initialize >*/
{
    sf_map str;
    
    str = (sf_map) sf_alloc (1, sizeof(*str));

    str->nt = n1; 
    str->t0 = o1; 
    str->dt = d1; 
    str->nd = nd; 
    str->eps = eps;
    str->narrow = narrow;
    
    str->x = sf_intalloc (nd);
    str->m = sf_boolalloc (nd);
    str->w = sf_floatalloc (nd);
    str->diag = sf_floatalloc (n1);
    str->offd = sf_floatalloc (n1-1);
  
    str->slv = sf_tridiagonal_init (n1);

    return str;
}

nfilter nallocate(int np   /* number of patches */, 
		  int nd   /* data size */, 
		  int *nh  /* filter size [np] */, 
		  int *pch /* patching [nd] */) 
/*< allocate >*/
{
    nfilter aa;
    int ip, id;
    
    aa = (nfilter) sf_alloc(1,sizeof(*aa));

    aa->np = np;
    aa->hlx = (sf_filter*) sf_alloc(np,sizeof(sf_filter));

    for (ip=0; ip < np; ip++) {
	aa->hlx[ip] = sf_allocatehelix(nh[ip]);
    }
    
    aa->pch = sf_intalloc(nd);
    for (id=0; id < nd; id++) {
	aa->pch[id] = pch[id];
    }

    aa->mis = NULL;
    return aa;
}

sf_upgrad sf_upgrad_init(int mdim        /* number of dimensions */,
			 const int *mm   /* [dim] data size */,
			 const float *d  /* [dim] data sampling */)
/*< initialize >*/
{
    sf_upgrad upg;
    int i;

    if (mdim > 3) sf_error("%s: dim=%d > 3",__FILE__,mdim);

    ndim = mdim;
    nn = mm;

    nt = 1;
    for (i=0; i < ndim; i++) {
	ss[i] = nt;
	nt *= nn[i];
	dd[i] = 1.0/(d[i]*d[i]);
    }

    upg = (sf_upgrad) sf_alloc(1,sizeof(*upg));

    upg->update = sf_ucharalloc2(2,nt);
    upg->ww = sf_floatalloc2(ndim+1,nt);
    upg->order = sf_intalloc(nt);

    return upg;
}

void fastmarch_init(int nx_init, int nz_init, int nt_init,
		    float hx_init, float hz_init, float ht_init,
		    float *x0_init, float *t0_init, float *v_init, float *s_init)
/*< initialize >*/
{
    int i,j,ind, nxz;

    nx = nx_init;
    nz = nz_init;
    nt = nt_init;
    nt1 = nt-1;
    nx1 = nx-1;

    hx = hx_init;
    hz = hz_init;
    ht = ht_init;

    x0 = x0_init;
    t0 = t0_init;

    v = v_init;
    s = s_init;

    nxz = nx*nz;

    pup= sf_intalloc(nxz);
    ms= sf_charalloc(nxz);
    r= sf_intalloc(nxz);
    p= sf_intalloc(nxz);

    count=0;

    for( i=0;i<nx;i++ ) {
	*(pup+i)=2;
	*(ms+i)='a';
	ind=i;
	for( j=1;j<nz;j++ ) {
	    ind+=nx;
	    *(pup+ind)=0;
	    *(ms+ind)='u';
	}
    }

}

void vp_plot_init(int n2)
{
    int  i, len;

    dashtype = sf_floatalloc(n2);
    if (!sf_getfloats ("dash",dashtype,n2)) {
	/*  line dash type	
	    0 continuos (default)
	    1 fine dash
	    2 fine dot
	    3 dash
	    4 large dash
	    5 dot dash
	    6 large dash small dash
	    7 double dot
	    8 double dash
	    9 loose dash  The part after the decimal point determines the pattern repetition interval */	   
	for (i = 0; i < n2; i++) 
	    dashtype[i] = 0.;
    }

    fat = sf_intalloc(n2);
    if (!sf_getints("plotfat",fat,n2)) {
	/* line fatness */
	for (i = 0; i < n2; i++)
	    fat[i] = 0;
    }
    
    col = sf_intalloc(n2);
    if (!sf_getints("plotcol",col,n2)) {
	/* line color 
	   7 white
	   6 yellow (default)
	   5 cyan
	   4 green
	   3 magenta
	   2 red
	   1 blue
	   0 black */
	for (i = 0; i < n2; i++)
	    col[i] = 6 - (i % 6);
    }
}

void* mflt_init(int m1, int mw)
/*< median filter initialize >*/
{
	struct tag_mflt *p;
	p=sf_alloc(1, sizeof(struct tag_mflt));

	p->n1 = m1;
	p->nw = mw;
	p->bi = sf_floatalloc(m1+2*mw);
	p->pi = sf_intalloc(2*mw+1);
	return p;
}

void medbin_init (float **coord       /* coordinates [nd][2] */, 
		  float o1, float o2, 
		  float d1, float d2,
		  int   n1, int   n2  /* axes */,
		  int nd_in           /* number of data points */) 
/*< initialize >*/
{
    int id, im, i1, i2, start;

    nd = nd_in; 
    nm = n1*n2;

    ct = sf_intalloc(nm);
    pt = sf_intalloc(nm);
    k = sf_intalloc(nd);
    m = sf_boolalloc(nd);
    bd = sf_floatalloc(nd);

    for (im=0; im < nm; im++) {
	ct[im] = 0.;
    }
    for (id = 0; id < nd; id++) {
	i1 = 0.5 + (coord[id][0] - o1)/d1;
	i2 = 0.5 + (coord[id][1] - o2)/d2;
	m[id] = (bool)
	    ((i1 >= 0) && (i1 < n1) &&    
	     (i2 >= 0) && (i2 < n2));
	if (m[id]) {
	    im = i2*n1+i1;
	    ct[im]++;
	    k[id] = im; 
	}
    }

    start = 0;
    for (im=0; im < nm; im++) {
	pt[im] = start;
	start += ct[im];
    }
}

sf_filter steep(int dim    /* number of dimensions */,
                int *n     /* data size [dim] */,
                int *a     /* filter size [dim] */,
                float *d   /* axis sampling [dim] */,
                float vel  /* velocity */,
                float tgap /* time gap */)
/*< define PEF >*/
{
    int *lag, c[SF_MAX_DIM], i, h, na, j;
    float x, t0;
    sf_filter aa;

    na = 1;
    for (j=0; j < dim; j++) {
        na *= a[j];
    }

    lag = sf_intalloc(na);

    for (h=i=0; i < na; i++) {
        sf_line2cart(dim, a, i, c);

        for (j=0; j < dim-1; j++) {
            c[j] -= (a[j]-1)/2;
        }

        t0 = 0.;
        for (j=0; j < dim-1; j++) {
            x = d[j]*c[j];
            t0 += x*x;
        }
        t0 = sqrtf(t0)/vel;
        if (t0 < tgap) t0 = tgap;

        x = d[dim-1]*c[dim-1];
        if(x >= t0) {
            lag[h] = sf_cart2line(dim, n, c);
            h++;
        }
    }

    aa = sf_allocatehelix(h);

    for (i=0; i < h; i++) {
        aa->lag[i] = lag[i];
        aa->flt[i] = 0.;
    }

    free(lag);

    return aa;
}

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;

}

void stretch_init (
		     int n1   /* trace length */, 
		     float o1 /* trace origin */, 
		     float d1 /* trace sampling */, 
		     int n2   /* number of data samples */)
/*< initialization >*/
{
    nt = n1; 
    t0 = o1; 
    dt = d1; 
    nd = n2;
    
    x = sf_intalloc(nd);
    m = sf_boolalloc(nd);
    w = sf_floatalloc(nd);
}

int main(int argc, char* argv[])
{
    int i1, n1, i2, n2, nv, two, *trace;
    float o1, o2, d1, d2, x, y, **vert;
    sf_file inp, out, poly;

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

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

    if (!sf_histfloat(inp,"d1",&d1)) d1=1.;
    if (!sf_histfloat(inp,"d2",&d2)) d2=1.;
    
    if (!sf_histfloat(inp,"o1",&o1)) o1=0.;
    if (!sf_histfloat(inp,"o2",&o2)) o2=0.;

    poly = sf_input("poly"); 
    /* list of polygon vertices */

    if (SF_FLOAT != sf_gettype(poly)) sf_error("Need float type in poly");
    if (!sf_histint(poly,"n1",&two) || 2 != two)
	sf_error("Need n1=2 in poly");
    if (!sf_histint(poly,"n2",&nv))
	sf_error("No n2= in poly");
    
    trace = sf_intalloc(n1);
    vert = sf_floatalloc2(2,nv);

    sf_floatread(vert[0],2*nv,poly);

    sf_settype(out,SF_INT);

    for (i2=0; i2 < n2; i2++) {
	y = o2+i2*d2;
	for (i1=0; i1 < n1; i1++) {
	    x = o1+i1*d1;
	    trace[i1] = pnpoly(nv, vert, x, y);
	}
	sf_intwrite(trace,n1,out);
    }

    exit(0);
}

/*------------------------------------------------------------*/
slo3d slow3_init(cub3d   cub,
		 sf_fslice  slice_,   /* slowness slice */
		 int     nrmax,   /* maximum number of references */
		 float   dsmax,
		 float  twoway
    )
/*< initialize slowness >*/
{
    int imz, jj;
    int ompith=0;

    /*------------------------------------------------------------*/
    slo3d slo;
    slo = (slo3d) sf_alloc(1,sizeof(*slo));

    slo->slice=slice_;
    slo->nrmax=nrmax;
    slo->dsmax=dsmax;

    if(twoway) { slo->twoway = 2;
    } else {     slo->twoway = 1;
    }

    slo->ss = sf_floatalloc3(cub->alx.n,cub->aly.n,cub->ompnth);  /* slowness */
    slo->so = sf_floatalloc3(cub->alx.n,cub->aly.n,cub->ompnth);  /* slowness */
    slo->sm = sf_floatalloc2(slo->nrmax,cub->amz.n);  /* ref slowness squared */
    slo->nr = sf_intalloc              (cub->amz.n);  /* nr of ref slownesses */
    
    for (imz=0; imz<cub->amz.n; imz++) {
	sf_fslice_get(slo->slice,imz,slo->ss[0][0]);
	slow3_twoway(cub,slo,slo->ss,ompith);

	slo->nr[imz] = slow3(slo->nrmax,
			     slo->dsmax,
			     cub->alx.n*cub->aly.n,
			     slo->ss[0][0],
			     slo->sm[imz]);
    }
    for (imz=0; imz<cub->amz.n-1; imz++) {
	for (jj=0; jj<slo->nr[imz]; jj++) {
	    slo->sm[imz][jj] = 0.5*(slo->sm[imz][jj]+slo->sm[imz+1][jj]);
	}
    }

    return slo;
}

void sf_trianglen_init (int ndim  /* number of dimensions */, 
			int *nbox /* triangle radius [ndim] */, 
			int *ndat /* data dimensions [ndim] */)
/*< initialize >*/
{
    int i;

    dim = ndim;
    n = sf_intalloc(dim);

    tr = (sf_triangle*) sf_alloc(dim,sizeof(sf_triangle));

    nd = 1;
    for (i=0; i < dim; i++) {
	tr[i] = (nbox[i] > 1)? sf_triangle_init (nbox[i],ndat[i],false): NULL;
	s[i] = nd;
	n[i] = ndat[i];
	nd *= ndat[i];
    }
    tmp = sf_floatalloc (nd);
}

void sf_psss_init(int nw0,int nx0,int nz,
	float dw,float dx, float dz0, float *v0)
/*< initialize >*/
{
	int i;
	nw=nw0;
	nx=nx0;

	dz=dz0*2.0*SF_PI/(dx*nx);

	vel=(float*)sf_floatalloc(nz);
	k2=(int*)sf_intalloc(nx);

	for(i=0; i<nx/2+1; i++)		k2[i] = i*i;
	for(i=nx/2+1; i<nx; i++) 	k2[i] = (i-nx)*(i-nx);

	for(i=0; i<nz; i++)
	{
		vel[i] = dw*dx*nx/v0[i];
		vel[i] *= vel[i];
	}
}