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

Revision 7903, 31.7 KB checked in by stefan, 4 months ago (diff)

cosmos putaway in progress

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