Titan2Reader: segy.c Source File

00001 /*
00002  * This file is part of Titan2Reader
00003  *
00004  * Copyright (C) 2006, Sebastien Judenherc <sebastien.judenherc@agecodagis.com>
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
00019  * USA
00020  *
00021  */
00022 #include <stdio.h>
00023 #include <ctype.h>
00024 #include <stdlib.h>
00025 #include <unistd.h>
00026 #include <time.h>
00027 #include <string.h>
00028 #include <assert.h>
00029 #include <sys/types.h>
00030 
00031 #include "common.h"
00032 
00033 
00034 /* PASSCAL SEGY header: http://www.mpl.ucsd.edu/obs/reports/tn061/p-segy.html */
00035 
00036 /* This is the header for the PASSCAL SEGY trace data.  
00037  *
00038  * The PASSCAL SEGY trace format is a modified form of the SEG-Y
00039  * trace format.  The modification comes is because we use some of
00040  * the unspecified header words to store information pertinent to 
00041  * the PASSCAL data.  The data values for each trace are preceded
00042  * by a 240 byte header. This format is given below.  All integer 
00043  * values are stored with the most significant byte first.  Data
00044  * values are either 16 0r 32 bit integers depending on byte 206
00045  * of the header, the field named "data_form".
00046  * 
00047  * SEGYHEAD is now typedef'ed 
00048  *
00049  * Reading bytes  directly into this header will allow access to
00050  * all of the fields.  The number in the comment is the byte offset
00051  * into the segy file.  An "X" in the comment indicates that field
00052  * is NEVER set.  Fields that are set to default values contain that
00053  * value and follow a ":" in the comment ("grep : segy.h" will spit
00054  * all default fields out).  Two pairs of fields exist to cover an
00055  * inherited limitation; sampleLength/num_samps and deltaSample/samp_rate.
00056  * When the value is too large to fit in a short, sampleLength or
00057  * deltaSample become flags and require their long counterparts,
00058  * num_samps and samp_rate, to contain that value.
00059  */
00060 //typedef struct SegyHead {   /*  Offset Description  */
00061 struct SegyHead {   /*  Offset Description  */
00062    int32_t  lineSeq;           /*   0 Sequence numbers within line */
00063    int32_t  reelSeq;           /*   4 Sequence numbers within reel */
00064    int32_t  event_number;      /*   8 Original field record number or trigger number */
00065    int32_t  channel_number;    /*  12 Trace channel number within the original field record */
00066    int32_t  energySourcePt;    /*  16 X */
00067    int32_t  cdpEns;            /*  20 X */
00068    int32_t  traceInEnsemble;   /*  24 X */
00069    int16_t traceID;           /*  28 Trace identification code: seismic data = 1 */
00070    int16_t vertSum;           /*  30 X */
00071    int16_t horSum;            /*  32 X */
00072    int16_t dataUse;           /*  34 X */
00073    int32_t  sourceToRecDist;   /*  36 X */
00074    int32_t  recElevation;      /*  40 X */
00075    int32_t  sourceSurfaceElevation; /*  44 X */
00076    int32_t  sourceDepth;       /*  48 X */
00077    int32_t  datumElevRec;      /*  52 X */
00078    int32_t  datumElevSource;   /*  56 X */
00079    int32_t  sourceWaterDepth;  /*  60 X */
00080    int32_t  recWaterDepth;     /*  64 X */
00081    int16_t elevationScale;    /*  68 Elevation Scaler: scale = 1 */
00082    int16_t coordScale;        /*  70 Coordinate Scaler: scale = 1 */
00083    int32_t  sourceLongOrX;     /*  72 X */
00084    int32_t  sourceLatOrY;      /*  76 X */
00085    int32_t  recLongOrX;        /*  80 X */
00086    int32_t  recLatOrY;         /*  84 X */
00087    int16_t coordUnits;        /*  88 Coordinate Units:  = 2 (Lat/Long) */
00088    int16_t weatheringVelocity;/*  90 X */
00089    int16_t subWeatheringVelocity; /*  92 X */
00090    int16_t sourceUpholeTime;  /*  94 X */
00091    int16_t recUpholeTime;     /*  96 X */
00092    int16_t sourceStaticCor;   /*  98 X */
00093    int16_t recStaticCor;      /* 100 X */
00094    int16_t totalStatic;       /* 102 Total Static in MILLISECS added to Trace Start Time (lower 2 bytes)*/
00095    int16_t lagTimeA;          /* 104 X */
00096    int16_t lagTimeB;          /* 106 X */
00097    int16_t delay;             /* 108 X */
00098    int16_t muteStart;         /* 110 X */
00099    int16_t muteEnd;           /* 112 X */
00100    int16_t sampleLength;      /* 114 Number of samples in this trace (unless == 32767) */
00101    int16_t deltaSample;       /* 116 Sampling interval in MICROSECONDS (unless == 1) */
00102    int16_t gainType;          /* 118 Gain Type: 1 = Fixed Gain */
00103    int16_t gainConst;         /* 120 Gain of amplifier */
00104    int16_t initialGain;       /* 122 X */
00105    int16_t correlated;        /* 124 X */
00106    int16_t sweepStart;        /* 126 X */
00107    int16_t sweepEnd;          /* 128 X */
00108    int16_t sweepLength;       /* 130 X */
00109    int16_t sweepType;         /* 132 X */
00110    int16_t sweepTaperAtStart; /* 134 X */
00111    int16_t sweepTaperAtEnd;   /* 136 X */
00112    int16_t taperType;         /* 138 X */
00113    int16_t aliasFreq;         /* 140 X */
00114    int16_t aliasSlope;        /* 142 X */
00115    int16_t notchFreq;         /* 144 X */
00116    int16_t notchSlope;        /* 146 X */
00117    int16_t lowCutFreq;        /* 148 X */
00118    int16_t hiCutFreq;         /* 150 X */
00119    int16_t lowCutSlope;       /* 152 X */
00120    int16_t hiCutSlope;        /* 154 X */
00121    int16_t year;              /* 156 year of Start of trace */
00122    int16_t day;               /* 158 day of year at Start of trace */
00123    int16_t hour;              /* 160 hour of day at Start of trace */
00124    int16_t minute;            /* 162 minute of hour at Start of trace */
00125    int16_t second;            /* 164 second of minute at Start of trace */
00126    int16_t timeBasisCode;     /* 166 Time basis code: 2 = GMT */
00127    int16_t traceWeightingFactor; /* 168 X */
00128    int16_t phoneRollPos1;     /* 170 X */
00129    int16_t phoneFirstTrace;   /* 172 X */
00130    int16_t phoneLastTrace;    /* 174 X */
00131    int16_t gapSize;           /* 176 X */
00132    int16_t taperOvertravel;   /* 178 X */
00133    char  station_name[6];   /* 180 Station Name code (5 chars + \0) */
00134    char  sensor_serial[8];  /* 186 Sensor Serial code (7 chars + \0) */
00135    char  channel_name[4];   /* 194 Channel Name code (3 chars + \0) */
00136    int16_t totalStaticHi;     /* 198 Total Static in MILLISECS added to Trace Start Time (high 2 bytes)*/
00137    int32_t  samp_rate;         /* 200 Sample interval in MICROSECS as a 32 bit integer */
00138    int16_t data_form;         /* 204 Data Format flag: 0=16 bit, 1=32 bit integer */
00139    int16_t m_secs;            /* 206 MILLISECONDS of seconds of Start of trace */
00140    int16_t trigyear;          /* 208 year of Trigger time */
00141    int16_t trigday;           /* 210 day of year at Trigger time */
00142    int16_t trighour;          /* 212 hour of day at Trigger time */
00143    int16_t trigminute;        /* 214 minute of hour at Trigger time */
00144    int16_t trigsecond;        /* 216 second of minute at Trigger time */
00145    int16_t trigmills;         /* 218 MILLISECONDS of seconds of Trigger time */
00146    float scale_fac;         /* 220 Scale Factor (IEEE 32 bit float) */
00147    int16_t inst_no;           /* 224 Instrument Serial Number */
00148    int16_t not_to_be_used;    /* 226 X */
00149    int32_t  num_samps;         /* 228 Number of Samples as a 32 bit integer
00150                              * (when sampleLength == 32767) */
00151    int32_t  max;               /* 232 Maximum value in Counts */
00152    int32_t  min;               /* 236 Minimum value in Counts */
00153 }; // SEGYHEAD;                 /* end of segy trace header */
00154 
00155 struct SegyHead segyHeader[MAX_CHANNEL];
00156 FILE *segyChannelFile[MAX_CHANNEL] = {
00157         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00158         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00159         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00160         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00161         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00162 };
00163 char segyChannelFileName[MAX_CHANNEL][128];
00164 
00165 
00166 void UpdateSegyHeader(int ch)
00167 {
00168         time_t T0;
00169         struct tm tm0;
00170         UpdateInfoHeader(ch);
00171         T0=(time_t)infoHeader[ch].correctedTime;
00172         memcpy(&tm0,gmtime(&T0),sizeof(struct tm));
00173 
00174 
00175         segyHeader[ch].year=    tm0.tm_year+1900;       /* year data recorded */
00176         segyHeader[ch].day=     tm0.tm_yday+1;          /* day of year */
00177         segyHeader[ch].hour=    tm0.tm_hour;            /* hour of day (24 hour clock) */
00178         segyHeader[ch].minute=  tm0.tm_min;             /* minute of hour */
00179         segyHeader[ch].second=  tm0.tm_sec;             /* second of minute */
00180         segyHeader[ch].m_secs=  (int32_t)(1000.*(infoHeader[ch].correctedTime-(double)T0));
00181 
00182    //??? //int16_t trigyear;          /* 208 year of Trigger time */
00183    //??? //int16_t trigday;           /* 210 day of year at Trigger time */
00184    //??? //int16_t trighour;          /* 212 hour of day at Trigger time */
00185    //??? //int16_t trigminute;        /* 214 minute of hour at Trigger time */
00186    //??? //int16_t trigsecond;        /* 216 second of minute at Trigger time */
00187    //??? //int16_t trigmills;         /* 218 MILLISECONDS of seconds of Trigger time */
00188 
00189         //int16_t sampleLength;      /* 114 Number of samples in this trace (unless == 32767) */
00190         //int32_t  num_samps;         /* 228 Number of Samples as a 32 bit integer (when sampleLength == 32767) */
00191         if (infoHeader[ch].npts<32767)
00192         {
00193                 segyHeader[ch].sampleLength=infoHeader[ch].npts;
00194                 segyHeader[ch].num_samps=infoHeader[ch].npts;
00195         } else
00196         {
00197                 segyHeader[ch].sampleLength=32767;
00198                 segyHeader[ch].num_samps=infoHeader[ch].npts;
00199         }
00200 }
00201 
00202 void InitSegyHeader(int ch)
00203 {
00204         memset(&(segyHeader[ch]),0,sizeof(struct SegyHead));
00205         segyHeader[ch].traceID=1;
00206         segyHeader[ch].elevationScale=1;
00207         segyHeader[ch].coordScale=1;
00208         segyHeader[ch].coordUnits=2;
00209         segyHeader[ch].gainType=1;
00210         segyHeader[ch].gainConst=1;
00211         segyHeader[ch].timeBasisCode=2;
00212         snprintf(segyHeader[ch].station_name,5,"%-5s",infoHeader[ch].station);
00213         segyHeader[ch].station_name[5]=0;
00214         segyHeader[ch].sensor_serial[0]=0;
00215         snprintf(segyHeader[ch].channel_name,3,"%-3s",infoHeader[ch].channel);
00216         segyHeader[ch].channel_name[3]=0;
00217         segyHeader[ch].samp_rate=
00218                 ((1000000*infoHeader[ch].sr.div)/infoHeader[ch].sr.base);
00219         // check for deltaSample width
00220         if (segyHeader[ch].samp_rate>32768)
00221         {
00222                 segyHeader[ch].deltaSample=1;
00223         }
00224         else
00225         {
00226                 segyHeader[ch].deltaSample=segyHeader[ch].samp_rate;
00227         }
00228         segyHeader[ch].data_form=1;
00229         segyHeader[ch].scale_fac=1.;
00230 }
00231 
00235 FILE* OpenSegyChannel(int channel)
00236 {
00237         if (segyChannelFile[channel]==NULL)
00238         {
00239                 int fd;
00240                 strcpy(segyChannelFileName[channel],"titanfileXXXXXX");
00241                 fd=OpenTempFile(segyChannelFileName[channel]);
00242                 assert((segyChannelFile[channel]=fdopen(fd,"wb"))!=NULL);
00243                 setvbuf(segyChannelFile[channel],NULL,_IOFBF,FILE_BUFFER_SIZE);
00244 
00245                 if (segyChannelFile[channel])
00246                 {
00247                         // initialized a new generic header
00248                         InitInfoHeader(channel);
00249                         // initiliaze a new segy header
00250                         InitSegyHeader(channel);
00251                         // and write the header to the file
00252                         fwrite(&segyHeader[channel],sizeof(struct SegyHead),1,segyChannelFile[channel]);
00253                 }
00254                 else
00255                 {
00256                         PrintError("unable to open %s: %m\n",segyChannelFileName[channel]);
00257                 }
00258         }
00259         return(segyChannelFile[channel]);
00260 }
00261 
00265 void CloseSegyFile(int channel)
00266 {
00267         char line[128];
00268         if (segyChannelFile[channel]==NULL)
00269         {
00270                 return;
00271         }
00272         // update the internal information header and the segy start time
00273         UpdateSegyHeader(channel);
00274         // then copy the appropriate fields to the SEGY header
00275 
00276         // rewrite the header
00277         fseek(segyChannelFile[channel],0,SEEK_SET);
00278         fwrite(&segyHeader[channel],sizeof(struct SegyHead),1,segyChannelFile[channel]);
00279         fclose(segyChannelFile[channel]);
00280         segyChannelFile[channel]=NULL;
00281         strcpy(line,FileName("", ".SEGY", channel, fileNameOptions));
00282         if (noWrite==0)
00283         {
00284                 Rename(segyChannelFileName[channel],line);
00285         }
00286         PrintLog("OUTPUT: %s %d samples (%.2fs)\n",line,
00287                         infoHeader[channel].npts,
00288                         (float)infoHeader[channel].npts/(double)infoHeader[channel].sr.npts);
00289         unlink(segyChannelFileName[channel]);
00290 }
00291 
00294 void AbortSegyFile(int channel)
00295 {
00296         PrintDebug("abort for %s\n",segyChannelFileName[channel]);
00297         if (segyChannelFile[channel])
00298         {
00299                 fclose(segyChannelFile[channel]);
00300                 segyChannelFile[channel]=NULL;
00301                 unlink(segyChannelFileName[channel]);
00302         }
00303 }
00304 
00306 void AddSegyData(int channel)
00307 {
00308         if ((HWConfig.srObs[channel].npts==0)||(HWConfig.srExp[channel].npts==0))
00309         {
00310                 return;
00311         }
00312         if (OpenSegyChannel(channel))
00313         {
00314                 if (fwrite(HWConfig.data[channel],4,  
00315                             HWConfig.srObs[channel].npts,segyChannelFile[channel])
00316                                 !=HWConfig.srObs[channel].npts)
00317                 {
00318                         PrintError("unable to write to %s: %m\n",segyChannelFileName[channel]);
00319                         AbortSegyFile(channel);
00320                 }
00321                 infoHeader[channel].npts+=HWConfig.srObs[channel].npts;
00322         }
00323 }