c语言中su,在C程序中读写SU格式文件

/* Copyright (c) Colorado School of Mines, 2011.*/

/* All rights reserved. */

/* segy.h - include file for SEGY traces** declarations for:*typedef struct {} segy - the trace identification header*typedef struct {} bhed - binary header** Note:*If header words are added, run the makefile in this directory*to recreate hdr.h.** Reference:*K. M. Barry, D. A. Cavers and C. W. Kneale, "Special Report:*Recommended Standards for Digital Tape Formats",*Geophysics, vol. 40, no. 2 (April 1975), P. 344-352.** $Author: john $* $Source: /usr/local/cwp/src/su/include/RCS/segy.h,v $* $Revision: 1.33 $ ; $Date: 2011/11/11 23:56:14 $*/

#include #include "par.h"

#ifndef SEGY_H #define SEGY_H #define TRCBYTES240

#define SU_NFLTS32767/* Arbitrary limit on data array size*/

/* TYPEDEFS */

typedef struct {/* segy - trace identification header */

int tracl;/* Trace sequence number within line--numbers continue to increase if thesame line continues across multipleSEG Y files.byte# 1-4*/

int tracr;/* Trace sequence number within SEG Y file---each file starts with trace sequenceonebyte# 5-8*/

int fldr;/* Original field record numberbyte# 9-12*/

int tracf;/* Trace number within original field recordbyte# 13-16*/

int ep;/* energy source point number---Used when more than one record occursat the same effective surface location.byte# 17-20*/

int cdp;/* Ensemble number (i.e. CDP, CMP, CRP,...)byte# 21-24*/

int cdpt;/* trace number within the ensemble---each ensemble starts with trace number one.byte# 25-28*/

short trid;/* trace identification code:-1 = Other0 = Unknown1 = Seismic data2 = Dead3 = Dummy4 = Time break5 = Uphole6 = Sweep7 = Timing8 = Water break9 = Near-field gun signature10 = Far-field gun signature11 = Seismic pressure sensor12 = Multicomponent seismic sensor- Vertical component13 = Multicomponent seismic sensor- Cross-line component14 = Multicomponent seismic sensor- in-line component15 = Rotated multicomponent seismic sensor- Vertical component16 = Rotated multicomponent seismic sensor- Transverse component17 = Rotated multicomponent seismic sensor- Radial component18 = Vibrator reaction mass19 = Vibrator baseplate20 = Vibrator estimated ground force21 = Vibrator reference22 = Time-velocity pairs23 ... N = optional use(maximum N = 32,767)Following are CWP id flags:109 = autocorrelation110 = Fourier transformed - no packingxr[0],xi[0], ..., xr[N-1],xi[N-1]111 = Fourier transformed - unpacked Nyquistxr[0],xi[0],...,xr[N/2],xi[N/2]112 = Fourier transformed - packed Nyquisteven N:xr[0],xr[N/2],xr[1],xi[1], ...,xr[N/2 -1],xi[N/2 -1](note the exceptional second entry)odd N:xr[0],xr[(N-1)/2],xr[1],xi[1], ...,xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2](note the exceptional second & last entries)113 = Complex signal in the time domainxr[0],xi[0], ..., xr[N-1],xi[N-1]114 = Fourier transformed - amplitude/phasea[0],p[0], ..., a[N-1],p[N-1]115 = Complex time signal - amplitude/phasea[0],p[0], ..., a[N-1],p[N-1]116 = Real part of complex trace from 0 to Nyquist117 = Imag part of complex trace from 0 to Nyquist118 = Amplitude of complex trace from 0 to Nyquist119 = Phase of complex trace from 0 to Nyquist121 = Wavenumber time domain (k-t)122 = Wavenumber frequency (k-omega)123 = Envelope of the complex time trace124 = Phase of the complex time trace125 = Frequency of the complex time trace130 = Depth-Range (z-x) traces201 = Seismic data packed to bytes (by supack1)202 = Seismic data packed to 2 bytes (by supack2)byte# 29-30*/

short nvs;/* Number of vertically summed traces yieldingthis trace. (1 is one trace,2 is two summed traces, etc.)byte# 31-32*/

short nhs;/* Number of horizontally summed traces yieldingthis trace. (1 is one trace2 is two summed traces, etc.)byte# 33-34*/

short duse;/* Data use:1 = Production2 = Testbyte# 35-36*/

int offset;/* Distance from the center of the source pointto the center of the receiver group(negative if opposite to direction in whichthe line was shot).byte# 37-40*/

int gelev;/* Receiver group elevation from sea level(all elevations above the Vertical datum arepositive and below are negative).byte# 41-44*/

int selev;/* Surface elevation at source.byte# 45-48*/

int sdepth;/* Source depth below surface (a positive number).byte# 49-52*/

int gdel;/* Datum elevation at receiver group.byte# 53-56*/

int sdel;/* Datum elevation at source.byte# 57-60*/

int swdep;/* Water depth at source.byte# 61-64*/

int gwdep;/* Water depth at receiver group.byte# 65-68*/

short scalel;/* Scalar to be applied to the previous 7 entriesto give the real value.Scalar = 1, +10, +100, +1000, +10000.If positive, scalar is used as a multiplier,if negative, scalar is used as a divisor.byte# 69-70*/

short scalco;/* Scalar to be applied to the next 4 entriesto give the real value.Scalar = 1, +10, +100, +1000, +10000.If positive, scalar is used as a multiplier,if negative, scalar is used as a divisor.byte# 71-72*/

int sx;/* Source coordinate - Xbyte# 73-76*/

int sy;/* Source coordinate - Ybyte# 77-80*/

int gx;/* Group coordinate - Xbyte# 81-84*/

int gy;/* Group coordinate - Ybyte# 85-88*/

short counit;/* Coordinate units: (for previous 4 entries andfor the 7 entries before scalel)1 = Length (meters or feet)2 = Seconds of arc3 = Decimal degrees4 = Degrees, minutes, seconds (DMS)In case 2, the X values are longitude andthe Y values are latitude, a positive value designatesthe number of seconds east of Greenwichor north of the equatorIn case 4, to encode +-DDDMMSScounit = +-DDD*10^4 + MM*10^2 + SS,with scalco = 1. To encode +-DDDMMSS.sscounit = +-DDD*10^6 + MM*10^4 + SS*10^2with scalco = -100.byte# 89-90*/

short wevel;/* Weathering velocity.byte# 91-92*/

short swevel;/* Subweathering velocity.byte# 93-94*/

short sut;/* Uphole time at source in milliseconds.byte# 95-96*/

short gut;/* Uphole time at receiver group in milliseconds.byte# 97-98*/

short sstat;/* Source static correction in milliseconds.byte# 99-100*/

short gstat;/* Group static correction in milliseconds.byte# 101-102*/

short tstat;/* Total static applied in milliseconds.(Zero if no static has been applied.)byte# 103-104*/

short laga;/* Lag time A, time in ms between end of 240-byte trace identification header and timebreak, positive if time break occurs afterend of header, time break is defined asthe initiation pulse which maybe recordedon an auxiliary trace or as otherwisespecified by the recording systembyte# 105-106*/

short lagb;/* lag time B, time in ms between the time breakand the initiation time of the energy source,may be positive or negativebyte# 107-108*/

short delrt;/* delay recording time, time in ms betweeninitiation time of energy source and timewhen recording of data samples begins(for deep water work if recording does notstart at zero time)byte# 109-110*/

short muts;/* mute time--startbyte# 111-112*/

short mute;/* mute time--endbyte# 113-114*/

unsigned short ns;/* number of samples in this tracebyte# 115-116*/

unsigned short dt;/* sample interval; in micro-secondsbyte# 117-118*/

short gain;/* gain type of field instruments code:1 = fixed2 = binary3 = floating point4 ---- N = optional usebyte# 119-120*/

short igc;/* instrument gain constantbyte# 121-122*/

short igi;/* instrument early or initial gainbyte# 123-124*/

short corr;/* correlated:1 = no2 = yesbyte# 125-126*/

short sfs;/* sweep frequency at startbyte# 127-128*/

short sfe;/* sweep frequency at endbyte# 129-130*/

short slen;/* sweep length in msbyte# 131-132*/

short styp;/* sweep type code:1 = linear2 = cos-squared3 = otherbyte# 133-134*/

short stas;/* sweep trace length at start in msbyte# 135-136*/

short stae;/* sweep trace length at end in msbyte# 137-138*/

short tatyp;/* taper type: 1=linear, 2=cos^2, 3=otherbyte# 139-140*/

short afilf;/* alias filter frequency if usedbyte# 141-142*/

short afils;/* alias filter slopebyte# 143-144*/

short nofilf;/* notch filter frequency if usedbyte# 145-146*/

short nofils;/* notch filter slopebyte# 147-148*/

short lcf;/* low cut frequency if usedbyte# 149-150*/

short hcf;/* high cut frequncy if usedbyte# 151-152*/

short lcs;/* low cut slopebyte# 153-154*/

short hcs;/* high cut slopebyte# 155-156*/

short year;/* year data recordedbyte# 157-158*/

short day;/* day of yearbyte# 159-160*/

short hour;/* hour of day (24 hour clock)byte# 161-162*/

short minute;/* minute of hourbyte# 163-164*/

short sec;/* second of minutebyte# 165-166*/

short timbas;/* time basis code:1 = local2 = GMT3 = otherbyte# 167-168*/

short trwf;/* trace weighting factor, defined as 1/2^Nvolts for the least sigificant bitbyte# 169-170*/

short grnors;/* geophone group number of roll switchposition onebyte# 171-172*/

short grnofr;/* geophone group number of trace one withinoriginal field recordbyte# 173-174*/

short grnlof;/* geophone group number of last trace withinoriginal field recordbyte# 175-176*/

short gaps;/* gap size (total number of groups dropped)byte# 177-178*/

short otrav;/* overtravel taper code:1 = down (or behind)2 = up (or ahead)byte# 179-180*/

#ifdef SLTSU_SEGY_H/* begin Unocal SU segy.h differences */

/* cwp local assignments */

float d1;/* sample spacing for non-seismic databyte# 181-184*/

float f1;/* first sample location for non-seismic databyte# 185-188*/

float d2;/* sample spacing between tracesbyte# 189-192*/

float f2;/* first trace locationbyte# 193-196*/

float ungpow;/* negative of power used for dynamicrange compressionbyte# 197-200*/

float unscale;/* reciprocal of scaling factor to normalizerangebyte# 201-204*/

short mark;/* mark selected tracesbyte# 205-206*/

/* SLTSU local assignments */

short mutb;/* mute time at bottom (start time)bottom mute ends at last samplebyte# 207-208*/

float dz;/* depth sampling interval in (m or ft)if =0.0, input are time samplesbyte# 209-212*/

float fz;/* depth of first sample in (m or ft)byte# 213-116*/

short n2;/* number of traces per cdp or per shotbyte# 217-218*/

short shortpad; /* alignment paddingbyte# 219-220*/

int ntr; /* number of tracesbyte# 221-224*/

/* SLTSU local assignments end */

short unass[8];/* unassignedbyte# 225-240*/

#else

/* cwp local assignments */

float d1;/* sample spacing for non-seismic databyte# 181-184*/

float f1;/* first sample location for non-seismic databyte# 185-188*/

float d2;/* sample spacing between tracesbyte# 189-192*/

float f2;/* first trace locationbyte# 193-196*/

float ungpow;/* negative of power used for dynamicrange compressionbyte# 197-200*/

float unscale;/* reciprocal of scaling factor to normalizerangebyte# 201-204*/

int ntr; /* number of tracesbyte# 205-208*/

short mark;/* mark selected tracesbyte# 209-210*/

short shortpad; /* alignment paddingbyte# 211-212*/

short unass[14];/* unassigned--NOTE: last entry causesa break in the word alignment, if we REALLYwant to maintain 240 bytes, the followingentry should be an odd number of short/UINT2OR do the insertion above the "mark" keywordentrybyte# 213-240*/

#endif

} segy;

