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