segdread.c

su 的源代码库

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

/* SEGDREAD: $Revision: 1.9 $ ; $Date: 2004/02/02 18:34:30 $	*/

#include "su.h"
#include "segy.h"
#include "segd.h"

#include "sfio.h"

/* Portable decoding utilities to retrieve char/u_char/short/u_short from SEG-D byte stream */
static char GET_C(Sfio_t *f)
{
 int n;
 char c;

 n = sfgetc(f);
 c = (char) ((n > 127) ? (n-256) : n);

 return (c);
}
 
static unsigned int GET_UC(Sfio_t *f)
{
 int n;
 unsigned char uc;

 n = sfgetc(f);
 uc = (unsigned char) n;

 return (uc);
}

static short GET_S(Sfio_t *f)
{
 int n1, n2;
 int n;
 short s;

 n1 = sfgetc(f);
 n2 = sfgetc(f);
 n = (n1<<8) | n2;
 s = (short) ((n > 32767)? n - 65536 : n);

 return (s);
}

static unsigned short GET_US(Sfio_t *f)
{
 int n1, n2;
 int n;
 unsigned short us;

 n1 = sfgetc(f);
 n2 = sfgetc(f);
 n = (n1<<8) | n2;
 us = (unsigned short) n;

 return (us);
}

#define REC_L (1*20/*demux header*/ + 255*32 /*up to 255 header extensions*/ + 4*SU_NFLTS/*4 bytes per sample max*/) /* expected max record size: 20 + 10*(SU_NFLTS)/4 bytes */

#define BCD_FF (15+10*(15))                   /* bcd interpretation of FF */
#define BCD_FFF (15+10*(15+10*(15)))          /* bcd interpretation of FFF */
#define BCD_FFFF (15+10*(15+10*(15+10*(15)))) /* bcd interpretation of FFFF */

/*********************** self documentation **********************/
char *sdoc[] = {
"                                                                      ",
" SEGDREAD - read an SEG-D tape                                        ",
"                                                                      ",
" segdread > stdout tape=                                              ",
"                                                                      ",
"                                                                      ",
" Required parameters:                                                 ",
"       tape=           input tape device                              ",
" 			tape=- to read from stdin                      ",
"                                                                      ",
" Optional parameters:                                                 ",
"       use_stdio=0     for record devices (9-track reel tape drive)   ",
"                       =1 for pipe, disk and fixed-block 8mm drives   ",
"       verbose=0       silent operation                               ",
"                       = 1 ; echo every 'vblock' traces               ",
"                       = 2 ; echo information about blocks            ",
"       vblock=50       echo every 'vblock' traces under verbose option",
"       ptmin=1         first shot to read                             ",
"       ptmax=INT_MAX   last shot to read                              ",
"       gain=0          no application of gain                         ",
"       aux=0           no recovery of auxiliary traces                ",
"       errmax=0        allowable number of consecutive tape IO errors ",
"	ns=0		override of computed ns to work around SEG-D   ",
"			flaws.  Ignored when use_stdio=0.              ",
"	pivot_year=30   Use current century for 2 digit yrs less than  ",
"	                pivot_year, previous century otherwise.        ",  
"                                                                      ",
"         type:   sudoc segdread   for further information             ",
NULL};

