Context Navigation

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

Visit:

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: