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

Revision 7211, 49.6 KB checked in by baker, 9 months 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
83static          int             K2EW_Debug;
84
85/* Functions in this source file
86 *******************************/
87static  int read_tag (FILE *fp, KFF_TAG *);
88static  int read_head (FILE *fp, MW_HEADER *, int tagLength);
89static  int read_frame (FILE *fp, FRAME_HEADER *, uint32_t *channels );
90
91int 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);
94int peak_ground(float *Data, int npts, int itype, float dt, SM_INFO *sm);
95void demean(float *A, int N);
96void locut(float *s, int nd, float fcut, float delt, int nroll, int icaus);
97void rdrvaa(float *acc, int na, float omega, float damp, float dt,
98            float *rd, float *rv, float *aa, int *maxtime);
99void 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
125int 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
354int 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
435int 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
772int 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
850int 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
915int 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
1075void 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
1117void 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
1199void 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
1265void 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.