/* Credits:
 *  for version 1:
 *    IPRA, Pau, France: Dominique Rousset, [e-mail]
 *  for version 2.2:
 *    EOST, Strasbourg, France: Celine Girard
 *  for versions 2.3:
 *    EOST, Strasbourg, France: Marc Schaming, [e-mail]
 *  for version 2.4:
 *    SEP, Stanford University: Stew Levin, [e-mail]
 *    a) Changed definitions of BCD_FF/BCD_FFFF
 *    b) Corrected decoding of general_header_1 in info_gh1.
 *    c) Changed buff=0 to use_stdio=1 to avoid confusion (stdio
 *       IS buffered I/O). Kept old buff= internally for backwards
 *       compatibility.
 *    d) Changed F8015 decoding of negative mantissas to avoid
 *       1 part in 2^14th decoding error on 2's complement platforms.
 *    e) Adapted F8015 to F0015 decoding routine. Unused, but now available.
 *    f) Use AT&T sfio package for tape read.
 *    g) Handle endian and wordsize dependencies portably (I think).
 *    g) Allow tape=- command line argument to accept input from stdin.
 *    h) Compute trace length explicitly from headers so that disk data
 *       input can work.
 *    i) Correct tape trace length calculation to account for demux
 *       trace header extensions.
 *    j) Fix a couple of typos in comments and selfdoc
 *    k) Added F8022 for 1 byte quaternary exponent demux.
 *    l) Added F8024 for 2 byte quaternary exponent demux.
 *    m) Added F8042 for 1 byte hexadecimal exponent demux.
 *    n) Added F8044 for 2 byte hexadecimal exponent demux.
 *    o) Added F8048 for 4 byte hexadecimal exponent demux.
 *    p) Added F8036 for 2 byte 2's complement integer demux.
 *    q) Added F8038 for 4 byte 2's complement integer demux.
 *    r) Added F8058 for 4 byte IEEE float demux.
 *    s) Added ns= parameter to work around bad SEG-D trace
 *       length specifications in headers.
 *  for version 2.5:
 *    SEP, Stanford University: Stew Levin, [e-mail]
 *    a) Added pivot_year to disambiguate decoding of 2-digit yrs
 *    b) Modified decode of 2-byte BCD to avoid endian problems
 *    c) Modified debug printout to fix endian BCD display problems
 *    d) Don't let dem_trace_header override ns specified on command line
 *    e) Removed extra factor of two in decoding of general_header_1.r
 *    f) Removed conditional disabling of sfio 
 *
 *--------------------------------------------------------------------
 * SEGDREAD: Version 2.1, 10/10/94
 *           Version 2.2, 17/08/95
 *           Version 2.3, 04/1997 Thu Apr 10 11:55:45 DFT 1997
 *           Version 2.4, 10/03/98 Tue Mar 10 1998
 *           Version 2.5, Feb 4, 2001
 *--------------------------------------------------------------------
 */
/**************** end self doc ***********************************/

/* subroutine prototypes */
static int bcd (unsigned char *ptr, int begin, int n) ;            
static void F0015_to_float (Sfio_t *from, float to[], int len);
static void F8015_to_float (Sfio_t *from, float to[], int len);
static void F8022_to_float (Sfio_t *from, float to[], int len);
static void F8024_to_float (Sfio_t *from, float to[], int len);
static void F8036_to_float (Sfio_t *from, float to[], int len);
static void F8038_to_float (Sfio_t *from, float to[], int len);
static void F8042_to_float (Sfio_t *from, float to[], int len);
static void F8044_to_float (Sfio_t *from, float to[], int len);
static void F8048_to_float (Sfio_t *from, float to[], int len);
static void F8058_to_float (Sfio_t *from, float to[], int len);

int get_gh1(general_header_1 * gh1, Sfio_t *tapeun) ;
int get_gh2(general_header_2 * gh2, Sfio_t *tapeun) ;
int get_ghn(general_header_n * ghn, Sfio_t *tapeun) ;
int get_gn_sn358(gen_head_sn358 * gh358, Sfio_t *tapeun) ;
int get_csh(channel_set_header * csh, Sfio_t *tapeun) ;
int get_ssh(sample_skew * ssh, Sfio_t *tapeun) ;
int get_ech(extended_header * ech, Sfio_t *tapeun) ;
int get_exh(external_header * exh, Sfio_t *tapeun) ;
int get_gt(general_trailer * gt, Sfio_t *tapeun) ;
int get_dth(dem_trace_header * dth, Sfio_t *tapeun) ;
int get_the(trace_header_ext * the, Sfio_t *tapeun) ;

void info_gh1(general_header_1 * gh1) ;
void info_gh2(general_header_2 * gh2) ;
void info_ghn(general_header_n * ghn) ;
void info_gn_sn358(gen_head_sn358 * gh358);
void info_csh(channel_set_header * csh) ;
void info_ssh(sample_skew * ssh) ;
void info_ech(extended_header * ech);
void info_exh(external_header * exh);
void info_gt(general_trailer * gt);
void info_dth(dem_trace_header * dth);

segy tr;

