Context Navigation

source: trunk/src/libsrc/util/k2evt2ew.c @ 7211

Visit:

Revision 7211, 49.6 KB checked in by baker, 2 years ago (diff)
fix printf()/logit() format mismatches, use inttypes.h formats
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`

Line
1	/*
2	* THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
3	* CHECKED IT OUT USING THE COMMAND CHECKOUT.
4	*
5	* $Id$
6	*
7	* Revision history:
8	* $Log$
9	* Revision 1.11 2008/02/23 21:19:29 luetgert
10	* Made sure reported time is not in the future.
11	* .
12	*
13	* Revision 1.10 2007/08/03 05:14:51 luetgert
14	* Gravity changed from 978.03 to 981.
15	* .
16	*
17	* Revision 1.9 2005/07/20 15:08:10 friberg
18	* cleaned up a warning about tags at end of #endif statements
19	*
20	* Revision 1.8 2004/03/10 23:17:48 luetgert
21	* Make sure both the site name and serial number are checked!
22	* Kinemetrics has the nasty habit of duplicating serial numbers!
23	*
24	* Revision 1.7 2004/02/23 19:12:46 luetgert
25	* Redefined sm.t as trigger time instead of stream start time.
26	*
27	* Revision 1.6 2003/01/09 00:45:29 davidk
28	* Improved handling of evt files that are greater that the size of
29	* MAXTRACELTH, which is the maximum number of samples per channel per
30	* EVT file.
31	* Improved error reporting and handling.
32	*
33	* Revision 1.5 2002/11/07 18:33:41 lucky
34	* The conversion from box/chan to SCNL was still wrong. Fixed it again,
35	* this time hopefully for good.
36	*
37	* Revision 1.4 2002/08/16 20:51:00 alex
38	* Mod to do checksums etc. Alex
39	*
40	* Revision 1.3 2002/02/25 21:28:27 lucky
41	* Fixed the translation from box to SCN
42	*
43	* Revision 1.1 2001/11/19 16:17:59 lucky
44	* Initial revision
45	*
46	*
47	*/
48
49	/* k2hder2ew.c
50	* Read a K2-format .evt data file, crack the header, and fill
51	* the return structure to be returned to the calling program.
52	* This code was cloned from nsmp2ew; Credit Jim Luetgert.
53	*/
54
55	/* k2evt2ew() and extract() modified by DK 20030108 in order to:
56	* 1) Prevent data from being overwritten(corrupted) when the
57	* length of the EVT file exceeds the length of the K2InfoStruct.
58	* 2) Return a warning to the k2evt2ew() caller when the data output
59	* is trunctated because the length of the EVT file exceeds the
60	* length of the K2InfoStruct.
61	*/
62
63	#include <inttypes.h> /* PRIu32 */
64	#include <math.h>
65	#include <stdint.h> /* int32_t, uint32_t */
66	#include <stdio.h>
67	#include <stdlib.h>
68	#include <string.h>
69	#include <sys/types.h>
70	#include <time.h>
71
72	#include "earthworm.h"
73	#include "chron3.h"
74	#include "kom.h"
75	#include "swap.h"
76	#include "k2evt2ew.h"
77
78	#define PACKET_MAX_SIZE 3003 /* 2992 + 8 + 3, same as QT */
79	#define GRAVITY 981.0 /* Gravity in cm/sec/sec */
80	#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp
81	#define DECADE 315532800 /* Number of seconds betweeen 1/1/1970 and 1/1/1980 */
82
83	static int K2EW_Debug;
84
85	/* Functions in this source file
86	*******************************/
87	static int read_tag (FILE fp, KFF_TAG );
88	static int read_head (FILE fp, MW_HEADER , int tagLength);
89	static int read_frame (FILE fp, FRAME_HEADER , uint32_t *channels );
90
91	int extract(int pktdatalen, const unsigned char *pktbuff,
92	float Data, int32_t Counts, float scale, int nchan, int jchan,
93	int *ind, int array_size);
94	int peak_ground(float Data, int npts, int itype, float dt, SM_INFO sm);
95	void demean(float *A, int N);
96	void locut(float *s, int nd, float fcut, float delt, int nroll, int icaus);
97	void rdrvaa(float *acc, int na, float omega, float damp, float dt,
98	float rd, float rv, float aa, int maxtime);
99	void amaxper(int npts, float dt, float fc, float amaxmm,
100	float aminmm, float pmax, int imin, int imax);
101
102	/******************************************************************************
103	* k2evt2ew() Read and crack one .evt file.
104	*
105	* fp: File pointer to the .evt file to crack; file must
106	* previuosly be opened for binary reads.
107	* fname: name of the file pointed to by fp (debug purposes)
108	* pk2info: pointer to a pre-allocated structure where info from
109	* the file will be stored.
110	* pChanName: pointer to the CHANNELNAME structure containing mappings
111	* from SM Box and channel number to SCN. This is optional;
112	* if it is NULL, then default SCN names will be assigned as
113	* follows: S=box, C=channel number, N=NetCode
114	* numChans: Number of channels defined in pChanName
115	* NetCode: Network Code to use if pChanName is null.
116	* Debug: Integer Debug level.
117	*
118	* Return Values: EW_SUCCESS and EW_FAILURE, with pk2info filled with
119	* other useful info.
120	* EW_WARNING when data output is truncated for
121	* one or more channels.
122	*
123	******************************************************************************/
124
125	int k2evt2ew (FILE fp, char fname, K2InfoStruct *pk2info,
126	CHANNELNAME pChanName, int numChans, char NetCode, int Debug)
127	{
128	unsigned char g_k2mi_buff[PACKET_MAX_SIZE-9]; /* received data buffer */
129	char serial[10];
130	int i, j, ci, chanind, nscans, ret, rc=EW_SUCCESS;
131	int flen, stat, nchans, itype[18], azim[18];
132	uint32_t decade, channels;
133	float dt, scale[18], sec1970;
134
135	if ((fp == NULL) \|\| (fname == NULL) \|\| (pk2info == NULL))
136	{
137	logit ("e", "k2evt2ew: Invalid arguments passed in.\n");
138	return EW_FAILURE;
139	}
140
141	if ((strstr (fname, ".EVT") == 0) && (strstr (fname, ".evt") == 0))
142	{
143	logit ("e", "k2evt2ew: invalid file name <%s>\n", fname);
144	return EW_FAILURE;
145	}
146
147	K2EW_Debug = Debug;
148	sec1970 = (float) 11676096000.00; /* # seconds between Carl Johnson's */
149	/* time 0 and 1970-01-01 00:00:00.0 GMT */
150
151	/* Initialize variables
152	**********************/
153	memset ( &pk2info->frame, 0, sizeof(FRAME_HEADER) );
154	memset ( &pk2info->head, 0, sizeof(K2_HEADER) );
155	memset ( &pk2info->tag, 0, sizeof(KFF_TAG) );
156	memset ( &pk2info->sm, 0, MAX_SM_PER_BOX*sizeof(SM_INFO) );
157
158
159	/* Kinemetric time starts 10 years after EW time */
160	decade = (36510+2)246060;
161
162	/* First, just read the header & first frame to get some basic info.
163	******************************************************************/
164	rewind(fp);
165
166	ret = read_tag(fp, &pk2info->tag);
167	if ( ret != EW_SUCCESS ) return EW_FAILURE;
168	stat = read_head(fp, &pk2info->head, pk2info->tag.length);
169
170	ret = read_tag(fp, &pk2info->tag);
171	if ( ret != EW_SUCCESS ) return EW_FAILURE;
172	flen = pk2info->tag.dataLength;
173	read_frame(fp, &pk2info->frame, &channels);
174
175	nscans = pk2info->head.roParms.stream.duration;
176	nchans = pk2info->head.rwParms.misc.nchannels;
177	dt = (float)(1.0/(10.0(flen/(3.0nchans))));
178
179	if(pk2info->head.roParms.stream.flags != 0)
180	{
181	logit("e", "%s: not data. flags = %d\n",
182	fname, pk2info->head.roParms.stream.flags);
183	return EW_FAILURE;
184	}
185
186	sprintf(serial, "%d", pk2info->head.rwParms.misc.serialNumber);
187
188
189	for (i=0; i<nchans; i++)
190	{
191	pk2info->sm[i].t = pk2info->head.roParms.stream.triggerTime +
192	(float)pk2info->head.roParms.stream.triggerTimeMsec/1000.0;
193	/* JHL 20040223 Set the sm time to the stream start time so that peak times
194	can be calculated. It will be reset to trigger time before
195	the sm message is delivered. */
196	pk2info->sm[i].t = pk2info->head.roParms.stream.startTime +
197	(float)pk2info->head.roParms.stream.startTimeMsec/1000.0;
198	pk2info->sm[i].t += decade;
199	pk2info->sm[i].talt = (float) time(NULL);
200	/* JHL 20080223 Make sure time is not in the future */
201	if(pk2info->sm[i].t > pk2info->sm[i].talt) pk2info->sm[i].t = pk2info->sm[i].talt;
202	pk2info->sm[i].altcode = SM_ALTCODE_RECEIVING_MODULE;
203	itype[i] = 0;
204	if( pk2info->head.rwParms.channel[i].sensorType==SENSOR_VELOCITY \|\|
205	pk2info->head.rwParms.channel[i].sensorType==SENSOR_MARKL22 \|\|
206	pk2info->head.rwParms.channel[i].sensorType==SENSOR_MARKL4C)
207	itype[i] = 2;
208	else
209	if((pk2info->head.rwParms.channel[i].sensorType>=SENSOR_FBA11 &&
210	pk2info->head.rwParms.channel[i].sensorType<=SENSOR_FBA23) \|\|
211	pk2info->head.rwParms.channel[i].sensorType==SENSOR_EPI \|\|
212	pk2info->head.rwParms.channel[i].sensorType==SENSOR_ACCELERATION)
213	itype[i] = 3;
214
215	scale[i] = (float)((2.0*pk2info->head.rwParms.channel[i].fullscale/(1<<24))
216	/ pk2info->head.rwParms.channel[i].sensitivity);
217
218	if (itype[i]==3)
219	scale[i] = scale[i] * (float) GRAVITY;
220
221	if (itype[i]==0)
222	{
223	logit("e", "Can't decode sensor type! = %d\n",
224	pk2info->head.rwParms.channel[i].sensorType);
225	}
226	azim[i] = pk2info->head.rwParms.channel[i].azimuth;
227
228	/* Find this channel in the list */
229	chanind = -1;
230	for (ci = 0; ci < numChans; ci++)
231	{
232	if ((strcmp(serial, pChanName[ci].box) == 0) &&
233	(pChanName[ci].chan == i) &&
234	(strcmp(pk2info->head.rwParms.misc.stnID, pChanName[ci].sta) == 0) )
235	{
236	chanind = ci;
237	break;
238	}
239	}
240
241
242	if (chanind < 0)
243	{
244	/* build default name */
245	strcpy (pk2info->sm[i].sta, pk2info->head.rwParms.misc.stnID);
246	sprintf (pk2info->sm[i].comp, "%d", i);
247	strcpy (pk2info->sm[i].net, NetCode);
248	strcpy (pk2info->sm[i].loc, "");
249	}
250	else
251	{
252	strcpy (pk2info->sm[i].sta, pChanName[chanind].sta);
253	strcpy (pk2info->sm[i].comp, pChanName[chanind].comp);
254	strcpy (pk2info->sm[i].net, pChanName[chanind].net);
255	strcpy (pk2info->sm[i].loc, pChanName[chanind].loc);
256	}
257	}
258
259	if(K2EW_Debug)
260	logit("e", "\n\|%s\| \|%d\| \|%f\| \|%s\| \|%s\|\n",
261	pk2info->head.rwParms.misc.stnID,
262	pk2info->head.rwParms.misc.serialNumber, dt,
263	fname, pk2info->head.rwParms.misc.comment);
264
265	/* Read and process each channel individually.
266	*********************************************/
267	for(i=0;i<nchans;i++)
268	{
269	rewind(fp);
270	ret = read_tag(fp, &pk2info->tag);
271	if ( ret != EW_SUCCESS ) return EW_FAILURE;
272	stat = read_head(fp, &pk2info->head, pk2info->tag.length);
273
274	/* initialize trace buffer */
275	for (j=0; j<MAXTRACELTH; j++)
276	{
277	pk2info->Databuf[i][j] = 0.0;
278	pk2info->Counts[i][j] = 0;
279	}
280	pk2info->numDataPoints[i] = 0;
281
282
283	for(j=0;j<nscans;j++)
284	{
285	ret = read_tag(fp, &pk2info->tag);
286	if ( ret != EW_SUCCESS ) return EW_FAILURE;
287
288	if(pk2info->numDataPoints[i] > MAXTRACELTH)
289	{
290	/* DK 20030108 BEGIN */
291	/* this isn't good, this means that we've written
292	passed the end of the data allocation area in
293	the struct, thus corrupting the next channel or
294	the number of datapoints per channel.
295	************************************************/
296	rc = EW_FAILURE;
297	logit("","k2evt2ew(): exceeded MAXTRACELTH (%d <= %d)\n",
298	MAXTRACELTH, pk2info->numDataPoints[i]);
299	/* DK 20030108 END */
300	break;
301	}
302	read_frame(fp, &pk2info->frame, &channels);
303	flen = pk2info->tag.dataLength;
304	stat = (int)fread(g_k2mi_buff, 1, flen, fp);
305	ret = extract(flen, g_k2mi_buff, pk2info->Databuf[i], pk2info->Counts[i],
306	scale[i], nchans, i, &pk2info->numDataPoints[i], MAXTRACELTH);
307	if(ret)
308	{
309	if(ret == 1)
310	{
311	rc = EW_WARNING;
312	logit("t","k2evt2ew(): extract(): returned warning while processing "
313	" channel (%s,%s,%s). Data was truncated after %d datapoints.\n",
314	pk2info->sm[i].sta, pk2info->sm[i].comp, pk2info->sm[i].net,
315	pk2info->numDataPoints[i]);
316	}
317	else
318	{
319	rc = EW_FAILURE;
320	logit("t","k2evt2ew(): extract(): returned ERROR while processing "
321	" channel (%s,%s,%s).\n",
322	pk2info->sm[i].sta, pk2info->sm[i].comp, pk2info->sm[i].net);
323	}
324
325	break; /* this moves us to the next channel */
326	}
327	/* DK 20030108 Added MAXTRACELTH param to extract() call */
328	}
329
330	if(K2EW_Debug)
331	logit("e", "%2d \|%s\| \|%s\| \|%s\| %3d\n",
332	i, pk2info->sm[i].sta,
333	pk2info->sm[i].comp, pk2info->sm[i].net, azim[i]);
334
335	dt = (float)(1.0/(10.0(flen/(3.0nchans))));
336
337	peak_ground(pk2info->Databuf[i], pk2info->numDataPoints[i],
338	itype[i], dt, &pk2info->sm[i]);
339
340	/* JHL 20040223 Set the sm time to the trigger time */
341	pk2info->sm[i].t = pk2info->head.roParms.stream.triggerTime +
342	(float)pk2info->head.roParms.stream.triggerTimeMsec/1000.0;
343	pk2info->sm[i].t += decade;
344	/* JHL 20101022 Set the length of time scanned. */
345	pk2info->sm[i].length = pk2info->numDataPoints[i]*dt;
346	} /* for loop over channels */
347	return rc;
348	}
349
350	/******************************************************************************
351	* read_tag(fp) Read the 16 byte tag, swap bytes if necessary, and print. *
352	******************************************************************************/
353
354	int read_tag (FILE fp, KFF_TAG tag)
355	{
356	if ((fp == NULL) \|\| (tag == NULL))
357	{
358	logit ("e", "read_tag: Invalid arguments passed in\n");
359	return EW_FAILURE;
360	}
361
362	if (fread(tag, 1, 16, fp) != 16)
363	{
364	logit ("e", "read_tag: read of file failed.\n");
365	return EW_FAILURE;
366	}
367
368	#ifdef _INTEL
369	SwapUint32 (&tag->type);
370	SwapShort ((short *)&tag->length);
371	SwapShort ((short *)&tag->dataLength);
372	SwapShort ((short *)&tag->id);
373	SwapShort ((short *)&tag->checksum);
374	#endif /* _INTEL */
375
376	if (K2EW_Debug > 0)
377	{
378	logit("e", "TAG: %c %d %d %d %d %d %d %d %d\n",
379	tag->sync,
380	(int)tag->byteOrder,
381	(int)tag->version,
382	(int)tag->instrumentType,
383	tag->type, tag->length, tag->dataLength,
384	tag->id, tag->checksum);
385	}
386
387	/* look ahead, and check on the upcoming record for sanity
388	**********************************************************/
389	{
390	long fpos;
391	char checkBuffer[MAX_REC];
392	unsigned short checksum=0;
393	int bytesToCheck;
394	int i;
395
396	fpos = ftell (fp); /* remember where we were */
397	bytesToCheck = tag->length + tag->dataLength;
398	if (bytesToCheck > MAX_REC)
399	{
400	logit ("e", "read_tag: record too long.\n");
401	logit ("e", "record + data length > MAX_REC.\n");
402	return EW_FAILURE;
403	}
404	if (fread(checkBuffer, 1, bytesToCheck, fp) != bytesToCheck)
405	{
406	logit ("e", "read_tag: read of file failed.\n");
407	return EW_FAILURE;
408	}
409	/* look at the synch character */
410	if( tag->sync != 'K')
411	{
412	logit ("e", "read_tag: bad synch character.\n");
413	return EW_FAILURE;
414	}
415	for ( i=0; i<bytesToCheck; i++)
416	checksum = checksum + (unsigned char) checkBuffer[i];
417	if (checksum != tag->checksum)
418	{
419	logit("","read_tag: checksum error\n");
420	return EW_FAILURE;
421	}
422
423	/* now put things back the way they were */
424	fseek(fp, fpos, SEEK_SET );
425	}
426
427
428	return EW_SUCCESS;
429	}
430
431	/******************************************************************************
432	* read_head(fp) Read the file header, swap bytes if necessary, and print. *
433	******************************************************************************/
434
435	int read_head (FILE fp, MW_HEADER head, int tagLength)
436	{
437	int i, maxchans, siz;
438
439	if ((fp == NULL) \|\| (head == NULL))
440	{
441	logit ("e", "read_head: Invalid arguments passed in\n");
442	return EW_FAILURE;
443	}
444
445	/* Read in the file header.
446	If a K2, there will be 2040 bytes,
447	otherwise assume a Mt Whitney.
448	************************************/
449	maxchans = tagLength==2040? MAX_K2_CHANNELS:MAX_MW_CHANNELS;
450
451	if (fread(head, 1, 8, fp) != 8)
452	{
453	logit ("e", "read_head: read of file failed.\n");
454	return EW_FAILURE;
455	}
456
457	siz = sizeof(struct MISC_RO_PARMS) + sizeof(struct TIMING_RO_PARMS);
458	if (fread(&head->roParms.misc, 1, siz, fp) != (unsigned int) siz)
459	{
460	logit ("e", "read_head: read of file failed.\n");
461	return EW_FAILURE;
462	}
463
464	siz = sizeof(struct CHANNEL_RO_PARMS)*maxchans;
465	if (fread(&head->roParms.channel[0], 1, siz, fp) != (unsigned int) siz)
466	{
467	logit ("e", "read_head: read of file failed.\n");
468	return EW_FAILURE;
469	}
470
471	siz = sizeof(struct STREAM_RO_PARMS);
472	if (fread(&head->roParms.stream, 1, siz, fp) != (unsigned int) siz)
473	{
474	logit ("e", "read_head: read of file failed.\n");
475	return EW_FAILURE;
476	}
477
478	siz = sizeof(struct MISC_RW_PARMS)+sizeof(struct TIMING_RW_PARMS);
479	if (fread(&head->rwParms.misc, 1, siz, fp) != (unsigned int) siz)
480	{
481	logit ("e", "read_head: read of file failed.\n");
482	return EW_FAILURE;
483	}
484
485	siz = sizeof(struct CHANNEL_RW_PARMS)*maxchans;
486	if (fread(&head->rwParms.channel[0], 1, siz, fp) != (unsigned int) siz)
487	{
488	logit ("e", "read_head: read of file failed.\n");
489	return EW_FAILURE;
490	}
491
492	if(tagLength==2040) {
493	siz = sizeof(struct STREAM_K2_RW_PARMS);
494	if (fread(&head->rwParms.stream, 1, siz, fp) != (unsigned int) siz)
495	{
496	logit ("e", "read_head: read of file failed.\n");
497	return EW_FAILURE;
498	}
499	} else {
500	siz = sizeof(struct STREAM_MW_RW_PARMS);
501	if (fread(&head->rwParms.stream, 1, siz, fp) != (unsigned int) siz)
502	{
503	logit ("e", "read_head: read of file failed.\n");
504	return EW_FAILURE;
505	}
506	}
507	siz = sizeof(struct MODEM_RW_PARMS);
508	if (fread(&head->rwParms.modem, 1, siz, fp) != (unsigned int) siz)
509	{
510	logit ("e", "read_head: read of file failed.\n");
511	return EW_FAILURE;
512	}
513
514	#ifdef _INTEL
515	/* Byte-swap values, if necessary */
516
517	/* roParms */
518	SwapShort ((short *)&head->roParms.headerVersion);
519	SwapShort ((short *)&head->roParms.headerBytes);
520
521	SwapShort ((short *)&head->roParms.misc.installedChan);
522	SwapShort ((short *)&head->roParms.misc.maxChannels);
523	SwapShort ((short *)&head->roParms.misc.sysBlkVersion);
524	SwapShort ((short *)&head->roParms.misc.bootBlkVersion);
525	SwapShort ((short *)&head->roParms.misc.appBlkVersion);
526	SwapShort ((short *)&head->roParms.misc.dspBlkVersion);
527	SwapShort (&head->roParms.misc.batteryVoltage);
528	SwapShort ((short *)&head->roParms.misc.crc);
529	SwapShort ((short *)&head->roParms.misc.flags);
530	SwapShort (&head->roParms.misc.temperature);
531
532	SwapShort ((short *)&head->roParms.timing.gpsLockFailCount);
533	SwapShort ((short *)&head->roParms.timing.gpsUpdateRTCCount);
534	SwapShort (&head->roParms.timing.acqDelay);
535	SwapShort (&head->roParms.timing.gpsLatitude);
536	SwapShort (&head->roParms.timing.gpsLongitude);
537	SwapShort (&head->roParms.timing.gpsAltitude);
538	SwapShort ((short *)&head->roParms.timing.dacCount);
539	SwapShort (&head->roParms.timing.gpsLastDrift[0]);
540	SwapShort (&head->roParms.timing.gpsLastDrift[1]);
541	SwapUint32 (&head->roParms.timing.gpsLastTurnOnTime[0]);
542	SwapUint32 (&head->roParms.timing.gpsLastTurnOnTime[1]);
543	SwapUint32 (&head->roParms.timing.gpsLastUpdateTime[0]);
544	SwapUint32 (&head->roParms.timing.gpsLastUpdateTime[1]);
545	SwapUint32 (&head->roParms.timing.gpsLastLockTime[0]);
546	SwapUint32 (&head->roParms.timing.gpsLastLockTime[1]);
547
548	SwapUint32 (&head->roParms.stream.startTime);
549	SwapUint32 (&head->roParms.stream.triggerTime);
550	SwapUint32 (&head->roParms.stream.duration);
551	SwapShort ((short *)&head->roParms.stream.errors);
552	SwapShort ((short *)&head->roParms.stream.flags);
553	SwapShort ((short *)&head->roParms.stream.startTimeMsec);
554	SwapShort ((short *)&head->roParms.stream.triggerTimeMsec);
555	SwapUint32 (&head->roParms.stream.nscans);
556
557
558	for(i=0;i<head->roParms.misc.maxChannels;i++)
559	{
560	SwapInt32 (&head->roParms.channel[i].maxPeak);
561	SwapUint32 (&head->roParms.channel[i].maxPeakOffset);
562	SwapInt32 (&head->roParms.channel[i].minPeak);
563	SwapUint32 (&head->roParms.channel[i].minPeakOffset);
564	SwapInt32 (&head->roParms.channel[i].mean);
565	}
566
567
568	/* rwParams */
569	SwapShort ((short *)&head->rwParms.misc.serialNumber);
570	SwapShort ((short *)&head->rwParms.misc.nchannels);
571	SwapShort (&head->rwParms.misc.elevation);
572	SwapFloat (&head->rwParms.misc.latitude);
573	SwapFloat (&head->rwParms.misc.longitude);
574	SwapShort (&head->rwParms.misc.userCodes[0]);
575	SwapShort (&head->rwParms.misc.userCodes[1]);
576	SwapShort (&head->rwParms.misc.userCodes[2]);
577	SwapShort (&head->rwParms.misc.userCodes[3]);
578	SwapUint32 (&head->rwParms.misc.cutler_bitmap);
579	SwapUint32 (&head->rwParms.misc.channel_bitmap);
580
581	SwapShort (&head->rwParms.timing.localOffset);
582
583	for(i=0;i<head->rwParms.misc.nchannels;i++)
584	{
585	SwapShort ((short *)&head->rwParms.channel[i].sensorSerialNumberExt);
586	SwapShort (&head->rwParms.channel[i].north);
587	SwapShort (&head->rwParms.channel[i].east);
588	SwapShort (&head->rwParms.channel[i].up);
589	SwapShort (&head->rwParms.channel[i].altitude);
590	SwapShort (&head->rwParms.channel[i].azimuth);
591	SwapShort ((short *)&head->rwParms.channel[i].sensorType);
592	SwapShort ((short *)&head->rwParms.channel[i].sensorSerialNumber);
593	SwapShort ((short *)&head->rwParms.channel[i].gain);
594	SwapShort ((short *)&head->rwParms.channel[i].StaLtaRatio);
595	SwapFloat (&head->rwParms.channel[i].fullscale);
596	SwapFloat (&head->rwParms.channel[i].sensitivity);
597	SwapFloat (&head->rwParms.channel[i].damping);
598	SwapFloat (&head->rwParms.channel[i].naturalFrequency);
599	SwapFloat (&head->rwParms.channel[i].triggerThreshold);
600	SwapFloat (&head->rwParms.channel[i].detriggerThreshold);
601	SwapFloat (&head->rwParms.channel[i].alarmTriggerThreshold);
602	}
603
604	SwapShort ((short *)&head->rwParms.stream.eventNumber);
605	SwapShort ((short *)&head->rwParms.stream.sps);
606	SwapShort ((short *)&head->rwParms.stream.apw);
607	SwapShort ((short *)&head->rwParms.stream.preEvent);
608	SwapShort ((short *)&head->rwParms.stream.postEvent);
609	SwapShort ((short *)&head->rwParms.stream.minRunTime);
610	SwapShort (&head->rwParms.stream.VotesToTrigger);
611	SwapShort (&head->rwParms.stream.VotesToDetrigger);
612	SwapShort (&head->rwParms.stream.Timeout);
613	SwapShort ((short *)&head->rwParms.stream.TxBlkSize);
614	SwapShort ((short *)&head->rwParms.stream.BufferSize);
615	SwapShort ((short *)&head->rwParms.stream.SampleRate);
616	SwapUint32 (&head->rwParms.stream.TxChanMap);
617	#endif /* _INTEL */
618
619	if(K2EW_Debug > 0)
620	{
621	logit("e", "HEADER: %c%c%c %d %hu %hu\n",
622	head->roParms.id[0], head->roParms.id[1], head->roParms.id[2],
623	(int)head->roParms.instrumentCode,
624	head->roParms.headerVersion,
625	head->roParms.headerBytes);
626
627	logit("e", "MISC_RO_PARMS: %d %d %d %hu %hu %hu %hu %hu %hu %d %hu %hu %d\n",
628	(int)head->roParms.misc.a2dBits,
629	(int)head->roParms.misc.sampleBytes,
630	(int)head->roParms.misc.restartSource,
631	head->roParms.misc.installedChan,
632	head->roParms.misc.maxChannels,
633	head->roParms.misc.sysBlkVersion,
634	head->roParms.misc.bootBlkVersion,
635	head->roParms.misc.appBlkVersion,
636	head->roParms.misc.dspBlkVersion,
637	head->roParms.misc.batteryVoltage,
638	head->roParms.misc.crc,
639	head->roParms.misc.flags,
640	head->roParms.misc.temperature );
641
642
643	logit("e", "TIMING_RO_PARMS: %d %d %d %hu %hu %d %d %d %d %hu %d %d"
644	" %" PRIu32 " %" PRIu32 " %" PRIu32 " %" PRIu32
645	" %" PRIu32 " %" PRIu32 "\n",
646	(int)head->roParms.timing.clockSource,
647	(int)head->roParms.timing.gpsStatus,
648	(int)head->roParms.timing.gpsSOH,
649	head->roParms.timing.gpsLockFailCount,
650	head->roParms.timing.gpsUpdateRTCCount,
651	head->roParms.timing.acqDelay,
652	head->roParms.timing.gpsLatitude,
653	head->roParms.timing.gpsLongitude,
654	head->roParms.timing.gpsAltitude,
655	head->roParms.timing.dacCount,
656	head->roParms.timing.gpsLastDrift[0],
657	head->roParms.timing.gpsLastDrift[1],
658	head->roParms.timing.gpsLastTurnOnTime[0],
659	head->roParms.timing.gpsLastTurnOnTime[1],
660	head->roParms.timing.gpsLastUpdateTime[0],
661	head->roParms.timing.gpsLastUpdateTime[1],
662	head->roParms.timing.gpsLastLockTime[0],
663	head->roParms.timing.gpsLastLockTime[1] );
664
665
666	for(i=0;i<head->roParms.misc.maxChannels;i++)
667	{
668	logit("e", "CHANNEL_RO_PARMS: %d %" PRIu32 " %d %" PRIu32 " %d %d\n",
669	head->roParms.channel[i].maxPeak,
670	head->roParms.channel[i].maxPeakOffset,
671	head->roParms.channel[i].minPeak,
672	head->roParms.channel[i].minPeakOffset,
673	head->roParms.channel[i].mean,
674	head->roParms.channel[i].aqOffset );
675	}
676
677	logit("e", "STREAM_RO_PARMS: %d %d %d %d %d %d %d %d\n",
678	head->roParms.stream.startTime,
679	head->roParms.stream.triggerTime,
680	head->roParms.stream.duration,
681	head->roParms.stream.errors,
682	head->roParms.stream.flags,
683	head->roParms.stream.startTimeMsec,
684	head->roParms.stream.triggerTimeMsec,
685	head->roParms.stream.nscans );
686
687	logit("e", "MISC_RW_PARMS: %hu %hu %.5s %.33s %d %f %f %d %d %d %d"
688	" %d %" PRIo32 " %d %.17s %d %d\n",
689	head->rwParms.misc.serialNumber,
690	head->rwParms.misc.nchannels,
691	head->rwParms.misc.stnID,
692	head->rwParms.misc.comment,
693	head->rwParms.misc.elevation,
694	head->rwParms.misc.latitude,
695	head->rwParms.misc.longitude,
696	(int)head->rwParms.misc.cutlerCode,
697	(int)head->rwParms.misc.minBatteryVoltage,
698	(int)head->rwParms.misc.cutler_decimation,
699	(int)head->rwParms.misc.cutler_irig_type,
700	head->rwParms.misc.cutler_bitmap,
701	head->rwParms.misc.channel_bitmap,
702	(int)head->rwParms.misc.cutler_protocol,
703	head->rwParms.misc.siteID,
704	(int)head->rwParms.misc.externalTrigger,
705	(int)head->rwParms.misc.networkFlag );
706
707	logit("e", "TIMING_RW_PARMS: %d %d %hu\n",
708	(int)head->rwParms.timing.gpsTurnOnInterval,
709	(int)head->rwParms.timing.gpsMaxTurnOnTime,
710	head->rwParms.timing.localOffset );
711
712	for(i=0;i<head->rwParms.misc.nchannels;i++)
713	{
714	logit("e", "CHANNEL_RW_PARMS: %s %d %d %d %d %d %d %d %d %d %d"
715	" %d %d %d %d %d %f %f %f %f %f %f %f\n",
716	head->rwParms.channel[i].id,
717	head->rwParms.channel[i].sensorSerialNumberExt,
718	head->rwParms.channel[i].north,
719	head->rwParms.channel[i].east,
720	head->rwParms.channel[i].up,
721	head->rwParms.channel[i].altitude,
722	head->rwParms.channel[i].azimuth,
723	head->rwParms.channel[i].sensorType,
724	head->rwParms.channel[i].sensorSerialNumber,
725	head->rwParms.channel[i].gain,
726	(int)head->rwParms.channel[i].triggerType,
727	(int)head->rwParms.channel[i].iirTriggerFilter,
728	(int)head->rwParms.channel[i].StaSeconds,
729	(int)head->rwParms.channel[i].LtaSeconds,
730	head->rwParms.channel[i].StaLtaRatio,
731	(int)head->rwParms.channel[i].StaLtaPercent,
732	head->rwParms.channel[i].fullscale,
733	head->rwParms.channel[i].sensitivity,
734	head->rwParms.channel[i].damping,
735	head->rwParms.channel[i].naturalFrequency,
736	head->rwParms.channel[i].triggerThreshold,
737	head->rwParms.channel[i].detriggerThreshold,
738	head->rwParms.channel[i].alarmTriggerThreshold);
739	}
740
741	logit("e", "STREAM_K2_RW_PARMS: \|%d\| \|%d\| \|%d\| %d %d %d %d %d %d %d"
742	" %d %d %d %d %d %d %d %" PRIu32 "\n",
743	(int)head->rwParms.stream.filterFlag,
744	(int)head->rwParms.stream.primaryStorage,
745	(int)head->rwParms.stream.secondaryStorage,
746	head->rwParms.stream.eventNumber,
747	head->rwParms.stream.sps,
748	head->rwParms.stream.apw,
749	head->rwParms.stream.preEvent,
750	head->rwParms.stream.postEvent,
751	head->rwParms.stream.minRunTime,
752	head->rwParms.stream.VotesToTrigger,
753	head->rwParms.stream.VotesToDetrigger,
754	(int)head->rwParms.stream.FilterType,
755	(int)head->rwParms.stream.DataFmt,
756	head->rwParms.stream.Timeout,
757	head->rwParms.stream.TxBlkSize,
758	head->rwParms.stream.BufferSize,
759	head->rwParms.stream.SampleRate,
760	head->rwParms.stream.TxChanMap);
761
762	} /* If debug */
763
764	return EW_SUCCESS;
765
766	} /* read_head */
767
768	/******************************************************************************
769	* read_frame(fp) Read the frame header, swap bytes if necessary. *
770	******************************************************************************/
771
772	int read_frame (FILE fp, FRAME_HEADER frame, uint32_t *channels )
773	{
774	unsigned short frameStatus, frameStatus2, samprate, streamnumber;
775	unsigned char BitMap[4];
776	uint32_t bmap;
777
778	if ((fp == NULL) \|\| (frame == NULL) \|\| (channels == NULL))
779	{
780	logit ("e", "read_frame: invalid arguments passed in.\n");
781	return EW_FAILURE;
782	}
783
784	if (fread(frame, 32, 1, fp) != 1)
785	{
786	logit ("e", "read_frame: read of file failed.\n");
787	return EW_FAILURE;
788	}
789
790	#ifdef _INTEL
791	SwapShort ((short *)&frame->recorderID);
792	SwapShort ((short *)&frame->frameSize);
793	SwapShort ((short *)&frame->blockTime1);
794	SwapShort ((short *)&frame->blockTime2);
795	SwapShort ((short *)&frame->channelBitMap);
796	SwapShort ((short *)&frame->streamPar);
797	SwapShort ((short *)&frame->msec);
798	#endif /* _INTEL */
799
800	BitMap[0] = frame->channelBitMap & 255;
801	BitMap[1] = frame->channelBitMap >> 8;
802	BitMap[2] = frame->channelBitMap1;
803	BitMap[3] = 0;
804
805	bmap = ((uint32_t )(BitMap));
806	frameStatus = frame->frameStatus;
807	frameStatus2 = frame->frameStatus2;
808	samprate = frame->streamPar & 4095;
809	streamnumber = frame->streamPar >> 12;
810
811	if(K2EW_Debug > 0)
812	logit("e", "FRAME: %d %d %d %d %" PRIu32 " X%0hX "
813	"%hu X%0hX %hu %hu X%0hX X%0hX %hu X%0hX\n",
814	(int)frame->frameType,
815	(int)frame->instrumentCode,
816	frame->recorderID,
817	frame->frameSize,
818	(uint32_t )&frame->blockTime1,
819	frame->channelBitMap,
820	frame->streamPar,
821	streamnumber, samprate, (unsigned short) ( samprate >> 8 ),
822	frameStatus,
823	frameStatus2,
824	frame->msec,
825	(unsigned short)frame->channelBitMap1);
826
827	*channels = bmap;
828
829	return EW_SUCCESS;
830	}
831
832	/******************************************************************************
833	* extract: Processes k2 file buffers to expand data from 24 bits to *
834	* 36 bits, sort channels and scale to cgs units. *
835	* pktdatalen - length of data in buffer *
836	* pktbuff - address of buffer containing the data bytes *
837	* Data - the destination array. *
838	* scale - the scaling factor (units/count). *
839	* nchan - the total number of channels in the buffer. *
840	* jchan - the channel to be extracted. *
841	* ind - the current index in the destination array. *
842	* array_size - the size of the destination array. *
843	* *
844	* return values: *
845	* 0 Success *
846	* 1 Extraction was truncated because packet *
847	* would have cause overrun in destination array. *
848	******************************************************************************/
849
850	int extract(int pktdatalen, const unsigned char pktbuff, float Data,
851	int32_t Counts, float scale, int nchan, int jchan, int ind,
852	int array_size)
853	{
854	int k, pktpos, ichan;
855	unsigned char btarr[4];
856	int rc = 0; /* DK 20030108 Added so that we could set return code in loop */
857
858
859	/* process and store serial data stream values */
860	k = *ind;
861	pktpos = ichan = 0;
862	while (pktpos < pktdatalen) {
863	/* loop through each serial data stream sample data value */
864	if(ichan == jchan) {
865	btarr[1] = pktbuff[pktpos++]; /* copy MSByte of 24-bit value */
866	btarr[2] = pktbuff[pktpos++]; /* copy middle-byte of 24-bit value */
867	btarr[3] = pktbuff[pktpos++]; /* copy LSByte of 24-bit value */
868	/* sign extend to 32-bits */
869
870	btarr[0] = ((btarr[1] & (unsigned char)0x80) ==
871	(unsigned char)0) ? (unsigned char)0 : (unsigned char)0xFF;
872
873	/* byte-swap and enter value */
874	#ifdef _INTEL
875	SwapInt32 ((int32_t *) btarr);
876	#endif
877	/* DK 20030108 Make sure we don't overrun the destination arrays */
878	if(k >= array_size)
879	{
880	/* Uh Oh! We are about to write past the end of the array.
881	Truncate operations here.
882	************************************************************/
883	rc = 1;
884	break;
885	}
886
887	Counts[k] = (((int32_t )(btarr)));
888	Data[k] = Counts[k]*scale;
889	k = k + 1;
890	} else {
891	pktpos += 3;
892	}
893	ichan += 1;
894	if(ichan >= nchan) ichan = 0;
895	}
896
897	*ind = k;
898	return rc; /* DK return rc which defaults to 0, but may be 1 */
899	}
900
901	/******************************************************************************
902	* subroutine for estimation of ground motion *
903	* *
904	* input: *
905	* Data - data array *
906	* npts - number of points in timeseries *
907	* itype - units for timeseries. 1=disp 2=vel 3=acc *
908	* dt - sample spacing in sec *
909	* return: *
910	* 0 - All OK *
911	* 1 - error *
912	* *
913	******************************************************************************/
914
915	int peak_ground(float Data, int npts, int itype, float dt, SM_INFO sm)
916	{
917	int imax_acc, imax_vel, imax_dsp, imax_raw;
918	int imin_acc, imin_vel, imin_dsp, imin_raw;
919	int ii, kk, kpts, id[4], npd[2], icaus, maxtime;
920	float totint, a, tpi, omega, damp, rd, rv, aa;
921	float amax_acc, amin_acc, pmax_acc;
922	float amax_vel, amin_vel, pmax_vel;
923	float amax_dsp, amin_dsp, pmax_dsp;
924	float amax_raw, amin_raw, pmax_raw, gd[4], sd[4];
925
926	/* Made these float arrays static because Solaris was Segfaulting on an
927	* allocation this big from the stack. These currently are 3.2 MB each
928	* DK 20030108
929	************************************************************************/
930	static float d1[MAXTRACELTH];
931	static float d2[MAXTRACELTH];
932	static float d3[MAXTRACELTH];
933
934	gd[0] = 0.05f;
935	gd[1] = 0.10f;
936	gd[2] = 0.20f;
937	gd[3] = 0.50f;
938	icaus = 1;
939
940	tpi = (float)(8.0*atan(1.0));
941
942	/* Find the raw maximum and its period
943	*************************************/
944
945	demean(Data, npts);
946	amaxper(npts, dt, Data, &amin_raw, &amax_raw, &pmax_raw, &imin_raw, &imax_raw);
947
948	if(itype == 1) { /* input data is displacement */
949	for(kk=0;kk<npts;kk++) d1[kk] = Data[kk];
950	locut(d1, npts, 0.17f, dt, 2, icaus);
951	for(kk=1;kk<npts;kk++) {
952	d2[kk] = (d1[kk] - d1[kk-1])/dt;
953	}
954	d2[0] = d2[1];
955	demean(d2, npts);
956	for(kk=1;kk<npts;kk++) {
957	d3[kk] = (d2[kk] - d2[kk-1])/dt;
958	}
959	d3[0] = d3[1];
960	demean(d3, npts);
961	} else
962	if(itype == 2) { /* input data is velocity */
963	for(kk=0;kk<npts;kk++) d2[kk] = Data[kk];
964	locut(d2, npts, 0.17f, dt, 2, icaus);
965	for(kk=1;kk<npts;kk++) {
966	d3[kk] = (d2[kk] - d2[kk-1])/dt;
967	}
968	d3[0] = d3[1];
969	demean(d3, npts);
970
971	totint = 0.0;
972	for(kk=0;kk<npts-1;kk++) {
973	totint = totint + (d2[kk] + d2[kk+1])0.5fdt;
974	d1[kk] = totint;
975	}
976	d1[npts-1] = d1[npts-2];
977	demean(d1, npts);
978
979	} else
980	if(itype == 3) { /* input data is acceleration */
981	for(kk=0;kk<npts;kk++) d3[kk] = Data[kk];
982	locut(d3, npts, 0.17f, dt, 2, icaus);
983
984	totint = 0.0;
985	for(kk=0;kk<npts-1;kk++) {
986	totint = totint + (d3[kk] + d3[kk+1])0.5fdt;
987	d2[kk] = totint;
988	}
989	d2[npts-1] = d2[npts-2];
990	demean(d2, npts);
991
992	totint = 0.0;
993	for(kk=0;kk<npts-1;kk++) {
994	totint = totint + (d2[kk] + d2[kk+1])0.5fdt;
995	d1[kk] = totint;
996	}
997	d1[npts-1] = d1[npts-2];
998	demean(d1, npts);
999	} else {
1000	return 1;
1001	}
1002
1003	/* Find the displacement(cm), velocity(cm/s), & acceleration(cm/s/s) maxima and their periods
1004	*********************************************************************************************/
1005
1006	amaxper(npts, dt, &d1[0], &amin_dsp, &amax_dsp, &pmax_dsp, &imin_dsp, &imax_dsp);
1007	amaxper(npts, dt, &d2[0], &amin_vel, &amax_vel, &pmax_vel, &imin_vel, &imax_vel);
1008	amaxper(npts, dt, &d3[0], &amin_acc, &amax_acc, &pmax_acc, &imin_acc, &imax_acc);
1009
1010	/* Find the spectral response
1011	****************************/
1012
1013	damp = 0.05f;
1014	kk = 0;
1015	sm->pdrsa[kk] = 0.3;
1016	omega = (float)(tpi/sm->pdrsa[kk]);
1017	rdrvaa(&d3[0], npts-1, omega, damp, dt, &rd, &rv, &aa, &maxtime);
1018	sm->rsa[kk] = aa;
1019	sm->trsa[kk] = sm->t + dt*maxtime;
1020	kk += 1;
1021
1022	sm->pdrsa[kk] = 1.0;
1023	omega = (float)(tpi/sm->pdrsa[kk]);
1024	rdrvaa(&d3[0], npts-1, omega, damp, dt, &rd, &rv, &aa, &maxtime);
1025	sm->rsa[kk] = aa;
1026	sm->trsa[kk] = sm->t + dt*maxtime;
1027	kk += 1;
1028
1029	sm->pdrsa[kk] = 3.0;
1030	omega = (float)(tpi/sm->pdrsa[kk]);
1031	rdrvaa(&d3[0], npts-1, omega, damp, dt, &rd, &rv, &aa, &maxtime);
1032	sm->rsa[kk] = aa;
1033	sm->trsa[kk] = sm->t + dt*maxtime;
1034	kk += 1;
1035
1036	sm->nrsa = kk;
1037
1038	/* Since we are here, determine the duration of strong shaking
1039	*************************************************************/
1040
1041	for(kk=0;kk<4;kk++) {
1042	id[kk] = 0;
1043	for(ii=1;ii<=npts-1;ii++) {
1044	a = (float)fabs(d3[ii]/GRAVITY);
1045	if (a >= gd[kk]) {
1046	id[kk]++;
1047	if (id[kk] == 1) npd[0] = ii;
1048	npd[1] = ii;
1049	}
1050	}
1051	if (id[kk] > 0) {
1052	kpts = npd[1] - npd[0] + 1;
1053	sd[kk] = kpts*dt;
1054	} else {
1055	sd[kk] = 0.0;
1056	}
1057	}
1058
1059	sm->pgd = fabs(amin_dsp)>fabs(amax_dsp)? fabs(amin_dsp):fabs(amax_dsp);
1060	sm->pgv = fabs(amin_vel)>fabs(amax_vel)? fabs(amin_vel):fabs(amax_vel);
1061	sm->pga = fabs(amin_acc)>fabs(amax_acc)? fabs(amin_acc):fabs(amax_acc);
1062
1063	sm->tpgd = fabs(amin_dsp)>fabs(amax_dsp)? sm->t + dtimin_dsp:sm->t + dtimax_dsp;
1064	sm->tpgv = fabs(amin_vel)>fabs(amax_vel)? sm->t + dtimin_vel:sm->t + dtimax_vel;
1065	sm->tpga = fabs(amin_acc)>fabs(amax_acc)? sm->t + dtimin_acc:sm->t + dtimax_acc;
1066
1067	return 0;
1068	}
1069
1070	/******************************************************************************
1071	* demean removes the mean from the n point series stored in array A. *
1072	* *
1073	******************************************************************************/
1074
1075	void demean(float *A, int n)
1076	{
1077	int i;
1078	float xm;
1079
1080	xm = 0.0;
1081	for(i=0;i<n;i++) xm = xm + A[i];
1082	xm = xm/n;
1083	for(i=0;i<n;i++) A[i] = A[i] - xm;
1084	}
1085
1086	/******************************************************************************
1087	* Butterworth locut filter order 2nroll (nroll<=8)
1088	* (see Kanasewich, Time Sequence Analysis in Geophysics, *
1089	* Third Edition, University of Alberta Press, 1981) *
1090	* written by W. B. Joyner 01/07/97 *
1091	* *
1092	* s[j] input = the time series to be filtered *
1093	* output = the filtered series *
1094	* dimension of s[j] must be at least as large as *
1095	* nd+3.0float(nroll)/(fcutdelt) *
1096	* nd = the number of points in the time series *
1097	* fcut = the cutoff frequency *
1098	* delt = the timestep *
1099	* nroll = filter order *
1100	* causal if icaus.eq.1 - zero phase shift otherwise *
1101	* *
1102	* The response is given by eq. 15.8-6 in Kanasewich: *
1103	* Y = sqrt((f/fcut)*(2n)/(1+(f/fcut)*(2n))), *
1104	* where n = 2nroll
1105	* *
1106	* Dates: 01/07/97 - Written by Bill Joyner *
1107	* 12/17/99 - D. Boore added check for fcut = 0.0, in which case *
1108	* no filter is applied. He also cleaned up the *
1109	* appearance of the code (indented statements in *
1110	* loops, etc.) *
1111	* 02/04/00 - Changed "n" to "nroll" to eliminate confusion with *
1112	* Kanesewich, who uses "n" as the order (=2nroll)
1113	* 03/01/00 - Ported to C by Jim Luetgert *
1114	* *
1115	******************************************************************************/
1116
1117	void locut(float *s, int nd, float fcut, float delt, int nroll, int icaus)
1118	{
1119	float fact[8], b1[8], b2[8];
1120	float pi, w0, w1, w2, w3, w4, w5, xp, yp, x1, x2, y1, y2;
1121	int j, k, np2, npad;
1122
1123	if (fcut == 0.0) return; /* Added by DMB */
1124
1125	pi = (float)(4.0*atan(1.0));
1126	w0 = (float)(2.0pifcut);
1127	w1 = (float)(2.0tan(w0delt/2.0));
1128	w2 = (float)((w1/2.0)*(w1/2.0));
1129	w3 = (float)((w1*w1)/2.0 - 2.0);
1130	w4 = (float)(0.25*pi/nroll);
1131
1132	for(k=0;k<nroll;k++) {
1133	w5 = (float)(w4(2.0k + 1.0));
1134	fact[k] = (float)(1.0/(1.0+sin(w5)*w1 + w2));
1135	b1[k] = w3*fact[k];
1136	b2[k] = (float)((1.0-sin(w5)w1 + w2)fact[k]);
1137	}
1138
1139	np2 = nd;
1140
1141	if(icaus != 1) {
1142	npad = (int)(3.0nroll/(fcutdelt));
1143	np2 = nd+npad;
1144	for(j=nd;j<np2;j++) s[j] = 0.0;
1145	}
1146
1147	for(k=0;k<nroll;k++) {
1148	x1 = x2 = y1 = y2 = 0.0;
1149	for(j=0;j<np2;j++) {
1150	xp = s[j];
1151	yp = (float)(fact[k](xp-2.0x1+x2) - b1[k]y1 - b2[k]y2);
1152	s[j] = yp;
1153	y2 = y1;
1154	y1 = yp;
1155	x2 = x1;
1156	x1 = xp;
1157	}
1158	}
1159
1160	if(icaus != 1) {
1161	for(k=0;k<nroll;k++) {
1162	x1 = x2 = y1 = y2 = 0.0;
1163	for(j=0;j<np2;j++) {
1164	xp = s[np2-j-1];
1165	yp = (float)(fact[k](xp-2.0x1+x2) - b1[k]y1 - b2[k]y2);
1166	s[np2-j-1] = yp;
1167	y2 = y1;
1168	y1 = yp;
1169	x2 = x1;
1170	x1 = xp;
1171	}
1172	}
1173	}
1174
1175	return;
1176	}
1177
1178	/******************************************************************************
1179	* rdrvaa *
1180	* *
1181	* This is a modified version of "Quake.For", originally *
1182	* written by J.M. Roesset in 1971 and modified by *
1183	* Stavros A. Anagnostopoulos, Oct. 1986. The formulation is *
1184	* that of Nigam and Jennings (BSSA, v. 59, 909-922, 1969). *
1185	* Dates: 02/11/00 - Modified by David M. Boore, based on RD_CALC *
1186	* 03/01/00 - Ported to C by Jim Luetgert *
1187	* *
1188	* acc = acceleration time series *
1189	* na = length of time series *
1190	* omega = 2pi/per
1191	* damp = fractional damping (e.g., 0.05) *
1192	* dt = time spacing of input *
1193	* rd = relative displacement of oscillator *
1194	* rv = relative velocity of oscillator *
1195	* aa = absolute acceleration of oscillator *
1196	* maxtime = time of max aa *
1197	******************************************************************************/
1198
1199	void rdrvaa(float *acc, int na, float omega, float damp, float dt,
1200	float rd, float rv, float aa, int maxtime)
1201	{
1202	float omt, d2, bom, d3, omd, om2, omdt, c1, c2, c3, c4, cc, ee;
1203	float s1, s2, s3, s4, s5, a11, a12, a21, a22, b11, b12, b21, b22;
1204	float y, ydot, y1, z, z1, z2, ra;
1205	int i;
1206
1207	omt = omega*dt;
1208	d2 = (float)(sqrt(1.0-damp*damp));
1209	bom = damp*omega;
1210	d3 = (float)(2.0*bom);
1211	omd = omega*d2;
1212	om2 = omega*omega;
1213	omdt = omd*dt;
1214	c1 = (float)(1.0/om2);
1215	c2 = (float)(2.0damp/(om2omt));
1216	c3 = c1+c2;
1217	c4 = (float)(1.0/(omega*omt));
1218	ee = (float)(exp(-damp*omt));
1219	cc = (float)(cos(omdt)*ee);
1220	s1 = (float)(sin(omdt)*ee/omd);
1221	s2 = s1*bom;
1222	s3 = s2 + cc;
1223	s4 = (float)(c4*(1.0-s3));
1224	s5 = s1*c4 + c2;
1225
1226	a11 = s3; a12 = s1;
1227	a21 = -om2*s1; a22 = cc - s2;
1228
1229	b11 = s3c3 - s5; b12 = -c2s3 + s5 - c1;
1230	b21 = -s1+s4; b22 = -s4;
1231
1232	y = ydot = rd = rv = *aa = 0.0;
1233	*maxtime = -1;
1234	for(i=0;i<na-1;i++) {
1235	y1 = a11y + a12ydot + b11acc[i] + b12acc[i+1];
1236	ydot = a21y + a22ydot + b21acc[i] + b22acc[i+1];
1237	y = y1;/* y is the oscillator output at time corresponding to index i */
1238	z = (float)fabs(y);
1239	if (z > rd) rd = z;
1240	z1 = (float)fabs(ydot);
1241	if (z1 > rv) rv = z1;
1242	ra = -d3ydot -om2y1;
1243	z2 = (float)fabs(ra);
1244	if (z2 > *aa) {
1245	*aa = z2;
1246	*maxtime = i;
1247	}
1248	}
1249	}
1250
1251	/******************************************************************************
1252	* compute maximum amplitude and its associated period *
1253	* *
1254	* input: *
1255	* npts - number of points in timeseries *
1256	* dt - sample spacing in sec *
1257	* fc - input timeseries *
1258	* output: *
1259	* amaxmm - raw maximum *
1260	* pmax - period of maximum *
1261	* imax - index of maxmimum point *
1262	* *
1263	******************************************************************************/
1264
1265	void amaxper(int npts, float dt, float fc, float aminmm, float *amaxmm,
1266	float pmax, int imin, int *imax)
1267	{
1268	float amin, amax, pp, pm, mean, frac;
1269	int i, j, jmin, jmax;
1270
1271	imax = jmax = imin = jmin = 0;
1272	amax = amin = amaxmm = aminmm = fc[0];
1273	aminmm = pmax = mean = 0.0;
1274	for(i=0;i<npts;i++) {
1275	mean = mean + fc[i]/npts;
1276	if (fc[i] > amax) { jmax = i; amax = fc[i]; }
1277	if (fc[i] < amin) { jmin = i; amin = fc[i]; }
1278	}
1279
1280	/* compute period of maximum */
1281
1282	pp = pm = 0.0;
1283	if (fc[jmax] > mean) {
1284	j = jmax+1;
1285	while(fc[j] > mean && j < npts) {
1286	pp += dt;
1287	j += 1;
1288	}
1289	frac = dt*(mean-fc[j-1])/(fc[j]-fc[j-1]);
1290	frac = 0.0;
1291	pp = pp + frac;
1292	j = jmax-1;
1293	while(fc[j] > mean && j >= 0) {
1294	pm += dt;
1295	j -= 1;
1296	}
1297	frac = dt*(mean-fc[j+1])/(fc[j]-fc[j+1]);
1298	frac = 0.0;
1299	pm = pm + frac;
1300	} else {
1301	j = jmax+1;
1302	if(fc[j] < mean && j < npts) {
1303	pp += dt;
1304	j += 1;
1305	}
1306	frac = dt*(mean-fc[j-1])/(fc[j]-fc[j-1]);
1307	frac = 0.0;
1308	pp = pp + frac;
1309	j = jmax-1;
1310	if(fc[j] < mean && j >= 0) {
1311	pm += dt;
1312	j -= 1;
1313	}
1314	frac = dt*(mean-fc[j+1])/(fc[j]-fc[j+1]);
1315	frac = 0.0;
1316	pm = pm + frac;
1317	}
1318
1319	*imin = jmin;
1320	*imax = jmax;
1321	pmax = (float)(2.0(pm+pp));
1322	*aminmm = amin;
1323	*amaxmm = amax;
1324
1325	return;
1326	}

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

Download in other formats: