sucoher.c

su 的源代码库

		t0   = (float)(trace[mid]->delrt )/    1000.0; /* milliseconds to seconds */
		x0   = (float)(trace[mid]->offset)           ; /* shotpoint to centre dist*/

/* allocate space for slow_n slowness sums for each sample of centre trace    */
/* and similarly, allocate space for slow_n counts of contributing traces for */
/* each sample of the centre trace                                            */

		if ( ns )                         /* If the centre trace has any data,    */
			if ( ns * slow_n > max_sums) {  /* and if insufficient space has been   */
				if ( max_sums ) {             /* allocated in the past, then free up  */
					free (sums)      ;          /* the space previously allocated and   */
					free (sum_counts); }        /* try to allocate a sufficient amount. */
				max_sums = ns * slow_n;
				if( (sums = (float*)calloc (max_sums, sizeof(float))) == (float *)NULL )
					err("couldn't allocate space for slant sums");
				if( (sum_counts = (int *)calloc(max_sums, sizeof(int)))== (int  *)NULL )
					err("couldn't allocate space for slant counts"); }

/* compute the total number of slowness sums to be held in memory. There is  */
/* one sum for each slowness at each sample. Then clear the sums and the     */
/* counts of traces contributing to each sum.                                */

		bzero ((char *)sums      , ns * slow_n * sizeof(float));
		bzero ((char *)sum_counts, ns * slow_n * sizeof(float));

/* For each trace in the gather, . . . */

		for (j = 0; j < nt; j++) {

			if ( verbose & LOG_TRACE )
				fprintf (stderr, "\nWorking on trace %d", j);

/* Get sample rate, time of first sample, number of samples, and distance from*/
/* the "centre" trace of the current trace. For our purposes, we will take the*/
/* distances to be negative for traces that occur "left" of the "centre" trace*/
/* with respect to trace index. So that if j < mid, then we multiply the mag- */
/* nitude of the difference in offsets by -1, and otherwise by +1.            */

			dt_x    = (float)trace[j]->dt    / 1000000.0;
			delrt_x = (float)trace[j]->delrt /    1000.0;
			ns_x    = trace[j]->ns;
			delta_x = abs ((float)trace[j]->offset - x0) * (float)(sign(j - mid));

/* Allocate enough space for the selected level of interpolation of the trace.*/
/* Suppose that the user wants each input sample interval to have ten interpo-*/
/* lation intervals. That is we want each interval to go from this:           */
/*                                                                            */
/* j-1                  j          j-1                  j                     */
/*  |-------------------| to this:  |-+-+-+-+-+-+-+-+-+-|                     */
/*                                  0 1 2 3 4 5 6 7 8 9                       */
/*                                                                            */
/* Then there will be ten output samples for each input sample except for the */
/* last input sample, which will be represented by itself alone on output.    */
/*                                                                            */
/* In general, then, if the user asks for n_itrp interpolation samples per in-*/
/* put sample, and if there are ns_x input samples, then there will be a total*/
/* of n_itrp * (ns_x - 1) + 1 interpolated output samples.                    */
/*                                                                            */
			if ( ns_x )
				if ( ni * (ns_x - 1) + 1 > max_out ) {    /* Need more space?         */
					if ( max_out ) {                        /* Space previously saved?  */
						free (xout);                          /* Free previously saved    */
						free (yout); }                        /* abscissae and ordinates. */
					max_out = ni * (ns_x - 1) + 1;          /* New maximum saved space. */

					if ((xout = (float *)calloc(max_out, sizeof(float))) == (float *)NULL)
						err("Unable to allocate interpolation abscissae array");
					if ((yout = (float *)calloc(max_out, sizeof(float))) == (float *)NULL)
						err("Unable to allocate interpolation array"); }

/* Compute total number of output interpolated samples, then clear them.      */

			nxout  = ni * (ns_x - 1) + 1;
			dt_out = dt / (float)ni;
			bzero ((char *)yout, nxout * sizeof(float));

/* Initialize the interpolation abscissae, in this case time.                */

			for (i = 0; i < nxout; i++)
				xout[i] = delrt + i * dt_out;

/* CWP canned routine to do sinc interpolation.                              */

			ints8r (ns_x          , dt_x             , delrt_x               ,
			        trace[j]->data, trace[j]->data[0], trace[j]->data[ns_x-1],
			        nxout         , xout             , yout                  );

/* Current trace now contributes to a sum for each sample of the "centre"    */
/* trace and for each slowness value. The interpolated current trace is used */
/* to supply amplitudes for those slowness lines that do not intersect the   */
/* current trace at an input sample.                                         */

/* For each sample of the "centre" trace . . .                               */

			for ( sample = 0; sample < ns; sample++) {

				if ( (verbose & LOG_SAMPLE) && sample % every_n == 0 )
					fprintf (stderr, "\nStarting sample %d of gather %d ... ", sample, k);

			  time = t0 + (float)sample * dt; /* Time corresponding to this sample */

/* For each slowness value, . . . */

				for ( i = 0; i < slow_n; i++) {

/* By definition lines with slopes equal to slownesses pass through each      */
/* sample of the "centre" trace. Thus, we only need to accumulate the "centre"*/
/* trace's sample values.                                                     */

					if ( j == mid ) {
						sums      [i + sample * slow_n] += trace[j]->data[sample];
						sum_counts[i + sample * slow_n]++;
						if ( verbose & LOG_INTERP )
							fprintf (stderr, "\ny[%d](%f) = %f +> sum[%d] = %f",
							                    j, time, trace[j]->data[sample], i, sums[i]);
						continue; }

					if ( verbose & LOG_SLOWNESS )
						fprintf (stderr, "\nWorking on slowness %d", i);

/* Compute the time (and thereby the sample) at which the slowness line     */
/* through the "centre" trace's current sample intercepts the current trace */

					slowness = slow_low + i * dslow;
					delta_t  = slowness * delta_x;
					time_j   = time + delta_t;
					tsample  = (time_j - delrt_x) * (float)(nxout - 1)
					         / (dt_x              * (float)(ns_x  - 1));

/* If the slowness line intercepts the current trace in its time range,   */
/* interpolate a value on that trace and accumulate it into the sum for   */
/* the current slowness through the current sample of the "centre" trace. */

					if ( 0.0 <= tsample && tsample < nxout ) {
						sums      [i + sample * slow_n] += yout[(int)tsample];
						sum_counts[i + sample * slow_n]++;
						if ( verbose & LOG_INTERP )
							fprintf (stderr, "\ny[%d](%f) = %f +> sum[%d] = %f",
							                    j, time_j, yout[(int)tsample], i,
							                    sums[i + sample * slow_n]);
						} /* ---------------------------------> End of tsample range test */
					}   /* ---------------------------------> End of slowness loop      */
				}     /* ---------------------------------> End of sample loop        */
			}       /* ---------------------------------> End of trace loop         */

