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

Revision 7976, 35.9 KB checked in by stefan, 5 months ago (diff)

read arc file

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