sunmo.c

su 的源代码库

			intlin(ntnmo,tnmo,vnmo,vnmo[0],vnmo[nvnmo-1],1,&tn,&v);
			ovv[icdp][it] = 1.0/(v*v);
		}
		for (it=0,tn=ft; it<nt; ++it,tn+=dt) {
			intlin(ntnmo,tnmo,anis1,anis1[0],anis1[nanis1-1],1,&tn,
			       &oa1[icdp][it]);
		}
		for (it=0,tn=ft; it<nt; ++it,tn+=dt) {
			intlin(ntnmo,tnmo,anis2,anis2[0],anis2[nanis2-1],1,&tn,
			       &oa2[icdp][it]);
		}
		free1float(vnmo);
		free1float(tnmo);
		free1float(anis1);
		free1float(anis2);
	}

	/* sort (by insertion) sloth and anis functions by increasing cdp */
	for (jcdp=1; jcdp<ncdp; ++jcdp) {
		acdp = cdp[jcdp];
		aovv = ovv[jcdp];
		aoa1 = oa1[jcdp];
		aoa2 = oa2[jcdp];
		for (icdp=jcdp-1; icdp>=0 && cdp[icdp]>acdp; --icdp) {
			cdp[icdp+1] = cdp[icdp];
			ovv[icdp+1] = ovv[icdp];
			oa1[icdp+1] = oa1[icdp];
			oa2[icdp+1] = oa2[icdp];
		}
		cdp[icdp+1] = acdp;
		ovv[icdp+1] = aovv;
		oa1[icdp+1] = aoa1;
		oa2[icdp+1] = aoa2;
	}

	/* get other optional parameters */
	if (!getparfloat("smute",&smute)) smute = 1.5;
	if (!getparint("ixoffset",&ixoffset)) ixoffset=0; 
	  if (ixoffset==0) sy = 0.0;
	if (smute<=0.0) err("smute must be greater than 0.0");
	if (!getparint("lmute",&lmute)) lmute = 25;
	if (!getparint("sscale",&sscale)) sscale = 1;
	if (!getparint("invert",&invert)) invert = 0;
	if (!getparint("upward",&upward)) upward = 0;

	/* allocate workspace */
	ovvt = ealloc1float(nt);
	oa1t = ealloc1float(nt);
	oa2t = ealloc1float(nt);
	ttn = ealloc1float(nt);
	atn = ealloc1float(nt);
	qtn = ealloc1float(nt);
	tnt = ealloc1float(nt);
	at = ealloc1float(nt);
	qt = ealloc1float(nt);

	/* interpolate sloth and anis function for first trace */
	interpovv(nt,ncdp,cdp,ovv,oa1,oa2,(float)tr.cdp,ovvt,oa1t,oa2t);

	/* set old cdp and old offset for first trace */
	oldcdp = tr.cdp;
	oldoffset = tr.offset-1;

	warn("sy = %f",sy);

	/* loop over traces */
	do {
		/* if necessary, compute new sloth and anis function */
		if (tr.cdp!=oldcdp && ncdp>1) {
			interpovv(nt,ncdp,cdp,ovv,oa1,oa2,(float)tr.cdp,
				  ovvt,oa1t,oa2t);
			newsloth = 1;
		} else {
			newsloth = 0;
		}

		/* if sloth and anis function or offset has changed */
		if (newsloth || tr.offset!=oldoffset) {
			/* compute time t(tn) (normalized) */
			temp = ((float) tr.offset*(float) tr.offset + sy*sy)/(dt*dt);
			for (it=0,tn=ft/dt; it<nt; ++it,tn+=1.0) {
				tsq = temp*ovvt[it] + \
				      oa1t[it]*temp*temp / (1.0+oa2t[it]*temp);
				if (tsq<0.0)
					err("negative moveout; check anis1, "
					    "anis2, or suwind far-offset "
					    "traces");
				if ((1.0+oa2t[it]*temp)<=0.0)
					err("anis2 negative and too small; "
					    "check anis2, or suwind far-offset"
					    " traces");
				ttn[it] = sqrt (tn*tn + tsq);
				}
			/* compute inverse of stretch factor a(tn) */
			atn[0] = ttn[1]-ttn[0];
			for (it=1; it<nt; ++it)
				atn[it] = ttn[it]-ttn[it-1];
			
			/* determine index of first sample to survive mute */
			osmute = 1.0/smute;
			if( !upward ) {
				for (it=0; it<nt-1 && atn[it]<osmute; ++it)
					;
			} else {
				/* scan samples from bottom to top */
				for (it=nt-1; it>0 && atn[it]>=osmute; --it)
					;
			}
			itmute = it;

			/* if inverse NMO will be performed */
			if (invert) {
							
				/* compute tn(t) from t(tn) */
				yxtoxy(nt-itmute,1.0,ft/dt+itmute,&ttn[itmute],
					nt-itmute,1.0,ft/dt+itmute,
					ft/dt-nt,ft/dt+nt,&tnt[itmute]);
			
				/* adjust mute time */
				itmute = 1.0+ttn[itmute]-ft/dt;
				itmute = MIN(nt-2,itmute);
								
				/* compute a(t) */
				if (sscale) {
					for (it=itmute+1; it<nt; ++it)
						at[it] = tnt[it]-tnt[it-1];
					at[itmute] = at[itmute+1];
				}
			}
		}
		
		/* if forward (not inverse) nmo */
		if (!invert) {
	
			/* do nmo via 8-point sinc interpolation */
			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;
			
			/* 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];
			
			/* copy NMO corrected trace to output trace */
			memcpy( (void *) tr.data,
					(const void *) qtn, nt*sizeof(float));
		
		/* else inverse nmo */
		} else {
	
			/* do inverse nmo via 8-point sinc interpolation */
			ints8r(nt,1.0,ft/dt,tr.data,0.0,0.0,
				nt-itmute,&tnt[itmute],&qt[itmute]);
			
			/* apply mute */
			for (it=0; it<itmute; ++it)
				qt[it] = 0.0;
			
			/* if specified, undo NMO stretch factor scaling */
			if (sscale)
				for (it=itmute; it<nt; ++it)
					qt[it] *= at[it];
			
			/* copy inverse NMO corrected trace to output trace */
			memcpy( (void *) tr.data,
					(const void *) qt,nt*sizeof(float));
		}

		/* write output trace */
		puttr(&tr);

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

	} while (gettr(&tr));

	return(CWP_Exit());
}


