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

Revision 7991, 34.0 KB checked in by stefan, 6 weeks ago (diff)

For PRISM changed line 8 to end in RcrdId?: (see comment)

  • 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); <- prism doesn't like this*/
405                        /* sprintf(tempstr2,    "RcrdId:%s.%lu.%s.%s.%s.%s     ",
406                                ptrReq->net, arcmsg.sum.qid, ptrReq->net, ptrReq->sta, ptrReq->chan, ptrReq->loc); <- prism likes this but it might go beyond legal 80 chars/line */
407                        sprintf(tempstr2, "RcrdId: (see comment)                   ");
408
409                        strncpy(tempstr + 51, tempstr2, 29);
410                        tempstr[80] = '\0'; /* Chop off any extra junk */
411                        fprintf(COSMOS0fp, "%s\n", tempstr);
412                        break;
413                case 10: /* 10 20-27 Length of raw record, as recorded (seconds); 45-54 Maximum of raw (uncorrected)
414                                 record in g (or other units);* 60-67 Time max occurs in raw record (secs after start).
415                                 Example:
416                                 Raw record length = 175.340  sec, Uncor max = 100000000, at   40.220 sec. */
417                        sprintf(tempstr, "Raw record length = %f", (ptrReq->reqEndtime - ptrReq->reqStarttime));
418                        strncpy(tempstr + 28, "sec, Uncor max =                                    ", 53);
419                        fprintf(COSMOS0fp, "%s\n", tempstr);
420                        break;
421                case 11: /*"Record Information; Line 11
422                                 11-40 Processing/release date, time and agency.
423                                 48-80 Maximum of data series in file, units, and time of maximum (seconds after start).*"
424                                 What should I put here? An example file looks like:
425                                 "Processed: 06/17/18 12:00 PDT UCB                                   "
426                                 */
427                        time(&rawtime);
428                        timeinfo = gmtime(&rawtime);
429                        sprintf(tempstr, "Processed: %d/%d/%d %d:%01d UTC                                       ",
430                                (timeinfo->tm_mon + 1),
431                                timeinfo->tm_mday,
432                                (timeinfo->tm_year + 1900),
433                                timeinfo->tm_hour,
434                                timeinfo->tm_min
435                        );
436                        fprintf(COSMOS0fp, "%s\n", tempstr);
437                        break;
438                /* Real Headers include Earthquake information similar to comments above, comments for humans, below for PRISM*/
439                case 27: /* Site geology stuff from the database are the first 4 params in this line, so read them in and leave
440                                    them be. But the 5th param we'll overwrite here with Earthquake latitude (decimal degrees, North positive)*/
441                        sprintf(fifteen, "%15.6f", arcmsg.sum.lat);
442                        strncpy(str + 60, fifteen, 15);
443                        fprintf(COSMOS0fp, "%s", str);
444                        break;
445                case 28: /* Earthquake longitude (decimal degrees, East positive), Earthquake depth (km)., Moment magnitude, M., Surface-wave magnitude, MS., Local magnitude, ML.*/
446                        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);
447                        fprintf(COSMOS0fp, "%s\n", tempstr);
448                        break;
449                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 */
450                        sprintf(fifteen, "%15.6f", arcmsg.sum.Md);
451                        strncpy(str, fifteen, 15);                     
452                        sprintf(fifteen, "%15.6f", -999.0);
453                        strncpy(str + 15, fifteen, 15);
454                        strncpy(str + 30, fifteen, 15);
455                        fprintf(COSMOS0fp, "%s", str);
456                        break;
457                case 47: /*AQMS normally writes this comment, but we're not using AQMS*/
458                                 /* Example: | RcrdId: NC.72282711.NC.C031.HNE.01 */
459                        sprintf(tempstr, "| RcrdId: %s.%lu.%s.%s.%s.%s ", 
460                                ptrReq->net, arcmsg.sum.qid, ptrReq->net, ptrReq->sta, ptrReq->chan, ptrReq->loc);
461                        fprintf(COSMOS0fp, "%s\n", tempstr);
462                        break;
463                case 48: /*AQMS normally writes this comment, but we're not using AQMS*/
464                                 /*Example:|<SCNL>C031.HNE.NC.01    <AUTH> 2015/03/01 16:46:25.000*/
465                        sprintf(tempstr, "|<SCNL> %s.%s.%s.%s    <AUTH> %d/%d/%d %d:%d:%d.000",
466                                ptrReq->sta, ptrReq->chan, ptrReq->net, ptrReq->loc,
467                                (timeinfo->tm_year + 1900), (timeinfo->tm_mon + 1), timeinfo->tm_mday,
468                                timeinfo->tm_hour, timeinfo->tm_min, timeinfo->tm_sec);
469                        fprintf(COSMOS0fp, "%s\n", tempstr);
470                        break;
471                default:
472                        fprintf(COSMOS0fp, "%s", str);
473                }
474
475        }
476        fclose(COSMOS0Libfp);
477        return (EW_SUCCESS);
478}
479
480/************************************************************************
481* This is the working entry point into the disposal system. This        *
482* routine gets called for each trace snippet which has been recovered.  *
483* It gets to see the corresponding SNIPPET structure, and the event id  *
484*                                                                       *
485* For COSMOS0, this routine writes to the COSMOS0 file, pointed to by   *
486* the COSMOS0fp pointer, all of the received trace data in COSMOS0      *
487* format:                                                               *
488*                                                                       *
489*      1. COSMOS0 tag - indicating what follows                         *
490*      2. COSMOS0_STATIONCOMP struct - describe the station             *
491*      3. COSMOS0 tag - indicating what follows                         *
492*      4. COSMOS0_DESCRIPTRACE struct - describe the trace data         *
493*      5. trace data                                                    *
494*                                                                       *
495*  One bit of complexity is that we need to write the files in the      *
496*  correct byte-order. Based on the OutputFormat parameter, determine   *
497*  whether or not to swap bytes before writing the COSMOS0 file.        *
498*                                                                       *
499* WARNING: we clip trace data to -2147483648 - +2147483647 so it will   *
500*  fit in a long int. Any incoming data that is longer than 32 bits     *
501*  will be CLIPPED. cjb 5/18/2001                                       *
502*************************************************************************/
503/* Process one channel of data */
504int COSMOS0PA_next(TRACE_REQ *getThis, double GapThresh,
505        long OutBufferLen, int debug,
506        int Cadence, double Multiplier)
507{
508        TRACE2_HEADER *wf;
509        char    datatype;
510        char    day_line[122];
511        char    fourdigits[5];
512        char    hour_line[81] = ";                                                                                "; 
513        char    sample[12];
514        char    elevendigits[12];
515        char    eightdigits[9];
516        char   *msg_p;        /* pointer into tracebuf data */
517        double  begintime = 0, starttime = 0, endtime = 0, currenttime = 0;
518        double  samprate, unrounded;
519        float  *f_data;
520        int     gap_count = 0;
521        int     i, j;
522        int     seconds = 0;
523        int     s_place = 0;
524        int     tabular_base;
525        int     total, raw_counts;
526        long    nfill_max = 0l;
527        long    nsamp, nfill;
528        long    nsamp_this_scn = 0l;
529        long    this_size;
530        long   *l_data;
531        short  *s_data;
532        int    current_int;
533        char    tempstr[COSMOSLINELEN];
534        char    tempstr2[COSMOSLINELEN];
535
536        /* Put this in the .d file once we know we want it. */
537        /*  double multiplier = 0.001; */
538
539        /* Check arguments */
540        if (getThis == NULL)
541        {
542                logit("e", "COSMOS0PA_next: invalid argument passed in.\n");
543                return EW_FAILURE;
544        }
545
546
547        COSMOS0PA_header(getThis, cosmos_info->OutDir, cosmos_info->LibDir, cosmos_info->EventDate, cosmos_info->EventTime, Cadence, debug, Multiplier);
548        /* Used for computing trace statistics */
549        total = 0;
550
551        if ((msg_p = getThis->pBuf) == NULL)   /* pointer to first message */
552        {
553                logit("e", "COSMOS0PA_next: Message buffer is NULL.\n");
554                return EW_FAILURE;
555        }
556        wf = (TRACE2_HEADER *)msg_p;
557
558        /* Look at the first TRACE2_HEADER and get set up of action */
559        if (WaveMsg2MakeLocal(wf) < 0)
560        {
561                logit("e", "COSMOS0PA_next: unknown trace data type: %s\n",
562                        wf->datatype);
563                return(EW_FAILURE);
564        }
565
566        nsamp = wf->nsamp;
567        starttime = wf->starttime;
568        endtime = wf->endtime;
569        samprate = wf->samprate;
570        if (samprate < 0.0001)
571        {
572                logit("et", "unreasonable samplerate (%f) for <%s.%s.%s.%s>\n",
573                        samprate, wf->sta, wf->chan, wf->net, wf->loc);
574                return(EW_FAILURE);
575        }
576        /* LAST header line now that we know sample rate*/
577        /* WRITE Dynamic header; this probably can't be here though because we need to calculate these results */
578        /* First we'll write a line that looks like this: */
579        /* 17770    raw accel.pts, approx  178 secs, units = counts(50), Format = (7I11) */
580        /*Line Cols
581        1st 1 - 8 Number of data points following; 10 - 21 Physical parameter of the data.
582        35 - 38 Approximate length of record(rounded to nearest sec; see Rhdr(66) for precise value).
583        52 - 58 Units of the data values(e.g., cm / sec2); 60 - 61 Index code for units
584        71-80 Format of the data values written on the following lines.
585        */
586        /*sprintf(hour_line, "17770    raw accel.pts, approx  178 secs, units=counts (50),Format=(10I8)\r\n"); /* we want leading spaces */
587
588/*      sprintf(tempstr, "17770    raw accel.pts,   approx  178  secs, units=counts (50),Format=(10I8)   \n");*/
589        sprintf(tempstr, "         raw accel.pts,   approx           , units=counts (50),Format=(10I8)    \n");
590        raw_counts = (getThis->reqEndtime - getThis->reqStarttime ) * samprate; 
591        sprintf(tempstr2, "%d", raw_counts);
592        strncpy(tempstr, tempstr2, strlen(tempstr2));
593        seconds = (getThis->reqEndtime - getThis->reqStarttime);
594        sprintf(tempstr2, "%d secs", seconds);
595        strncpy(tempstr + 34, tempstr2, strlen(tempstr2));
596        if (!COSMOS0fp) {
597                logit("et", "Tried to write to a null file pointer, exiting. Does your library match your waveserver request?\n");
598                exit(0);
599        }
600        if (fwrite(tempstr, 81, 1, COSMOS0fp) != 1)
601        {
602                logit("et", "COSMOS0PA_next: error writing COSMOS0 dynamic header line. \n");
603                return EW_FAILURE;
604        }
605        begintime = starttime;
606        /* datatype i4 = intel byte order 4 bytes, s2 = sparc byte order; 2 bytes */
607        datatype = 'n';
608        if (wf->datatype[0] == 's' || wf->datatype[0] == 'i')
609        {
610                if (wf->datatype[1] == '2') datatype = 's';
611                else if (wf->datatype[1] == '4') datatype = 'l';
612        }
613        else if (wf->datatype[0] == 't' || wf->datatype[0] == 'f')
614        {
615                if (wf->datatype[1] == '4') datatype = 'f';
616        }
617        if (datatype == 'n')
618        {
619                logit( "et", "COSMOS0PA_next: unsupported datatype: %s\n", wf->datatype );
620                return(EW_FAILURE);
621        }
622
623        if (debug == 1)
624                logit("et", "COSMOS0PA_next: working on <%s/%s/%s/%s> datatype: %c \n",
625                        wf->sta, wf->chan, wf->net, wf->loc, datatype);
626
627        /********************** loop through all the messages for this s-c-n **********************/
628        while (1)
629        {
630                /* advance message pointer to the data */
631                msg_p += sizeof(TRACE2_HEADER);
632
633                /* check for sufficient memory in output buffer */
634                this_size = (nsamp_this_scn + nsamp) * sizeof(long);
635                if (OutBufferLen < this_size)
636                {
637                        logit("e", "out of space for <%s.%s.%s.%s>; saving long trace.\n",
638                                wf->sta, wf->chan, wf->net, wf->loc);
639                        break;
640                }
641
642                /* if data are floats, clip to longs cjb 5/18/2001 */
643                switch (datatype)
644                {
645                case 's':
646                        s_data = (short *)msg_p;
647                        for (j = 0; j < nsamp; j++, nsamp_this_scn++)
648                                COSMOS0Buffer[nsamp_this_scn] = (long)s_data[j];
649                        msg_p += sizeof(short) * nsamp;
650                        break;
651                case 'l':
652                        l_data = (long *)msg_p;
653                        for (j = 0; j < nsamp; j++, nsamp_this_scn++)
654                                COSMOS0Buffer[nsamp_this_scn] = l_data[j];
655                        msg_p += sizeof(long) * nsamp;
656                        break;
657                case 'f':
658                        f_data = (float *)msg_p;
659                        /* CLIP the data to long int */
660                        for (j = 0; j < nsamp; j++, nsamp_this_scn++)
661                        {
662                                if (l_data[j] < (float)LONG_MIN)
663                                        COSMOS0Buffer[nsamp_this_scn] = LONG_MIN;
664                                else if (l_data[j] > (float)LONG_MAX)
665                                        COSMOS0Buffer[nsamp_this_scn] = LONG_MAX;
666                                else
667                                        COSMOS0Buffer[nsamp_this_scn] = (long)l_data[j];
668                        }
669                        msg_p += sizeof(float) * nsamp;
670                        break;
671                }
672
673                /* msg_p has been advanced to the next TRACE_BUF; localize bytes *
674                * and check for gaps.                                            */
675                wf = (TRACE2_HEADER *)msg_p;
676                if (WaveMsg2MakeLocal(wf) < 0)
677                {
678                        if (debug == 1)
679                                logit("e", "COSMOS0PA_next: unknown trace data type or unexpected end of data: %s\n",
680                                        wf->datatype);
681                        else
682                                logit("e", "COSMOS0PA_next: unknown trace data type or unexpected end of data.\n");
683                        break;
684                        //return(EW_FAILURE);
685                }
686                nsamp = wf->nsamp;
687                starttime = wf->starttime;
688                /* starttime is set for new packet; endtime is still set for old packet */
689                if (endtime + (1.0 / samprate) * GapThresh < starttime)
690                {
691                        /* there's a gap, so fill it */
692                        if (debug == 1)
693                                logit("e", "gap in %s.%s.%s.%s: %lf: %lf\n", wf->sta, wf->chan, wf->net,
694                                        wf->loc, endtime, starttime - endtime);
695                        nfill = (long)(samprate * (starttime - endtime) - 1);
696                        if ((nsamp_this_scn + nfill) * (long)sizeof(long) > OutBufferLen)
697                        {
698                                logit( "e", "bogus gap (%ld); skipping\n", nfill );
699                                return(EW_FAILURE);
700                        }
701                        /* do the filling */
702                        for (j = 0; j < nfill; j++, nsamp_this_scn++)
703                                COSMOS0Buffer[nsamp_this_scn] = getThis->fill; // changed from local variable fill swl
704                                                                                                                   /* keep track of how many gaps and the largest one */
705                        gap_count++;
706                        if (nfill_max < nfill)
707                                nfill_max = nfill;
708                }
709                /* Advance endtime to the new packet;        *
710                * process this packet in the next iteration */
711                endtime = wf->endtime;
712        } /******************** end while(1) ***************************************************/
713
714          /* If the sample rate is 1 sample per minute then we'll have a sample rate of .016 */
715          /* For minute data we want 24 rows of 60 samples each, 1440 samples in a day. */
716          /* A single file fills the month of October for SJG, and includes four trace types */
717          /* FYXZ, so there are a total of 2976 lines in this file. */
718
719          /* Match our metadata with our waveserver request, if possible.
720          * observatory's value is -1 if there isn't a match. */
721
722        currenttime = begintime;
723        j = 0;
724
725
726        while ((j < nsamp_this_scn) && (currenttime < getThis->reqEndtime) && (j < raw_counts))
727        {
728                /* Only give them what they asked for, not each sample we got back.
729                Tracebufs contain multiple samples, and we may need to request
730                an earlier one to get the start sample we need, or a later one
731                for the last sample*/
732                while (currenttime < getThis->reqStarttime) {
733                        currenttime = currenttime + 1 / samprate;
734                        j++;
735                }
736
737                s_place = 0;
738
739                /* 35-394 60I6   60 6-digit 1-minute values for the given element for that data hour.
740                *               The values are in tenth-minutes for D and I, and in
741                *               nanoTeslas for the intensity elements.
742                */
743                total = 0;
744                LineMean = 0;
745
746                i = 0;
747
748                /* WRITE DATA */
749                /* Even if there are more samples in this scn, we shouldn't have more than the raw counts the user asked for*/
750                while (i < 10 && j < nsamp_this_scn && j < raw_counts)
751                {
752                        sprintf(eightdigits, "        "); /* we want leading spaces */
753                        if (COSMOS0Buffer[j] != getThis->fill) {
754                                if (((int)(COSMOS0Buffer[j] * Multiplier) > 999999) || ((int)(COSMOS0Buffer[j] * Multiplier) < -99999)) {
755                                        sprintf(sample, GAP_FILL);
756                                        /* prevent out of range string */
757                                }
758                                else {
759                                        sprintf(sample, "%d", (int)(COSMOS0Buffer[j] * Multiplier));
760                                }
761                                strcpy(eightdigits + 8 - strlen(sample), sample);
762                                strcpy(hour_line + s_place, eightdigits);
763                        }
764                        else {
765                                /* We have a gap, this is where gap data is written */
766                                //strcpy(hour_line + s_place, "  9999");
767                                sprintf(sample, GAP_FILL);
768                                strcpy(eightdigits + 8 - strlen(sample), sample);
769                                strcpy(hour_line + s_place, eightdigits);                               
770                        }
771                        s_place = s_place + 8;
772
773                        total += (int)(COSMOS0Buffer[j] * Multiplier);
774
775                        j++; i++;
776                }
777
778
779                /* 401-402       Record end marker.
780                *               Two chars 'cr'= 13 and 'nl'= 10.
781                */
782                /*hour_line[77] = ' '; /*Replace that null that sprintf got us*/
783                hour_line[80] = '\n';
784
785                /* Write out line */
786                if (fwrite(&hour_line, sizeof(hour_line), 1, COSMOS0fp) != 1)
787                {
788                        logit("et", "COSMOS0PA_next: error writing COSMOS0 line. \n");
789                        return EW_FAILURE;
790                }
791        }
792        /* After we process all of the data, we have to write:
793                        End-of-Data Flag Line:
794                N+1 1 -11 End of data flag string, “End-of-data for..”.
795                21-36 Station channel number and physical parameter of data (a checksum may optionally
796                be included on the remainder of this line).  */
797        sprintf(tempstr, "                                                                                ");/*clear it*/
798        sprintf(tempstr, "End-of-data for %s.%s.%s.%s acceleration", getThis->sta, getThis->chan, getThis->net, getThis->loc );
799
800        if (fwrite(tempstr, 81, 1, COSMOS0fp) != 1)
801        {
802                logit("et", "COSMOS0PA_next: error writing COSMOS0 dynamic header line. \n");
803                return EW_FAILURE;
804        }
805        return EW_SUCCESS;
806}
807
808
809
810/************************************************************************
811* This is the Put Away end event routine. It's called after we've       *
812* finished processing one event                                         *
813*                                                                       *
814* For PC-COSMOS0 - close the COSMOS0 file, pointed to by COSMOS0fp      *
815*               free COSMOS0Buffer memory to be nice to everyone else   *
816*************************************************************************/
817int COSMOS0PA_end_ev(int debug)
818{
819        /* Actually we opened multiple files, one for each SCNL, and they should
820          all be closed by the time we get here; we don't know who they all are to
821          close them here. */
822/*      fclose(COSMOS0fp);
823
824        if (debug == 1)
825                logit("t", "Closing COSMOS0 file \n"); */
826
827        return(EW_SUCCESS);
828}
829
830/************************************************************************
831*       This is the Put Away close routine. It's called after when      *
832*       we're being shut down.                                          *
833*************************************************************************/
834int COSMOS0PA_close(int debug)
835{
836
837        free((char *)COSMOS0Buffer);
838        return(EW_SUCCESS);
839}
840
841/*
842*
843*  Byte swapping functions
844*/
845
846/* swap bytes if necessary, depending on what machine */
847/* we have been compiled, and what machine the data   */
848/* is being read from/being written to                */
849
850static int StructMakeLocal(void *ptr, int struct_type, char data_type,
851        int debug)
852{
853        if (ptr == NULL)
854        {
855                logit("e", "StructMakeLocal: NULL pointer passed in\n");
856                return EW_FAILURE;
857        }
858
859#if defined (_INTEL)
860        if (data_type != '6')
861        {
862                if (debug == 1)
863                        logit("et", "Hoping we don't really have to do any swapping from Intel to Sun \n");
864
865        }
866
867#elif defined (_SPARC)
868        if (data_type == '6')
869        {
870                if (debug == 1)
871                        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");
872        }
873
874#endif
875        return EW_SUCCESS;
876}
877
878long Median(int number_of_array_elements, long *Array)
879{
880        long *LocalArray;
881
882        LocalArray = (long*)malloc(number_of_array_elements * sizeof(long));
883
884        qsort(&LocalArray[0], number_of_array_elements, sizeof(long), longCompareForQsort);
885        /* Get Median */
886        return(LocalArray[(long)(number_of_array_elements / 2)]);
887}
888/*****************************************************************/
889/* Just a simple compare function so that Q sort does it's thing */
890/*****************************************************************/
891int longCompareForQsort(const void *x, const void *y)
892{
893        /* compare must have const void ptr params for qsort to work */
894        const long   *ix = (const long *)x;
895        const long   *iy = (const long *)y;
896        /* returns 1, -1 or 0 */
897        return (*ix > *iy) - (*ix < *iy);
898}
Note: See TracBrowser for help on using the repository browser.