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

Revision 7974, 34.4 KB checked in by stefan, 8 weeks ago (diff)

read arc file

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