/* linearly interpolate/extrapolate sloth and anis between cdps */
static void interpovv (int nt, int ncdp, float *cdp, float **ovv, float **oa1, 
	float **oa2, float cdpt, float *ovvt, float *oa1t, float *oa2t)
{
	static int indx=0;
	int it;
	float a1,a2;

	/* if before first cdp, constant extrapolate */
	if (cdpt<=cdp[0]) {
		for (it=0; it<nt; ++it) {
			ovvt[it] = ovv[0][it];
			oa1t[it] = oa1[0][it];
			oa2t[it] = oa2[0][it];
		      };
	
	/* else if beyond last cdp, constant extrapolate */
	} else if (cdpt>=cdp[ncdp-1]) {
		for (it=0; it<nt; ++it) {
			ovvt[it] = ovv[ncdp-1][it];
			oa1t[it] = oa1[ncdp-1][it];
			oa2t[it] = oa2[ncdp-1][it];
		      };
	
	/* else, linearly interpolate */
	} else {
		xindex(ncdp,cdp,cdpt,&indx);
		a1 = (cdp[indx+1]-cdpt)/(cdp[indx+1]-cdp[indx]);
		a2 = (cdpt-cdp[indx])/(cdp[indx+1]-cdp[indx]);
		for (it=0; it<nt; ++it) {
			ovvt[it] = a1*ovv[indx][it]+a2*ovv[indx+1][it];
			oa1t[it] = a1*oa1[indx][it]+a2*oa1[indx+1][it];
			oa2t[it] = a1*oa2[indx][it]+a2*oa2[indx+1][it];
		      };
	}
}