/* For each sample in the output trace, . . . */
			for (i = 0; i < ns; i++) {

/* Point to the first slowness as having given the "best" sum, whether it did */
				j = 0;

/* Then for each slowness at this sample . . . */
				for (k = 0; k < slow_n; k++) {

/* Calculate the average based on the number of traces that contributed to it.*/

					sums[k + i * slow_n] /= (float)sum_counts[k + i * slow_n];

/* If this average is of greater magnitude than previous best, point to it.*/
					if (abs (sums[k + i * slow_n]) > abs(sums[j + i * slow_n]) ) j = k; }

/* Set the output sample to the "best" of the average sums at the input sample*/
				out_trace.data[i] = sums[j + i * slow_n]; }

/* Write out the trace. */
		puttr (&out_trace); }   /* End of GetReceiverGather loop */

	if ( verbose & LOG_TRACES )
		fprintf (stderr, "\ntotal input traces = %d\n", total_traces);

	(void)gettimeofday (&tp, (struct timezone *)NULL);
	if ( verbose & LOG_TIME )
		fprintf (stderr, "\nelapsed time in seconds = %d\n", tp.tv_sec - start_sec);

	exit(0); }

int
GetReceiverGather (segy **trace, int centre, int *n, int max_traces) {

	                            /* trace      = array of pointers to traces     */
	                            /* centre     = index of "centre" trace         */
	                            /* n          = number traces allocated so far  */
	                            /* max_traces = maximum # traces permitted      */
	int         i, j          ; /* counters of traces in next receiver gather   */
	long        gx            ; /* receiver group x- and y-coordinates defining */
	long        gy            ; /*   the next input receiver group              */
	Bool        done          ; /* Goes true on EOF or new receiver group       */
	static segy previous_trace; /* holds 1st trace of subsequent receiver gather*/

	i = 0;

	if ( feof (stdin) ) return(i);

	if ( *n ) {
		bcopy ((char *)&previous_trace, (char *)trace[0], sizeof(segy));
		i = 1; }
	else if ( (trace[0] = (segy *)calloc (1, sizeof(segy))) == (segy *)NULL ) {
		err   ("can't allocate space for first trace.");
		*n = 1; }
	else if ( !gettr (trace[0]) )
		err   ("can't get first trace of gather.");

	gx = trace[0]->gx;
	gy = trace[0]->gy;

	done = false;
	for (j = 1; j < *n; j++ ) {
		if ( !gettr (trace[j]) ) {
			done = true;
			break; }
		if ( trace[j]->gx != gx || trace[j]->gy != gy ) {
			bcopy ( (char *)trace[j], (char *)&previous_trace, sizeof(segy) );
			done = true;
			break; }
		i++; }

	if ( ! done ) {
		for (*n = j; gettr (&previous_trace); *n = *n + 1) {
			if ( previous_trace.gx != gx || previous_trace.gy != gy ) break;
			if ( *n >= max_traces ) err ("can't allocate any more traces.");
			if ( (trace[*n] = (segy *)calloc (1, sizeof(segy))) == (segy *)NULL )
				err ("can't allocate space for a new trace.");
			bcopy ((char *)&previous_trace, (char *)trace[*n], sizeof(segy));
			i++; } }

	return (i); }