sucoher.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* SUCOHER: $Revision: 1.0 $ ; $Date: 93/05/10 13:39:15 $		*/

/*----------------------------------------------------------------------
 * Copyright (c) Colorado School of Mines
 * All rights reserved.
 *
 * This code is part of SU.  SU stands for Seismic Unix, a processing line
 * developed at the Colorado School of Mines, partially based on Stanford
 * Exploration Project (SEP) software.  Inquiries should be addressed to:
 *
 *  Jack K. Cohen, Center for Wave Phenomena, Colorado School of Mines,
 *  Golden, CO 80401  ([e-mail])
 *
 * Modified from $CWPROOT/su/main/suvelan.c to replace each sample of a
 * given "central" trace of a receiver gather with the "best" summation
 * over all the traces in the gather through each of a set of slowness
 * lines through the given sample.
 *
 * -----------------------------------------------------------------------------
 * Tony Kocurko
 * Seismological Systems Manager                e-mail: [e-mail]
 * Departmant of Earth Sciences                         [e-mail]
 * Alexander Murray Building - Room ER-4063     office: [phone]
 * Memorial University of Newfoundland          fax   : [phone]
 * St. John's, Newfoundland                     TELEX : 016-4101 MEMORIAL SNF
 * Canada      A1B 3X5
 *----------------------------------------------------------------------
 */

#include "su.h"
#include "segy.h"
#include <sys/time.h>

/*********************** self documentation ******************************/
char *sdoc[] = {
"									",
" SUCOHER - replace each sample of centre trace with best slowness sum",
"									",
" sucoher <stdin >stdout [optional parameters]				",
"									",
" Optional Parameters:							",
" ns=12                   number of slowness intervals",
" ni=10                   number of interpolation intervals",
" fv=1500.0               first velocity (inverse of slowness)",
" lv=3000.0               last velocity  (inverse of slowness)",
" g1=0                    first receiver gather to treat",
" gn=-1                   receiver gather on which to end (< 0 => all)",
" m=-1                    index of \"centre\" trace",
"                         set to center trace of each group by default",
" verbose=0               &   1 track gathers",
"                         &   2 track samples",
"                         &   4 track slowness",
"                         &   8 track traces",
"                         &  16 track interpolation",
"                         &  32 track summation",
"                         &  64 track elapsed time",
"                         & 128 track total input traces time",
" every_n=1               if (verbose > 1) report every_n samples",
"									",
"									",
" Notes:",
" Input traces should be sorted by receiver group co-ordinates -",
" sucoher outputs a single trace every time the receiver group",
" co-ordinates change. Therefore, the output will ",
" be useful only if receiver gathers are input.				",
"									",
" Trace header fields accessed:  ns, dt, delrt, offset, gx, gy",
" Trace header fields modified:  none",
NULL};

/**************** end self doc *******************************************/

/* Credits:
 * -----------------------------------------------------------------------------
 * Tony Kocurko
 * Seismological Systems Manager                e-mail: [e-mail]
 * Departmant of Earth Sciences                         [e-mail]
 * Alexander Murray Building - Room ER-4063     office: [phone]
 * Memorial University of Newfoundland          fax   : [phone]
 * St. John's, Newfoundland                     TELEX : 016-4101 MEMORIAL SNF
 * Canada      A1B 3X5
 *----------------------------------------------------------------------
 * 
 * Modified from Dave Hale's suvelan to replace each sample of a 
 * given "central" trace of a receiver gather with the "best" summation
 * over all the traces in the gather through each of a set of slowness
 * lines through the given sample.
 *
 */

#define MAX_TRACES   1000

#define QUIET           0
#define LOG_GATHER      1
#define LOG_SAMPLE      2
#define LOG_SLOWNESS    4
#define LOG_TRACE       8
#define LOG_INTERP     16
#define LOG_SUM        32
#define LOG_TIME       64
#define LOG_TRACES    128

#define min(x,y)   ( (x) < (y) ?  (x) : (y) )
#define max(x,y)   ( (x) > (y) ?  (x) : (y) )
#define abs(x)     ( (x) <  0  ? -(x) : (x) )
#define sign(x)    ( (x) <  0  ?  -1  :  1  )

