sutzconv.c

seismic software,very useful

#include <math.h>

#include "usgrid.h"
#include "su.h"
#include "segy.h"
#include "par.h"

static segytrace tr, tro;

char    *sdoc =
      "SUTZCONV - time-depth conversion                                 \n"
      "\n"
      "sutzconv vgrid= [parameters] <input-data >output-data            \n"
      "\n"
      "Required parameters:                                                     \n"
      "vgrid=                 name of interval or average velocity grid \n"
      "fout=                  time (ms) or depth of first sample to output      \n"
      "dout=                  time (ms) or depth interval to output             \n"
      "nout=                  number of samples per trace to output             \n"
      "\n"
      "Optional parameters:                                                     \n"
      "tracekey=tracl         segy key word defining trace number within line \n"
      "linekey=tracr          segy key word defining line number \n"
      "vgtype=1               velocity type (1=interval 0=time-average)         \n"
      "torz=0                 0 --- time to depth; input is time data \n"
      "                       1 --- depth to time; input is depth data \n"
      "vtorz=0                0 --- velocity is in time-domain \n"
      "                       1 --- velocity is in depth-domain \n"
      "ocdp2=from_vgrid       inline trace number of first velocity grid trace \n"
      "dcdp2=from_vgrid       inline trace number increment of velocity grid \n"
      "oline3=from_vgrid      line number of first velocity grid trace \n"
      "dline3=from_vgrid      line number increment of velocity grid \n"
      "orient=1               velocity grid orientation                         \n"
      "                       1=(time/depth,inline,xline)                       \n"
      "                       4=(inline,xline,time/depth)                       \n"
      "print=1                print velocity grid location (1=yes 0=no)       \n"
      "\n"
      "Note:                                                                    \n"
      "    1. the fold of the input 3D data must be constant and specified at \n"
      "       the binary header (the number of traces per cdp must be the same) \n"
      "       if fold=0 in the binary header, 1 is assumed (stack data) \n"
      "    2. ocdp2,dcdp2 in velocity grid must be specified correctly to       \n"
      "       represent the first trace number within line and trace number     \n"
      "       increment (inline) of velocity grid                               \n"
      "    3. oline3,dline3 in velocity grid must be specified correctly to     \n"
      "       represent the first line number and line number   \n"
      "       increment of velocity grid                                \n"
      "    4. bilinear interpolation is used to extra the velocity function     \n"
      "       at input trace location from velocity grid                        \n"
      "    5. when vtorz=0, velocity grid is in time with time unit of ms;\n"
      "       dt in the trace haeader will be used as time sampling interval of\n"
      "       input data                                                        \n"
      "    6. when vtorz=1, velocity grid is in depth with depth unit of m\n"
      "       or ft; dz (or dt/1000 if dz=0) in the trace haeader will be used \n"
      "       as depth sampling interval of input data. fz (or delrt if fz=0) \n"
      "       in the trace header will be used as first sample depth            \n"
      "    7. vgrid consists of nx by ny traces, where nx is the number of      \n"
      "       inline traces and ny is number of lines                           \n"
      "\n"

      "AUTHOR:          Zhiming Li,       ,     12/13/94                \n";

void     bilint_(int *n1, int *nx, int *ny, float *x0, float *y0,
                 float *dx, float *dy, float *x, float *y, float *vs,

                 float *v);

void     lin1dn_(float *xin, float *yin, int *nin, float *xo, float *yo,
                 int *nout, int *indx, float *dydx0, float *dydxn);

