gridsalt.c

seismic software,very useful

	if(!getparint("nc",&nc)) nc = 0;
	if(vscalbot!=1.0 || vscalzmx!=1.0) iscal = 1;
	if(ibot2==0) iscal = 0; 

	if(iscal==1) {
		one =1;
		scale1 = alloc3float(n2,n3,one);
		scale2 = alloc3float(n2,n3,one);
		for(i3=0;i3<n3;i3++) {
		for(i2=0;i2<n2;i2++) {
			if(zbot2[i3*n2+i2]<zmax && zbot2[i3*n2+i2]!=0.) {
				scale1[0][i3][i2] = vscalbot;
				scale2[0][i3][i2] = vscalzmx;
			} else {
				scale1[0][i3][i2] = 1.0;
				scale2[0][i3][i2] = 1.0;
			}
		}
		}
		r1 = 0.;
        	r2 = (d2>0)?dscal/d2:0;
        	r3 = (d3>0)?dscal/d3:0;
        	r1 = 0.5*r1*r1 ;
        	r2 = 0.5*r2*r2 ;
        	r3 = 0.5*r3*r3 ;
		lam = 1.0;
		slowness=0;
		niter = 2;
		depth = 3;
		vsm3d(scale1,one,n3,n2,niter,depth,r1,r3,r2,lam,slowness);
		vsm3d(scale2,one,n3,n2,niter,depth,r1,r3,r2,lam,slowness);
	}


	efseek(infp,0,0);
        if(check==1) {
		cube= (float*) emalloc(n1*n2*n3*sizeof(float));
		vs = (float*) emalloc(8*sizeof(float));
	}

	for(i3=0;i3<n3;i3++) {
		for(i2=0;i2<n2;i2++) {

			efread(grid,sizeof(float),n1,infp);

			top = ztop[i3*n2+i2];
			bot = zbot[i3*n2+i2];
			top2 = ztop2[i3*n2+i2];
			bot2 = zbot2[i3*n2+i2];

			tmp = (top - o1)/d1 + 0.5; 
			i1top = tmp;
			tmp = (bot - o1)/d1 + 0.5; 
			i1bot = tmp;
			tmp = (top2 - o1)/d1 + 0.5; 
			i1top2 = tmp;
			tmp = (bot2 - o1)/d1 + 0.5; 
			i1bot2 = tmp;

			if(top==depthnull) i1top = n1;
			if(bot==depthnull) i1bot = 0;
			if(top2==depthnull) i1top2 = n1;
			if(bot2==depthnull) i1bot2 = 0;

/*
	fprintf(stderr, "top=%g bot=%g top2=%g bot2=%g line=%d trace=%d\n",
			top,bot,top2,bot2,(int)(i3*d3+o3),(int)(i2*d2+o2)); 
	fprintf(stderr, "i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n",
			i1top,i1bot,i1top2,i1bot2);
*/
			if(ivtop==1) {
				vt = vtop[i3*n2+i2];
			} else {
				if(i1top>0) {
					vt = grid[i1top];
				} else {
					vt = grid[0];
				}
			}
			if(ivbot2==1) {
				vb = vbot2[i3*n2+i2];
			} else {
				vb = vtop[i3*n2+i2];
			}


			/* salt between top and bot2 */
			if( (itop2==0 && ibot==0 && top<=bot2) || 
	    ( (i1top2==n1 || i1bot==0) && top<=bot2 && bot2!=depthnull) || 
			  (bot>=top2 && bot!=depthnull && top2!=depthnull) 
			  || (ibot2==1 && itop2==0 && ibot==0) ) {

/*
if(i2==381 && i3==169) fprintf(stderr," top to bot2 i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n", i1top,i1bot,i1top2,i1bot2);
if(i2==381 && i3==169) fprintf(stderr," top to bot2 top=%g bot=%g top2=%g bot2=%g \n", top,bot,top2,bot2);
*/

		if(i1top2==n1&&i1bot!=0&&i1bot2>i1bot) i1bot2=i1bot;  

	/*
	if(i1top>i1bot2)
	fprintf(stderr, "i1top=%d i1bot2=%d line=%d trace=%d\n",
			i1top,i1bot2,(int)(i3*d3+o3),(int)(i2*d2+o2)); 
	*/
				if(saltvelo!=-9898) {
					for(i1=i1top;i1<=i1bot2;i1++) {
						if(i1>=0 && i1<=n1-1) grid[i1] = saltvelo;
					}
				} else {
					for(i1=i1top;i1<=i1bot2;i1++) {
						if(i1>=0 && i1<=n1-1) {
							tmp = i1*d1+o1;
							tmp = tmp - top;
							if(bot2>top && tmp>0.) {
								grid[i1] = vt + tmp*(vb-vt)/(bot2-top);
							} else {
								grid[i1] = vt;
							}
						}
					}
				}
			/* overhangs */
			} else if(top<=bot && bot<=top2 && top2<=bot2) {
			/*
if(i2==178 && i3==271) fprintf(stderr," overhang i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n", i1top,i1bot,i1top2,i1bot2);
*/
				if(saltvelo!=-9898) {
					for(i1=i1top;i1<=i1bot;i1++) {
						if(i1>=0 && i1<=n1-1) grid[i1] = saltvelo;
					}
					for(i1=i1top2;i1<=i1bot2;i1++) {
						if(i1<n1 && i1>=0) grid[i1] = saltvelo;
					}
					if(i1bot2<i1top2) i1bot2 = n1;
				} else {
					for(i1=i1top;i1<=i1bot;i1++) {
						if(i1>=0 && i1<=n1-1) {
							tmp = i1*d1+o1;
							tmp = tmp - top;
							grid[i1] = vt + tmp * (vb-vt) / (bot2-top);
						}
					}
					for(i1=i1top2;i1<=i1bot2;i1++) {
						if(i1<n1 && i1>=0) {
							tmp = i1*d1+o1;
							tmp = tmp - top;
							grid[i1] = vt + tmp * (vb-vt) / (bot2-top);
						}
					}
					if(i1bot2<i1top2) i1bot2 = n1;
				}
			/* salt from top2 to bot2 */
			} else if(top2<=bot2 && bot<=top && bot2!=depthnull ) {
			/*
if(i2==178 && i3==271) fprintf(stderr," top2 to bot2 i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n", i1top,i1bot,i1top2,i1bot2);
*/
				if(saltvelo!=-9898) {
					for(i1=i1top2;i1<=i1bot2;i1++) {
						if(i1>=0 && i1<=n1-1) 
							grid[i1] = saltvelo;
					}
				} else {
					for(i1=i1top2;i1<=i1bot2;i1++) {
						if(i1>=0 && i1<=n1-1) {
							tmp = i1*d1+o1;
							tmp = tmp - top;
							grid[i1] = vt + tmp * (vb-vt) /(bot2-top);
						}
					}					
				}
			/* salt below top */
			} else if(top2<=bot) {
			/*
if(i2==178 && i3==271) fprintf(stderr," below top i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n", i1top,i1bot,i1top2,i1bot2);
*/
				if(i1bot2>i1bot || bot2==depthnull) i1bot2=i1bot;
				for(i1=i1top;i1<i1bot2;i1++) {
					if(i1<n1 && i1>=0) grid[i1] = saltvelo;
				}
				if(i1bot2<i1top) i1bot2 = n1;
			/* top to bot salt */
			} else if((top2<=top || top2>bot2) &&  top<= bot) {
			/*
if(i2==178 && i3==271) fprintf(stderr," top to bot i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n", i1top,i1bot,i1top2,i1bot2);
*/
				/*
				if(i1bot>i1bot2) i1bot=i1bot2;
				*/
				for(i1=i1top;i1<=i1bot;i1++) {
					if(i1>=0 && i1<=n1-1) 
						grid[i1] = saltvelo;
				}
				i1bot2 = i1bot;	
				if(i1bot2<i1top) i1bot2 = n1;
			/* salt below top2 */
			} else if(bot<=top) {
			/*
if(i2==178 && i3==271) fprintf(stderr," below top2 i1top=%d i1bot=%d i1top2=%d i1bot2=%d \n", i1top,i1bot,i1top2,i1bot2);
*/
				if(i1top2<i1top) i1top2 = i1top;
				for(i1=i1top2;i1<=i1bot2;i1++) {
					if(i1>=0 && i1<n1) grid[i1] = saltvelo;
				}
				if(i1bot2<i1top2) i1bot2 = n1;
			}

			if(iscal==1) {
				r1 = scale1[0][i3][i2];
				r2 = scale2[0][i3][i2];
				tmp = r2 - r1;
				if((n1-1)>i1bot2) {
					tmp = tmp/(n1-1-i1bot2);
					for(i1=i1bot2;i1<n1;i1++)
						grid[i1] *=(r1+tmp*(i1-i1bot2));
				}
			}
			
			if(check==1) { 
				bcopy(grid,cube+(i3*n2+i2)*n1,
					n1*sizeof(float)); 
			} else {
				efwrite(grid,sizeof(float),n1,outfp);
				for(i1=0;i1<n1;i1++) {
					if(gmin>grid[i1]) gmin = grid[i1]; 
					if(gmax<grid[i1]) gmax = grid[i1]; 
				}
			}
		}
	}


	if(check==1) {
	  for(jter=0;jter<iter;jter++) {
	   for(i3=1;i3<n3-1;i3++) {
	      for(i2=1;i2<n2-1;i2++) {
	         for(i1=0;i1<n1;i1++) {
			v0 = cube[i1+i2*n1+i3*n1*n2];
			vs[0] = cube[i1+i2*n1+(i3-1)*n1*n2];
			vs[1] = cube[i1+i2*n1+(i3+1)*n1*n2];
			vs[2] = cube[i1+(i2-1)*n1+i3*n1*n2];
			vs[3] = cube[i1+(i2+1)*n1+i3*n1*n2];
			vs[4] = cube[i1+(i2-1)*n1+(i3-1)*n1*n2];
			vs[5] = cube[i1+(i2+1)*n1+(i3-1)*n1*n2];
			vs[6] = cube[i1+(i2-1)*n1+(i3+1)*n1*n2];
			vs[7] = cube[i1+(i2+1)*n1+(i3+1)*n1*n2];
			nsalt = 0;
			velo = 0.;
			for(ii=0;ii<8;ii++) {
				if(vs[ii]==saltvelo) {
					nsalt += 1;
				} else {
					velo = velo + vs[ii];
				}
			}
			if(nsalt<=nc && v0==saltvelo) {
				cube[i1+i2*n1+i3*n1*n2] = velo/(8.-nsalt);
			} else if(nsalt>=(8-nc) && v0!=saltvelo) {
				cube[i1+i2*n1+i3*n1*n2] = saltvelo;
			}	
		 }
	      }	
	   }	
          }
	   for(i3=0;i3<n3;i3++) {
	   	for(i2=0;i2<n2;i2++) {
			bcopy(cube+i3*n1*n2+i2*n1,grid,n1*sizeof(float));
			efwrite(grid,sizeof(float),n1,outfp);
			for(i1=0;i1<n1;i1++) {
				if(gmin>grid[i1]) gmin = grid[i1]; 
				if(gmax<grid[i1]) gmax = grid[i1]; 
			}
	   	}
	   }
	}

	usghin.gmin = gmin;
	usghin.gmax = gmax;

	ierr = fputusghdr(outfp,&usghin);
	
	free(ztop);
	free(zbot);
	free(ztop2);
	free(grid);
	if(iscal==1) free3float(scale1);	
	if(iscal==1) free3float(scale2);	
	if(check==1) free(cube);
	if(check==1) free(vs);
	exit(0);
}