typedef struct {/* bhed - binary header */

int jobid;/* job identification number */

int lino;/* line number (only one line per reel) */

int reno;/* reel number */

short ntrpr;/* number of data traces per record */

short nart;/* number of auxiliary traces per record */

unsigned short hdt; /* sample interval in micro secs for this reel */

unsigned short dto; /* same for original field recording */

unsigned short hns; /* number of samples per trace for this reel */

unsigned short nso; /* same for original field recording */

short format;/* data sample format code:1 = floating point, 4 byte (32 bits)2 = fixed point, 4 byte (32 bits)3 = fixed point, 2 byte (16 bits)4 = fixed point w/gain code, 4 byte (32 bits)5 = IEEE floating point, 4 byte (32 bits)8 = two's complement integer, 1 byte (8 bits)*/

short fold;/* CDP fold expected per CDP ensemble */

short tsort;/* trace sorting code:1 = as recorded (no sorting)2 = CDP ensemble3 = single fold continuous profile4 = horizontally stacked */

short vscode;/* vertical sum code:1 = no sum2 = two sum ...N = N sum (N = 32,767) */

short hsfs;/* sweep frequency at start */

short hsfe;/* sweep frequency at end */

short hslen;/* sweep length (ms) */

short hstyp;/* sweep type code:1 = linear2 = parabolic3 = exponential4 = other */

short schn;/* trace number of sweep channel */

short hstas;/* sweep trace taper length at start iftapered (the taper starts at zero timeand is effective for this length) */

short hstae;/* sweep trace taper length at end (the endingtaper starts at sweep length minus the taperlength at end) */

short htatyp;/* sweep trace taper type code:1 = linear2 = cos-squared3 = other */

short hcorr;/* correlated data traces code:1 = no2 = yes */

short bgrcv;/* binary gain recovered code:1 = yes2 = no */

short rcvm;/* amplitude recovery method code:1 = none2 = spherical divergence3 = AGC4 = other */

short mfeet;/* measurement system code:1 = meters2 = feet */

short polyt;/* impulse signal polarity code:1 = increase in pressure or upwardgeophone case movement givesnegative number on tape2 = increase in pressure or upwardgeophone case movement givespositive number on tape */

short vpol;/* vibratory polarity code:codeseismic signal lags pilot by1337.5 to 22.5 degrees2 22.5 to 67.5 degrees3 67.5 to 112.5 degrees4112.5 to 157.5 degrees5157.5 to 202.5 degrees6202.5 to 247.5 degrees7247.5 to 292.5 degrees8293.5 to 337.5 degrees */

short hunass[170];/* unassigned */

} bhed;

