gridsalt.new.c

seismic software,very useful

		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];
			}

			ivvt = 0;
			ivvb = 0;
			ivvt2 = 0;
			ivvb2  0;

			/* 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(i1top2==n1&&i1bot!=0&&i1bot2>i1bot) i1bot2=i1bot;  

			   if(saltvelo!=-9898) {
				
				/* boundary value at the top */
				tmp = ztop[i3*n2+i2];
				tmp = (tmp - o1)/d1; 
				itmp = tmp;
				dzz1 = ztop[i3*n2+i2]-(o1+itmp*d1);
				dzz2 = dzz - dzz1;
				if(itmp>=0 && itmp<n1) { 
					v11 = grid[itmp];
					v22 = saltvelo;
					calcv(v11,v22,dzz,dzz1,dzz2,ivatbn,&vvt);
					ivvt = 1;
				}

				/* boundary value at the bottom */
				tmp = zbot2[i3*n2+i2];
				tmp = (tmp - o1)/d1;
				itmp = tmp;
				dzz1 = zbot2[i3*n2+i2]-(o1+itmp*d1);
				dzz2 = dzz - dzz1;
				if(itmp+1>=0 && itmp+1<n1) { 
					v11 = saltvelo;
					v22 = grid[itmp+1];
					calcv(v11,v22,dzz,dzz1,dzz2,ivatbn,&vvb2);
					ivvb2 = 1;
				}

				for(i1=i1top;i1<=i1bot2;i1++) {
					if(i1>=0 && i1<=n1-1) grid[i1] = saltvelo;
				}

				if(i1top>=0&&i1top<=n1-1&&ivvt==1) grid[i1top] = vvt;
				if(i1bot2>=0&&i1bot2<=n1-1&&ivvb2==1) grid[i1bot2] = vvb2;

			   } 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);
}

/* velocity calulation at the nearest grid point of a defined horizon */
/* input:  v1  --- velocity above the horizon 
	   v2  --- velocity below the horizon
	   dz  --- vertical grid interval 
	   dz1 --- distance of horizon to the grid point above
	   dz2 --- distance of horizon to the grid point below
	  ivatbn - method of calculating the velocity 
		   1=minimize the deeper travel time error
		   2=distance-weighted slowness average
   output:
	   v  --- velocity at the grid point nearest to the horizon
*/

void calcv(float v1,float v2, float dz, float dz1, float dz2, 
	int ivatbn, float *v) {

	float a, b, c, tmp1, tmp2, tmp, vv;

  if(ivatbn==1) {
	if(dz1>=dz2) {
		a = (dz1/v1+dz2/v2+dz/v2)/2./dz;
	} else { 
		a = (dz1/v1+dz2/v2+dz/v1)/2./dz;
	}	
	a = 1/.a;
	b = v1 + v2 - 2 * a;
	c = v1*v2 - a*(v1+v2);

	tmp = b*b - 4.*c;
	if (tmp>=0.) {
		tmp = sqrt(tmp);
		tmp1 = (-b+tmp)/2.;
		tmp2 = (-b-tmp)/2.;
		vv = tmp1;
		if(tmp2>tmp1) v = tmp2;
	} else {
		tmp = (dz-dz1)/dz/v1 + (dz-dz2)/dz/v2;
		vv = 1./tmp;
	}

  } else {
		tmp = (dz-dz1)/dz/v1 + (dz-dz2)/dz/v2;
		vv = 1./tmp;
  }

  *v = vv;
}