ktmig.c,v

seismic software,very useful

      if (jtrace > 0) {
         fprintf(jpfp,
                 "Number of Traces Processed for this input: %d\n",
                 jtrace);
      }
   } else {
      hffp = efopen(hisfile, "w");
      ntrpre = 0;
   }


   /* compute memory requirement for migration */

   ntmp = nl * ndl;
   nl = ntmp;
   ntmp = ns * nds;
   ns = ntmp;

   nlcore = nl;

   incdxw = 10;
   nindxw = ntau / incdxw + 1;
   ires = 4;
   ntres = nt * ires + 1;


   if (ntrend == 0) {
      nsout = ns;
      nlout = nl;
   } else {
      nsout = ntrend;
      nlout = 1;
      nlcore = 1;
      dsout = (send - sstart) / (nsout - 1);
      dlout = (lend - lstart) / (nsout - 1);
      dend = sqrt((send - sstart) * (send - sstart)
                  + (lend - lstart) * (lend - lstart)) / (nsout - 1);
   }

   lwork = nsout * nlcore;
   lwork = lwork * ntau * nofo * sizeof(float);

   ltmp = nsout * nlout * nofo + nt + ntau + ntl + nt * mtrace;
   ltmp =
         ltmp + 2 * ntau + 2 * ntau * mtrace + 2 * nsout * nlout + ntv +
         mtrace * 4;
   ltmp = ltmp + ntau + nfft * (nf + nsout * nlout) + ntl * nf + nsave;
   ltmp = ltmp + nsout * nofo + ntres * nf + ntau;
   ltmp = ltmp + ntau * nsout * nlout * 2 + ntau * mtrace + ntau * 3;
   ltmp = ltmp * sizeof(float);

   ltmp =
         ltmp + (nsout * nlout * nofo + mtrace * 2 + nf * 2 +
                 nindxw * nsout * nlout) * sizeof(int);

   lwork = lwork + ltmp;
   ltmp = ntv * nsv * nlv * sizeof(float);

   lwork = lwork + ltmp;

   ivdisk = 0;
   if (lwork > mlimit) {
      warn("Velocity Disk I/O Used to Reduce Memory \n");
      lwork = lwork - ltmp;
      ivdisk = 1;
   }

   if (lwork > mlimit)
      err("mlimit too small ");

   fprintf(jpfp, " Total Memory used (in Byte) = %g \n", (float) lwork);
   fprintf(jpfp, "\n");

   /* allocation of memory */
   mig =

         (float *) emalloc(nsout * nlcore * nofo * ntau *
                           sizeof(float));
   fold = (float *) emalloc(nsout * nlout * nofo * sizeof(float));
   off = (int *) emalloc(nsout * nlout * nofo * sizeof(int));
   tgain = (float *) emalloc(nt * sizeof(float));
   tugain = (float *) emalloc(ntau * sizeof(float));
   trace = (float *) emalloc(ntl * sizeof(float));
   tras = (float *) emalloc(nt * mtrace * sizeof(float));
   imutel = (int *) emalloc(mtrace * sizeof(int));
   tm = (float *) emalloc(ntau * sizeof(float));
   tmig = (float *) emalloc(ntau * sizeof(float));
   fovt2 = (float *) emalloc(ntau * mtrace * sizeof(float));
   s = (float *) emalloc(nsout * nlout * sizeof(float));
   l = (float *) emalloc(nsout * nlout * sizeof(float));
   vel = (float *) emalloc(ntv * sizeof(float));
   ss = (float *) emalloc(mtrace * sizeof(float));
   sl = (float *) emalloc(mtrace * sizeof(float));
   gs = (float *) emalloc(mtrace * sizeof(float));
   gl = (float *) emalloc(mtrace * sizeof(float));
   iofs = (int *) emalloc(mtrace * sizeof(int));
   w1 = (float *) emalloc(ntau * sizeof(float));
   work1 = (float *) emalloc(nfft * nf * sizeof(float));
   tracef = (float *) emalloc(ntl * nf * sizeof(float));
   work2 = (float *) emalloc(nfft * nsout * nlout * sizeof(float));
   wsave = (float *) emalloc(nsave * sizeof(float));
   ifcut = (int *) emalloc(nf * sizeof(int));
   ltaper = (int *) emalloc(nf * sizeof(int));
   disend = (float *) emalloc(nsout * nofo * sizeof(float));
   resamp = (float *) emalloc(ntres * nf * sizeof(float));
   indxw = (int *) emalloc(nindxw * nsout * nlout * sizeof(int));

   ss2 = (float *) emalloc(ntau * sizeof(float));
   scg = (float *) emalloc(ntau * nsout * nlout * sizeof(float));
   scs = (float *) emalloc(ntau * nsout * nlout * sizeof(float));
   vr2 = (float *) emalloc(ntau * mtrace * sizeof(float));
   vi2 = (float *) emalloc(ntau * sizeof(float));
   tau = (float *) emalloc(ntau * sizeof(float));
   vq4 = (float *) emalloc(ntau * sizeof(float));

   if (ivdisk == 0)
      velos = (float *) emalloc(ntv * nsv * nlv * sizeof(float));

   itrace = 0;

   bzero(fold, nofo * nsout * nlout * sizeof(float));
   bzero(mig, nofo * nsout * nlcore * ntau * sizeof(float));

   if (ntrend > 0) {
      for (io = 0; io < nsout * nofo; io++)
         disend[io] = 99999999999.;
   }

   if (ntrend == 0) {
      for (il = 0; il < nlout; il++) {
         for (is = 0; is < nsout; is++) {
            s[is + il * nsout] =
                  sstart + is / nds * ds + is % nds * dds;
         }
      }
      for (il = 0; il < nlout; il++) {
         for (is = 0; is < nsout; is++) {
            l[is + il * nsout] =
                  lstart + il / ndl * dl + il % ndl * ddl;
         }
      }
   } else {
      fprintf(jpfp, " Output s and l Locations of Target Line\n");
      for (is = 0; is < ntrend; is++) {
         s[is] = sstart + is * (send - sstart) / (ntrend - 1.);
         l[is] = lstart + is * (lend - lstart) / (ntrend - 1.);
         fprintf(jpfp,
                 "   s=%g l=%g trace=%d \n", s[is], l[is], is + 1);
      }
   }

   sinmin = s[0];
   sinmax = s[0];
   linmin = l[0];
   linmax = l[0];

   for (is = 0; is < nsout * nlout; is++) {
      if (sinmin > s[is])
         sinmin = s[is];
      if (sinmax < s[is])
         sinmax = s[is];
      if (linmin > l[is])
         linmin = l[is];
      if (linmax < l[is])
         linmax = l[is];
   }

   sinmin = sinmin - apers;
   sinmax = sinmax + apers;
   linmin = linmin - aperl;
   linmax = linmax + aperl;

   /* allocate disk space for image */
   if (ntrpre > 0) {

      fprintf(jpfp, " open diskimg file... \n");
      imgfp = efopen(diskimg, "r+");
      fseek2g(imgfp, 0, 1);
      ltmp = ntau * nofo * nlout * nsout;
      efread(mig, sizeof(float), ltmp, imgfp);

      fprintf(jpfp, " diskimg file opened  \n");

      fprintf(jpfp, " open diskfld file... \n");
      fldfp = efopen(diskfld, "r+");
      fseek2g(fldfp, 0, 1);
      efread(fold, sizeof(float), nofo * nsout * nlout, fldfp);

      efclose(fldfp);
      fprintf(jpfp, " diskfld file opened  \n");

      fprintf(jpfp, " open diskhdr file... \n");
      hdrfp = efopen(diskhdr, "r+");
      fseek2g(hdrfp, 0, 1);
      fprintf(jpfp, " diskhdr file opened  \n");

      for (ix = 0; ix < nsout * nlout * nofo; ix++) {
         efread(&tr, sizeof(char), 240, hdrfp);

         off[ix] = tr.offset;
         if (ntrend > 0)
            disend[ix] = tr.dz;
      }
      ipre = 1;
   } else if (isave == 1) {
      imgfp = efopen(diskimg, "w+r");
      fseek2g(imgfp, 0, 1);
      bzero(tm, ntau * sizeof(float));

      fprintf(jpfp, " initialize diskimg file... \n");
      for (ix = 0; ix < nsout * nlout * nofo; ix++)
         efwrite(tm, sizeof(float), ntau, imgfp);

      fprintf(jpfp, " initialization of diskimg file done \n");
      if (nlcore == nlout)
         efclose(imgfp);

      fprintf(jpfp, " initialize diskfld file... \n");
      fldfp = efopen(diskfld, "w+r");
      fseek2g(fldfp, 0, 1);
      bzero(fold, nofo * nsout * nlout * sizeof(float));
      efwrite(fold, sizeof(float), nofo * nsout * nlout, fldfp);

      efclose(fldfp);
      fprintf(jpfp, " initialization of diskfld file done  \n");
   }

   if (ntrpre == 0) {
      fprintf(jpfp, " initializing diskhdr file... \n");
      hdrfp = efopen(diskhdr, "w");
      bzero(&tr, 240);
      tr.offset = 99999999;
      tr.ns = tra.ns;
      tr.dt = tra.dt;
      tr.trid = 2;
      if (ntrend > 0)
         tr.dz = 99999999999.;
      for (iy = 0; iy < nlout; iy++) {
         tr.tracr = iy + 1;
         for (ix = 0; ix < nsout; ix++) {
            tmp = (l[iy * nsout + ix] - l1) / ddl + 0.5;
            iiy = tmp;
            tmp = (s[iy * nsout + ix] - s1) / dds + 0.5;
            iix = tmp;
            if (cdpnum == 0) {
               tr.cdp = iiy * ncdppl + iix + cdp1;
            } else {
               tr.cdp = iix * nlines + iiy + cdp1;
            }
            tr.tracl = ix + 1;
            ms = s[iy * nsout + ix];
            ml = l[iy * nsout + ix];
            sl2xy(s1, l1, x1, y1, s2, l2, x2, y2, s3, l3, x3, y3, ms,
                  ml, &mx, &my);
            if (ikey == 1) {
               fline = (ml - l1) * (ln3 - ln1) / (l3 - l1) + ln1 + 0.5;
               ftrace = (ms - s1) * (tr2 - tr1) / (s2 - s1) + tr1 + 0.5;
               intline = fline;
               inttrace = ftrace;
               itov(trktype, &trkval, inttrace);
               itov(lnktype, &lnkval, intline);
               puthval(&tr, indxtrk, &trkval);
               puthval(&tr, indxlnk, &lnkval);
            }
            for (io = 0; io < nofo; io++) {
               tr.cdpt = io + 1;
               ofs = ofo[io];
               tr.sx = mx - ofs / 2;
               tr.gx = mx + ofs / 2;
               tr.sy = my;
               tr.gy = my;
               efwrite(&tr, sizeof(char), 240, hdrfp);
            }
         }
      }
      fprintf(jpfp, " initializing diskhdr file done \n");
      for (ix = 0; ix < nofo * nsout * nlout; ix++)
         off[ix] = 99999999;
      fflush(hdrfp);

   }


   ktrace = 0;
   cdppre = tra.cdp - 1;
   ispre = -9999;
   ilpre = -9999;
   m1000 = 0;
   for (it = 0; it < ntau; it++) {
      tmig[it] = tau0 + it * dtau;
      tm[it] = tmig[it] * 0.5 / dtl;
   }

   for (it = 0; it < ntau; it++)
      tau[it] = tau0 + (it + 1) * dtau;

   if (itgain == 1) {
      for (it = 0; it < nt; it++) {
         tmp = t0 + it * dt;
         if (tmp != 0.0) {
            tgain[it] = pow(tmp, tpow);
         } else {
            tgain[it] = 0.;
         }
      }
      for (it = 0; it < ntau; it++) {
         tmp = tau0 + it * dtau;
         if (tmp != 0.0) {
            tugain[it] = pow(tmp, tpowaf);
         } else {
            tugain[it] = 0.;
         }
      }
   }

   fprintf(jpfp, " \n");
   fprintf(jpfp, " Migration Parameters \n");
   fprintf(jpfp, " -------------------- \n");
   fprintf(jpfp, " x1=%g y1=%g s1=%g l1=%g cdp1=%d\n", x1, y1, s1, l1,
           cdp1);
   fprintf(jpfp, " x2=%g y2=%g s2=%g l2=%g cdp2=%d\n", x2, y2, s2, l2,
           cdp2);
   fprintf(jpfp, " x3=%g y3=%g s3=%g l3=%g cdp3=%d\n", x3, y3, s3, l3,
           cdp3);
   fprintf(jpfp, " sstart=%g lstart=%g ns=%d nl=%d ds=%g dl=%g \n",
           sstart, lstart, nsread, nlread, ds, dl);
   if (ntrend > 0)
      fprintf(jpfp, " send=%g lend=%g ntrend=%d dend=%g \n", send, lend,
              ntrend, dend);
   fprintf(jpfp, " nds=%d ndl=%d dds=%g ddl=%g \n", nds, ndl, dds, ddl);
   fprintf(jpfp, " nss=%d nll=%d is0=%d il0=%d \n", nss, nll, is0, il0);
   fprintf(jpfp,
           " cdppdds=%d cdppds=%d cdppddl=%d cdppdl=%d cdpnum=%d\n",
           cdppdds, cdppds, cdppddl, cdppdl, cdpnum);
   fprintf(jpfp, " ntrpre=%d nt=%d ntau=%d ntl=%d nfft=%d\n", ntrpre,
           nt, ntau, ntl, nfft);
   fprintf(jpfp, " dt=%f dtau=%f t0=%f tau0=%f \n", dt * 1000.,
           dtau * 1000., t0 * 1000., tau0 * 1000.);
   fprintf(jpfp, " tpow=%g ofomin=%g dofo=%g nofo=%d isave=%d \n", tpow,
           ofomin, dofo, nofo, isave);
   fprintf(jpfp, " apers=%g aperl=%g angmax=%g tapecntl=%d\n", apers,