/* DEFINES */

#define gettr(x)fgettr(stdin, (x)) #define vgettr(x)fvgettr(stdin, (x)) #define puttr(x)fputtr(stdout, (x)) #define vputtr(x)fvputtr(stdout, (x)) #define gettra(x, y) fgettra(stdin, (x), (y))

/* TOTHER represents "other"*/

#defineTOTHER-1 /* TUNK represents time traces of an unknown type*/

#defineTUNK0 /* TREAL represents real time traces */

#defineTREAL1 /* TDEAD represents dead time traces */

#defineTDEAD2 /* TDUMMY represents dummy time traces */

#defineTDUMMY3 /* TBREAK represents time break traces */

#defineTBREAK4 /* UPHOLE represents uphole traces */

#defineUPHOLE5 /* SWEEP represents sweep traces */

#defineSWEEP6 /* TIMING represents timing traces */

#defineTIMING7 /* WBREAK represents timing traces */

#defineWBREAK8 /* NFGUNSIG represents near field gun signature */

#defineNFGUNSIG9 /* FFGUNSIG represents far field gun signature */

#defineFFGUNSIG10 /* SPSENSOR represents seismic pressure sensor */

#defineSPSENSOR11 /* TVERT represents multicomponent seismic sensor- vertical component */

#defineTVERT12 /* TXLIN represents multicomponent seismic sensor- cross-line component */

