supstack.c

seismic software,very useful

	/* loop over traces */
	do {

		if(cdp1==0) {
			gethval(&tr,indxtrk,&trkval);
			ix = vtoi(trktype,trkval);
			gethval(&tr,indxlnk,&lnkval);
			iy = vtoi(lnktype,lnkval);
			x = (ix - fx)/dx + 0.5;
			i2 = x;
			y = (iy - fy)/dy + 0.5;
			i3 = y;
			newcdp = iy*10000+ix;
			x = ix;
			y = iy;
		} else {
			y = (tr.cdp - cdp1)/ncdppl;
			i3 = y;
			x = ( (tr.cdp - cdp1) - i3*ncdppl )/dcdpx + 0.5;
			i2 = x;
			newcdp = tr.cdp;
			y = y / dline;
			i3 = y;
			x = (tr.cdp - cdp1) - i3*ncdppl;
			fx = 0.;
			dx = dcdpx;
			y = (tr.cdp - cdp1)/ncdppl;
			fy = 0.;
			dy = dline;
		}

		if ( newcdp!=oldcdp ) {

			/* output partial stacked gather */
			if(inmo==0)
			pstkout(gather,headers,nofo,fold,nt,ofomin,dofo,outfp);
			bzero(fold,nofo*sizeof(int));
			bzero(gather,nofo*nt*sizeof(float));

			/* compute new sloth function ovv(t) */
			if(i2<0) i2=0; 
			if(i2>n2-1) i2=n2-1;
			if(i3<0) i3=0; 
			if(i3>n3-1) i3=n3-1;
			if(iprint==1) {
				if(cdp1>0) {
					fprintf(stderr," vgrid location i2=%d i3=%d for cdp=%d \n", 
						i2+1,i3+1,tr.cdp);
				} else {
		  fprintf(stderr," vgrid location i2=%d i3=%d for trace=%d line=%d\n", 
						i2+1,i3+1,ix,iy);
				}
			}
			readovv(vgfp,ntvgrid,n2,n3,i2,i3,dtvgrid,ftvgrid,ntc,dtc,ft,ovvt,
					idisk,vs,fx,fy,dx,dy,x,y);
			newsloth = 1;
			mcdp = mcdp + 1;
		} else {
			newsloth = 0;
		}

		if(inmo==0) {
			if(key==1) {
				gethdval(&tr, binkey, &valkey);
				keyvalue = vtof(typekey,valkey);
				ofo = (keyvalue - minkey)/dkey + 0.5;
				iof = ofo;
			} else {
				ofo = tr.offset;
				ofo = (fabs(ofo) - ofomin)/dofo + 0.5;
				iof = ofo;
			}
		} else { 
			iof = 0;
		}

/*
fprintf(stderr,"keyvalue=%g iof=%d \n",keyvalue,iof);
*/
		
		if(retain==1) {
			if(iof<0) iof = 0;
			if(iof>nofo-1) iof = nofo - 1;
		}
		if(iof<0 || iof > nofo-1 ) continue;

		/* if sloth function or offset has changed */
		if (newsloth || tr.offset!=oldoffset) {
		
			if(inmo==0) {
				ofo = ofomin + iof * dofo;
				ofo2 = ofo*ofo; 
				temp = (tr.offset*tr.offset-ofo2)*odtt2;
				temp2 = ofo2 * odtt2;
			} else {
				temp = (tr.offset*tr.offset)*odtt2;
			}

			/* compute time t(tn) (normalized) */
			/* nmo */
			if(inmo==1) {
				for (it=0,tn=ftc; it<ntc; ++it,tn+=fratio)
					ttnc[it] = sqrt(tn*tn+temp*ovvt[it]);
			/* inverse nmo */
			} else if(inmo==-1) {
				/* compute inverse nmo time */
				for (it=0,tn=ftc; it<ntc; ++it,tn+=fratio)
					ttnc[it] = sqrt(tn*tn+temp*ovvt[it]);

				/* interpolate sloth */
				finds(ttnc,ovvt,ntc,tntc,aovv,ntc);

				for (it=0,tn=ftc; it<ntc; ++it,tn+=fratio) {
					tmp=tn*tn - temp*aovv[it];
					if(tmp>0.) {
						ttnc[it] = sqrt(tmp);
					} else {
						ttnc[it] = nt*3;
					}
				}
			/* partial nmo */
			} else {
				/* compute inverse nmo time */
				for (it=0,tn=ftc; it<ntc; ++it,tn+=fratio)
					ttnc[it] = sqrt(tn*tn+temp2*ovvt[it]);
				/* interpolate sloth */
				finds(ttnc,ovvt,ntc,tntc,aovv,ntc);

				for (it=0,tn=ftc; it<ntc; ++it,tn+=fratio) {
					tmp=temp*aovv[it]+tn*tn;
					if(tmp>0.) {
						ttnc[it] = sqrt(tmp);
					} else {
						ttnc[it] = nt*3;
					}
				}
			}
			
			/* linear interpolate ttnc to ttn */
			if(iratio==1) {
				bcopy((char*)ttnc,(char*)ttn,nt*sizeof(float));
			} else {
				intsln(ntc,fratio,ftc,ttnc,
					tmin,tmax,nt,tnt,ttn);
			}
			
			/* compute inverse of stretch factor a(tn) */
			if(ttn[0]<nt && ttn[1]<nt) {
				atn[0] = ttn[1]-ttn[0];
			} else {
				atn[0] = 0.;
			}
			for (it=1; it<nt-1; ++it) {
				if(ttn[it]<nt && ttn[it]<nt) {
					atn[it] = ttn[it]-ttn[it-1];
				} else {
					atn[it] = 0.; 
				}
			}
			if(ttn[nt-1]<nt && ttn[nt-2]<nt) {
				atn[nt-1] = ttn[nt-1]-ttn[nt-2];
			} else {
				atn[nt-1] = 0.;
			}

			for(it=1;it<nt;it++) {
				if(atn[it]<0.) {
					atn[it-1] = 0.;
					ttn[it-1] = nt + 1;
				}
			}

			if(atn[1]==0) {
				atn[0] = 0.;
				ttn[0] = 0.;	
			}
			for(it=0;it<nt;it++) {
				if(atn[it]>smute*2) {
					atn[it] = 0.;
					ttn[it] = nt+1;
				}
			}
			/* determine index of first sample to survive mute */
			for (it=0,itmute=0; it<nt && atn[it]<osmute; ++it) 
				itmute++;

		}

		
		/* intype=0 do partial nmo via linear interpolation */
		/* intype=1 do partial nmo via 8-point sinc interpolation */

		if(intype==0) {
			intsln(nt,1.0,ft/dt,tr.data,0.0,0.0,
				nt-itmute,&ttn[itmute],&qtn[itmute]);
		} else {
			ints8r(nt,1.0,ft/dt,tr.data,0.0,0.0,
				nt-itmute,&ttn[itmute],&qtn[itmute]);
		}
			
		/* apply mute */
		for (it=0; it<itmute; ++it)
			qtn[it] = 0.0;

		/* update mute time */
		mutime = (tr.mute*0.001-ft)/dt;
		imut = mutime;

		if(imut>0) {
			for(it=itmute;it<nt;it++) {
				if(imut<ttn[it]) {
					mutime = it;
					break;
				}
			}
		}
		imut = mutime;


		if(itmute > imut) imut = itmute;
		mutime = (imut*dt+ft)*1000.;
		tr.mute = (int) mutime;

		/* apply linear ramp */
		for (it=itmute; it<itmute+lmute && it<nt; ++it)
			qtn[it] *= (float)(it-itmute+1)/(float)lmute;

			
		/* if specified, scale by the NMO stretch factor */
		if (sscale)
			for (it=itmute; it<nt; ++it)
				qtn[it] *= atn[it];
			
		/* stack the partial NMO corrected trace to output gather */
		if(inmo==0) {
			for(it=0;it<nt;it++) 
				gather[it+iof*nt] += qtn[it]; 
			fold[iof] += 1;
			bcopy((char*)&tr,headers+iof*HDRBYTES,HDRBYTES);

		} else {
			for(it=0;it<nt;it++) tr.data[it] = qtn[it]; 
			fputtr(outfp,&tr);
		}

		/* remember offset and cdp */
		oldoffset = tr.offset;
		oldcdp = newcdp;

	} while (fgettr(infp,&tr));

	/* output last gather */

	if(inmo==0) pstkout(gather,headers,nofo,fold,nt,ofomin,dofo,outfp);
	
	if(cdp1!=0) {
		fprintf(stderr,
		" Total %d cdp gathers starting at cdp=%d processed \n",mcdp,ocdp); 
	} else {
		fprintf(stderr,
		" Total %d gathers starting at cdplbl=%d processed \n",mcdp,ocdp); 
	}

	free(gather);
	free(headers);
	free(fold);

	return EXIT_SUCCESS;
}

