gridsalt.c

seismic software,very useful

#include "velo.h"
#include "usgrid.h"
#include "par.h"

char *sdoc =
"GRIDSALT - 3D attribute grid building around a salt body \n"
"\n"
"gridsalt <infile >outfile [optional parameters]\n"
"\n"
"Required Parameters:\n"
"infile=        name of the input 3D attribute grid file	(velocity)	\n"
"outfile=       name of the output 3D attribute grid file (velocity) \n"
"saltvelo=      salt velocity (attribute) value				\n"
"depthnull=     depth value to be used where there is no salt in the 	\n"
"               depth grid files when the grid value is greater or equal \n"
"               to depthnull 					\n"
"salttop=       name of the 2D depth grid defining the top of the salt \n"
 "Optional Parameters:\n"
"salttop2=      name of the 2D depth grid defining the top of the salt \n"
"               below the overhangh 					\n"
"               (if not given, only salttop is used to define salt top)\n"
"saltbot=       name of the 2D depth grid defining the overhang bottom of \n"
"               the salt (if not given,  default to below the bottom of the  \n"
"               model --- no salt overhand bottom )		\n"
"saltbot2=      name of the 2D depth grid defining the bottom of the salt \n"
"               (if not given,  default to below the bottom of the  \n"
"               model --- no salt bottom )		\n"
"vscalbot=1.0   velocity at the salt bottom (saltbot2) will be scaled \n"
"               by vscalbot \n"
"vscalzmx=1.0   velocity at the maximum depth will be scaled by vscalzmx \n"
"               (linearly interpolated scale will be applied for velocities \n"
"               at depths between base of salt and the maximum depth) \n" 
"               when vscalbot and vscalzmx are both 1.0, no scaling applied \n"
"dscal=500.     distance from salt edge where the full scaling applied \n"
"               (linear tapering will be appled to the salt edge)    \n"
"               dscal must have the same unit as d2 and d3 in the velocity \n"
"               grid file						\n"
"check=0        check grided model				\n"
"               0=no						\n"
"               1=yes  						\n"
"               remove salt if <=nc of 8 adjacent lateral points are salt \n"  
"               fill salt if >=(8-nc) of 8 adjacent lateral points are salt \n"
"nc=0           number of adjacent points to have the same property as \n"
"               that of the checked point		\n"
"iter=1         number of iterations to check the grided volume	\n"
"vtopgrid=      grid value file  at the top of the salt top (salttop)\n" 
"vbot2grid=     grid value file at the bottom of the salt bottom (saltbot2)\n" 
"               if both vtopgrid and vbot2grid are not specified, saltvelo \n"
"               is used within the salt \n"
"               if vtopgrid is given and vbot2grid is not, vtopgrid value \n"
"               will be used below salttop \n"
"               if vbot2grid is given and vtopgrid is not, grid values \n"
"               within the salt body will be interpolated from the grid \n"
"               value at the salttop to the the vbot2grid value at the \n"
"               saltbot2 \n"
" Notes:						\n"
" 1. when there are overhungs, use four horizons: \n"
"       salttop, saltbot, salttop2 and saltbot2 as showen below \n"
"       						\n"
" ---- salttop						\n"
" **** saltbot					\n"
" ++++ salttop2 					\n"
" &&&& saltbottom2 					\n"
" SED  sediment region					\n"
" SALT salt region					\n"
"       						\n"
"       						\n"
"           SED						\n"
"         ---------- 					\n"
"  SED   .         *.      				\n"
"       .  SALT   *  .					\n"
" ------         *    .    	  			\n"
"  SALT          * SED .  SED				\n"
"                +      .				\n"
"                  +     .				\n"
"&&&&&&&&&&          + + +------------ 			\n"
"          & 						\n"
"           &						\n"
"            &&&&&&&&&&&&&&&&&&&&&&&&&& 	        \n"
"       						\n"
" 2. when there are no overhungs, use two horizons: \n"
"       salttop and saltbot2 to define the salt body \n"
"\n"
"AUTHOR:  Zhiming Li,         12/05/95			\n"
"\n";