#defineTXLIN13 /* TINLIN represents multicomponent seismic sensor- in-line component */

#defineTINLIN14 /* ROTVERT represents rotated multicomponent seismic sensor- vertical component */

#defineROTVERT15 /* TTRANS represents rotated multicomponent seismic sensor- transverse component */

#defineTTRANS16 /* TRADIAL represents rotated multicomponent seismic sensor- radial component */

#defineTRADIAL17 /* VRMASS represents vibrator reaction mass */

#defineVRMASS18 /* VBASS represents vibrator baseplate */

#defineVBASS19 /* VEGF represents vibrator estimated ground force */

#defineVEGF20 /* VREF represents vibrator reference */

#defineVREF21

/*** CWP trid assignments ***/

/* ACOR represents autocorrelation */

#defineACOR109 /* FCMPLX represents fourier transformed - no packingxr[0],xi[0], ..., xr[N-1],xi[N-1] */

#defineFCMPLX110 /* FUNPACKNYQ represents fourier transformed - unpacked Nyquistxr[0],xi[0],...,xr[N/2],xi[N/2] */

#defineFUNPACKNYQ111 /* FTPACK represents fourier transformed - packed Nyquisteven N: xr[0],xr[N/2],xr[1],xi[1], ...,xr[N/2 -1],xi[N/2 -1](note the exceptional second entry)odd N:xr[0],xr[(N-1)/2],xr[1],xi[1], ...,xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2](note the exceptional second & last entries)*/

