void sfFloatWrite1d(const char *fn, const float *dat, int n1, float d1, float o1) {
  sf_file f = sf_output(fn);
  sf_putint(f, "n1", n1);
  sf_putfloat(f, "d1", d1);
  sf_putfloat(f, "o1", o1);
  sf_floatwrite(const_cast<float *>(dat), n1, f);
}

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 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 nz, nx, ix, iz;
    float *trace, x0, z0, x, z, dx, dz, v0, g, v, a, d;
    sf_file mod;

    sf_init (argc,argv);
    mod = sf_output("out");

    if (!sf_getint("nz",&nz)) nz=601; 
    /* vertical dimension */
    if (!sf_getint("nx",&nx)) nx=401; 
    /* horizontal dimension */

    dz = 4./(nz-1);
    dx = 6./(nx-1);
    
    sf_putint   (mod,"n1",nz); 
    sf_putfloat (mod,"d1",dz); 
    sf_putfloat (mod,"o1",0.);
    sf_putint   (mod,"n2",nx); 
    sf_putfloat (mod,"d2",dx); 
    sf_putfloat (mod,"o2",0.);
    sf_setformat (mod,"native_float");

    if(!sf_getfloat ("x0",&x0)) x0=4.;
    /* anomaly center in x */
    if(!sf_getfloat ("z0",&z0)) z0=1.5;
    /* anomaly center in z */
    if(!sf_getfloat ("v0",&v0)) v0=1.5;
    /* surface velocity */
    if(!sf_getfloat ("g",&g)) g=0.6;
    /* velocity gradient */
    if(!sf_getfloat ("a",&a)) a=1.;
    /* anomaly magnitude */
    if(!sf_getfloat ("d",&d)) d=1.;
    /* anomaly radius */
    d *= d;

    trace = sf_floatalloc(nz);

    for (ix=0; ix <  nx; ix++) {
	x = ix*dx - x0;
	x *= x;
	for (iz=0; iz < nz; iz++) {
	    z = iz*dz;
	    v = v0 + g*z;
	    z = z - z0;
	    z = z*z + x;
	    trace[iz] = v + a*exp(-z/d);
	}
	sf_floatwrite (trace,nz,mod);
    }

    exit (0);
}