/* interpolate sloth */ 
void finds(float *t, float *s, int nt, float *tnew, float *snew, int ntnew){

	float *tsort, *ssort;	
	int *indx, it, itnew;
	float tmp, dt;

	tsort = emalloc(nt*sizeof(float));
	ssort = emalloc(nt*sizeof(float));
	indx = emalloc(nt*sizeof(int));
	for(it=0;it<nt;it++) indx[it] = it;
	qkisort(nt,t,indx);
	for(it=0;it<nt;it++) {
		tsort[it] = t[indx[it]];
		ssort[it] = s[indx[it]];
	}	

	for(itnew=0;itnew<ntnew;itnew++) {
		tmp = tnew[itnew];
		if(tmp<=tsort[0]) {
			snew[itnew] = ssort[0];
		} else if(tmp>=tsort[nt-1]) {
			snew[itnew] = ssort[nt-1];
		} else {
			for(it=0;it<nt-1;it++) {
				if(tmp>=tsort[it] && tmp<=tsort[it+1]) {
					dt = tsort[it+1]-tsort[it]; 
					if(dt!=0.) {
						snew[itnew] = ssort[it] + 
							(tmp-tsort[it])*
							(ssort[it+1]-ssort[it])
							/dt;
					} else {
						snew[itnew] = ssort[it]; 
					}
					break;
				}
			}
		}
	}

	free(tsort);
	free(ssort);
	free(indx);
}



