gridsediment.c

seismic software,very useful

		if(usghin.n3!=usghbot.n2) err("check layerbot header n2");
		if(usghin.o2!=usghbot.o1) err("check layerbot header o1");
		if(usghin.o3!=usghbot.o2) err("check layerbot header o2");
		if(usghin.d2!=usghbot.d1) err("check layerbot header d1");
		if(usghin.d3!=usghbot.d2) err("check layerbot header d2");
		efseek(botfp,0,0);
		efread(zbot,sizeof(float),n2*n3,botfp);

		if(ivtop==1) {
		if(usghin.n2!=usghvtop.n1) err("check gtopgrid header n1");
		if(usghin.n3!=usghvtop.n2) err("check gtopgrid header n2");
		if(usghin.o2!=usghvtop.o1) err("check gtopgrid header o1");
		if(usghin.o3!=usghvtop.o2) err("check gtopgrid header o2");
		if(usghin.d2!=usghvtop.d1) err("check gtopgrid header d1");
		if(usghin.d3!=usghvtop.d2) err("check gtopgrid header d2");
		efseek(vtopfp,0,0);
		efread(vtop,sizeof(float),n2*n3,vtopfp);
		} 

		if(ivbot==1) {
		if(usghin.n2!=usghvbot.n1) err("check gbotgrid header n1");
		if(usghin.n3!=usghvbot.n2) err("check gbotgrid header n2");
		if(usghin.o2!=usghvbot.o1) err("check gbotgrid header o1");
		if(usghin.o3!=usghvbot.o2) err("check gbotgrid header o2");
		if(usghin.d2!=usghvbot.d1) err("check gbotgrid header d1");
		if(usghin.d3!=usghvbot.d2) err("check gbotgrid header d2");
		efseek(vbotfp,0,0);
		efread(vbot,sizeof(float),n2*n3,vbotfp);
		}
	} else {
		for(iz=0;iz<nlayer;iz++) {
			getnparstring(iz+1,"layers",&layers);
			layerfp=efopen(layers,"r");
			ierr = fgetusghdr(layerfp,&usghlayer);
			if(ierr!=0) err(" error open layers=%s \n",layers);
			if(usghin.n2!=usghlayer.n1) err("check %s header n1",layers);
			if(usghin.n3!=usghlayer.n2) err("check %s header n2",layers);
			if(usghin.o2!=usghlayer.o1) err("check %s header o1",layers);
			if(usghin.o3!=usghlayer.o2) err("check %s header o2",layers);
			if(usghin.d2!=usghlayer.d1) err("check %s header d1",layers);
			if(usghin.d3!=usghlayer.d2) err("check %s header d2",layers);
			efseek(layerfp,0,0);
			efread(zs+iz*n2*n3,sizeof(float),n2*n3,layerfp);
			efclose(layerfp);
		}
		for(iz=0;iz<nglayer;iz++) {
			getnparstring(iz+1,"glayers",&glayers);
			glayerfp=efopen(glayers,"r");
			ierr = fgetusghdr(glayerfp,&usghglayer);
			if(ierr!=0) err(" error open layers=%s \n",layers);
			if(usghin.n2!=usghglayer.n1) err("check %s header n1",glayers);
			if(usghin.n3!=usghglayer.n2) err("check %s header n2",glayers);
			if(usghin.o2!=usghglayer.o1) err("check %s header o1",glayers);
			if(usghin.o3!=usghglayer.o2) err("check %s header o2",glayers);
			if(usghin.d2!=usghglayer.d1) err("check %s header d1",glayers);
			if(usghin.d3!=usghglayer.d2) err("check %s header d2",glayers);
			efseek(glayerfp,0,0);
			efread(vs+iz*n2*n3,sizeof(float),n2*n3,glayerfp);
			efclose(glayerfp);
		}
		nz = nlayer - 1;
	}

/* compute mini layer depth and grid values */
	if(nlayer==0) {
		for(i2=0;i2<n2*n3;i2++) {
			zs[i2] = ztop[i2];
			vs[i2] = vtop[i2];
		}
		tmp = 0.;
		for(iz=0;iz<nz;iz++) tmp = tmp + dzrl[iz];
		vscale = 0.;
		zscale = 0.;
		for(iz=1;iz<nz;iz++) {
			zscale += dzrl[iz-1]/tmp;
			if(ginterp==0) {
				vscale += dzrl[iz-1]/tmp;
			} else {
				vscale = (float)iz / nz;
			}

			for(i2=0;i2<n2*n3;i2++) {
				zs[i2+iz*n2*n3] = ztop[i2] +
					zscale*(zbot[i2]-ztop[i2]);
				vs[i2+iz*n2*n3] = vtop[i2] + 
					vscale*(vbot[i2]-vtop[i2]);
			}
		}
		for(i2=0;i2<n2*n3;i2++) {
			zs[i2+nz*n2*n3] = zbot[i2];
			vs[i2+nz*n2*n3] = vbot[i2];
		}
	}

/* compute smoothing window for mini layers */
	tmp = 0.;
	for(iz=0;iz<nz;iz++) tmp = tmp + dzrl[iz];
	scale = 0.;
	sm2s[0] = sm2top;
	sm3s[0] = sm3top;
	for(iz=1;iz<nz+1;iz++) {
		scale += dzrl[iz-1]/tmp;
		sm2s[iz] = sm2top + scale*(sm2bot-sm2top);
		sm3s[iz] = sm3top + scale*(sm3bot-sm3top);
	} 

