vgrid3d.c

seismic software,very useful

				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.e-6;
	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;

}