suread.c

su 的源代码库

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

/***************************************************************************
SUREAD:  read SU input stream
***************************************************************************
code by J.E. Anderson, Mobil/CWP (June, 1993)
modified by J.E. Anderson, April, 1995, updated trace header words 
	and internal documentation
***************************************************************************/
#include     "cwpsegy.h"
#include     <cpromax.h>
#include     <pwd.h>

#if defined SOLARIS

#else
extern FILE *popen(char *cmnd,char *type);
#endif

void init_suread_(int *lensav, int *itooltype);
void exec_suread_(float *trace, int *ihead, float *rhead);
static FILE *suread_open(int ichoose,char *cmnd,char *path,char *prog); 
static void suread_inith(int *ind);
static void suread_puth(cwpsegy *tr, int *ihead,float *rhead,int *ind);
static char *suread_stringPrep(char *s);
#ifndef ABS
#define ABS(x) ((x) < 0 ? -(x) : (x))
#endif

BEGINPARMS
FILE *fp;
int ichoose;
char *cmnd;
char *datafile;
char *cmndfile;
int ipsort;
int idtyp;
int nread;
int nbytes;
int maxdtr;
int ipkey;
int iskey;
float samprat;
long numsmp;
long iens;
long ens_value;
long tr_in_ens;
int *ind;
cwpsegy *tr;
ENDPARMS(parms)

void init_suread_(int *lensav, int *itooltype)
/***************************************************************************
initialization phase for SUREAD
***************************************************************************/
{
	char *temp;
	exParGetInt("CHOOSE1",&parms->ichoose);
	exParGetString("CMND",&parms->cmnd);
	exParGetString("CMNDFILE",&temp);
	parms->cmndfile = suread_stringPrep(temp);
	exParGetString("DATAFILE",&temp);
	parms->datafile = suread_stringPrep(temp);
	exParGetInt("SORT",&parms->ipsort);
	exParGetInt("TYPE",&parms->idtyp);
	exParGetInt("NPENS",&parms->maxdtr);
	suread_stringPrep(parms->cmndfile);
	suread_stringPrep(parms->datafile);

/***************************************************************************
Allocate save buffers and open the input stream
***************************************************************************/
	parms->tr=(cwpsegy *)malloc(sizeof(cwpsegy));
	parms->ind=(int *)malloc(120*sizeof(int));
	parms->fp=suread_open(parms->ichoose,parms->cmnd,
				parms->datafile,parms->cmndfile);

/***************************************************************************
Read the first trace header.  Then read the trace data for the first trace.
***************************************************************************/
	parms->nread=fread(parms->tr,1,240,parms->fp);
	if(!parms->nread)
		exErrFatal("suread: Could not get first trace header.");
	parms->nbytes=parms->tr->ns*sizeof(float);
	parms->nread=fread(&parms->tr->data[0],1,parms->nbytes,parms->fp);
	if(!parms->nread)
		exErrFatal("suread: Could not get data for first trace.");
	if(parms->nread!=parms->nbytes)
		exErrFatal("suread: Read wrong number of bytes on first trace");
	parms->nbytes+=240;
	parms->nread=parms->nbytes;
	fprintf(stderr,
"Successfully read first trace. Sample rate = %d usec.  Number of samples = %d.\n",
parms->tr->dt,parms->tr->ns);
	parms->samprat=0.001*parms->tr->dt;
	parms->numsmp=parms->tr->ns;

/***************************************************************************
Initialize ensemble, common block and trace header information
***************************************************************************/
	parms->iens=1;
	parms->tr_in_ens=1;
	suread_inith(parms->ind);
	parms->iskey=1+hdrIndex("SEQNO");
	switch(parms->ipsort){
	case ICDP:
		parms->ipkey=1+hdrIndex("CDP");
		parms->ens_value=parms->tr->cdp;
		fprintf(stderr,"suread: Input ensembles based on cdp\n");
		break;
	case ISIN:
		parms->ipkey=1+hdrIndex("SIN");
		parms->ens_value=parms->tr->fldr;
		fprintf(stderr,"suread: Input ensembles based on fldr\n");
		break;
	case IRECSLOC:
		parms->ipkey=1+hdrIndex("REC_SLOC");
		parms->ens_value=parms->tr->wevel;
		fprintf(stderr,"suread: Input ensembles based on wevel\n");
		break;
	case IOFFSET:
		parms->ipkey=1+hdrIndex("OFFSET");
		parms->ens_value=parms->tr->offset;
		fprintf(stderr,"suread: Input ensembles based on offset\n");
		break;
	case ICHAN:
		parms->ipkey=1+hdrIndex("CHAN");
		parms->ens_value=parms->tr->tracf;
		fprintf(stderr,"suread: Input ensembles based on tracf\n");
		break;
	case IUNKNOWN:
		parms->ipkey=1+hdrIndex("TRACENO");
		parms->ens_value=parms->tr->tracl;
		fprintf(stderr,"suread: Input ensembles based on tracl\n");
		break;
	}
	globalRuntime->numsmp     	= parms->numsmp;
	globalRuntime->samprat     		= parms->samprat;
	globalRuntime->igeom_match 	= FALSE;
	globalRuntime->itrno_valid 		= FALSE;
	globalRuntime->idomain     		= ITX;
	globalRuntime->maxdtr      		= parms->maxdtr;
	globalRuntime->ipsort  		= parms->ipsort;
	globalRuntime->idtyp   		= parms->idtyp;
	globalRuntime->iskey		= parms->iskey;
	globalRuntime->ipkey		= parms->ipkey;

/***************************************************************************
Return the length of the saved common block and type of module.
***************************************************************************/
	*lensav = NPARMS(parms);
	*itooltype = INPUT;
	return;
}