int main(int argc, char* argv[])
{
	int i3, n3, ix, nx, iz, nz;
	int nx2, nz2, scalex, scalez;
	float dx, dz;
	float **a, **b;
	sf_file in, out;

	sf_init(argc, argv);

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

	if(!sf_getint("scalex", &scalex)) sf_error("Need scalex=");
	if(!sf_getint("scalez", &scalez)) sf_error("Need scalez=");

	if(!sf_histint(in, "n1", &nz)) sf_error("No n1= in input");
	if(!sf_histfloat(in, "d1", &dz)) sf_error("No d1= in input");
	if(!sf_histint(in, "n2", &nx)) sf_error("No n2= in input");
	if(!sf_histfloat(in, "d2", &dx)) sf_error("No d2= in input");

	nx2=(nx-1)*scalex+1;
	nz2=(nz-1)*scalez+1;
	n3=sf_leftsize(in, 2);

	sf_putint(out, "n1", nz2);
	sf_putfloat(out, "d1", dz/scalez);
	sf_putint(out, "n2", nx2);
	sf_putfloat(out, "d2", dx/scalex);

	a=sf_floatalloc2(nz, nx);
	b=sf_floatalloc2(nz2, nx2);

	for (i3=0; i3<n3; i3++){
		sf_floatread(a[0], nz*nx, in);

#pragma omp parallel for private(ix, iz)
		for(ix=0; ix<nx2; ix++){
			for(iz=0; iz<nz2; iz++){
				b[ix][iz]=0.;
			}
		}

#pragma omp parallel for private(ix, iz)
		for(ix=0; ix<nx; ix++){
			for(iz=0; iz<nz; iz++){
				b[ix*scalex][iz*scalez]=a[ix][iz];
			}
		}

		sf_floatwrite(b[0], nz2*nx2, 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;
}

void sfFloatWrite2d(const char *fn, const float *dat, int n1, int n2,
    float d1, float d2, float o1, float o2) {
  sf_file f = sf_output(fn);
  sf_putint(f, "n1", n1);
  sf_putint(f, "n2", n2);
  sf_putfloat(f, "d1", d1);
  sf_putfloat(f, "d2", d2);
  sf_putfloat(f, "o1", o1);
  sf_putfloat(f, "o2", o2);
  sf_floatwrite(const_cast<float *>(dat), n1 * n2, f);
}

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

void sfDoubleWrite2d(const char *fn, const double *dat, int n1, int n2,
    float d1, float d2, float o1, float o2) {
  sf_file f = sf_output(fn);
  sf_putint(f, "n1", n1);
  sf_putint(f, "n2", n2);
  sf_putfloat(f, "d1", d1);
  sf_putfloat(f, "d2", d2);
  sf_putfloat(f, "o1", o1);
  sf_putfloat(f, "o2", o2);

  std::vector<float> tmp(dat, dat + n1 * n2);
  sf_floatwrite(&tmp[0], n1 * n2, f);
}

void puthead4x(sf_file Fo, int n1, int n2, int n3, int n4, float d1, float d2, float d3, float d4, float o1, float o2, float o3, float o4)
/*<  put head for (y,x,z,t) domain float-type 3D data sets>*/
{
        /* Read/Write axes */
        sf_putint(Fo,"n1",n1);
        sf_putint(Fo,"n2",n2);
        sf_putint(Fo,"n3",n3);
        sf_putint(Fo,"n4",n4);
        sf_putfloat(Fo,"d1",d1);
        sf_putfloat(Fo,"d2",d2);
        sf_putfloat(Fo,"d3",d3);
        sf_putfloat(Fo,"d4",d4);
        sf_putfloat(Fo,"o1",o1);
        sf_putfloat(Fo,"o2",o2);
        sf_putfloat(Fo,"o3",o3);
        sf_putfloat(Fo,"o4",o4);
        sf_putstring(Fo,"label1","z");
        sf_putstring(Fo,"label2","x");
        sf_putstring(Fo,"label3","y");
        sf_putstring(Fo,"label4","t");
        sf_putstring(Fo,"unit1","km");
        sf_putstring(Fo,"unit2","km");
        sf_putstring(Fo,"unit3","km");
        sf_putstring(Fo,"unit4","s");
}

void FwiParams::putOutputParams() {
  sf_putint(vupdates, "n1", nz);
  sf_putint(vupdates, "n2", nx);
  sf_putfloat(vupdates, "d1", dz);
  sf_putfloat(vupdates, "d2", dx);
  sf_putstring(vupdates, "label1", "Depth");
  sf_putstring(vupdates, "label2", "Distance");
  sf_putstring(vupdates, "label3", "Iteration");
  sf_putint(vupdates, "n3", niter);
  sf_putint(vupdates, "d3", 1);
  sf_putint(vupdates, "o3", 1);
  sf_putint(grads, "n1", nz);
  sf_putint(grads, "n2", nx);
  sf_putint(grads, "n3", niter);
  sf_putfloat(grads, "d1", dz);
  sf_putfloat(grads, "d2", dx);
  sf_putint(grads, "d3", 1);
  sf_putint(grads, "o3", 1);
  sf_putstring(grads, "label1", "Depth");
  sf_putstring(grads, "label2", "Distance");
  sf_putstring(grads, "label3", "Iteration");
  sf_putint(illums, "n1", nz);
  sf_putint(illums, "n2", nx);
  sf_putfloat(illums, "d1", dz);
  sf_putfloat(illums, "d2", dx);
  sf_putint(illums, "n3", niter);
  sf_putint(illums, "d3", 1);
  sf_putint(illums, "o3", 1);
  sf_putint(objs, "n1", niter);
  sf_putint(objs, "n2", 1);
  sf_putfloat(objs, "d1", 1);
  sf_putfloat(objs, "o1", 1);
}

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[])
{
    float x;	    /* trace value       */
    int ns=512;	    /* number of samples */
    int ntr=8;	    /* number of traces  */
    int dtime=16;   /* one way time in samples through "ocean" layer */
    float rbot=0.8; /* reflection strength of "ocean" bottom         */
    int h=100;	    /* location in samples of two way reverb train 	 */
    float amp=0.2;  /* strength of reflector */
    int loc=170;    /* location of reflector on trace 1 in time samples */
    int dip=12;	    /* dip of reflector in time samples */
    int i, j, k, sgn;
    float *trace;
    sf_file traces;

    sf_init(argc,argv);
    traces = sf_output("out");
    sf_setformat(traces,"native_float");
    sf_putint(traces,"n1",ns);
    sf_putint(traces,"n2",ntr);

    sf_putfloat(traces,"o1",0);
    sf_putfloat(traces,"d1",0.004);
    sf_putfloat(traces,"o2",1);
    sf_putfloat(traces,"d2",1);

    sf_putstring(traces,"label1","Time");
    sf_putstring(traces,"unit1","s");
    sf_putstring(traces,"label2","Trace");

    trace = sf_floatalloc(ns);

    for (j = 0; j < ntr; j++) {
        for (i = 0; i < ns; i++) {
            if (i >= h && ((i-h) % dtime == 0)) {
                k = (i-h)/dtime;
                sgn = (ISODD(k) ? -1 : 1);
                x = sgn * (k+1) * powf(rbot, k);
            } else {
                x = 0.0f;
            }
            if (i == loc + j*dip) x += amp;
            trace[i] = x;
        }
        sf_floatwrite(trace,ns,traces);
    }

    exit(0);
}

void puthead2kx(sf_file Fo, int n1, int n2, float d1, float d2, float o1, float o2)
/*<  put head for (kx,kz) domain float-type 2D data sets>*/
{
        /* Read/Write axes */
        sf_putint(Fo,"n1",n1);
        sf_putint(Fo,"n2",n2);
        sf_putfloat(Fo,"d1",d1);
        sf_putfloat(Fo,"d2",d2);
        sf_putfloat(Fo,"o1",o1);
        sf_putfloat(Fo,"o2",o2);
        sf_putstring(Fo,"label1","kz");
        sf_putstring(Fo,"label2","kx");
        sf_putstring(Fo,"unit1","2*pi/m");
        sf_putstring(Fo,"unit2","2*pi/m");
}

int main(int argc, char* argv[])
{
    int  i1, i2, n1,n2,n3,n12;
    float  s=0, *pp1, *pp2;
    sf_file inp1, inp2,  dif;

    sf_init(argc,argv);
    inp1 = sf_input("in");
    inp2 = sf_input("match");
    dif = sf_output("out");
    
    sf_putfloat(dif,"o2",0);
    sf_putfloat(dif,"o1",0);
    sf_putint(dif,"n1",1);  
    sf_putint(dif,"n2",1);    
	
    if (!sf_histint(inp1,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(inp1,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(inp1,"n3",&n3)) n3=1;

    n2 = sf_leftsize(inp1,1); /* left dimensions after the first one */
    n12 = n1*n2;
    
    if (n3!=1) {sf_putint(dif,"n3",1); sf_putint(dif,"d3",1); }
    
    pp1 = sf_floatalloc(n12);
    pp2 = sf_floatalloc(n12);
    sf_floatread(pp1,n12,inp1);
    sf_floatread(pp2,n12,inp2);

    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    s = s + pow((pp1[i2*n1+i1]-pp2[i2*n1+i1]),2);
	}
    }
    
    sf_warning("The difference is %f", s );
    sf_floatwrite(&s,1,dif);
    exit(0);
}

int main(int argc,char*argv[])
{
	sf_file in, out;
	int i1, i2, j2, i3, n1, n2, n3, nf;
	float **u1, **u2, sigma, d2;

	sf_init(argc,argv);

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

	n3 = sf_leftsize(in, 2);

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

	if (!sf_getfloat("sigma", &sigma)) sigma=1.0;
	/* sigma */
	if (!sf_getint("nf", &nf)) nf=100;
	/* frequency samples [0, 0.5] */

	d2 = SF_PI/(nf-1);

	sf_putint(out, "n2", nf);
	sf_putfloat(out, "d2", 0.5/(nf-1));
	dgauss_init(n1, 1.0/sigma);

	u1 = sf_floatalloc2(n1, n2);
	u2 = sf_floatalloc2(n1, nf);

	for (i3=0; i3<n3; i3++)
	{
		sf_floatread(u1[0], n1*n2, in);
		for(i1=0; i1<n1; i1++)
		{
			for (j2=0; j2<nf; j2++)
			for (i2=0; i2<n2; i2++)
			{
				u2[j2][i1] = u1[i2][i1] * dgauss_val(i2, d2*j2) 
					* powf(0.5/SF_PI, i2);
			}
		}
		sf_floatwrite(u2[0], n1*nf, out);
	}

	free(u1[0]);
	free(u2[0]);
	free(u1);
	free(u2);
	dgauss_close();
	exit(0);
}

int main(int argc, char* argv[])
{
    int nx, nz, ix, iz;
    float* trace;
    float dx, dz, x[2];
    sf_file mod;

    sf_init (argc,argv);
    mod = sf_output("out");

    if (!sf_getint("nx",&nx)) nx=400; 
    /* horizontal dimension */
    if (!sf_getint("nz",&nz)) nz=800; 
    /* vertical dimension */

    dx = 1./(nx-1);
    dz = 2./(nz-1);

    sf_putint   (mod,"n1",nz); 
    sf_putfloat (mod,"d1",dz); 
    sf_putfloat (mod,"o1",0.);
    sf_putint   (mod,"n2",nx); 
    sf_putfloat (mod,"d2",dx); 
    sf_putfloat (mod,"o2",0.);
    sf_setformat (mod,"native_float");

    trace = sf_floatalloc(nz);

    for (ix = 0; ix < nx; ix++) {
	x[1] = ix*dx;
	for (iz = 0; iz < nz; iz++) {
	    x[0] = iz*dz;
	    trace[iz] = 1./sqrtf(symes_vel(NULL,x));
	}
	sf_floatwrite(trace,nz,mod);
    }

    exit (0);
}

void write3dfdm(const char *file, float ***dat, fdm3d fdm) 
{
  sf_file f = sf_output(file);
  sf_putint(f, "n1", fdm->nzpad); sf_putfloat(f, "o1", fdm->ozpad); sf_putfloat(f, "d1", fdm->dz);
  sf_putint(f, "n2", fdm->nxpad); sf_putfloat(f, "o2", fdm->oxpad); sf_putfloat(f, "d2", fdm->dx);
  sf_putint(f, "n3", fdm->nypad); sf_putfloat(f, "o3", fdm->oypad); sf_putfloat(f, "d3", fdm->dy);
  sf_floatwrite(dat[0][0], fdm->nzpad*fdm->nxpad*fdm->nypad, f);
}

void write3df(const char *file, float ***dat, int n1, int n2, int n3) 
{
  sf_file f = sf_output(file);
  sf_putint(f, "n1", n1); sf_putfloat(f, "o1", 0); sf_putfloat(f, "d1", 0);
  sf_putint(f, "n2", n2); sf_putfloat(f, "o2", 0); sf_putfloat(f, "d2", 0);
  sf_putint(f, "n3", n3); sf_putfloat(f, "o3", 0); sf_putfloat(f, "d3", 0);
  sf_floatwrite(dat[0][0], n1*n2*n3, f);
}