/* compute grid values at the mini layers if not specified */
	fseek2g(infp,0,0);
	if( ( (ivtop==0 || ivbot==0) && nlayer==0 ) || (nglayer==0 && nlayer>0) ) {
		for(i3=0;i3<n3;i3++) {
			for(i2=0;i2<n2;i2++) {
				efread(grid,sizeof(float),n1,infp);
				for(iz=0;iz<nz+1;iz++) {
					tmp = (zs[i2+i3*n2+iz*n2*n3] - o1)/d1 + 0.5; 
					i1 = tmp;
					if(i1<0) {
						vs[i2+i3*n2+iz*n2*n3] = grid[0];
					} else if(i1>=n1-1) {
						vs[i2+i3*n2+iz*n2*n3] = grid[n1-1];
					} else {
						vs[i2+i3*n2+iz*n2*n3] = grid[i1]+
								(tmp-i1)*(grid[i1+1]-grid[i1]);
					}
				}
			}
		}
	} else if(ivtop==-1 && ivbot==-1) {
		for(i3=0;i3<n3;i3++) {
			for(i2=0;i2<n2;i2++) {
				for(iz=0;iz<nz+1;iz++) {
					tmp = (zs[i2+i3*n2+iz*n2*n3] - ztop[i3*n2+i2]);
					vs[i2+i3*n2+iz*n2*n3] = g0 + r0 * tmp;
				}	
			}
		}
	}
/* smooth mini layer grids */
	for(iz=0;iz<nz+1;iz++) {
		tmp = sm2s[iz]/d2;
		ismx = tmp + 1.5;
		tmp = sm3s[iz]/d3;
		ismy = tmp + 1.5;
		fsmx = (float*) emalloc(ismx*sizeof(float));
		fsmy = (float*) emalloc(ismy*sizeof(float));
		/* 2d smoothing */
		smth2d_(vs+iz*n2*n3,work,fsmx,fsmy,&n2,&n3,&ismx,&ismy);
		/*
		dump2xplot(vs+iz*n2*n3,n2,n3,0,"vsmoth");
		*/
		free(fsmx);
		free(fsmy);
	}


/* output grid */
	fseek2g(infp,0,0);
	for(i3=0;i3<n3;i3++) {
		for(i2=0;i2<n2;i2++) {

			efread(grid,sizeof(float),n1,infp);
			itmp = i2+i3*n2;
			if(nlayer==0) {
				top = ztop[itmp];
				bot = zbot[itmp];
			} else {
				top = zs[itmp];
				bot = zs[itmp+(nlayer-1)*n2*n3];
			}
			tmp = (top - o1)/d1 + 0.5; 
			i1top = tmp;
			tmp = (bot - o1)/d1 + 0.5; 
			i1bot = tmp;
			if(i1top<0) i1top = 0;
			if(i1bot>n1-1) i1bot = n1-1;
			for(i1=i1top;i1<=i1bot;i1++) {
				z = o1 + i1*d1;
                                if(z>=zs[itmp] && z<=zs[itmp+nz*n2*n3]) {
					for(iz=0;iz<nz;iz++) {
						if( z == zs[itmp+iz*n2*n3] && z == zs[itmp+(iz+1)*n2*n3] ){
							grid[i1] = vs[itmp+iz*n2*n3];

						}else if(z>=zs[itmp+iz*n2*n3] && z<zs[itmp+(iz+1)*n2*n3]) {
							tmp = (z-zs[itmp+iz*n2*n3]) /
					  			(zs[itmp+(iz+1)*n2*n3]-zs[itmp+iz*n2*n3]);
							grid[i1] = vs[itmp+iz*n2*n3] + 
								tmp*(vs[itmp+(iz+1)*n2*n3]-vs[itmp+iz*n2*n3]);
                                                	if( grid[i1] != grid[i1] ){
                                                   		fprintf( stderr ,"NaN at i3=%d " ,i3 );
                                                   		fprintf( stderr ,"i2=%d i1=%d\n" ,i2, i3);
                                                	}
							break;
						}
					}
				}
			}

			if(igabovetop==1) {
				for(i1=0;i1<i1top;i1++) {
					grid[i1] = gabovetop;
				}
			}
			if(igbelowbot==1) {
				for(i1=i1bot;i1<n1;i1++) {
					grid[i1] = gbelowbot;
				}
			}

			if(i2==0 && i3==0) {
				gmin = grid[0];
				gmax = grid[0];
			}
			for(i1=0;i1<n1;i1++) {
				if(gmin>grid[i1]) gmin = grid[i1];
				if(gmax<grid[i1]) gmax = grid[i1];
			}

                        if( grid[0] != grid[0] ){
                           fprintf( stderr ,"NaN\n" );
                        }
			efwrite(grid,sizeof(float),n1,outfp);

		}
	}

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

	ierr = fputusghdr(outfp,&usghin);
	
	free(ztop);
	free(zbot);
	free(vtop);
	free(vbot);
	free(zs);
	free(vs);
	free(work);
	free(grid);
	exit(0);
}