void exec_suread_(float *trace, int *ihead, float *rhead)
/***************************************************************************
execution phase for SUREAD
***************************************************************************/
{
	int j; 
	long ens;

/***************************************************************************
Check to see if this is the last pass or the cleanup pass and close files, etc.
***************************************************************************/
	if (globalRuntime->cleanup) {
		(void)free(parms->ind);
		return;
	}
	if(!parms->nread){
		if(parms->ichoose==0) {
			j=pclose(parms->fp);
		}else{
			j=fclose(parms->fp);
		}
		exPipeMode();
		return;
	}

/***************************************************************************
Fix up the output trace header and output trace from information in save buffer
***************************************************************************/
	suread_puth(parms->tr,ihead,rhead,parms->ind);
	ihead[STDHDR(isource)]  = parms->iens;
	ihead[STDHDR(iseqno)]   = parms->tr_in_ens;
	rhead[STDHDR(itlive_s)] = 0.;
	rhead[STDHDR(itfull_s)] = 0.;
	rhead[STDHDR(itlive_e)] = 
		(globalRuntime->numsmp-1)*globalRuntime->samprat;
	rhead[STDHDR(itfull_e)] = 
		(globalRuntime->numsmp-1)*globalRuntime->samprat;
	for(j=0;j<globalRuntime->numsmp;j++)
		trace[j]=parms->tr->data[j];

readnext:
/***************************************************************************
Read one trace ahead for next pass.  Always need to know ahead of time 
about changes to a new ensemble or if an end of file will be encountered to
set end of ensemble and end of data trace header flags for current trace
so have to always check values for one trace ahead of where current trace 
is being output.  Zero buffer prior to read so traces shorter than first
can be read ok.  Have to bomb if encounter a trace longer than first trace.
***************************************************************************/
	for(j=0;j<globalRuntime->numsmp;j++) parms->tr->data[j]=0.;
	parms->nread = fread(parms->tr,1,parms->nbytes,parms->fp);
	if(parms->nread==0) {
		fprintf(stderr,"suread: Hit EOF on input stream.\n");
		fprintf(stderr,"suread: Read %ld ensembles.\n",parms->iens);
		fprintf(stderr,"suread: Key for last ensemble equals %ld\n",parms->ens_value);
	}else if(parms->nread!=parms->nbytes) {
		fprintf(stderr,"suread: Read wrong number of bytes from input stream.\n");
		fprintf(stderr,"suread: Can not handle variable length records!\n");
		if(parms->nread>parms->nbytes) 
			exErrFatal("suread: encountered trace longer than first trace read");
	}
	if(parms->nread) {
		switch(parms->ipsort){
		case ICDP: 
			ens = parms->tr->cdp;
			break;
		case ISIN:
			ens = parms->tr->fldr;
			break;
		case IRECSLOC:
			ens = parms->tr->wevel;
			break;
		case IOFFSET:
			ens = parms->tr->offset;
			break;
		case ICHAN:
			ens = parms->tr->tracf;
			break;
		case IUNKNOWN:
			ens = parms->tr->tracl;
			break;
		}
		if(parms->ens_value!=ens) {
/***************************************************************************
			Encountered new ensemble.  
			Current trace is at end of ensemble.
***************************************************************************/
			if(parms->tr_in_ens>globalRuntime->maxdtr) 
				fprintf(stderr,