ktmig.c,v

seismic software,very useful

           aperl, angmax, tapecntl);
   fprintf(jpfp, " nlcore=%d mtrace=%d ntras=%d mlimit=%d \n", nlcore,
           mtrace, ntras, mlimit / (1024 * 1024));
   fprintf(jpfp, " f0=%g df=%g nf=%d ftaper=%g ksmax=%g klmax=%g\n", f0,
           df, nf, ftaper, ksmax, klmax);
   fprintf(jpfp, " t0v=%g dtv=%g ntv=%d \n", t0v, dtv, ntv);
   fprintf(jpfp, " s0v=%g dsv=%g nsv=%d \n", s0v, dsv, nsv);
   fprintf(jpfp, " l0v=%g dlv=%g nlv=%d \n", l0v, dlv, nlv);
   fprintf(jpfp, " diskimg=%s\n", diskimg);
   fprintf(jpfp, " diskfld=%s\n", diskfld);
   fprintf(jpfp, " diskhdr=%s\n", diskhdr);
   fprintf(jpfp, " backupi=%s\n", backupi);
   fprintf(jpfp, " backupo=%s\n", backupo);
   fprintf(jpfp, " hisfile=%s\n", hisfile);
   fprintf(jpfp, " ntrdsk=%d tpowaf=%g \n", ntrdsk, tpowaf);
   fprintf(jpfp, " ncpu=%d traceout=%d \n", ncpu, traceout);
   if (on2trace == 1)
      fprintf(jpfp, " total_number_of_input_trace=%d \n", ntotal);
   if (ikey == 1) {
      fprintf(jpfp, " tr1=%g tr2=%g tr3=%g \n", tr1, tr2, tr3);
      fprintf(jpfp, " ln1=%g ln2=%g ln3=%g \n", ln1, ln2, ln3);
      fprintf(jpfp, " tracekey=%s linekey=%s \n", tracekey, linekey);
      fprintf(jpfp, " strace=%d jthtrace=%d sline=%d jthline=%d \n",
              strace, jthtrace, sline, jthline);
      fprintf(jpfp, " mintrace=%d maxtrace=%d \n", mintrace, maxtrace);
      fprintf(jpfp, " minline=%d maxline=%d \n", minline, maxline);
   }



   fprintf(jpfp, " \n");


   fprintf(jpfp, " start migration ... \n");
   fflush(jpfp);

   t0v = t0v * 0.001;
   dtv = dtv * 0.001;

   m1000 = ntrpre / 1000;
   ntrace = 0;

   if (tapecntl == 1) {
      if (jtrace > 0) {
         lpos = nt * sizeof(float) + 240;

         lpos = lpos * jtrace + 3600;
         fseek2g(infp, lpos, 0);
         fgettr(infp, &tra);
      }
      ltrace = jtrace;
      jtrace = ntrpre;
   } else if (tapecntl == -1) {
      ltrace = 0;
      jtrace = ntrpre;
      lpos = nt * sizeof(float) + 240;

      lpos = lpos * ntrpre + 3600;
      fseek2g(infp, lpos, 0);
      fgettr(infp, &tra);
   } else if (tapecntl == 0) {
      ltrace = jtrace;
      jtrace = ntrpre - jtrace;
   } else if (tapecntl == 2) {
      ltrace = jtrace;
      jtrace = ntrpre;
   }

   i3d = 3;
   if (nlines < 10)
      i3d = 2;

   /* initialization of fft and taper for antialiasing filter */
   fhigint_(&f0, &df, &nf, &ftaper, &dtl, ifcut, ltaper, wsave, &nsave,
            &nfft);

   /* read in velocity file if ivdisk=0 */
   if (ivdisk == 0) {
      fseek2g(velfp, 0, 0);
      efread(velos, sizeof(float), ntv * nsv * nlv, velfp);
   }

   do {
      jtrace = jtrace + 1;
      if (jtrace <= ntrpre)
         continue;
      ntrace = ntrace + 1;

      /* dead trace */
      if (tra.trid != 1) {
         fprintf(btfp, " trid not 1, at trace=%d cdp=%d cdpt=%d\n",
                 jtrace, tra.cdp, tra.cdpt);
         fflush(btfp);
         continue;
      }

      /* trace select */
      if (ikey == 1) {
         gethval(&tra, indxtrk, &trkval);
         inttrace = vtoi(trktype, trkval);
         if ((inttrace - strace) % jthtrace != 0
             || inttrace < mintrace || inttrace > maxtrace)
            continue;

         gethval(&tra, indxlnk, &lnkval);
         intline = vtoi(lnktype, lnkval);
         if ((intline - sline) % jthline != 0
             || intline < minline || intline > maxline)
            continue;
      }

      /* compute rms and max of trace amplitudes */
      itmp = (tra.mute - tra.delrt);
      itmp = itmp * 1000 / (int) tra.dt;

      ampmax = 0.;
      amprms = 0.;
      if (itmp < 0)
         itmp = 0;
      it0 = 0;

      for (it = itmp; it < nt; it++) {
         tmp = fabs(tra.data[it]);
         if (it0 == 0 && tmp > 0.)
            it0 = it;
         amprms = amprms + tmp * tmp;
         if (tmp > ampmax)
            ampmax = tmp;
      }
      itmp = it0;

      it0 = tra.delrt + itmp * (int) tra.dt / 1000;
      if (it0 > tra.mute)
         tra.mute = it0;
      /* set the mute time to be 10 samples shallower */
      /* just to be save */
      tra.mute = tra.mute - 10 * tra.dt / 1000;
      if (tra.mute < 0)
         tra.mute = 0;

      if (nt > itmp)
         amprms = sqrt(amprms / (nt - itmp));
      if (amprms == 0.) {
         fprintf(btfp, " zero trace, at trace=%d cdp=%d cdpt=%d\n",
                 jtrace, tra.cdp, tra.cdpt);
         fflush(btfp);
         continue;
      } else {
         if (ampmax > amprms * rmsmar) {
            fprintf(btfp,
                    " bad trace with ampmax=%g  rms=%g, at trace=%d cdp=%d cdpt=%d\n",
                    ampmax, amprms, jtrace, tra.cdp, tra.cdpt);
            fflush(btfp);
            continue;
         }
      }


      /* source and group x, y locations */
      xs = tra.sx;
      ys = tra.sy;
      xg = tra.gx;
      yg = tra.gy;

      if (tra.scalco > 1) {
         xs = xs * tra.scalco;
         ys = ys * tra.scalco;
         xg = xg * tra.scalco;
         yg = yg * tra.scalco;
      } else if (tra.scalco < 0) {
         xs = xs / (-tra.scalco);
         ys = ys / (-tra.scalco);
         xg = xg / (-tra.scalco);
         yg = yg / (-tra.scalco);
      }

      /* source and group s, l locations */
      xy2sl(x1, y1, s1, l1, x2, y2, s2, l2, x3, y3, s3, l3, xs, ys,
            &ss[ktrace], &sl[ktrace]);

      xy2sl(x1, y1, s1, l1, x2, y2, s2, l2, x3, y3, s3, l3, xg, yg,
            &gs[ktrace], &gl[ktrace]);


      ms = (ss[ktrace] + gs[ktrace]) / 2.;
      ml = (sl[ktrace] + gl[ktrace]) / 2.;

      fprintf(stderr, "ml=%g ms=%g cdp=%d offset=%d \n",
              ml, ms, tra.cdp, tra.offset);

      if (ntrend == 0) {

         /* 
            itmp = tra.cdp - cdp1; if(cdpnum==0) { itmp = itmp/ncdppl; iiy =
            itmp; iix = tra.cdp - cdp1 - iiy*ncdppl; } else { itmp =
            itmp/nlines; iix = itmp; iiy = tra.cdp - cdp1 - iix*nlines; } */

         /* replaced by the following */
         tmp = (ml - l1) / ddl + 0.5;
         iiy = tmp;
         tmp = (ms - s1) / dds + 0.5;
         iix = tmp;

         is = (iix - is0) / nss;
         ids = iix - is * nss - is0;

         if (ids >= 0 && ids < nds) {
            is = is * nds + ids;
         } else {
            is = -1;
         }

         il = (iiy - il0) / nll;
         idl = iiy - il * nll - il0;

         if (idl >= 0 && idl < ndl) {
            il = il * ndl + idl;
         } else {
            il = -1;
         }

      } else {
         is = 0;
         il = 0;
      }

      tmp = tra.offset;
      if (iofoin == 0 || nofo == 1) {
         if (tmp < 0.)
            tmp = -tmp;
         tmp = (tmp - ofomin) / dofo + .5;
         io = tmp;
      } else {
         if (tmp < ofol) {
            io = -1;
         } else if (tmp > ofor) {
            io = nofo;
         } else {
            bisear_(&nofo, &one, ofo, &tmp, &io);
            io = io - 1;
            if (io < nofo - 1) {
               if (abs(tmp - ofo[io]) > abs(tmp - ofo[io + 1]))
                  io = io + 1;
            }
         }
      }

      if (io >= 0 && io <= nofo - 1) {
         tmp = tra.offset;
         if (iofoin == 0 || nofo == 1) {
            if (tmp < 0.)
               tmp = -tmp;
         }
         tmp = tmp - ofo[io];
         if (tmp < 0.)
            tmp = -tmp;
         if (tmp > obtol * 1.001)
            io = -1;
      }

      fprintf(stderr,
              "il=%d is=%d ml=%g ms=%g cdp=%d offset=%d io=%d\n", il,
              is, ml, ms, tra.cdp, tra.offset, io);

      if (io < 0 || io >= nofo)
         continue;

      if (ms < sinmin || ms > sinmax)
         continue;
      if (ml < linmin || ml > linmax)
         continue;

      iofs[ktrace] = io + 1;

      /* apply rho filter */
      f2p5n_(tra.data, &nt, &i3d, &mm, &ll);

      /* apply tpow */
      if (itgain == 1)
         for (it = 0; it < nt; it++)
            tra.data[it] *= tgain[it];
      for (it = 0; it < nt; it++)
         tras[it + ktrace * nt] = tra.data[it];


      itmp = (tra.mute - tra.delrt);
      itmp = itmp * 1000 / (int) tra.dt * ifact + 1;
      if (itmp < 1)
         itmp = 1;
      if (itmp > ntl)
         itmp = ntl;
      imutel[ktrace] = itmp;



      if (il >= 0 && il < nlout && is >= 0 && is < nsout) {
         if (ntrend == 0) {
            if (abs(tra.offset) <
                off[io + is * nofo + il * nofo * nsout]) {
               off[io + is * nofo + il * nofo * nsout] =
                     abs(tra.offset);
               lpos = io + is * nofo + il * nofo * nsout;
               lpos = lpos * 240;
               fseek2g(hdrfp, lpos, 0);
               tmp = tra.mute * 0.001;
               tmp = tmp * tmp - tra.offset * tra.offset / gmin / gmin;
               if (tmp >= 0.) {
                  tmp = sqrt(tmp);
                  tmp = tmp * 1000.;
                  itmp = tra.mute;
                  tra.mute = tmp;
               }
               efwrite(&tra, sizeof(char), 240, hdrfp);

               if (tmp >= 0.)
                  tra.mute = itmp;
            }
         } else {
            if (lend == lstart) {
               ps = ms;
               pl = lstart;
            } else {
               tmp = (send - sstart) / (lend - lstart);
               pl = ml + tmp * (lstart * tmp + ms - sstart);
               pl = pl / (1. + tmp * tmp);
               ps = tmp * (pl - lstart) + sstart;
            }
            if (dsout != 0.) {
               tmp = (ps - sstart) / dsout + 0.5;
               isend = tmp;
            } else {
               tmp = (pl - lstart) / dlout + 0.5;
               isend = tmp;
            }
            if (isend >= 0 && isend < ntrend) {
               if (abs(tra.offset) < off[io + isend * nofo]) {
                  tmp = sqrt((ms - s[isend]) * (ms - s[isend])
                             + (ml - l[isend]) * (ml - l[isend]));
                  if (tmp < disend[io + isend * nofo]) {
                     disend[io + isend * nofo] = tmp;