main (int argc, char **argv) {

	segy  *trace[MAX_TRACES] ; /* array of pointers to input traces           */
	segy   out_trace         ; /* output trace                                */
	struct timeval tp        ; /* timer structure                             */

	int    i                 ; /* loop counter over slownesses                */
	int    j                 ; /* loop counter over traces                    */
	int    k                 ; /* loop counter over gathers                   */
	int    m                 ; /* nearest index to crossing time of slowness  */
	int    g1                ; /* first receiver gather to treat              */
	int    gn                ; /* last receiver gather to treat               */
	int    nt                ; /* number of traces in current gather          */
	int    ni                ; /* number of interpolation intervals           */
	int    ns                ; /* number of samples in  "centre" trace        */
	int    ns_x              ; /* number of samples in  current trace         */
	int    mid               ; /* index of "centre" trace of receiver gather  */
	int    centre            ; /* sets index of "centre": < 0 => (mid = nt/2) */
	int    sample            ; /* index of "current" sample being computed    */
	int    slow_n            ; /* number of slownesses over which to sum      */
	int    max_traces        ; /* maximum number of traces / gather so far    */
	int    max_sums          ; /* maximum number of sums allocated so far     */
	int    max_out           ; /* maximum number of interpolated samples so.. */
	int    verbose           ; /* verboseness flag                            */
	int    every_n           ; /* if (verbose > 1) report every_n samples     */
	int   *sum_counts        ; /* counts of traces contributing to each sum   */
	int    nxout             ; /* number of interpolation values              */
	long   start_sec         ; /* seconds corresponding to start of program   */
	long   total_traces      ; /* count of total input traces                 */
	float  slowness          ; /* slowness over which to compute sample       */
	float  slow_low          ; /* first slowness to use                       */
	float  slow_high         ; /* last slowness to use                        */
	float  dslow             ; /* slowness increment                          */
	float  delrt             ; /* time of first sample of "centre" trace      */
	float  delrt_x           ; /* time of first sample of current trace       */
	float  dt                ; /* sample rate for "centre" trace              */
	float  dt_x              ; /* sample rate for current trace               */
	float  dt_out            ; /* sample rate of interpolated data            */
	float  delta_x           ; /* distance of jth trace from "centre" trace   */
	float  delta_t           ; /* time between "centre" sample and jth trace's*/
	float  x0                ; /* offset of "centre" trace from its shot point*/
	float  t0                ; /* time of first sample of "centre" trace      */
	float  time              ; /* time corresponding to "centre" trace sample */
	float  time_j            ; /* slowness crossing time on current trace     */
	float  tsample           ; /* floating point index of sample to accumulate*/
	float  d                 ; /* distance of crossing slowness to lower index*/
	float *sums              ; /* slowness sums                               */
	float *xout              ; /* interpolation abscissae                     */
	float *yout              ; /* interpolation ordinates                     */

	int    GetReceiverGather (segy **, int  , int   *, int           );
	void   ints8r            (int    , float, float  , float *, float,
	                          float  , int  , float *, float *       );

	/* get the current time in seconds */

	(void)gettimeofday (&tp, (struct timezone *)NULL);
	start_sec = tp.tv_sec;

	/* hook up getpar */

	initargs(argc,argv);
	requestdoc(0);

	/* get optional parameters */

	if ( !getparint   ("g1"     , &g1       )) g1        =    0  ;
	if ( !getparint   ("gn"     , &gn       )) gn        =   -1  ;
	if ( !getparint   ("ns"     , &slow_n   )) slow_n    =   12  ;
	if ( !getparint   ("ni"     , &ni       )) ni        =   10  ;
	if ( !getparfloat ("fv"     , &slow_high)) slow_high = 1500.0;
	if ( !getparfloat ("lv"     , &slow_low )) slow_low  = 3000.0;
	if ( !getparint   ("m"      , &centre   )) centre    =   -1  ;
	if ( !getparint   ("every_n", &every_n  )) every_n   =  100  ;
	if ( !getparint   ("verbose", &verbose  )) verbose   = QUIET ;

	if ( slow_low             == 0.0 ) err ("first velocity must be non-zero.");
	if ( slow_high            == 0.0 ) err ("final velocity must be non-zero.");
	if ( slow_n               <  1   ) err ("must have at least one velocity.");
	if ( slow_low * slow_high <= 0.0 ) err ("velocities must have same sign." );

	slow_low     =  1.0 / slow_low ;   /* higher velocity to lower  slowness */
	slow_high    =  1.0 / slow_high;   /* lower  velocity to higher slowness */

	dslow        =  (slow_high - slow_low) / (float)slow_n;
	max_traces   = 0;
	max_sums     = 0;
	max_out      = 0;
	total_traces = 0;

/* Read in a receiver gather */


	for (k=0; (nt=GetReceiverGather (trace,centre,&max_traces,MAX_TRACES)); k++) {

		if ( k <  g1            ) continue; /* haven't reached first gather yet  */
		if ( k >= gn && gn >= 0 ) break   ; /* finished final gather             */

		total_traces += nt;

		if ( verbose & LOG_GATHER )
			fprintf (stderr, "\ngather %d has %d traces", k, nt);

/* get index of "centre" trace */

		mid  = (centre < 0) ? nt / 2 : ( (centre >= nt) ? nt - 1 : centre);

/* copy "centre" trace to output trace */

		bcopy ((char *)trace[mid], (char *)&out_trace, sizeof(segy));

/* clear output trace data */

		bzero ((char *)&(out_trace.data), 4 * SU_NFLTS);

		ns   = trace[mid]->ns                        ; /* samples in centre trace */
		dt   = (float)(trace[mid]->dt    )/ 1000000.0; /* microseconds to seconds */