main(int argc, char **argv)
{
   segychdr ch;
   segybhdr bh;
   usghed   usgh;
   FILE    *infp = stdin, *outfp = stdout;
   FILE    *vgfp;
   char    *vgrid;

   float    x0, y0, dx, dy, x, y;
   int      nx, ny;
   float    fout, dout;
   int      nout;
   int      n1, torz, vtorz, orient, idisk;
   float    o1, d1;
   float    tmp;
   int      ix, iy, change;
   float    xpre, ypre;
   long     lpos;
   float    tzmax;

   int      i1, it, iz, iti, izi, nt, nz;
   float    dt, dz, ft, fz, zi, ti;
   float    res, dzdt0, dzdtn, dtdz0, dtdzn;
   float   *tv, *zv, *vs, *v, *to, *zo, *vread;
   int     *indx;

   int      vgtype;
   int      iprint;

   String   tracekey = "tracl", linekey = "tracr", trktype, lnktype;
   Value    trkval, lnkval;
   int      indxtrk, indxlnk, trstart, lnstart;


   /* get parameters */
   initargs(argc, argv);
   askdoc(1);

   if (!getparstring("vgrid", &vgrid))
      err("vgrid missing");
   if (!getparfloat("fout", &fout))
      err("fout missing");
   if (!getparfloat("dout", &dout))
      err("dout missing");
   if (!getparint("nout", &nout))
      err("nout missing");
   if (!getparint("orient", &orient))
      orient = 1;
   if (!getparint("torz", &torz))
      torz = 0;
   if (!getparint("vtorz", &vtorz))
      vtorz = 0;
   if (!getparint("vgtype", &vgtype))
      vgtype = 1;
   if (!getparint("print", &iprint))
      iprint = 1;
   if (torz == 1 && vgtype != 1)
      err(" depth to time requires vgtype=1");

   vgfp = efopen(vgrid, "r");
   file2g(vgfp);
   if (fgetusghdr(vgfp, &usgh) != 0)
      err("error fgetusghdr");

   if (!getparfloat("oline3", &y0))
      y0 = usgh.oline3;
   if (y0 == 0.)
      err("vghd.oline3 zero");
   if (!getparfloat("dline3", &dy))
      dy = usgh.dline3;
   if (dy == 0.)
      err("vghd.dline3 zero");
   if (!getparfloat("ocdp2", &x0))
      x0 = usgh.ocdp2;
   if (x0 == 0.)
      err("vghd.ocdp2 zero");
   if (!getparfloat("dcdp2", &dx))
      dx = usgh.dcdp2;
   if (dx == 0.)
      err("vghd.dcdp2 zero");
   if (orient == 1) {
      ny = usgh.n3;
      if (ny == 0) {
         warn(" usgh.n3 = 0; use 1 instead ");
         ny = 1;
      }
      nx = usgh.n2;
      n1 = usgh.n1;
      d1 = usgh.d1;
      o1 = usgh.o1;
   } else {
      ny = usgh.n2;
      nx = usgh.n1;
      n1 = usgh.n3;
      if (n1 == 0)
         err(" usgh.n1 = 0 ");
      d1 = usgh.d3;
      o1 = usgh.o3;
   }

   getparstring("tracekey", &tracekey);
   getparstring("linekey", &linekey);
   trktype = hdtype(tracekey);
   lnktype = hdtype(linekey);
   indxtrk = getindex(tracekey);
   indxlnk = getindex(linekey);

   if (orient != usgh.orient)
      warn(" orient in vgrid different: %d \n", usgh.orient);

   /* make file size to be able to exceed 2 G on convex */
   file2g(infp);
   file2g(outfp);

   fgethdr(infp, &ch, &bh);
   if (bh.fold == 0)
      warn("fold in binary is 0; reset to 1 \n");
   if (bh.fold == 0)
      bh.fold = 1;
   bh.hns = nout;
   bh.hdt = dout * 1000.;
   // fputhdr(outfp, &ch, &bh);

   fseeko(infp, 0, SEEK_SET);

   /* read in first trace for nt and dt */
   if (!fgettr(infp, &tr))
      err("can't get first trace");

   idisk = 0;
   /* if(orient==1 && dx==1. && dy==1. ) idisk = 1; */
   if (orient == 1)
      idisk = 1;
   idisk = 0;

   /* memory allocations */
   if (orient == 4)
      vread = (float *) emalloc(nx * ny * sizeof(float));

   if (idisk == 0)
      vs = (float *) emalloc(n1 * nx * ny * sizeof(float));
   v = (float *) emalloc(n1 * sizeof(float));
   tv = (float *) emalloc(n1 * sizeof(float));
   zv = (float *) emalloc(n1 * sizeof(float));
   to = (float *) emalloc(nout * sizeof(float));
   zo = (float *) emalloc(nout * sizeof(float));
   indx = (int *) emalloc(nout * sizeof(float));

   if (torz == 0) {
      nt = tr.ns;
      dt = tr.dt * 0.001;
      ft = tr.delrt;
      nz = nout;
      dz = dout;
      fz = fout;
      for (it = 0; it < n1; it++)
         tv[it] = o1 + it * d1;
      for (iz = 0; iz < nout; iz++)
         zo[iz] = fout + iz * dout;
   } else {
      nz = tr.ns;
      dz = tr.dz;
      fz = tr.fz;
      if (tr.fz == 0.) {
         fz = tr.delrt;