vgrid3d.c

seismic software,very useful

		  vi[j0 + j] = sqrt(work[j]);
	       }
	    }
	 } else {
	    if (vitype == 0) {
	       vv[0] = vpicks[i * n1];
	       for (j = 1; j < np; j++) {
		  vv[j] = (vpicks[i * n1 + j] * vpicks[i * n1 + j] *
			   tpicks[i * n1 + j] -
			   vpicks[i * n1 + j - 1] * vpicks[i * n1 + j -
							   1] *
			   tpicks[i * n1 + j - 1]) / (tpicks[i * n1 +
							     j] -
						      tpicks[i * n1 +
							     j - 1]);
		  if (vv[j] < 0.)
		     err
			   ("check interval velocity at t=%g x=%g y=%g location \n",
			    tpicks[i * n1 + j], xs[i], ys[i]);

		  vv[j] = sqrt(vv[j]);
	       }
	       bisear_(&np, &nvs, tpicks + i * n1, zs, indx);
	       for (j = 0; j < nvs; j++) {
		  jj = indx[j] - 1;
		  if (jj < 0 || zs[j] <= tpicks[i * n1]) {
		     vi[j0 + j] = vv[0];
		  } else if (jj > np - 2) {
		     vi[j0 + j] = vv[np - 1];
		  } else if (abs(zs[j] - tpicks[i * n1 + jj])
			     < 0.001 * dvs) {
		     vi[j0 + j] = vv[jj];
		  } else {
		     vi[j0 + j] = vv[jj + 1];
		  }
	       }
	    } else {
	       lin1d_(tpicks + i * n1, vpicks + i * n1, &np, zs,
		      vi + j0, &nvs, indx);
	    }
	 }
      } else {
	 /* output in depth */
	 if (intype == 1) {
	    /* interpolate vs3d to nvs times */
	    lin1d_(tpicks + i * n1, vpicks + i * n1, &np, time,
		   vi + j0, &nvs, indx);
	    tmax = tpicks[i * n1 + np - 1];
	    /* compute interval velocities using dix formula */

	    if (vitype == 0) {
	       work[0] = vi[j0] * vi[j0];
	       for (j = 1; j < nvs; j++) {
		  if (time[j] <= tmax) {
		     work[j] = (vi[j + j0] * vi[j + j0] *
				time[j] -
				vi[j - 1 + j0] * vi[j - 1 + j0] *
				time[j - 1]) / dt;
		  } else {
		     work[j] = work[j - 1];
		  }
	       }
	       for (j = 0; j < nvs; j++) {
		  if (work[j] < 0.)
		     err
			   ("check interval velocity at t=%g x=%g y=%g location \n",
			    j * dvs, xs[i], ys[i]);
		  work[j] = sqrt(work[j]);
	       }
	    } else {
	       for (j = 0; j < nvs; j++)
		  work[j] = vi[j + j0];
	    }
	 } else {
	    if (vitype == 0) {
	       vv[0] = vpicks[i * n1];
	       for (j = 1; j < np; j++) {
		  vv[j] = (vpicks[i * n1 + j] * vpicks[i * n1 + j] *
			   tpicks[i * n1 + j] -
			   vpicks[i * n1 + j - 1] * vpicks[i * n1 + j -
							   1] *
			   tpicks[i * n1 + j - 1]) / (tpicks[i * n1 +
							     j] -
						      tpicks[i * n1 +
							     j - 1]);

		  if (vv[j] < 0.)
		     err
			   ("check interval velocity at t=%g x=%g y=%g location \n",
			    tpicks[i * n1 + j], xs[i], ys[i]);

		  vv[j] = sqrt(vv[j]);
	       }
	    } else {
	       for (j = 1; j < np; j++) {
		  vv[j] = vpicks[i * n1 + j];
	       }
	    }

	    bisear_(&np, &nvs, tpicks + i * n1, time, indx);
	    for (j = 0; j < nvs; j++) {
	       jj = indx[j] - 1;
	       if (jj < 0 || time[j] <= tpicks[i * n1]) {
		  work[j] = vv[0];
	       } else if (jj > np - 2) {
		  work[j] = vv[np - 1];
	       } else if (abs(time[j] - tpicks[i * n1 + jj])
			  < 0.001 * dt) {
		  work[j] = vv[jj];
	       } else {
		  work[j] = vv[jj + 1];
	       }
	    }

	 }
	 depth[0] = time[0] * work[0] * 0.5 * 0.001;
	 for (j = 1; j < nvs; j++) {
	    depth[j] = depth[j - 1] +
		  (time[j] - time[j - 1]) * work[j] * 0.5 * 0.001;
	 }
	 lin1d_(depth, work, &nvs, zs, vi + j0, &nvs, indx);
      }
      if (vintpr == 1) {
	 for (j = 0; j < nvs; j++) {
	    fprintf(stderr, "  time/depth=%g   vint=%g \n",
		    zs[j], vi[j + j0]);
	 }
      }
      /* convert to rms velocity if needed */
      if (votype == 0) {
	 tmp = 0.;
	 for (j = 1; j < nvs; j++) {
	    tmp = tmp + vi[j0 + j] * vi[j0 + j] * dvs;
	    vi[j0 + j] = sqrt(tmp / zs[j]);
	 }
	 /* convert to average velocity if needed */
      } else if (votype == 2) {
	 tmp = 0.;
	 for (j = 1; j < nvs; j++) {
	    tmp = tmp + vi[j0 + j] * dvs;
	    vi[j0 + j] = tmp / zs[j];
	 }
      }
      /* check velocity range and apply scale */
      for (j = 0; j < nvs; j++) {
	 tmp = vi[j0 + j] * vscale;
	 if (tmp < vmin)
	    tmp = vmin;
	 if (tmp > vmax)
	    tmp = vmax;
	 vi[j0 + j] = tmp;
      }
      /* smooth if needed */
      if (sms > 1)
	 smth2d_(vi + j0, work, fsms, fsmx, &nvs, &one, &sms, &one);

   }

   free(tpicks);
   free(vpicks);
   free(nps);
   free(zs);
   free(depth);
   free(time);
   free(fsms);


   /* output time/depth slices */
   free(work);
   free(indx);
   work = (float *) malloc(nxy * sizeof(float));
   vtx = (float *) malloc(nvs * nx * sizeof(float));
   vx = (float *) malloc(nx * sizeof(float));
   indx = (int *) malloc(nxy * sizeof(int));

   if (smx == 1 && smy == 1) {
      datafp = stdout;
   } else {
      datafp = etempfile(NULL);
   }

   np = nf;

   fprintf(stderr, " Start x-y plane interpolation \n");
   /* 2d interpolation */
   for (iy = 0; iy < ny; iy++) {
      y = fy + iy * dy;
      for (ix = 0; ix < nx; ix++) {
	 x = fx + ix * dx;
	 if (noxyintp == 0) {
	    if (nf == 0) {
	       plint_(xs, ys, vi, &nvs, &nxy, &x, &y,
		      vtx + ix * nvs, work, &dismax, indx);
	    } else {
	       intp2d_(xs, ys, vi, &nvs, &nxy, &x, &y,
		       vtx + ix * nvs, &np, indx, work);
	    }
	 } else {
	    bcopy(vi + (iy * nx + ix) * nvs, vtx + ix * nvs,
		  nvs * sizeof(float));
	 }
      }
      if (tslice == 1) {
	 for (is = 0; is < nvs; is++) {
	    for (ix = 0; ix < nx; ix++)
	       vx[ix] = vtx[is + ix * nvs];
	    ip = (is * nx * ny + iy * nx) * sizeof(float);

	    efseek(datafp, ip, 0);
	    efwrite(vx, sizeof(float), nx, datafp);
	 }
      } else {
	 bcopy(vtx, vgrid + iy * nx * nvs, nx * nvs * sizeof(float));
      }
      if (smx == 1 && smy == 1) {
	 fminmax(vtx, nvs * nx, &gmin, &gmax);
	 if (iy == 0) {
	    vmin = gmin;
	    vmax = gmax;
	 } else {
	    if (vmin > gmin)
	       vmin = gmin;
	    if (vmax < gmax)
	       vmax = gmax;
	 }
      }
   }

   fprintf(stderr, " x-y plane interpolation done \n");
   free(work);
   free(indx);
   free(vtx);
   free(xs);
   free(ys);
   free(vi);
   free(vx);

   if (smx > 1 || smy > 1) {
      if (tslice == 1)
	 rewind(datafp);
      outfp = stdout;
      nxny = nx * ny;
      vxy = (float *) emalloc(nxny * sizeof(float));
      work = (float *) emalloc(nxny * sizeof(float));

      for (is = 0; is < nvs; is++) {
	 if (tslice == 1) {
	    efread(vxy, sizeof(float), nxny, datafp);
	 } else {
	    for (iy = 0; iy < ny; iy++)
	       for (ix = 0; ix < nx; ix++)
		  vxy[ix + iy * nx] = vgrid[(iy * nx + ix) * nvs + is];
	 }
	 smth2d_(vxy, work, fsmx, fsmy, &nx, &ny, &smx, &smy);
	 if (tslice == 1) {
	    efwrite(vxy, sizeof(float), nxny, outfp);
	 } else {
	    for (iy = 0; iy < ny; iy++)
	       for (ix = 0; ix < nx; ix++)
		  vgrid[(iy * nx + ix) * nvs + is] = vxy[ix + iy * nx];
	 }
	 fminmax(vxy, ny * nx, &gmin, &gmax);
	 if (is == 0) {
	    vmin = gmin;
	    vmax = gmax;
	 } else {
	    if (vmin > gmin)
	       vmin = gmin;
	    if (vmax < gmin)
	       vmax = gmax;
	 }
      }
      free(vxy);
      free(work);
   }
   free(fsmx);
   free(fsmy);

   if (tslice == 0) {
      if (smx == 1 && smy == 1)
	 outfp = datafp;
      efwrite(vgrid, sizeof(float), nx * ny * nvs, outfp);
   } else {
      if (smx == 1 && smy == 1)
	 outfp = datafp;
   }

   fflush(outfp);


   bzero((char *) &gh, GHDRBYTES);
   gh.scale = 1.;
   if (tslice == 1) {
      putgval(&gh, "n1", (float) nx);
      putgval(&gh, "n2", (float) ny);
      putgval(&gh, "n3", (float) nvs);
      putgval(&gh, "o1", fx);
      putgval(&gh, "o2", fy);
      putgval(&gh, "o3", 0.);
      putgval(&gh, "d1", dx);
      putgval(&gh, "d2", dy);
      putgval(&gh, "d3", dvs);
   } else {
      putgval(&gh, "n1", (float) nvs);
      putgval(&gh, "n2", (float) nx);
      putgval(&gh, "n3", (float) ny);
      putgval(&gh, "o1", 0.);
      putgval(&gh, "o2", fx);
      putgval(&gh, "o3", fy);
      putgval(&gh, "d1", dvs);
      putgval(&gh, "d2", dx);
      putgval(&gh, "d3", dy);
   }


   putgval(&gh, "dtype", 4.);
   putgval(&gh, "gmin", vmin);
   putgval(&gh, "gmax", vmax);
   putgval(&gh, "gtype", gtype);
   putgval(&gh, "orient", orient);
   putgval(&gh, "ocdp2", ocdp2);
   putgval(&gh, "dcdp2", dcdp2);
   putgval(&gh, "oline3", oline3);
   putgval(&gh, "dline3", dline3);


   ierr = fputghdr(outfp, &gh);
   if (ierr != 0)
      err("error output grid header ");
   return 0;

}