sutzconv.c

seismic software,very useful

         warn(" tr.fz equals 0; reset to tr.delrt=%f \n", fz);
      }
      if (tr.dz == 0.) {
         dz = tr.dt * 0.001;
         warn(" tr.dz equals 0; reset to tr.dt/1000=%f \n", dz);
      }
      nt = nout;
      dt = dout;
      ft = fout;
      for (iz = 0; iz < n1; iz++)
         zv[iz] = o1 + iz * d1;
      for (it = 0; it < nout; it++)
         to[it] = fout + it * dout;
   }


   fprintf(stderr, " \n");
   fprintf(stderr, " sutzconv parameters \n");
   fprintf(stderr, " =================== \n");
   fprintf(stderr, " tracekey=%s linekey=%s \n", tracekey, linekey);
   fprintf(stderr, " vgrid=%s \n", vgrid);
   fprintf(stderr, " torz=%d vtorz=%d \n", torz, vtorz);
   fprintf(stderr, " fout=%g dout=%g nout=%d \n", fout, dout, nout);
   fprintf(stderr, " orient=%d \n", orient);
   fprintf(stderr, " vgrid: ocdp2=%g dcdp2=%g nx=%d \n", x0, dx, nx);
   fprintf(stderr, " vgrid: oline3=%g dline3=%g ny=%d \n", y0, dy, ny);
   fprintf(stderr, " vgrid: o1=%g d1=%g nt/nz=%d \n", o1, d1, n1);

   fseek2g(vgfp, 0, 0);
   if (orient == 1) {
      if (idisk == 0)
         efread(vs, sizeof(float), n1 * nx * ny, vgfp);
   } else {
      for (i1 = 0; i1 < n1; i1++) {
         efread(vread, sizeof(float), nx * ny, vgfp);

         for (it = 0; it < nx * ny; it++)
            vs[i1 + it * n1] = vread[it];
      }
   }

   if (orient == 4)
      free(vread);

   if (idisk == 0)
      fprintf(stderr, " velocity function read \n");

   gethval(&tr, indxtrk, &trkval);
   ix = vtoi(trktype, trkval);
   gethval(&tr, indxlnk, &lnkval);
   iy = vtoi(lnktype, lnkval);

   xpre = ix - 100;
   ypre = iy - 100;
   change = 0;
   tzmax = 9999999;

   /* loop over output traces */
   do {
      /* compute velocity at output location */
      /* via bilinear interpolation */
      gethval(&tr, indxtrk, &trkval);
      ix = vtoi(trktype, trkval);
      gethval(&tr, indxlnk, &lnkval);
      iy = vtoi(lnktype, lnkval);

      x = ix;
      y = iy;
      if (x != xpre || y != ypre) {
         change = 1;
         xpre = x;
         ypre = y;
      } else {
         change = 0;
      }
      if (change == 1) {

         if (idisk == 0) {
            bilint_(&n1, &nx, &ny, &x0, &y0, &dx, &dy, &x, &y, vs, v);
         } else {
            tmp = (x - x0) / dx + 0.5;
            ix = tmp;
            if (ix < 0)
               ix = 0;
            if (ix > nx - 1)
               ix = nx - 1;
            tmp = (y - y0) / dy + 0.5;
            iy = tmp;
            if (iy < 0)
               iy = 0;
            if (iy > ny - 1)
               iy = ny - 1;
            lpos = iy * nx + ix;
            lpos = lpos * n1 * sizeof(float);

            fseek2g(vgfp, lpos, 0);
            efread(v, sizeof(float), n1, vgfp);

            if (iprint == 1)
               fprintf(stderr,
                       "line=%g trace=%g i3vgrid=%d i2vgrid=%d \n", y,
                       x, iy + 1, ix + 1);
         }
         if (torz == 0) {
            if (vgtype == 1) {
               if (vtorz == 0) {
                  zv[0] = o1 * v[0] * 0.0005;
                  for (i1 = 1; i1 < n1; i1++)
                     zv[i1] =
                           zv[i1 - 1] + v[i1] * d1 * 0.0005;
               } else if (vtorz == 1) {
                  for (i1 = 0; i1 < n1; i1++) {
                     zv[i1] = o1 + i1 * d1;
                  }
                  tv[0] = o1 / v[0] * 2000.;
                  for (i1 = 1; i1 < n1; i1++) {
                     if( v[i1] != v[i1-1] ){
                        tv[i1] = tv[i1 - 1] + log(v[i1]/v[i1-1])*d1 / (v[i1]-v[i1-1]) * 2000.;
                     }else{
                        tv[i1] = tv[i1-1] + d1/v[i1]*2000.0;
                     }
                  }
               }
            } else {
               for (i1 = 0; i1 < n1; i1++)
                  zv[i1] = (o1 + i1 * d1) * v[i1] * 0.0005;
            }
            dzdt0 = (zv[1] - zv[0]) / (tv[1] - tv[0]);
            dzdtn =
                  (zv[n1 - 1] - zv[n1 - 2]) / (tv[n1 - 1] - tv[n1 - 2]);
            lin1dn_(zv, tv, &n1, zo, to, &nz, indx, &dzdt0, &dzdtn);
         } else {
            if (vgtype == 1) {
               if (vtorz == 1) {
                  tv[0] = o1 / v[0] * 2000.;
                  for (i1 = 1; i1 < n1; i1++)
                     tv[i1] =
                           tv[i1 - 1] + d1 / v[i1] * 2000.;
               } else if (vtorz == 0) {
                  for (i1 = 0; i1 < n1; i1++)
                     tv[i1] = o1 + d1 * i1;
                  zv[0] = o1 * v[0] * 0.0005;
                  for (i1 = 1; i1 < n1; i1++)
                     zv[i1] =
                           zv[i1 - 1] + v[i1] * d1 * 0.0005;
               }
            } else {
               for (i1 = 0; i1 < n1; i1++)
                  tv[i1] = o1 + d1 * i1;
               for (i1 = 0; i1 < n1; i1++)
                  zv[i1] = (o1 + i1 * d1) * v[i1] * 0.0005;
            }
            dtdz0 = (tv[1] - tv[0]) / (zv[1] - zv[0]);
            dtdzn =
                  (tv[n1 - 1] - tv[n1 - 2]) / (zv[n1 - 1] - zv[n1 - 2]);
            lin1dn_(tv, zv, &n1, to, zo, &nt, indx, &dtdz0, &dtdzn);
         }
      }

      if (torz == 0) {
         /* linear interpolate trace */
         for (iz = 0; iz < nz; iz++) {
            ti = (to[iz] - ft) / dt;
            iti = ti;
            res = ti - iti;
            if (iti < 0) {
               tro.data[iz] = 0.;
            } else if (iti > nt - 1) {
               tro.data[iz] = 0.;
               if (tzmax > iz)
                  tzmax = iz;
            } else if (iti == nt - 1) {
               tro.data[iz] = tr.data[iti];
            } else {
               tro.data[iz] = tr.data[iti] * (1. - res) +
                     res * tr.data[iti + 1];
            }
         }
         bcopy(&tr, &tro, 240);
         tro.ns = nz;
         tro.dz = dz;
         tro.fz = fz;
         tro.dt = dz * 1000.;
         tro.delrt = fz;
      } else {
         /* linear interpolate trace */
         for (it = 0; it < nt; it++) {
            zi = (zo[it] - fz) / dz;
            izi = zi;
            res = zi - izi;
            if (izi < 0) {
               tro.data[it] = 0.;
            } else if (izi > nz - 1) {
               tro.data[it] = 0.;
               if (tzmax > it)
                  tzmax = it;
            } else if (izi == nz - 1) {
               tro.data[it] = tr.data[izi];
            } else {
               tro.data[it] = tr.data[izi] * (1. - res) +
                     res * tr.data[izi + 1];
            }
         }
         bcopy(&tr, &tro, 240);
         tro.ns = nt;
         tro.dt = dt * 1000.;
         tro.delrt = ft;
         tro.dz = 0.;
         tro.fz = 0.;
      }

      fputtr(outfp, &tro);

   } while (fgettr(infp, &tr));

   if (tzmax == 9999999.)
      tzmax = (nout - 1);
   tzmax = fout + tzmax * dout;
   fprintf(stderr, " ==========================\n");
   fprintf(stderr, "   minimum length of output live sample ===%g\n",
           tzmax);
   fprintf(stderr, " ==========================\n");

   free(to);
   free(zo);
   free(tv);
   free(zv);
   free(indx);
   free(vs);
   if (idisk == 0)
      free(v);

   return 0;

}