source: trunk/src/libsrc/util/cosmos0putaway.c @ 7978

Revision 7978, 31.4 KB checked in by stefan, 9 months ago (diff)

gapwork

  • Property svn:executable set to *
RevLine 
[7683]1/*
2*
3*    Based on SUDS putaway and Geomag WDC putaway.
[7876]4*    For Erol Kalcan
[7683]5*
6*/
7
8/* cosmos0putaway.c
9
10Routines for writing trace data in COSMOS V0 format.
11
12Volume 0 files have raw data values in digital counts, obtained directly from the native binary files of
13the recording instrument. (Records obtained by digitizing analog film generally do not have a V0
14file.) V0 files may have almost no quality control or checking, and so they may be treated as internal
15files by networks. V0 files have adequate information in the headers to be converted to files with
16physical values. There is one component per file, so one record obtained by a triaxial recorder will
17result in three files.
18
19Header information is populated by a database of files. We don't need to know
20about it after we use the SCNL to find the right file and prepend it.
21
22Count data is extracted from a wave server:
23
24Data Section:
25• The first line of the data section includes the type of data (accel., velocity, displ., etc.) and its
26units, the number of data points, and the format to be used in reading the data. The format
27will type be for integer or real values, as appropriate; the format for real values can be
28floating point (e.g., 8f10.5) or exponential (e.g., 5e16.7). Although not required, 80-character
29line lengths are most convenient for many engineering strong motion data users.
30• The first line is followed by as many lines of data as are needed to accommodate the number
31of points, given the format in which the data values are written.
32
33In our case we'll alway use units=counts, and interger Format=(10I8):
34First line example followed by a data line.
3517629    raw accel.   pts, approx  176 secs, units=counts (50),Format=(10I8)
36-28596  -28611  -28630  -28617  -28609  -28550  -28543  -28654  -28698  -28661
37
38*/
39
[7699]40#include <errno.h>
[7683]41#include <limits.h>
42#include <stdio.h>
43#include <stdlib.h>
44#include <string.h>
45#include <time.h>
46
[7699]47#include "earthworm.h"
48#include "trace_buf.h"
49#include "swap.h"
50#include "ws_clientII.h"
51#include "cosmos0head.h"
52#include "cosmos0putaway.h"
53#include "pa_subs.h"
54#include "earthworm_simple_funcs.h"
55#include "chron3.h"
[7974]56#include "read_arc.h"
57 
[7699]58
[7683]59#define TAG_FILE        '.tag'        /* file containing the last known file tag */
[7978]60#define GAP_FILL        'NaN'
[7683]61#define MAXTXT           150
62
63FILE    *COSMOS0fp;                      /* file pointer for the COSMOS0 file */
[7974]64FILE    *COSMOS0Libfp;                   /* file pointer for the COSMOS0 library file */
65FILE    *EVENTARCfp;                     /* file pointer to the seismic event file */
[7683]66
67static  long   *COSMOS0Buffer;           /* write out COSMOS0 data as long integers */
68static  char    COSMOS0OutputFormat[MAXTXT];
69static  long    LineMean;            /* average of the 60 samples in the line */
70
[7978]71
72
[7699]73/* Internal Function Prototypes */
[7683]74static int StructMakeLocal(void *, int, char, int);
[7978]75char global_outputdir[80];
76sprintf(global_outputdir, "e:\earthworm\memphis\data\database_v0");
[7683]77/************************************************************************
78* Initialization function,                                              *
79*       This is the Put Away startup intializer. This is called when    *
80*       the system first comes up. Here is a chance to look around      *
81*       and see if it's possible to do business, and to complain        *
82*       if not ,BEFORE an event has to be processed.                    *
83*                                                                       *
[7978]84*       For COSMOS0, all we want to do is to make sure that the       *
[7683]85*       directory where files should be written exists.                 *
86*************************************************************************/
87int COSMOS0PA_init(int OutBufferLen, char *OutDir, char *OutputFormat,
88        int debug)
89{
90        /** Allocate COSMOS0Buffer and COSMOS0BufferShort
91        We waste RAM by allocating both the long and short buffers here
92        at the beginning of the code because some fluke (feature?) of NT
93        which I don't understand becomes unhappy if we do the allocation
94        later. Win2000, of course, behaves differently, and is quite happy
95        with buffer allocation after we have determined the format of the
96        incoming data */
97        if ((COSMOS0Buffer = (long *)malloc(OutBufferLen * sizeof(char))) == NULL)
98        {
99                logit("et", "COSMOS0PA_init: couldn't malloc COSMOS0Buffer\n");
100                return EW_FAILURE;
101        }
[7876]102
[7683]103        /* Make sure that the top level output directory exists */
104        if (RecursiveCreateDir(OutDir) != EW_SUCCESS)
105        {
106                logit("e", "COSMOS0PA_init: Call to RecursiveCreateDir failed\n");
107                return EW_FAILURE;
108        }
109
110        if (strlen(OutputFormat) >= sizeof(COSMOS0OutputFormat))
111        {
112                logit("", "COSMOS0PA_init: Error: OutputFormat(%s) is too long! Quitting!\n",
113                        OutputFormat);
114                return(EW_FAILURE);
115        }
116        else
117        {
118                strcpy(COSMOS0OutputFormat, OutputFormat);
119        }
120        return EW_SUCCESS;
121}
122
123/************************************************************************
124*   This is the Put Away event initializer. It's called when a snippet  *
125*   has been received, and is about to be processed.                    *
126*   It gets to see the pointer to the TraceRequest array,               *
127*   and the number of loaded trace structures.                          *
128*                                                                       *
129*   We need to make sure that the target directory                      *
130*   exists, create it if it does not, then open the COSMOS0 file        *
[7699]131*   for writing.                                                        *
132*                                                                       *
133*   This also opens the library file, and writes its contents to our    *
134*   output file                                                         *
[7683]135*************************************************************************/
136int COSMOS0PA_next_ev(TRACE_REQ *ptrReq,
[7899]137        char *OutDir, char *LibDir, char *EventDate, char *EventTime,
[7683]138        int cadence, int debug, double Multiplier)
139
140{
141        char    COSMOS0File[4 * MAXTXT];
[7678]142        char    COSMOS0LibFile[4 * MAXTXT];
[7974]143        char    EventArcFile[4 * MAXTXT];
[7876]144/*      char    c; */
[7683]145        char    year[5];
146        char    yr[3];
147        char    mo[3];
148        char    dy[3];
[7899]149        char    hr[3];
150        char    mn[3];
151        char    sec[7];
[7976]152        char    cos_date[35]; /* Start date of data requested */
[7899]153        char    record_id[35];
[7876]154        char    str[COSMOSLINELEN];
[7899]155        char    tempstr[COSMOSLINELEN];
[7876]156        int     LineNumber = 0;
[7974]157        time_t  rawtime;
158        struct tm * timeinfo;
159        size_t  nfbuffer;               /* Read bytes */
160        char    *fbuffer;               /* File read buffer */
161        static long      MaxMessageSize = 100000;  /* size (bytes) of largest msg */
162        HypoArc                 arcmsg;                 /* ARC message */
163        char            timestr[80];                                    /* Holds time messages */
164        time_t          ot;
165        struct tm * (*timefunc)(const time_t *);
166        char            time_type[30] = { 0 };                                  /* Time type UTC or local */
[7683]167
[7974]168
[7683]169        /* Grab the date-related substrings that we need for filenames. */
170        strncpy(year, EventDate, 4);
171        year[4] = '\0';
172        strncpy(yr, EventDate + 2, 2);
173        yr[2] = '\0';
174        strncpy(mo, EventDate + 4, 2);
175        mo[2] = '\0';
176        strncpy(dy, EventDate + 6, 2);
177        dy[2] = '\0';
[7899]178
179        strncpy(hr, EventTime, 2);
180        hr[2] = '\0';
181        strncpy(mn, EventTime + 2, 2);
182        mn[2] = '\0';
183        strncpy(sec, EventTime + 4, 5);
184        sec[5] = '\0';
185
[7876]186        cos_date[34] = '\0';
[7899]187        record_id[34] = '\0';   
188        tempstr[0]='\0';
[7683]189
190
[7899]191        sprintf(COSMOS0File, "%s/%s%s%s%s.v0", OutDir,
192                ptrReq->sta, yr, mo, dy);
[7683]193
[7899]194        /* cos_date "06/17/2018, 18:34:38.004 GMT (Q=5) "<-35 chars; Q stands for Quality, we have no way to know that */
[7976]195        sprintf(cos_date, "%s/%s/%s, %s:%s:%s  UTC", mo, dy, year, hr, mn, sec);
[7899]196
[7683]197        RecursiveCreateDir(OutDir);
[7899]198        sprintf(COSMOS0LibFile, "%s/%s_%s_%s_%s.dlv0", 
199                LibDir, ptrReq->net, ptrReq->sta, ptrReq->loc, ptrReq->chan);
[7683]200
[7678]201        if (debug == 1)
202                logit("t", "Opening COSMOS0 library file header %s\n", COSMOS0LibFile);
203
204        /* open library file just for reading */
205        if ((COSMOS0Libfp = fopen(COSMOS0LibFile, "r")) == NULL)
206        {
207                logit("e", "COSMOS0PA_next_ev: unable to open file %s: %s\n",
208                        COSMOS0LibFile, strerror(errno));
209                return EW_FAILURE;
210        }
[7683]211
212        if (debug == 1)
213                logit("et", "Opening COSMOS0 file %s\n", COSMOS0File);
214
215        /* open file */
216        if ((COSMOS0fp = fopen(COSMOS0File, "wt")) == NULL)
217        {
218                logit("e", "COSMOS0PA_next_ev: unable to open file %s: %s\n",
219                        COSMOS0File, strerror(errno));
220                return EW_FAILURE;
221        }
[7974]222        /* If there's enough time to implement, this filename should be able to be specified in the dot D file. */
223        sprintf(EventArcFile, "%s/event.arc", LibDir);
224        if (debug == 1)
225                logit("t", "Attempting to open event file which always must be named event.arc and located in the LibraryDirectory. \n LibraryDirectory is currently configured as: %s\n", LibDir);
[7683]226
[7974]227        /* open Event file just for reading */
228        if ((EVENTARCfp = fopen(EventArcFile, "r")) == NULL)
229        {
230                logit("e", "COSMOS0PA_next_ev: unable to open file %s: %s\nAn .arc file is required to set event information in the COSMOS file format.",
231                        EventArcFile, strerror(errno));
232                return EW_FAILURE;
233        }
[7899]234
[7974]235        /* Read file to buffer
236        *********************/
237        fbuffer = (char*)malloc(sizeof(char) * MaxMessageSize); /*check logic*/
238        if (fbuffer == NULL)
239        {
240                logit("et", "Unable to allocate memory for filter buffer\n");
241                return -1;
242        }
243
244        nfbuffer = fread(fbuffer, sizeof(char),
245                (size_t)MaxMessageSize, EVENTARCfp);
246        fclose(EVENTARCfp); // Done reading
247        if (nfbuffer == 0)
248        {
249                logit( "eto", "COSMOS0PA_next_ev: No data read from %s\n\nAn .arc file is required to set event information in the COSMOS file format.",
250                        EventArcFile);
251                return EW_FAILURE;
252        }
253        /* We should be able to do something like               origintime = arcmsg.sum.ot - GSEC1970; after parse_arc */
254        if (debug == 1) logit("et", "COSMOS0PA_next_ev: Debug: Arc parsing %s\n", EventArcFile);
255        if (parse_arc(fbuffer, &arcmsg) != 0 )
256        {
257                logit("et", "COSMOS0PA_next_ev: Error parsing %s\n",
258                        EventArcFile);
259        }
260        timefunc = localtime;
[7976]261        ot = (time_t)(arcmsg.sum.ot - GSEC1970); 
262        timeinfo = timefunc(&ot);
[7678]263        /* Copy the library file to the putaway file*/
[7876]264/*      c = fgetc(COSMOS0Libfp);
[7678]265        while (c != EOF)
266        {
267                fputc(c, COSMOS0fp);
268                c = fgetc(COSMOS0Libfp);
269        }
[7876]270*/
271        while (fgets(str, COSMOSLINELEN, COSMOS0Libfp) != NULL) {
272                LineNumber++ ;
273                switch (LineNumber) {
[7899]274                case 2: /*  2 1-40 Earthquake name (before a name has been assigned, may appear as “Record of”; test
275                                        records may use “Test Record of”, etc.).
276                                        41-80 Earthquake date and time (including time zone). */
[7976]277                        strftime(timestr, 80, "Record of                 Earthquake of %B %d, %Y, %H:%M:%S   UTC", timeinfo);
278                        fprintf(COSMOS0fp, "%s\n", timestr);
[7876]279                        break;
[7976]280                case 3: /*3 12-46 Hypocenter: 12-20 Latitude (positive North); 22-31 Longitude (positive East); and
281                                        35-37 Depth (km); 40-46 Source agency for hypocenter information.
282                                        47-80 Magnitude(s), including source agency for values (may be entered manually, so
283                                        spacing may be variable).**/
284                        sprintf(tempstr, "Hypocenter: %9.4f %9.4f H=%3.0fkm M=%3.1f, Md=%3.1f",
285                                arcmsg.sum.lat, arcmsg.sum.lon, arcmsg.sum.z, arcmsg.sum.Mpref, arcmsg.sum.Md);
286                        fprintf(COSMOS0fp, "%s\n", tempstr);
287                        break;
[7974]288                case 4: /*      4 9-42 Earthquake origin date and time, in UTC (i.e., GMT), with source agency; 43-80,
289                                        For use by preparing agency; may include note, e.g., “Preliminary Processing”, etc.*/
[7976]290
291                        /* Prepare origin time. Second can have tenths place, so leaving two spaces here to indicate
292                           we're not going to that level of precision, the timestruct doesn't have tenths. Probably worth revisiting.*/
293                        strftime(timestr, 80, "%m/%d/%Y, %H:%M:%S  ", timeinfo); 
[7974]294                        sprintf(tempstr, "Origin: %s UTC", timestr);
295                        fprintf(COSMOS0fp, "%s\n", tempstr);
296                        break;
[7899]297                case 8: /*  Record information:
298                                        8 17-50 Record start time (date and time of first sample; may be approximate - precise value
299                                        should be obtained from real header) and time quality (0 through 5 correspond to
300                                        lowest through highest quality).
301                                        59-80 Record identification, assigned by preparing agency.
302                                    "06/17/2018, 18:34:38.004 GMT (Q=5) "<- 35 chars   " (Q=5) " "38199368.SB.WLA.00.HN " <- 22char*/ 
303                        strncpy(record_id, str + 45, 34);
304                        fprintf(COSMOS0fp, "Rcrd start time:%s %s\n", cos_date, record_id);
305                        break;
306                case 10: /* 10 20-27 Length of raw record, as recorded (seconds); 45-54 Maximum of raw (uncorrected)
307                                        record in g (or other units);* 60-67 Time max occurs in raw record (secs after start).
308                                        Example:
309                                        Raw record length = 175.340  sec, Uncor max = 100000000, at   40.220 sec. */
310                        sprintf(tempstr, "Raw record length = %f", (ptrReq->reqEndtime - ptrReq->reqStarttime));
311                        strncpy(tempstr + 28, "sec, Uncor max =                                    ", 53);
312                        fprintf(COSMOS0fp, "%s\n", tempstr);
313                        break;
[7974]314                case 11: /*"Record Information; Line 11
315                                        11-40 Processing/release date, time and agency.
316                                        48-80 Maximum of data series in file, units, and time of maximum (seconds after start).*"
317                                        What should I put here? An example file looks like:
318                                        "Processed: 06/17/18 12:00 PDT UCB                                   "
319                                */
320                        time(&rawtime);
321                        timeinfo = gmtime(&rawtime);
[7976]322                        sprintf(tempstr, "Processed: %d/%d/%d %d:%d UTC                                       ", 
[7974]323                                (timeinfo->tm_mon + 1),
324                                timeinfo->tm_mday,
325                                (timeinfo->tm_year + 1900),
326                                timeinfo->tm_hour, 
327                                timeinfo->tm_min
328                        );
329                        fprintf(COSMOS0fp, "%s\n", tempstr);
330                        break;
[7976]331                case 47: /*AQMS normally writes this comment, but we're not using AQMS*/
332                        /* Example: | RcrdId: NC.72282711.NC.C031.HNE.01 */
333                        sprintf(tempstr, "| RcrdId: %s.72282711.%s.%s.%s.%s ", /* FIX: If anyone knows what I can put here instead of 72282711, I'll do it!*/
334                                ptrReq->net, ptrReq->net, ptrReq->sta, ptrReq->chan, ptrReq->loc);
335                        fprintf(COSMOS0fp, "%s\n", tempstr);
336                        break;
337                case 48: /*AQMS normally writes this comment, but we're not using AQMS*/
338                        /*Example:|<SCNL>C031.HNE.NC.01    <AUTH> 2015/03/01 16:46:25.000*/
339                        sprintf(tempstr, "|<SCNL> %s.%s.%s.%s    <AUTH> %d/%d/%d %d:%d:%d.000",
340                                ptrReq->sta, ptrReq->chan, ptrReq->net, ptrReq->loc,
341                                (timeinfo->tm_year + 1900),     (timeinfo->tm_mon + 1), timeinfo->tm_mday,
342                                timeinfo->tm_hour, timeinfo->tm_min, timeinfo->tm_sec);
343                        fprintf(COSMOS0fp, "%s", tempstr);
344                        break;
[7876]345                default: 
346                        fprintf(COSMOS0fp, "%s", str);         
347                }
348
349        }
[7683]350        fclose(COSMOS0Libfp);
351        return (EW_SUCCESS);
352}
353
354/************************************************************************
355* This is the working entry point into the disposal system. This        *
356* routine gets called for each trace snippet which has been recovered.  *
357* It gets to see the corresponding SNIPPET structure, and the event id  *
358*                                                                       *
[7978]359* For COSMOS0, this routine writes to the COSMOS0 file, pointed to by   *
[7699]360* the COSMOS0fp pointer, all of the received trace data in COSMOS0      *
361* format:                                                               *
[7683]362*                                                                       *
[7699]363*      1. COSMOS0 tag - indicating what follows                         *
364*      2. COSMOS0_STATIONCOMP struct - describe the station             *
365*      3. COSMOS0 tag - indicating what follows                         *
366*      4. COSMOS0_DESCRIPTRACE struct - describe the trace data         *
[7683]367*      5. trace data                                                    *
368*                                                                       *
369*  One bit of complexity is that we need to write the files in the      *
370*  correct byte-order. Based on the OutputFormat parameter, determine   *
[7699]371*  whether or not to swap bytes before writing the COSMOS0 file.        *
[7683]372*                                                                       *
373* WARNING: we clip trace data to -2147483648 - +2147483647 so it will   *
374*  fit in a long int. Any incoming data that is longer than 32 bits     *
375*  will be CLIPPED. cjb 5/18/2001                                       *
376*************************************************************************/
377/* Process one channel of data */
378int COSMOS0PA_next(TRACE_REQ *getThis, double GapThresh,
379        long OutBufferLen, int debug,
380        int Cadence, double Multiplier)
381{
382        TRACE2_HEADER *wf;
383        char    datatype;
384        char    day_line[122];
385        char    fourdigits[5];
[7978]386        char    hour_line[81] = ";                                                                                "; 
[7683]387        char    sample[12];
388        char    elevendigits[12];
389        char   *msg_p;        /* pointer into tracebuf data */
[7903]390        double  begintime = 0, starttime = 0, endtime = 0, currenttime = 0;
[7683]391        double  samprate, unrounded;
392        float  *f_data;
393        int     gap_count = 0;
394        int     i, j;
[7974]395        int     seconds = 0;
[7683]396        int     s_place = 0;
397        int     tabular_base;
[7899]398        int     total, raw_counts;
[7683]399        long    nfill_max = 0l;
400        long    nsamp, nfill;
401        long    nsamp_this_scn = 0l;
402        long    this_size;
403        long   *l_data;
404        short  *s_data;
405        int    current_int;
[7974]406        char    tempstr[COSMOSLINELEN];
407        char    tempstr2[COSMOSLINELEN];
[7683]408
409        /* Put this in the .d file once we know we want it. */
410        /*  double multiplier = 0.001; */
411
412        /* Check arguments */
413        if (getThis == NULL)
414        {
415                logit("e", "COSMOS0PA_next: invalid argument passed in.\n");
416                return EW_FAILURE;
417        }
418        /* Used for computing trace statistics */
419        total = 0;
420
421        if ((msg_p = getThis->pBuf) == NULL)   /* pointer to first message */
422        {
423                logit("e", "COSMOS0PA_next: Message buffer is NULL.\n");
424                return EW_FAILURE;
425        }
426        wf = (TRACE2_HEADER *)msg_p;
427
428        /* Look at the first TRACE2_HEADER and get set up of action */
429        if (WaveMsg2MakeLocal(wf) < 0)
430        {
431                logit("e", "COSMOS0PA_next: unknown trace data type: %s\n",
432                        wf->datatype);
433                return(EW_FAILURE);
434        }
435
436        nsamp = wf->nsamp;
437        starttime = wf->starttime;
438        endtime = wf->endtime;
439        samprate = wf->samprate;
440        if (samprate < 0.0001)
441        {
442                logit("et", "unreasonable samplerate (%f) for <%s.%s.%s.%s>\n",
443                        samprate, wf->sta, wf->chan, wf->net, wf->loc);
444                return(EW_FAILURE);
445        }
[7899]446        /* LAST header line now that we know sample rate*/
447        /* WRITE Dynamic header; this probably can't be here though because we need to calculate these results */
448        /* First we'll write a line that looks like this: */
449        /* 17770    raw accel.pts, approx  178 secs, units = counts(50), Format = (7I11) */
450        /*Line Cols
451        1st 1 - 8 Number of data points following; 10 - 21 Physical parameter of the data.
452        35 - 38 Approximate length of record(rounded to nearest sec; see Rhdr(66) for precise value).
453        52 - 58 Units of the data values(e.g., cm / sec2); 60 - 61 Index code for units
[7974]454        71-80 Format of the data values written on the following lines.
[7899]455        */
456        /*sprintf(hour_line, "17770    raw accel.pts, approx  178 secs, units=counts (50),Format=(10I8)\r\n"); /* we want leading spaces */
457
[7974]458/*      sprintf(tempstr, "17770    raw accel.pts,   approx  178  secs, units=counts (50),Format=(10I8)   \n");*/
459        sprintf(tempstr, "         raw accel.pts,   approx           , units=counts (50),Format=(10I8)    \n");
460        raw_counts = (getThis->reqEndtime - getThis->reqStarttime ) * samprate;
461        sprintf(tempstr2, "%d", raw_counts);
462        strncpy(tempstr, tempstr2, strlen(tempstr2));
463        seconds = (getThis->reqEndtime - getThis->reqStarttime);
464        sprintf(tempstr2, "%d secs", seconds);
465        strncpy(tempstr + 34, tempstr2, strlen(tempstr2));
466        if (fwrite(tempstr, 81, 1, COSMOS0fp) != 1)
467        {
468                logit("et", "COSMOS0PA_next: error writing COSMOS0 dynamic header line. \n");
469                return EW_FAILURE;
470        }
[7683]471        begintime = starttime;
472        /* datatype i4 = intel byte order 4 bytes, s2 = sparc byte order; 2 bytes */
473        datatype = 'n';
474        if (wf->datatype[0] == 's' || wf->datatype[0] == 'i')
475        {
476                if (wf->datatype[1] == '2') datatype = 's';
477                else if (wf->datatype[1] == '4') datatype = 'l';
478        }
479        else if (wf->datatype[0] == 't' || wf->datatype[0] == 'f')
480        {
481                if (wf->datatype[1] == '4') datatype = 'f';
482        }
483        if (datatype == 'n')
484        {
[7699]485                logit( "et", "COSMOS0PA_next: unsupported datatype: %s\n", wf->datatype );
[7683]486                return(EW_FAILURE);
487        }
488
489        if (debug == 1)
490                logit("et", "COSMOS0PA_next: working on <%s/%s/%s/%s> datatype: %c \n",
491                        wf->sta, wf->chan, wf->net, wf->loc, datatype);
492
493        /********************** loop through all the messages for this s-c-n **********************/
494        while (1)
495        {
496                /* advance message pointer to the data */
497                msg_p += sizeof(TRACE2_HEADER);
498
499                /* check for sufficient memory in output buffer */
500                this_size = (nsamp_this_scn + nsamp) * sizeof(long);
501                if (OutBufferLen < this_size)
502                {
503                        logit("e", "out of space for <%s.%s.%s.%s>; saving long trace.\n",
504                                wf->sta, wf->chan, wf->net, wf->loc);
505                        break;
506                }
507
508                /* if data are floats, clip to longs cjb 5/18/2001 */
509                switch (datatype)
510                {
511                case 's':
512                        s_data = (short *)msg_p;
513                        for (j = 0; j < nsamp; j++, nsamp_this_scn++)
514                                COSMOS0Buffer[nsamp_this_scn] = (long)s_data[j];
515                        msg_p += sizeof(short) * nsamp;
516                        break;
517                case 'l':
518                        l_data = (long *)msg_p;
519                        for (j = 0; j < nsamp; j++, nsamp_this_scn++)
520                                COSMOS0Buffer[nsamp_this_scn] = l_data[j];
521                        msg_p += sizeof(long) * nsamp;
522                        break;
523                case 'f':
524                        f_data = (float *)msg_p;
525                        /* CLIP the data to long int */
526                        for (j = 0; j < nsamp; j++, nsamp_this_scn++)
527                        {
528                                if (l_data[j] < (float)LONG_MIN)
529                                        COSMOS0Buffer[nsamp_this_scn] = LONG_MIN;
530                                else if (l_data[j] > (float)LONG_MAX)
531                                        COSMOS0Buffer[nsamp_this_scn] = LONG_MAX;
532                                else
533                                        COSMOS0Buffer[nsamp_this_scn] = (long)l_data[j];
534                        }
535                        msg_p += sizeof(float) * nsamp;
536                        break;
537                }
538
539                /* msg_p has been advanced to the next TRACE_BUF; localize bytes *
[7699]540                * and check for gaps.                                            */
[7683]541                wf = (TRACE2_HEADER *)msg_p;
542                if (WaveMsg2MakeLocal(wf) < 0)
543                {
[7876]544                        if (debug == 1)
545                                logit("e", "COSMOS0PA_next: unknown trace data type or unexpected end of data: %s\n",
546                                        wf->datatype);
547                        else
548                                logit("e", "COSMOS0PA_next: unknown trace data type or unexpected end of data.\n");
[7724]549                        break;
550                        //return(EW_FAILURE);
[7683]551                }
552                nsamp = wf->nsamp;
553                starttime = wf->starttime;
554                /* starttime is set for new packet; endtime is still set for old packet */
555                if (endtime + (1.0 / samprate) * GapThresh < starttime)
556                {
557                        /* there's a gap, so fill it */
558                        if (debug == 1)
559                                logit("e", "gap in %s.%s.%s.%s: %lf: %lf\n", wf->sta, wf->chan, wf->net,
560                                        wf->loc, endtime, starttime - endtime);
561                        nfill = (long)(samprate * (starttime - endtime) - 1);
562                        if ((nsamp_this_scn + nfill) * (long)sizeof(long) > OutBufferLen)
563                        {
[7699]564                                logit( "e", "bogus gap (%ld); skipping\n", nfill );
[7683]565                                return(EW_FAILURE);
566                        }
567                        /* do the filling */
568                        for (j = 0; j < nfill; j++, nsamp_this_scn++)
569                                COSMOS0Buffer[nsamp_this_scn] = getThis->fill; // changed from local variable fill swl
570                                                                                                                   /* keep track of how many gaps and the largest one */
571                        gap_count++;
572                        if (nfill_max < nfill)
573                                nfill_max = nfill;
574                }
575                /* Advance endtime to the new packet;        *
576                * process this packet in the next iteration */
577                endtime = wf->endtime;
578        } /******************** end while(1) ***************************************************/
579
580          /* If the sample rate is 1 sample per minute then we'll have a sample rate of .016 */
581          /* For minute data we want 24 rows of 60 samples each, 1440 samples in a day. */
582          /* A single file fills the month of October for SJG, and includes four trace types */
583          /* FYXZ, so there are a total of 2976 lines in this file. */
584
585          /* Match our metadata with our waveserver request, if possible.
586          * observatory's value is -1 if there isn't a match. */
587
588        currenttime = begintime;
589        j = 0;
590
591
[7978]592        while ((j < nsamp_this_scn) && (currenttime < getThis->reqEndtime))
593        {
594                /* Only give them what they asked for, not each sample we got back.
595                Tracebufs contain multiple samples, and we may need to request
596                an earlier one to get the start sample we need, or a later one
597                for the last sample*/
598                while (currenttime < getThis->reqStarttime) {
599                        currenttime = currenttime + 1 / samprate;
[7683]600                        j++;
601                }
602
[7978]603                s_place = 0;
[7683]604
[7978]605                /* 35-394 60I6   60 6-digit 1-minute values for the given element for that data hour.
606                *               The values are in tenth-minutes for D and I, and in
607                *               nanoTeslas for the intensity elements.
608                */
609                total = 0;
610                LineMean = 0;
[7683]611
[7978]612                i = 0;
[7683]613
[7978]614                /* WRITE DATA */
[7683]615
[7978]616                while (i < 7 && j < nsamp_this_scn)
617                {
618                        if (COSMOS0Buffer[j] != getThis->fill) {
619                                sprintf(elevendigits, "           "); /* we want leading spaces */
620                                if (((int)(COSMOS0Buffer[j] * Multiplier) > 999999) || ((int)(COSMOS0Buffer[j] * Multiplier) < -99999)) {
621                                        sprintf(sample, "9999");
622                                        /* prevent out of range string */
[7683]623                                }
624                                else {
[7978]625                                        sprintf(sample, "%d", (int)(COSMOS0Buffer[j] * Multiplier));
[7683]626                                }
[7978]627                                strcpy(elevendigits + 11 - strlen(sample), sample);
628                                strcpy(hour_line + s_place, elevendigits);
[7683]629                        }
[7978]630                        else {
631                                /* We have a gap, this is where gap data is written */
632                                //strcpy(hour_line + s_place, "  9999");
633                                sprintf(sample, " 0");
634                                strcpy(elevendigits + 11 - strlen(sample), sample);
635                                strcpy(hour_line + s_place, elevendigits);                             
[7683]636                        }
[7978]637                        s_place = s_place + 11;
[7683]638
[7978]639                        total += (int)(COSMOS0Buffer[j] * Multiplier);
[7683]640
[7978]641                        j++; i++;
642                }
[7683]643
[7978]644
645                /* 401-402       Record end marker.
646                *               Two chars 'cr'= 13 and 'nl'= 10.
647                */
648                hour_line[77] = ' '; /*Replace that null that sprintf got us*/
649                hour_line[80] = '\n';
650
651                /* Write out line */
652                if (fwrite(&hour_line, sizeof(hour_line), 1, COSMOS0fp) != 1)
653                {
654                        logit("et", "COSMOS0PA_next: error writing COSMOS0 line. \n");
655                        return EW_FAILURE;
656                }
657        }
[7683]658        return EW_SUCCESS;
659}
660
[7978]661
662
[7683]663/************************************************************************
664* This is the Put Away end event routine. It's called after we've       *
665* finished processing one event                                         *
666*                                                                       *
[7699]667* For PC-COSMOS0 - close the COSMOS0 file, pointed to by COSMOS0fp      *
668*               free COSMOS0Buffer memory to be nice to everyone else   *
[7683]669*************************************************************************/
670int COSMOS0PA_end_ev(int debug)
671{
672        fclose(COSMOS0fp);
673
674        if (debug == 1)
675                logit("t", "Closing COSMOS0 file \n");
676
677        return(EW_SUCCESS);
678}
679
680/************************************************************************
681*       This is the Put Away close routine. It's called after when      *
682*       we're being shut down.                                          *
683*************************************************************************/
684int COSMOS0PA_close(int debug)
685{
686
687        free((char *)COSMOS0Buffer);
688        return(EW_SUCCESS);
689}
690
691/************************************************************************
[7699]692* We first try and find a tabular base that will apply to all data      *
693* that has been returned from our waveserver request. If there is       *
694* no base that will fit all the data and be within -999 to 9998, we     *
695* will keep removing low and high values 'till we get one. This base    *
696* will work for most days, but not spike days. Spike days will          *
697* recalculate their own tabular base                                    *
[7683]698*************************************************************************/
699int FindTabularBase(int number_samples, long *SampleBuffer, TRACE_REQ *getThis, double Multiplier, int skip)
700{
701        int daymax;
702        int daymin;
703        int tabular_base = 999;
704        int range;
705        int j;
706        int sample_COSMOS0_int;
707
708        daymax = (int)getThis->fill;
709        daymin = (int)getThis->fill;
710
711        for (j = skip; j<number_samples; j++) {
712                /* From the definition above: For the intensity elements we have that
713                * hourly value (nanoTeslas) = tab.base*100 + tab.value
714                * The Multiplier is specific to the USGS waveservers so they can get
715                * decimal values out of the integer type...                             */
716                /* Won't divide by 100 'till we're done */
717
718                /* Skip ahead past all gaps */
719                if (SampleBuffer[j] == getThis->fill) {
720                        continue;
721                }
722                sample_COSMOS0_int = (int)((SampleBuffer[j] * Multiplier));
723                /* Keep track of the line's max and min for the tabular base generation */
724                if ((sample_COSMOS0_int > daymax) || (daymax == getThis->fill)) {
725                        daymax = sample_COSMOS0_int;
726                }
727                if ((sample_COSMOS0_int < daymin) || (daymin == getThis->fill)) {
728                        daymin = sample_COSMOS0_int;
729                }
730                range = daymax - daymin;
731        }
732        /* Really we want our values ideally to be less than 1000 and > than 1.
733        * second most ideally less than 1000 and > than -999.
734        * after that, not ideal, anything that fits between 9998 and -999
735        * after that, grab the median and make that 5500;
736        * everything out of range becomes a spike = 9999 */
737        if (range < 899) {
738                /* let's have the lowest number in our range = 0 *
739                * so we will subtract that number from all numbers. */
740                tabular_base = daymin / 100;
741        }
742        else if (range < 1899) {
743                /* let's have the highest number in our range = 999
744                * so we have as few negative numbers as possible */
745                tabular_base = (daymax - 900) / 100;
746        }
747        else if (range < (9898 + 999 + 1)) { /* the 1 is to account for 0; 999 is for to -999*/
748                                                                                 /* let's have the lowest number in our range = -999
749                                                                                 * so our numbers are as low as possible, and are more
750                                                                                 * likely to be less than 3 digits . */
751                tabular_base = (daymin + 999) / 100;
752        }
753        else if (number_samples < 25) {
754                /* out of range: if we're out of range then let's make the
755                * median value equal to 5500. We'll have to drop (ie: make 9999)
756                * spike numbers on this line that are > 9998 and < -999  */
757                tabular_base = (Median(number_samples, SampleBuffer) + 5500) / 100;
758        }
759        else {
760                /* If we're finding the tabular base for a larger range than a day
761                * and we can't find one, don't use the median, just return failure
762                * We might be doing this we we want to find a common base for an
763                * entire COSMOS0 file for a particular element. So if this fails,
764                * First line will determine it's own tabular base. Each subsequent
765                * line will try to use the previous line's base. If that fails the
766                * subsequent line will calculate a new base. */
767                tabular_base = 9999;
768        }
769        return(tabular_base);
770}
771
772/*
773*
774*  Byte swapping functions
775*/
776
777/* swap bytes if necessary, depending on what machine */
778/* we have been compiled, and what machine the data   */
779/* is being read from/being written to                */
780
781static int StructMakeLocal(void *ptr, int struct_type, char data_type,
782        int debug)
783{
784        if (ptr == NULL)
785        {
786                logit("e", "StructMakeLocal: NULL pointer passed in\n");
787                return EW_FAILURE;
788        }
789
790#if defined (_INTEL)
791        if (data_type != '6')
792        {
793                if (debug == 1)
794                        logit("et", "Hoping we don't really have to do any swapping from Intel to Sun \n");
795
796        }
797
798#elif defined (_SPARC)
799        if (data_type == '6')
800        {
801                if (debug == 1)
802                        logit("et", "Hoping we don't really have to do any swapping from Sun to Intel because we've deleted from COSMOS0 the SwapDo function that suds used. \n");
803        }
804
805#endif
806        return EW_SUCCESS;
807}
808
809long Median(int number_of_array_elements, long *Array)
810{
811        long *LocalArray;
812
813        LocalArray = (long*)malloc(number_of_array_elements * sizeof(long));
814
815        qsort(&LocalArray[0], number_of_array_elements, sizeof(long), longCompareForQsort);
816        /* Get Median */
817        return(LocalArray[(long)(number_of_array_elements / 2)]);
818}
819/*****************************************************************/
820/* Just a simple compare function so that Q sort does it's thing */
821/*****************************************************************/
822int longCompareForQsort(const void *x, const void *y)
823{
824        /* compare must have const void ptr params for qsort to work */
825        const long   *ix = (const long *)x;
826        const long   *iy = (const long *)y;
827        /* returns 1, -1 or 0 */
828        return (*ix > *iy) - (*ix < *iy);
[7784]829}
Note: See TracBrowser for help on using the repository browser.