/* output partial stacked gather  */
void pstkout(float *gather, char *headers, int nofo, int *fold, int nt,
		float ofomin, float dofo, FILE *outfp) {  

	segytrace tro;
	int it, iof, mx, my;
	float scale, dd, ofo;

	for(iof=0;iof<nofo;iof++) {
		if(fold[iof]>=1) {
			scale=1./fold[iof];
			for(it=0;it<nt;it++) 
				tro.data[it] = gather[it+iof*nt]*scale;
			bcopy(headers+iof*HDRBYTES,(char*)&tro, HDRBYTES);
			mx = (tro.gx+tro.sx)/2;
			my = (tro.gy+tro.sy)/2;
			dd = tro.offset;
			dd = fabs(dd);
			ofo = ofomin + iof*dofo;
			if(dd>0.) {
				tro.sx = mx+ofo*(tro.sx-mx)/dd; 
				tro.gx = mx+ofo*(tro.gx-mx)/dd; 
				tro.sy = my+ofo*(tro.sy-my)/dd; 
				tro.gy = my+ofo*(tro.gy-my)/dd; 
			} else {
				tro.sx = mx-ofo/2; 
				tro.gx = mx-ofo/2; 
				tro.sy = my; 
				tro.gy = my; 
			}
			tro.offset = ofo;
			fputtr(outfp,&tro);
		}
	}
}


/* read and interpolate sloth function from an Vnmo grid file */

void readovv(FILE *vgfp, int ntvgrid, int nx, int ny, int ix, int iy,
	float dtvgrid, float ftvgrid, int nt, float dt, float ft, float *ovvt,
	int idisk, float *vs, float x0, float y0, float dx, float dy, float x, float y) {
	
	float *vread, ratio, t, ftr;
	int icdp, it, itread;
	float tmp;
	long long lpos;
	off_t lofset;

	icdp = ix + iy*nx;

	lpos = icdp;
	lpos = lpos*ntvgrid*sizeof(float);

	vread = (float*) emalloc(ntvgrid*sizeof(float));

	if(idisk==1) {
		bcopy(&lpos,&lofset,8);
		fseek2g(vgfp,lofset,0);	
		fread(vread,sizeof(float),ntvgrid,vgfp);
	} else {
		bilint_(&ntvgrid,&nx,&ny,&x0,&y0,&dx,&dy,&x,&y,vs,vread);
	}


	ratio = dt/dtvgrid;
	ftr = (ft - ftvgrid)/dtvgrid;

	for(it=0;it<nt;it++) {
		t = it*ratio + ftr;
		itread = t;
		if(itread < 0) {
			tmp = vread[0];
		} else if(itread >= ntvgrid-1) {
			tmp = vread[ntvgrid-1];
		} else {
			tmp = vread[itread] + (t-itread)*
				(vread[itread+1]-vread[itread]);
		}
		ovvt[it] = 1./(tmp*tmp);
	}

	free(vread);
}