#defineFTPACK112 /* TCMPLX represents complex time traces */

#defineTCMPLX113 /* FAMPH represents freq domain data in amplitude/phase form*/

#defineFAMPH114 /* TAMPH represents time domain data in amplitude/phase form*/

#defineTAMPH115 /* REALPART represents the real part of a trace to Nyquist*/

#defineREALPART116 /* IMAGPART represents the real part of a trace to Nyquist*/

#defineIMAGPART117 /* AMPLITUDE represents the amplitude of a trace to Nyquist*/

#defineAMPLITUDE118 /* PHASE represents the phase of a trace to Nyquist*/

#definePHASE119 /* KT represents wavenumber-time domain data */

#defineKT121 /* KOMEGA represents wavenumber-frequency domain data*/

#defineKOMEGA122 /* ENVELOPE represents the envelope of the complex time trace*/

#defineENVELOPE123 /* INSTPHASE represents the phase of the complex time trace*/

#defineINSTPHASE124 /* INSTFREQ represents the frequency of the complex time trace*/

#defineINSTFREQ125 /* DEPTH represents traces in depth-range (z-x)*/

#defineTRID_DEPTH130 /* 3C data... v,h1,h2=(11,12,13)+32 so a bitmask will convert */

/* between conventions */

/* CHARPACK represents byte packed seismic data from supack1*/

#defineCHARPACK201 /* SHORTPACK represents 2 byte packed seismic data from supack2*/

#defineSHORTPACK202

#define ISSEISMIC(id) (( (id)==TUNK || (id)==TREAL || (id)==TDEAD || (id)==TDUMMY || (id)==TBREAK || (id)==UPHOLE || (id)==SWEEP || (id)==TIMING || (id)==WBREAK || (id)==NFGUNSIG || ( id)==FFGUNSIG || (id)==SPSENSOR || (id)==TVERT || (id)==TXLIN || (id)==TINLIN || (id)==ROTVERT || (id)==TTRANS || (id)==TRADIAL || (id)==ACOR ) ? cwp_true : cwp_false )

/* FUNCTION PROTOTYPES */

#ifdef __cplusplus/* if C++, specify external linkage to C functions */ extern "C" {

#endif

/* get trace and put trace */

int fgettr(FILE *fp, segy *tp);

int fvgettr(FILE *fp, segy *tp);

void fputtr(FILE *fp, segy *tp);

void fvputtr(FILE *fp, segy *tp);

int fgettra(FILE *fp, segy *tp, int itr);

/* get gather and put gather */

segy **fget_gather(FILE *fp, cwp_String *key,cwp_String *type,Value *n_val,

int *nt,int *ntr, float *dt,int *first);

segy **get_gather(cwp_String *key, cwp_String *type, Value *n_val,

int *nt, int *ntr, float *dt, int *first);

segy **fput_gather(FILE *fp, segy **rec,int *nt, int *ntr);

segy **put_gather(segy **rec,int *nt, int *ntr);

/* hdrpkge */

void gethval(const segy *tp, int index, Value *valp);

void puthval(segy *tp, int index, Value *valp);

void getbhval(const bhed *bhp, int index, Value *valp);

void putbhval(bhed *bhp, int index, Value *valp);

void gethdval(const segy *tp, char *key, Value *valp);

void puthdval(segy *tp, char *key, Value *valp);

char *hdtype(const char *key);

char *getkey(const int index);

int getindex(const char *key);

void swaphval(segy *tp, int index);

void swapbhval(bhed *bhp, int index);

void printheader(const segy *tp);

void tabplot(segy *tp, int itmin, int itmax);

#ifdef __cplusplus/* if C++, end external linkage specification */ }

#endif

#endif