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 = (365*10+2)*24*60*60; |
---|
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.0*nchans)))); |
---|
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.0*nchans)))); |
---|
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.5f*dt; |
---|
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.5f*dt; |
---|
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.5f*dt; |
---|
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 + dt*imin_dsp:sm->t + dt*imax_dsp; |
---|
1064 | sm->tpgv = fabs(amin_vel)>fabs(amax_vel)? sm->t + dt*imin_vel:sm->t + dt*imax_vel; |
---|
1065 | sm->tpga = fabs(amin_acc)>fabs(amax_acc)? sm->t + dt*imin_acc:sm->t + dt*imax_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 2*nroll (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.0*float(nroll)/(fcut*delt) * |
---|
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)**(2*n)/(1+(f/fcut)**(2*n))), * |
---|
1104 | * where n = 2*nroll * |
---|
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 (=2*nroll) * |
---|
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.0*pi*fcut); |
---|
1127 | w1 = (float)(2.0*tan(w0*delt/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.0*k + 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.0*nroll/(fcut*delt)); |
---|
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.0*x1+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.0*x1+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 = 2*pi/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.0*damp/(om2*omt)); |
---|
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 = s3*c3 - s5; b12 = -c2*s3 + 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 = a11*y + a12*ydot + b11*acc[i] + b12*acc[i+1]; |
---|
1236 | ydot = a21*y + a22*ydot + b21*acc[i] + b22*acc[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 = -d3*ydot -om2*y1; |
---|
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 | } |
---|