int
main(int argc, char **argv)
{
 general_header_1   segd_general_header_1;
 general_header_2   segd_general_header_2;
 general_header_n   segd_general_header_n;
 gen_head_sn358     segd_gen_head_sn358;
 channel_set_header *segd_channel_set_header;
 sample_skew        segd_sample_skew;
 extended_header    segd_extended_header;
 external_header    segd_external_header;
 general_trailer    segd_general_trailer;
 dem_trace_header   segd_dem_trace_header;
 trace_header_ext   segd_trace_header_ext;

 channel_set_header **csd = NULL;
			/* array[n_str][n_cs] of channel_set_header */

 char *tape=NULL;        /* name of raw tape device */
 int tapefd=0;           /* file descriptor for tape */
 Sfio_t  *tapeun;		/* input for Sfio_t reads. May be memory or stdio */
 Sfoff_t startpos;		/* for rewind on memory Sfio_t */

 short scan_type;        /* scan type number */
 short chan_set;         /* channel set number */

 register int i, j;
 register int i_scan, i_cs, i_tr;
 int i_ss;
 int nread;               /* bytes read */
 int ns;                                  /* number of data samples */
 int n_gh;         /* number of additional blocks in general header */
 int n_str;                      /* number of scan types per record */
 int n_cs;                  /* number of channel sets per scan type */ 
 int n_sk;               /* number of 32 byte field for sample skew */
 int n_ec;                                /* extended header length */
 int n_ex;                                /* external header length */
 int n_gt=0;                 /* number of blocks of general trailer */
 int n_chan;                            /* total number of channels */
 int itr;                                   /* current trace number */
 int ipt;                                    /* current shot number */
 int ptmin;                                  /* first trace to read */
 int ptmax;                                   /* last trace to read */
 int verbose;                   /* echo every ...                   */
 int vblock;                    /* ... vblock traces with verbose=1 */
 int buff;                   /* flag for buffered/unbuffered device */
 int gain;                                /* flag for applying gain */
 int aux;                      /* flag for keeping auxiliary traces */
 int errmax;                      /* max consecutive tape io errors */
 int pivot_year;		  /* for choosing correct century   */
 int errcount = 0;                    /* counter for tape io errors */
 int hdr1_i, hdr1_r;                   /* i and r decoded from hdr1 */
 int ns_override;				/* for trace length fudging */
 unsigned int nsamp_hdr1;                   /* number of samples/trace using */
                            /* only general header 1 i and r fields */
 unsigned int nsamp_hdr2;            /* number of samples/trace from general */
                      /* header 2 erl and general header 1 i fields */
 unsigned int nsamp_hdr358;   /* Sercel header rec_length plus hdr1_i fields */
 unsigned int nsamp_cs;      /* nsamp from channel set named in demux header */
 unsigned int nsamp_the;         /* nsamp from demux trace header extensions */

 float  mp;                                   /* descaling exponent */  
 float **mmp = NULL;    /* array[n_cs][n_str] for escaling exponent */

 char *bloc1;                              /* pointer on data block */
 

 /* Initialize */
 initargs(argc, argv);
 requestdoc(0); /* stdin not used */

 /* Make sure stdout is a file or pipe */
 switch(filestat(STDOUT)) {
 case TTY:
 	err("stdout can't be tty");
 break;
 case DIRECTORY:
 	err("stdout must be a file, not a directory");
 break;
 case BADFILETYPE:
 	err("stdout is illegal filetype");
 break;
 default: /* Others OK */
 break;
 }

 /* Set filenames */
 MUSTGETPARSTRING("tape",  &tape);

 /* Set parameters */
 if (!getparint("ptmin", &ptmin))        ptmin = 1;
 if (!getparint("ptmax", &ptmax))        ptmax = INT_MAX;
 if (!getparint("verbose", &verbose))    verbose = 0;
 if (verbose==2) {ptmax=ptmin; warn("ptmax set to ptmin for verbose=2");}
 if (!getparint("vblock", &vblock))      vblock = 50;
 if (verbose==2) vblock = 1;
 if (!getparint("use_stdio",&buff)) {
	if (!getparint("buff", &buff))          buff = 1;
 } else {
	buff = !buff;
 }
 if (!getparint("gain", &gain))          gain = 0;
 if (!getparint("aux", &aux))            aux = 0;
 if (!getparint("errmax", &errmax))      errmax = 0;
 if (!getparint("pivot_year", &pivot_year))      pivot_year = 30;
 if (!getparint("ns", &ns_override))      	 ns_override = 0;

 /* Allocate space for the record block */
 if ((bloc1 = alloc1(REC_L, sizeof(char))) == NULL) err("error at bloc1 allocation");
 if ((unsigned long) bloc1 % 2)
	warn("there may be a problem since bloc1 is not on a short boundary (%ul)",(unsigned long) bloc1);

 /* Open the tape */
 if ( STREQ(tape,"-") ) {
    if (buff) tapefd = fileno(stdin);
    else tapeun = sfstdin;
 } else {
	 if (buff) tapefd = eopen(tape, O_RDONLY, 0444);
	 else      {
			tapeun = sfopen((Sfio_t *) 0, tape, "rb");
			if (tapeun == ((Sfio_t *) 0)) err("Unable to open tape %s\n");
		     }
 }
 if (verbose) warn("tape opened successfully");

 /* Create Sfio_t input stream */
 if(buff)
	tapeun = sfnew((Sfio_t *) 0, (Void_t *) bloc1, (size_t) REC_L, 0, SF_STRING|SF_READ);
 startpos = sftell(tapeun);

 /* Read the traces */
 ipt = 0;    /*current shot number */
 itr = 0;    /*current trace number */
 while (ipt < ptmax) {

 	/************************
 	* Read the Header Block *
 	************************/
 
 	if (buff) {
 		if (-1 == (nread = (int) read(tapefd, (void *) bloc1, (size_t) REC_L))){
 			if (verbose)
 				warn("tape read error on header block from shot %d", (ipt+1));
 			if (++errcount > errmax)
 				err("exceeded maximum io errors");
 		} else { /* Reset counter on successful tape IO */
 			errcount = 0;
 		}
	(void) sfseek(tapeun,startpos,SEEK_SET);  /* reset Sfio_t pointer to start of block */
 	if (!nread) break; /* middle exit loop instead of mile-long while */
	}

 	/* General Header #1 */

	if ( EXIT_FAILURE == get_gh1(&segd_general_header_1, tapeun) ) break;

      ns = 0;
      nsamp_hdr1 = 0; nsamp_hdr2 = 0; nsamp_hdr358 = 0;
	hdr1_r = bcd(&segd_general_header_1.z_r1,1,3);
        hdr1_i = segd_general_header_1.i;
        if(hdr1_r != BCD_FFF && hdr1_i != 0) {
              hdr1_r *= 2;  /* range is 10 to 1990 */
              nsamp_hdr1 = (hdr1_r*512*16)/(10*hdr1_i)  + 1;
                         /* 20*r*512 msec /10*(i/16) msec */
	      if(ns == 0) ns = nsamp_hdr1;
          if(verbose) warn("nsamp_hdr1=%d\n",nsamp_hdr1);
	}
  

 	tr.fldr =   bcd ((unsigned char *) &(segd_general_header_1.f[0]), 0, 4);
 	tr.year =   bcd ((unsigned char *) &segd_general_header_1.yr, 0, 2);
        if(tr.year < pivot_year) tr.year += 2000;
        else tr.year += 1900;
 	n_gh =      bcd ((unsigned char *) &segd_general_header_1.gh_dy1, 0, 1);
 	tr.day =    bcd ((unsigned char *) &segd_general_header_1.gh_dy1, 1, 3);
 	tr.hour =   bcd ((unsigned char *) &segd_general_header_1.h, 0, 2);
 	tr.minute = bcd ((unsigned char *) &segd_general_header_1.mi, 0, 2);
 	tr.sec =    bcd ((unsigned char *) &segd_general_header_1.se, 0, 2);
 	tr.dt =     (segd_general_header_1.i*1000) >> 4;
 	n_str =     bcd ((unsigned char *) &segd_general_header_1.str, 0, 2);
 	n_cs =      bcd ((unsigned char *) &segd_general_header_1.cs, 0, 2);