int main(int argc, char **argv)
{
	usghed usghin, usghtop, usghbot, usghtop2, usghbot2;
	usghed usghvtop, usghvbot2;
	FILE *infp,*outfp,*topfp,*botfp,*top2fp,*bot2fp;
	FILE *vtopfp,*vbot2fp;
	char *infile,*outfile,*salttop, *saltbot,*salttop2,*saltbot2;
	char *vtopgrid, *vbot2grid;
	float saltvelo, depthnull, zmax;
	int ibot, itop2, ibot2;
	int ierr;

	int n1,n2,n3;
	int i1,i2,i3;
	float d1,o1;

	float *grid, *ztop, *zbot, *ztop2, *zbot2, gmin, gmax;
	float *vtop, *vbot2;
	float top, bot, top2, bot2;
	float tmp;
	int i1top, i1bot, i1top2, i1bot2;

	float dscal, vscalbot, vscalzmx;
	float ***scale1, ***scale2;
	float r1, r2, r3, d2, d3, lam, o2, o3;
	int iscal=0,depth,slowness,iter,one,jter;

	int check,nc,ii,nsalt,niter;
	float *cube,*vs,v0,velo;
	int ivtop=0, ivbot2=0;
	float vt, vb;


	/* hook up getpar */
	initargs(argc,argv);
	askdoc(0);

	/* get parameters */
	if(getparstring("infile",&infile)) {
		infp = efopen(infile,"r");
	} else {
		infp = stdin;
	}
	ierr = fgetusghdr(infp,&usghin);
    	if(ierr!=0) err(" input grid header error ");
	if(getparstring("outfile",&outfile)) {
		outfp = efopen(outfile,"w");
	} else {
		outfp = stdout;
	}
	if(!getparfloat("saltvelo",&saltvelo)) err(" saltvelo missing ");
	if(!getparfloat("depthnull",&depthnull)) err(" depthnull missing ");


	if (getparstring("salttop",&salttop)) {
		topfp = efopen(salttop,"r");
		ierr = fgetusghdr(topfp,&usghtop);
      		if(ierr!=0) err(" salttop grid header error ");
	} else {
		err(" salttop missing ");
	}

	if (getparstring("saltbot",&saltbot)) {
		botfp = efopen(saltbot,"r");
		ierr = fgetusghdr(botfp,&usghbot);
      	if(ierr!=0) err(" saltbot grid header error ");
		ibot = 1;	
	} else {
		ibot = 0;
	}

	if (getparstring("salttop2",&salttop2)) {
		top2fp = efopen(salttop2,"r");
		ierr = fgetusghdr(top2fp,&usghtop2);
      	if(ierr!=0) err(" salttop2 grid header error ");
		itop2 = 1;	
	} else {
		itop2 = 0;
	}

	if (getparstring("saltbot2",&saltbot2)) {
		bot2fp = efopen(saltbot2,"r");
		ierr = fgetusghdr(bot2fp,&usghbot2);
      	if(ierr!=0) err(" saltbot2 grid header error ");
		ibot2 = 1;	
	} else {
		ibot2 = 0;
	}

	if (getparstring("vtopgrid",&vtopgrid)) {
		vtopfp = efopen(vtopgrid,"r");
		ierr = fgetusghdr(vtopfp,&usghvtop);
      	if(ierr!=0) err(" vtopgrid header error ");
		ivtop = 1;
		saltvelo = -9898;
	}
	if(getparstring("vbot2grid",&vbot2grid) ) {
		vbot2fp = efopen(vbot2grid,"r");
		ierr = fgetusghdr(vbot2fp,&usghvbot2);
      	if(ierr!=0) err(" vbot2grid header error ");
		ivbot2 = 1;
		saltvelo = -9898;
	}

	n1 = usghin.n1;
	n2 = usghin.n2;
	n3 = usghin.n3;
	o1 = usghin.o1;
	d1 = usghin.d1;
	d2 = usghin.d2;
	d3 = usghin.d3;
	o2 = usghin.o2;
	o3 = usghin.o3;
	gmin = usghin.gmin;
	gmax = usghin.gmax;

	zmax = o1 + n1*d1;

	/* memory allocations */
	ztop = (float*) emalloc(n2*n3*sizeof(float));
	zbot = (float*) emalloc(n2*n3*sizeof(float));
	ztop2 = (float*) emalloc(n2*n3*sizeof(float));
	zbot2 = (float*) emalloc(n2*n3*sizeof(float));
	grid = (float*) emalloc(n1*sizeof(float));
	vtop = (float*) emalloc(n2*n3*sizeof(float));
	vbot2 = (float*) emalloc(n2*n3*sizeof(float));

	
	if(usghin.n2!=usghtop.n1) err("check top grid header n1");
	if(usghin.n3!=usghtop.n2) err("check top grid header n2");
	if(usghin.o2!=usghtop.o1) err("check top grid header o1");
	if(usghin.o3!=usghtop.o2) err("check top grid header o2");
	if(usghin.d2!=usghtop.d1) err("check top grid header d1");
	if(usghin.d3!=usghtop.d2) err("check top grid header d2");
	efseek(topfp,0,0);
	efread(ztop,sizeof(float),n2*n3,topfp);


	if(ibot==1) {
		if(usghin.n2!=usghbot.n1) err("check bot grid header n1");
		if(usghin.n3!=usghbot.n2) err("check bot grid header n2");
		if(usghin.o2!=usghbot.o1) err("check bot grid header o1");
		if(usghin.o3!=usghbot.o2) err("check bot grid header o2");
		if(usghin.d2!=usghbot.d1) err("check bot grid header d1");
		if(usghin.d3!=usghbot.d2) err("check bot grid header d2");
		efseek(botfp,0,0);
		efread(zbot,sizeof(float),n2*n3,botfp);
	} else {
		for(i1=0;i1<n2*n3;i1++) zbot[i1] = o1+(n1-1)*d1;
	}

	if(saltvelo==-9898) {
		if(ivtop==1) {
		if(usghin.n2!=usghvtop.n1) err("check vtopgrid header n1");
		if(usghin.n3!=usghvtop.n2) err("check vtopgrid header n2");
		if(usghin.o2!=usghvtop.o1) err("check vtopgrid header o1");
		if(usghin.o3!=usghvtop.o2) err("check vtopgrid header o2");
		if(usghin.d2!=usghvtop.d1) err("check vtopgrid header d1");
		if(usghin.d3!=usghvtop.d2) err("check vtopgrid header d2");
		efseek(vtopfp,0,0);
		efread(vtop,sizeof(float),n2*n3,vtopfp);
		} 
		if(ivbot2==1) {
		if(usghin.n2!=usghvbot2.n1) err("check vbot2grid header n1");
		if(usghin.n3!=usghvbot2.n2) err("check vbot2grid header n2");
		if(usghin.o2!=usghvbot2.o1) err("check vbot2grid header o1");
		if(usghin.o3!=usghvbot2.o2) err("check vbot2grid header o2");
		if(usghin.d2!=usghvbot2.d1) err("check vbot2grid header d1");
		if(usghin.d3!=usghvbot2.d2) err("check vbot2grid header d2");
		efseek(vbot2fp,0,0);
		efread(vbot2,sizeof(float),n2*n3,vbot2fp);
		}
	}

	if(itop2==1) {
		if(usghin.n2!=usghtop2.n1) err("check top2 grid header n1");
		if(usghin.n3!=usghtop2.n2) err("check top2 grid header n2");
		if(usghin.o2!=usghtop2.o1) err("check top2 grid header o1");
		if(usghin.o3!=usghtop2.o2) err("check top2 grid header o2");
		if(usghin.d2!=usghtop2.d1) err("check top2 grid header d1");
		if(usghin.d3!=usghtop2.d2) err("check top2 grid header d2");
		efseek(top2fp,0,0);
		efread(ztop2,sizeof(float),n2*n3,top2fp);
	} else {
		for(i1=0;i1<n2*n3;i1++) ztop2[i1] = depthnull;
	}

	if(ibot2==1) {
		if(usghin.n2!=usghbot2.n1) err("check bot2 grid header n1");
		if(usghin.n3!=usghbot2.n2) err("check bot2 grid header n2");
		if(usghin.o2!=usghbot2.o1) err("check bot2 grid header o1");
		if(usghin.o3!=usghbot2.o2) err("check bot2 grid header o2");
		if(usghin.d2!=usghbot2.d1) err("check bot2 grid header d1");
		if(usghin.d3!=usghbot2.d2) err("check bot2 grid header d2");
		efseek(bot2fp,0,0);
		efread(zbot2,sizeof(float),n2*n3,bot2fp);
	} else {
		for(i1=0;i1<n2*n3;i1++) zbot2[i1] = o1+(n1-1)*d1;
	}

/* compute scales for velocity below the salt */
	if(!getparfloat("vscalbot",&vscalbot)) vscalbot = 1.0;
	if(!getparfloat("vscalzmx",&vscalzmx)) vscalzmx = 1.0;
	if(!getparfloat("dscal",&dscal)) dscal = 500.0;
	if(!getparint("check",&check)) check = 0;
	if(!getparint("iter",&iter)) iter = 1;