source: branches/cosmos/src/libsrc/util/cosmos0putaway.c @ 7990

Revision 7990, 33.8 KB checked in by stefan, 3 months ago (diff)

Fixed extra second bug and read arc issue

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