Main Page Class Hierarchy Compound List File List Compound Members File Members

# k2evt2ew.c

Go to the documentation of this file.

    00001 /*
    00002  *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
    00003  *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
    00004  *
    00005  *    $Id: k2evt2ew_8c-source.html 2161 2006-05-19 16:55:03Z paulf $
    00006  *
    00007  *    Revision history:
    00008  *     $Log$
    00008  *     Revision 1.1  2006/05/19 16:55:02  paulf
    00008  *     first inclusion
    00008  *
    00009  *     Revision 1.3  2002/02/25 21:28:27  lucky
    00010  *     Fixed the translation from box to SCN
    00011  *
    00012  *     Revision 1.1  2001/11/19 16:17:59  lucky
    00013  *     Initial revision
    00014  *
    00015  *
    00016  */
    00017
    00018 /*  k2hder2ew.c
    00019  *  Read a K2-format .evt data file, crack the header, and fill
    00020  *  the return structure to be returned to the calling program.
    00021  *   This code was cloned from nsmp2ew; Credit Jim Luetgert.
    00022  */
    00023
    00024 #include <stdio.h>
    00025 #include <stdlib.h>
    00026 #include <string.h>
    00027 #include <math.h>
    00028 #include <sys/types.h>
    00029 #include <time.h>
    00030 #include <earthworm.h>
    00031 #include <chron3.h>
    00032 #include <kom.h>
    00033 #include <swap.h>
    00034 #include <k2evt2ew.h>
    00035
    00036 #define PACKET_MAX_SIZE  3003      /* 2992 + 8 + 3, same as QT */
    00037 #define GRAVITY 978.03      /* Gravity in cm/sec/sec */
    00038 #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp
    00039 #define DECADE   315532800  /* Number of seconds betweeen 1/1/1970 and 1/1/1980 */
    00040
    00041
    00042
    00043 static          int             K2EW_Debug;
    00044
    00045
    00046 /* Functions in this source file
    00047  *******************************/
    00048 static  int read_tag (FILE *fp, KFF_TAG *);
    00049 static  int read_head (FILE *fp, MW_HEADER *, int tagLength);
    00050 static  int read_frame (FILE *fp, FRAME_HEADER *, unsigned long *channels );
    00051
    00052 int extract(int pktdatalen, const unsigned char *pktbuff,
    00053             float *Data, long *Counts, float scale, int nchan, int jchan, int *ind);
    00054 int peak_ground(float *Data, int npts, int itype, float dt, SM_INFO *sm);
    00055 void demean(float *A, int N);
    00056 void locut(float *s, int nd, float fcut, float delt, int nroll, int icaus);
    00057 void rdrvaa(float *acc, int na, float omega, float damp, float dt,
    00058             float *rd, float *rv, float *aa);
    00059 void amaxper(int npts, float dt, float *fc, float *amaxmm,
    00060                                         float *aminmm, float *pmax, int *imin, int *imax);
    00061
    00062
    00063
    00064
    00065 /******************************************************************************
    00066  * k2evt2ew()   Read and crack one .evt file.
    00067  *
    00068  *  fp:                 File pointer to the .evt file to crack; file must
    00069  *              previuosly be opened for binary reads.
    00070  *  fname:      name of the file pointed to by fp (debug purposes)
    00071  *  pk2info:    pointer to a pre-allocated structure where info from
    00072  *              the file will be stored.
    00073  *  pChanName:  pointer to the CHANNELNAME structure containing mappings
    00074  *              from SM Box and channel number to SCN. This is optional;
    00075  *              if it is NULL, then default SCN names will be assigned as
    00076  *              follows: S=box, C=channel number, N=NetCode
    00077  *  numChans:   Number of channels defined in pChanName
    00078  *  NetCode:    Network Code to use if pChanName is null.
    00079  *  Debug:      Integer Debug level.
    00080  *
    00081  *  Return Values: EW_SUCCESS and EW_FAILURE, with pk2info filled with
    00082  *   other useful info.
    00083  *
    00084  ******************************************************************************/
    00085 int k2evt2ew (FILE *fp, char *fname, K2InfoStruct *pk2info,
    00086                 CHANNELNAME *pChanName, int numChans, char *NetCode, int Debug)
    00087 {
    00088         unsigned char g_k2mi_buff[PACKET_MAX_SIZE-9]; /* received data buffer   */
    00089     char     serial[10];
    00090     int      i, j, ci, chanind, nscans;
    00091     int      tlen, flen,  stat, nchans, itype[18], azim[18];
    00092     unsigned long  decade, channels;
    00093     float   dt, scale[18], sec1970;
    00094
    00095         if ((fp == NULL) || (fname == NULL) || (pk2info == NULL))
    00096         {
    00097                 logit ("e", "k2evt2ew: Invalid arguments passed in.\n");
    00098                 return EW_FAILURE;
    00099         }
    00100
    00101     if ((strstr (fname, ".EVT") == 0) && (strstr (fname, ".evt") == 0))
    00102         {
    00103                 logit ("e", "k2evt2ew: invalid file name <%s>\n", fname);
    00104                 return EW_FAILURE;
    00105         }
    00106
    00107         K2EW_Debug = Debug;
    00108         sec1970 = (float) 11676096000.00;   /* # seconds between Carl Johnson's        */
    00109                                         /* time 0 and 1970-01-01 00:00:00.0 GMT    */
    00110
    00111 /* Initialize variables
    00112  **********************/
    00113     memset ( &pk2info->frame,  0, sizeof(FRAME_HEADER) );
    00114     memset ( &pk2info->head,  0, sizeof(K2_HEADER) );
    00115     memset ( &pk2info->tag,  0, sizeof(KFF_TAG) );
    00116     memset ( &pk2info->sm,  0, MAX_SM_PER_BOX*sizeof(SM_INFO) );
    00117
    00118
    00119         /*  Kinemetric time starts 10 years after EW time  */
    00120     decade = (365*10+2)*24*60*60;
    00121
    00122 /* First, just read the header & first frame to get some basic info.
    00123  ******************************************************************/
    00124     rewind(fp);
    00125
    00126     tlen = read_tag(fp, &pk2info->tag);
    00127     stat = read_head(fp, &pk2info->head, pk2info->tag.length);
    00128
    00129     tlen = read_tag(fp, &pk2info->tag);
    00130     flen = pk2info->tag.dataLength;
    00131     read_frame(fp, &pk2info->frame, &channels);
    00132
    00133     nscans = pk2info->head.roParms.stream.duration;
    00134     nchans = pk2info->head.rwParms.misc.nchannels;
    00135     dt  = 1.0/(10.0*(flen/(3.0*nchans)));
    00136
    00137     if(pk2info->head.roParms.stream.flags != 0)
    00138         {
    00139         logit("e", "%s: not data. flags = %d \n",
    00140                     fname, pk2info->head.roParms.stream.flags);
    00141         return EW_FAILURE;
    00142     }
    00143
    00144          sprintf(serial, "%d", pk2info->head.rwParms.misc.serialNumber);
    00145
    00146         /* Find this instrument in the list */
    00147         chanind = -1;
    00148         for (ci = 0; ci < numChans; ci++)
    00149         {
    00150                 if (strcmp(serial, pChanName[ci].box) == 0)
    00151                         chanind = ci;
    00152     }
    00153
    00154     for (i=0; i<nchans; i++)
    00155         {
    00156          pk2info->sm[i].t = pk2info->head.roParms.stream.triggerTime +
    00157                      (float)pk2info->head.roParms.stream.triggerTimeMsec/1000.0;
    00158          pk2info->sm[i].t = pk2info->head.roParms.stream.startTime   +
    00159                      (float)pk2info->head.roParms.stream.startTimeMsec/1000.0;
    00160          pk2info->sm[i].t += decade;
    00161          pk2info->sm[i].talt    = (float) time(NULL);
    00162          pk2info->sm[i].altcode = SM_ALTCODE_RECEIVING_MODULE;
    00163          itype[i] = 0;
    00164          if( pk2info->head.rwParms.channel[i].sensorType==SENSOR_VELOCITY ||
    00165                         pk2info->head.rwParms.channel[i].sensorType==SENSOR_MARKL22 ||
    00166                         pk2info->head.rwParms.channel[i].sensorType==SENSOR_MARKL4C)
    00167                                 itype[i] = 2;
    00168          else
    00169          if((pk2info->head.rwParms.channel[i].sensorType>=SENSOR_FBA11 &&
    00170                         pk2info->head.rwParms.channel[i].sensorType<=SENSOR_FBA23) ||
    00171                         pk2info->head.rwParms.channel[i].sensorType==SENSOR_EPI ||
    00172                         pk2info->head.rwParms.channel[i].sensorType==SENSOR_ACCELERATION)
    00173                                 itype[i] = 3;
    00174
    00175          scale[i] = (2.0*pk2info->head.rwParms.channel[i].fullscale/(1<<24))
    00176                                                                 / pk2info->head.rwParms.channel[i].sensitivity;
    00177
    00178          if (itype[i]==3)
    00179                                 scale[i] = scale[i] * (float) GRAVITY;
    00180
    00181          if (itype[i]==0)
    00182                  {
    00183                 logit("e", "Can't decode sensor type! = %d \n",
    00184                                             pk2info->head.rwParms.channel[i].sensorType);
    00185          }
    00186          azim[i] = pk2info->head.rwParms.channel[i].azimuth;
    00187
    00188
    00189                 if (chanind < 0)
    00190                 {
    00191                         /* build default name */
    00192                         strcpy (pk2info->sm[i].sta,  pk2info->head.rwParms.misc.stnID);
    00193                         sprintf (pk2info->sm[i].comp, "%d", i);
    00194                         strcpy (pk2info->sm[i].net, NetCode);
    00195                         strcpy (pk2info->sm[i].loc,  "");
    00196                 }
    00197                 else
    00198                 {
    00199                         strcpy (pk2info->sm[i].sta,  pChanName[chanind].sta);
    00200                         strcpy (pk2info->sm[i].comp, pChanName[chanind].comp);
    00201                         strcpy (pk2info->sm[i].net,  pChanName[chanind].net);
    00202                         strcpy (pk2info->sm[i].loc,  pChanName[chanind].loc);
    00203                 }
    00204     }
    00205
    00206     if(K2EW_Debug)
    00207         logit("e", "\n|%s| |%d| |%f| |%s| |%s| \n",
    00208             pk2info->head.rwParms.misc.stnID,
    00209                         pk2info->head.rwParms.misc.serialNumber, dt,
    00210             fname, pk2info->head.rwParms.misc.comment);
    00211
    00212
    00213
    00214
    00215 /* Read and process each channel individually.
    00216  *********************************************/
    00217     for(i=0;i<nchans;i++)
    00218         {
    00219            rewind(fp);
    00220            tlen = read_tag(fp, &pk2info->tag);
    00221            stat = read_head(fp, &pk2info->head, pk2info->tag.length);
    00222
    00223                 /* initialize trace buffer */
    00224                 for (j=0; j<MAXTRACELTH; j++)
    00225                 {
    00226                         pk2info->Databuf[i][j] = 0.0;
    00227                         pk2info->Counts[i][j] = 0;
    00228                 }
    00229                 pk2info->numDataPoints[i] = 0;
    00230
    00231
    00232                 for(j=0;j<nscans;j++)
    00233                 {
    00234                tlen = read_tag(fp, &pk2info->tag);
    00235 /*
    00236 logit ("", "tlen = %d\n", tlen);
    00237                if(tlen != 16) break;
    00238 */
    00239                if(pk2info->numDataPoints[i] >= MAXTRACELTH) break;
    00240                read_frame(fp, &pk2info->frame, &channels);
    00241                flen = pk2info->tag.dataLength;
    00242                stat = fread(g_k2mi_buff, 1, flen, fp);
    00243                extract(flen, g_k2mi_buff, pk2info->Databuf[i], pk2info->Counts[i],
    00244                                                         scale[i], nchans, i, &pk2info->numDataPoints[i]);
    00245                 }
    00246
    00247             if(K2EW_Debug)
    00248                 logit("e", "%2d |%s| |%s| |%s| %3d\n",
    00249                     i, pk2info->sm[i].sta,
    00250                                         pk2info->sm[i].comp, pk2info->sm[i].net, azim[i]);
    00251
    00252             dt  = 1.0/(10.0*(flen/(3.0*nchans)));
    00253
    00254             peak_ground(pk2info->Databuf[i], pk2info->numDataPoints[i],
    00255                                                                         itype[i], dt, &pk2info->sm[i]);
    00256
    00257     } /* for loop over channels */
    00258     return EW_SUCCESS;
    00259 }
    00260
    00261
    00262
    00263 /******************************************************************************
    00264  * read_tag(fp)  Read the 16 byte tag, swap bytes if necessary, and print.    *
    00265  ******************************************************************************/
    00266 int read_tag (FILE *fp, KFF_TAG *tag)
    00267 {
    00268         if ((fp == NULL) || (tag == NULL))
    00269         {
    00270                 logit ("e", "read_tag: Invalid arguments passed in\n");
    00271                 return EW_FAILURE;
    00272         }
    00273
    00274     if (fread(tag, 1, 16, fp) != 16)
    00275         {
    00276                 logit ("e", "read_tag: read of file failed.\n");
    00277                 return EW_FAILURE;
    00278         }
    00279
    00280 #ifdef _INTEL
    00281     SwapLong (&tag->type);
    00282     SwapShort (&tag->length);
    00283     SwapShort (&tag->dataLength);
    00284     SwapShort (&tag->id);
    00285     SwapShort (&tag->checksum);
    00286 #endif _INTEL
    00287
    00288     if (K2EW_Debug > 0)
    00289         {
    00290         logit("e", "TAG: %c %d %d %d %d %d %d %d %d  \n",
    00291                 tag->sync,
    00292                 (int)tag->byteOrder,
    00293                 (int)tag->version,
    00294                 (int)tag->instrumentType,
    00295                 tag->type, tag->length, tag->dataLength,
    00296                 tag->id, tag->checksum);
    00297     }
    00298
    00299     return EW_SUCCESS;
    00300 }
    00301
    00302
    00303 /******************************************************************************
    00304  * read_head(fp)  Read the file header, swap bytes if necessary, and print.   *
    00305  ******************************************************************************/
    00306 int read_head (FILE *fp, MW_HEADER *head, int tagLength)
    00307 {
    00308    int        i, maxchans, siz;
    00309
    00310         if ((fp == NULL) || (head == NULL))
    00311         {
    00312                 logit ("e", "read_head: Invalid arguments passed in\n");
    00313                 return EW_FAILURE;
    00314         }
    00315
    00316 /* Read in the file header.
    00317    If a K2, there will be 2040 bytes,
    00318    otherwise assume a Mt Whitney.
    00319  ************************************/
    00320     maxchans = tagLength==2040? MAX_K2_CHANNELS:MAX_MW_CHANNELS;
    00321
    00322     if (fread(head, 1, 8, fp) != 8)
    00323         {
    00324                 logit ("e", "read_head: read of file failed.\n");
    00325                 return EW_FAILURE;
    00326         }
    00327
    00328         siz = sizeof(struct MISC_RO_PARMS) + sizeof(struct TIMING_RO_PARMS);
    00329     if (fread(&head->roParms.misc, 1, siz, fp) != (unsigned int) siz)
    00330         {
    00331                 logit ("e", "read_head: read of file failed.\n");
    00332                 return EW_FAILURE;
    00333         }
    00334
    00335         siz = sizeof(struct CHANNEL_RO_PARMS)*maxchans;
    00336     if (fread(&head->roParms.channel[0], 1, siz, fp) != (unsigned int) siz)
    00337         {
    00338                 logit ("e", "read_head: read of file failed.\n");
    00339                 return EW_FAILURE;
    00340         }
    00341
    00342         siz = sizeof(struct STREAM_RO_PARMS);
    00343     if (fread(&head->roParms.stream, 1, siz, fp) != (unsigned int) siz)
    00344         {
    00345                 logit ("e", "read_head: read of file failed.\n");
    00346                 return EW_FAILURE;
    00347         }
    00348
    00349         siz = sizeof(struct MISC_RW_PARMS)+sizeof(struct TIMING_RW_PARMS);
    00350     if (fread(&head->rwParms.misc, 1, siz, fp) != (unsigned int) siz)
    00351         {
    00352                 logit ("e", "read_head: read of file failed.\n");
    00353                 return EW_FAILURE;
    00354         }
    00355
    00356     siz = sizeof(struct CHANNEL_RW_PARMS)*maxchans;
    00357     if (fread(&head->rwParms.channel[0], 1, siz, fp) != (unsigned int) siz)
    00358         {
    00359                 logit ("e", "read_head: read of file failed.\n");
    00360                 return EW_FAILURE;
    00361         }
    00362
    00363     if(tagLength==2040) {
    00364         siz = sizeof(struct STREAM_K2_RW_PARMS);
    00365         if (fread(&head->rwParms.stream, 1, siz, fp) != (unsigned int) siz)
    00366                 {
    00367                         logit ("e", "read_head: read of file failed.\n");
    00368                         return EW_FAILURE;
    00369                 }
    00370     } else {
    00371         siz = sizeof(struct STREAM_MW_RW_PARMS);
    00372         if (fread(&head->rwParms.stream, 1, siz, fp) != (unsigned int) siz)
    00373                 {
    00374                         logit ("e", "read_head: read of file failed.\n");
    00375                         return EW_FAILURE;
    00376                 }
    00377     }
    00378     siz = sizeof(struct MODEM_RW_PARMS);
    00379     if (fread(&head->rwParms.modem, 1, siz, fp) != (unsigned int) siz)
    00380         {
    00381                 logit ("e", "read_head: read of file failed.\n");
    00382                 return EW_FAILURE;
    00383         }
    00384
    00385
    00386 #ifdef _INTEL
    00387         /* Byte-swap values, if necessary */
    00388
    00389         /* roParms */
    00390     SwapShort (&head->roParms.headerVersion);
    00391     SwapShort (&head->roParms.headerBytes);
    00392
    00393     SwapShort (&head->roParms.misc.installedChan);
    00394     SwapShort (&head->roParms.misc.maxChannels);
    00395     SwapShort (&head->roParms.misc.sysBlkVersion);
    00396     SwapShort (&head->roParms.misc.bootBlkVersion);
    00397     SwapShort (&head->roParms.misc.appBlkVersion);
    00398     SwapShort (&head->roParms.misc.dspBlkVersion);
    00399     SwapShort (&head->roParms.misc.batteryVoltage);
    00400     SwapShort (&head->roParms.misc.crc);
    00401     SwapShort (&head->roParms.misc.flags);
    00402     SwapShort (&head->roParms.misc.temperature);
    00403
    00404     SwapShort (&head->roParms.timing.gpsLockFailCount);
    00405     SwapShort (&head->roParms.timing.gpsUpdateRTCCount);
    00406     SwapShort (&head->roParms.timing.acqDelay);
    00407     SwapShort (&head->roParms.timing.gpsLatitude);
    00408     SwapShort (&head->roParms.timing.gpsLongitude);
    00409     SwapShort (&head->roParms.timing.gpsAltitude);
    00410     SwapShort (&head->roParms.timing.dacCount);
    00411     SwapShort (&head->roParms.timing.gpsLastDrift[0]);
    00412     SwapShort (&head->roParms.timing.gpsLastDrift[1]);
    00413     SwapLong (&head->roParms.timing.gpsLastTurnOnTime[0]);
    00414     SwapLong (&head->roParms.timing.gpsLastTurnOnTime[1]);
    00415     SwapLong (&head->roParms.timing.gpsLastUpdateTime[0]);
    00416     SwapLong (&head->roParms.timing.gpsLastUpdateTime[1]);
    00417     SwapLong (&head->roParms.timing.gpsLastLockTime[0]);
    00418     SwapLong (&head->roParms.timing.gpsLastLockTime[1]);
    00419
    00420     SwapLong (&head->roParms.stream.startTime);
    00421     SwapLong (&head->roParms.stream.triggerTime);
    00422     SwapLong (&head->roParms.stream.duration);
    00423     SwapShort (&head->roParms.stream.errors);
    00424     SwapShort (&head->roParms.stream.flags);
    00425     SwapShort (&head->roParms.stream.startTimeMsec);
    00426     SwapShort (&head->roParms.stream.triggerTimeMsec);
    00427     SwapLong (&head->roParms.stream.nscans);
    00428
    00429
    00430     for(i=0;i<head->roParms.misc.maxChannels;i++)
    00431         {
    00432         SwapLong (&head->roParms.channel[i].maxPeak);
    00433         SwapLong (&head->roParms.channel[i].maxPeakOffset);
    00434         SwapLong (&head->roParms.channel[i].minPeak);
    00435         SwapLong (&head->roParms.channel[i].minPeakOffset);
    00436         SwapLong (&head->roParms.channel[i].mean);
    00437         }
    00438
    00439
    00440         /* rwParams */
    00441     SwapShort (&head->rwParms.misc.serialNumber);
    00442     SwapShort (&head->rwParms.misc.nchannels);
    00443     SwapShort (&head->rwParms.misc.elevation);
    00444     SwapFloat (&head->rwParms.misc.latitude);
    00445     SwapFloat (&head->rwParms.misc.longitude);
    00446     SwapShort (&head->rwParms.misc.userCodes[0]);
    00447     SwapShort (&head->rwParms.misc.userCodes[1]);
    00448     SwapShort (&head->rwParms.misc.userCodes[2]);
    00449     SwapShort (&head->rwParms.misc.userCodes[3]);
    00450     SwapLong (&head->rwParms.misc.cutler_bitmap);
    00451     SwapLong (&head->rwParms.misc.channel_bitmap);
    00452
    00453     SwapShort (&head->rwParms.timing.localOffset);
    00454
    00455     for(i=0;i<head->rwParms.misc.nchannels;i++)
    00456         {
    00457         SwapShort (&head->rwParms.channel[i].sensorSerialNumberExt);
    00458         SwapShort (&head->rwParms.channel[i].north);
    00459         SwapShort (&head->rwParms.channel[i].east);
    00460         SwapShort (&head->rwParms.channel[i].up);
    00461         SwapShort (&head->rwParms.channel[i].altitude);
    00462         SwapShort (&head->rwParms.channel[i].azimuth);
    00463         SwapShort (&head->rwParms.channel[i].sensorType);
    00464         SwapShort (&head->rwParms.channel[i].sensorSerialNumber);
    00465         SwapShort (&head->rwParms.channel[i].gain);
    00466         SwapShort (&head->rwParms.channel[i].StaLtaRatio);
    00467         SwapFloat (&head->rwParms.channel[i].fullscale);
    00468         SwapFloat (&head->rwParms.channel[i].sensitivity);
    00469         SwapFloat (&head->rwParms.channel[i].damping);
    00470         SwapFloat (&head->rwParms.channel[i].naturalFrequency);
    00471         SwapFloat (&head->rwParms.channel[i].triggerThreshold);
    00472         SwapFloat (&head->rwParms.channel[i].detriggerThreshold);
    00473         SwapFloat (&head->rwParms.channel[i].alarmTriggerThreshold);
    00474     }
    00475
    00476     SwapShort (&head->rwParms.stream.eventNumber);
    00477     SwapShort (&head->rwParms.stream.sps);
    00478     SwapShort (&head->rwParms.stream.apw);
    00479     SwapShort (&head->rwParms.stream.preEvent);
    00480     SwapShort (&head->rwParms.stream.postEvent);
    00481     SwapShort (&head->rwParms.stream.minRunTime);
    00482     SwapShort (&head->rwParms.stream.VotesToTrigger);
    00483     SwapShort (&head->rwParms.stream.VotesToDetrigger);
    00484     SwapShort (&head->rwParms.stream.Timeout);
    00485     SwapShort (&head->rwParms.stream.TxBlkSize);
    00486     SwapShort (&head->rwParms.stream.BufferSize);
    00487     SwapShort (&head->rwParms.stream.SampleRate);
    00488     SwapLong (&head->rwParms.stream.TxChanMap);
    00489 #endif _INTEL
    00490
    00491
    00492     if(K2EW_Debug > 0)
    00493         {
    00494             logit("e", "HEADER: %c%c%c %d %hu %hu \n",
    00495                  head->roParms.id[0], head->roParms.id[1], head->roParms.id[2],
    00496                  (int)head->roParms.instrumentCode,
    00497              head->roParms.headerVersion,
    00498              head->roParms.headerBytes);
    00499
    00500             logit("e", "MISC_RO_PARMS: %d %d %d %hu %hu %hu %hu %hu %hu %d %hu %hu %d \n",
    00501             (int)head->roParms.misc.a2dBits,
    00502             (int)head->roParms.misc.sampleBytes,
    00503             (int)head->roParms.misc.restartSource,
    00504             head->roParms.misc.installedChan,
    00505             head->roParms.misc.maxChannels,
    00506             head->roParms.misc.sysBlkVersion,
    00507             head->roParms.misc.bootBlkVersion,
    00508             head->roParms.misc.appBlkVersion,
    00509             head->roParms.misc.dspBlkVersion,
    00510             head->roParms.misc.batteryVoltage,
    00511             head->roParms.misc.crc,
    00512             head->roParms.misc.flags,
    00513             head->roParms.misc.temperature );
    00514
    00515
    00516             logit("e", "TIMING_RO_PARMS: %d %d %d %hu %hu %d %d %d %d %hu %d %d %lu %lu %lu %lu %lu %lu \n",
    00517             (int)head->roParms.timing.clockSource,
    00518             (int)head->roParms.timing.gpsStatus,
    00519             (int)head->roParms.timing.gpsSOH,
    00520             head->roParms.timing.gpsLockFailCount,
    00521             head->roParms.timing.gpsUpdateRTCCount,
    00522             head->roParms.timing.acqDelay,
    00523             head->roParms.timing.gpsLatitude,
    00524             head->roParms.timing.gpsLongitude,
    00525             head->roParms.timing.gpsAltitude,
    00526             head->roParms.timing.dacCount,
    00527             head->roParms.timing.gpsLastDrift[0],
    00528             head->roParms.timing.gpsLastDrift[1],
    00529             head->roParms.timing.gpsLastTurnOnTime[0],
    00530             head->roParms.timing.gpsLastTurnOnTime[1],
    00531             head->roParms.timing.gpsLastUpdateTime[0],
    00532             head->roParms.timing.gpsLastUpdateTime[1],
    00533             head->roParms.timing.gpsLastLockTime[0],
    00534             head->roParms.timing.gpsLastLockTime[1] );
    00535
    00536
    00537             for(i=0;i<head->roParms.misc.maxChannels;i++)
    00538                 {
    00539                 logit("e", "CHANNEL_RO_PARMS: %d %lu %d %lu %d %d \n",
    00540                     head->roParms.channel[i].maxPeak,
    00541                 head->roParms.channel[i].maxPeakOffset,
    00542                     head->roParms.channel[i].minPeak,
    00543                 head->roParms.channel[i].minPeakOffset,
    00544                 head->roParms.channel[i].mean,
    00545                 head->roParms.channel[i].aqOffset );
    00546         }
    00547
    00548         logit("e", "STREAM_RO_PARMS: %d %d %d %d %d %d %d %d \n",
    00549             head->roParms.stream.startTime,
    00550             head->roParms.stream.triggerTime,
    00551             head->roParms.stream.duration,
    00552             head->roParms.stream.errors,
    00553             head->roParms.stream.flags,
    00554             head->roParms.stream.startTimeMsec,
    00555             head->roParms.stream.triggerTimeMsec,
    00556             head->roParms.stream.nscans  );
    00557
    00558
    00559
    00560             logit("e", "MISC_RW_PARMS: %hu %hu %.5s %.33s %d %f %f %d %d %d %d %d %lo %d %.17s %d %d\n",
    00561             head->rwParms.misc.serialNumber,
    00562             head->rwParms.misc.nchannels,
    00563             head->rwParms.misc.stnID,
    00564             head->rwParms.misc.comment,
    00565             head->rwParms.misc.elevation,
    00566             head->rwParms.misc.latitude,
    00567             head->rwParms.misc.longitude,
    00568            (int)head->rwParms.misc.cutlerCode,
    00569            (int)head->rwParms.misc.minBatteryVoltage,
    00570            (int)head->rwParms.misc.cutler_decimation,
    00571            (int)head->rwParms.misc.cutler_irig_type,
    00572             head->rwParms.misc.cutler_bitmap,
    00573             head->rwParms.misc.channel_bitmap,
    00574            (int)head->rwParms.misc.cutler_protocol,
    00575             head->rwParms.misc.siteID,
    00576            (int)head->rwParms.misc.externalTrigger,
    00577            (int)head->rwParms.misc.networkFlag );
    00578
    00579             logit("e", "TIMING_RW_PARMS: %d %d %hu \n",
    00580            (int)head->rwParms.timing.gpsTurnOnInterval,
    00581            (int)head->rwParms.timing.gpsMaxTurnOnTime,
    00582                head->rwParms.timing.localOffset  );
    00583
    00584             for(i=0;i<head->rwParms.misc.nchannels;i++)
    00585                 {
    00586                 logit("e",
    00587                                 "CHANNEL_RW_PARMS: %s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %f %f %f %f %f %f %f \n",
    00588                    head->rwParms.channel[i].id,
    00589                    head->rwParms.channel[i].sensorSerialNumberExt,
    00590                    head->rwParms.channel[i].north,
    00591                    head->rwParms.channel[i].east,
    00592                    head->rwParms.channel[i].up,
    00593                    head->rwParms.channel[i].altitude,
    00594                    head->rwParms.channel[i].azimuth,
    00595                    head->rwParms.channel[i].sensorType,
    00596                    head->rwParms.channel[i].sensorSerialNumber,
    00597                    head->rwParms.channel[i].gain,
    00598                (int)head->rwParms.channel[i].triggerType,
    00599                (int)head->rwParms.channel[i].iirTriggerFilter,
    00600                (int)head->rwParms.channel[i].StaSeconds,
    00601                (int)head->rwParms.channel[i].LtaSeconds,
    00602                    head->rwParms.channel[i].StaLtaRatio,
    00603                (int)head->rwParms.channel[i].StaLtaPercent,
    00604                    head->rwParms.channel[i].fullscale,
    00605                    head->rwParms.channel[i].sensitivity,
    00606                    head->rwParms.channel[i].damping,
    00607                    head->rwParms.channel[i].naturalFrequency,
    00608                    head->rwParms.channel[i].triggerThreshold,
    00609                    head->rwParms.channel[i].detriggerThreshold,
    00610                    head->rwParms.channel[i].alarmTriggerThreshold);
    00611                 }
    00612
    00613             logit("e", "STREAM_K2_RW_PARMS: |%d| |%d| |%d| %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lu\n",
    00614             (int)head->rwParms.stream.filterFlag,
    00615             (int)head->rwParms.stream.primaryStorage,
    00616             (int)head->rwParms.stream.secondaryStorage,
    00617                 head->rwParms.stream.eventNumber,
    00618                 head->rwParms.stream.sps,
    00619                 head->rwParms.stream.apw,
    00620                 head->rwParms.stream.preEvent,
    00621                 head->rwParms.stream.postEvent,
    00622                 head->rwParms.stream.minRunTime,
    00623                 head->rwParms.stream.VotesToTrigger,
    00624                 head->rwParms.stream.VotesToDetrigger,
    00625             (int)head->rwParms.stream.FilterType,
    00626             (int)head->rwParms.stream.DataFmt,
    00627                 head->rwParms.stream.Timeout,
    00628                 head->rwParms.stream.TxBlkSize,
    00629                 head->rwParms.stream.BufferSize,
    00630                 head->rwParms.stream.SampleRate,
    00631                 head->rwParms.stream.TxChanMap);
    00632
    00633         } /* If debug */
    00634
    00635     return EW_SUCCESS;
    00636
    00637 } /* read_head */
    00638
    00639
    00640 /******************************************************************************
    00641  * read_frame(fp)  Read the frame header, swap bytes if necessary.            *
    00642  ******************************************************************************/
    00643 int read_frame (FILE *fp, FRAME_HEADER *frame, unsigned long *channels )
    00644 {
    00645     unsigned short   frameStatus, frameStatus2, samprate, streamnumber;
    00646     unsigned char    BitMap[4];
    00647     unsigned long    bmap;
    00648
    00649         if ((fp == NULL) || (frame == NULL) || (channels == NULL))
    00650         {
    00651                 logit ("e", "read_frame: invalid arguments passed in.\n");
    00652                 return EW_FAILURE;
    00653         }
    00654
    00655
    00656     if (fread(frame, 32, 1, fp) != 1)
    00657         {
    00658                 logit ("e", "read_frame: read of file failed.\n");
    00659                 return EW_FAILURE;
    00660         }
    00661
    00662
    00663 #ifdef _INTEL
    00664     SwapShort (&frame->recorderID);
    00665     SwapShort (&frame->frameSize);
    00666     SwapShort (&frame->blockTime);
    00667     SwapShort (&frame->blockTime2);
    00668     SwapShort (&frame->channelBitMap);
    00669     SwapShort (&frame->streamPar);
    00670     SwapShort (&frame->msec);
    00671 #endif _INTEL
    00672
    00673     BitMap[0] = frame->channelBitMap & 255;
    00674     BitMap[1] = frame->channelBitMap >> 8;
    00675     BitMap[2] = frame->channelBitMap1;
    00676     BitMap[3] = 0;
    00677
    00678     bmap = *((long *)(BitMap));
    00679     frameStatus      = frame->frameStatus;
    00680     frameStatus2     = frame->frameStatus2;
    00681     samprate         = frame->streamPar & 4095;
    00682     streamnumber     = frame->streamPar >> 12;
    00683
    00684     if(K2EW_Debug > 0)
    00685             logit("e", "FRAME: %d %d %d %d   %lu X%ho   %hu X%ho %hu %hu     X%ho X%ho %hu X%ho \n",
    00686             (int)frame->frameType,
    00687             (int)frame->instrumentCode,
    00688             frame->recorderID,
    00689             frame->frameSize,
    00690             frame->blockTime,
    00691             frame->channelBitMap,
    00692             frame->streamPar, streamnumber, samprate, samprate>>8,
    00693             frameStatus,
    00694             frameStatus2,
    00695             frame->msec,
    00696             (int)frame->channelBitMap1);
    00697
    00698     *channels = bmap;
    00699
    00700     return EW_SUCCESS;
    00701 }
    00702
    00703
    00704 /******************************************************************************
    00705  * extract: Processes k2 file buffers to expand data from 24 bits to          *
    00706  *         36 bits, sort channels and scale to cgs units.                     *
    00707  *         pktdatalen - length of data in buffer                              *
    00708  *         pktbuff    - address of buffer containing the data bytes           *
    00709  *         Data       - the destination array.                                *
    00710  *         scale      - the scaling factor (units/count).                     *
    00711  *         nchan      - the total number of channels in the buffer.           *
    00712  *         jchan      - the channel to be extracted.                          *
    00713  *         ind        - the current index in the destination array.           *
    00714  ******************************************************************************/
    00715
    00716 int extract(int pktdatalen, const unsigned char *pktbuff, float *Data,
    00717                                         long *Counts, float scale, int nchan, int jchan, int *ind)
    00718 {
    00719     int     k, pktpos, ichan;
    00720     unsigned char btarr[4];
    00721
    00722   /* process and store serial data stream values */
    00723     k = *ind;
    00724     pktpos = ichan = 0;
    00725     while (pktpos < pktdatalen) {
    00726         /* loop through each serial data stream sample data value */
    00727         if(ichan == jchan) {
    00728             btarr[1] = pktbuff[pktpos++];    /* copy MSByte of 24-bit value */
    00729             btarr[2] = pktbuff[pktpos++];    /* copy middle-byte of 24-bit value */
    00730             btarr[3] = pktbuff[pktpos++];    /* copy LSByte of 24-bit value */
    00731             /* sign extend to 32-bits */
    00732
    00733             btarr[0] = ((btarr[1] & (unsigned char)0x80) ==
    00734                       (unsigned char)0) ? (unsigned char)0 : (unsigned char)0xFF;
    00735
    00736             /* byte-swap and enter value */
    00737 #ifdef _INTEL
    00738                         SwapLong ((long *) btarr);
    00739 #endif
    00740             Counts[k] = (*((long *)(btarr)));
    00741             Data[k] = Counts[k]*scale;
    00742                         k = k + 1;
    00743         } else {
    00744             pktpos += 3;
    00745         }
    00746         ichan += 1;
    00747         if(ichan >= nchan) ichan = 0;
    00748     }
    00749
    00750     *ind = k;
    00751     return 0;
    00752 }
    00753
    00754 /******************************************************************************
    00755  *   subroutine for estimation of ground motion                               *
    00756  *                                                                            *
    00757  *   input:                                                                   *
    00758  *           Data   - data array                                              *
    00759  *           npts   - number of points in timeseries                          *
    00760  *           itype  - units for timeseries. 1=disp 2=vel 3=acc                *
    00761  *           dt     - sample spacing in sec                                   *
    00762  *   return:                                                                  *
    00763  *           0      - All OK                                                  *
    00764  *           1      - error                                                   *
    00765  *                                                                            *
    00766  ******************************************************************************/
    00767 int peak_ground(float *Data, int npts, int itype, float dt, SM_INFO *sm)
    00768 {
    00769     int     imax_acc, imax_vel, imax_dsp, imax_raw;
    00770     int     imin_acc, imin_vel, imin_dsp, imin_raw;
    00771     int     ii, kk, kpts, id[4], npd[4][2], icaus;
    00772     float  totint, a, tpi, omega, damp, rd, rv, aa;
    00773     float  amax_acc, amin_acc, pmax_acc;
    00774     float  amax_vel, amin_vel, pmax_vel;
    00775     float  amax_dsp, amin_dsp, pmax_dsp;
    00776     float  amax_raw, amin_raw, pmax_raw, gd[4], sd[4];
    00777         float    d1[MAXTRACELTH];
    00778         float    d2[MAXTRACELTH];
    00779         float    d3[MAXTRACELTH];
    00780
    00781     gd[0] = 0.05;
    00782     gd[1] = 0.10;
    00783     gd[2] = 0.20;
    00784     gd[3] = 0.50;
    00785     icaus = 1;
    00786
    00787     tpi  = 8.0*atan(1.0);
    00788
    00789 /* Find the raw maximum and its period
    00790  *************************************/
    00791
    00792     demean(Data, npts);
    00793     amaxper(npts, dt, Data, &amin_raw, &amax_raw, &pmax_raw, &imin_raw, &imax_raw);
    00794
    00795     if(itype == 1) {  /* input data is displacement  */
    00796         for(kk=0;kk<npts;kk++) d1[kk] = Data[kk];
    00797         locut(d1, npts, 0.17, dt, 2, icaus);
    00798         for(kk=1;kk<npts;kk++) {
    00799             d2[kk] = (d1[kk] - d1[kk-1])/dt;
    00800         }
    00801         d2[0] = d2[1];
    00802         demean(d2, npts);
    00803         for(kk=1;kk<npts;kk++) {
    00804             d3[kk] = (d2[kk] - d2[kk-1])/dt;
    00805         }
    00806         d3[0] = d3[1];
    00807         demean(d3, npts);
    00808     } else
    00809     if(itype == 2) {  /* input data is velocity      */
    00810         for(kk=0;kk<npts;kk++) d2[kk] = Data[kk];
    00811         locut(d2, npts, 0.17, dt, 2, icaus);
    00812         for(kk=1;kk<npts;kk++) {
    00813             d3[kk] = (d2[kk] - d2[kk-1])/dt;
    00814         }
    00815         d3[0] = d3[1];
    00816         demean(d3, npts);
    00817
    00818         totint = 0.0;
    00819         for(kk=0;kk<npts-1;kk++) {
    00820             totint = totint + (d2[kk] + d2[kk+1])*0.5*dt;
    00821             d1[kk] = totint;
    00822         }
    00823         d1[npts-1] = d1[npts-2];
    00824         demean(d1, npts);
    00825
    00826     } else
    00827     if(itype == 3) {  /* input data is acceleration  */
    00828         for(kk=0;kk<npts;kk++) d3[kk] = Data[kk];
    00829         locut(d3, npts, 0.17, dt, 2, icaus);
    00830
    00831         totint = 0.0;
    00832         for(kk=0;kk<npts-1;kk++) {
    00833             totint = totint + (d3[kk] + d3[kk+1])*0.5*dt;
    00834             d2[kk] = totint;
    00835         }
    00836         d2[npts-1] = d2[npts-2];
    00837         demean(d2, npts);
    00838
    00839         totint = 0.0;
    00840         for(kk=0;kk<npts-1;kk++) {
    00841             totint = totint + (d2[kk] + d2[kk+1])*0.5*dt;
    00842             d1[kk] = totint;
    00843         }
    00844         d1[npts-1] = d1[npts-2];
    00845         demean(d1, npts);
    00846     } else {
    00847         return 1;
    00848     }
    00849
    00850 /* Find the displacement(cm), velocity(cm/s), & acceleration(cm/s/s) maxima  and their periods
    00851  *********************************************************************************************/
    00852
    00853     amaxper(npts, dt, &d1[0], &amin_dsp, &amax_dsp, &pmax_dsp, &imin_dsp, &imax_dsp);
    00854     amaxper(npts, dt, &d2[0], &amin_vel, &amax_vel, &pmax_vel, &imin_vel, &imax_vel);
    00855     amaxper(npts, dt, &d3[0], &amin_acc, &amax_acc, &pmax_acc, &imin_acc, &imax_acc);
    00856
    00857 /* Find the spectral response
    00858  ****************************/
    00859
    00860     damp = 0.05;
    00861     kk = 0;
    00862     sm->pdrsa[kk] = 0.3;
    00863     omega = tpi/sm->pdrsa[kk];
    00864     rdrvaa(&d3[0], npts-1, omega, damp, dt, &rd, &rv, &aa);
    00865     sm->rsa[kk] = aa;
    00866     kk += 1;
    00867
    00868     sm->pdrsa[kk] = 1.0;
    00869     omega = tpi/sm->pdrsa[kk];
    00870     rdrvaa(&d3[0], npts-1, omega, damp, dt, &rd, &rv, &aa);
    00871     sm->rsa[kk] = aa;
    00872     kk += 1;
    00873
    00874     sm->pdrsa[kk] = 3.0;
    00875     omega = tpi/sm->pdrsa[kk];
    00876     rdrvaa(&d3[0], npts-1, omega, damp, dt, &rd, &rv, &aa);
    00877     sm->rsa[kk] = aa;
    00878     kk += 1;
    00879
    00880     sm->nrsa = kk;
    00881
    00882 /* Since we are here, determine the duration of strong shaking
    00883  *************************************************************/
    00884
    00885     for(kk=0;kk<4;kk++) {
    00886         id[kk] = npd[kk][1] = npd[kk][2] = 0;
    00887         for(ii=1;ii<=npts-1;ii++) {
    00888             a = fabs(d3[ii]/GRAVITY);
    00889             if (a >= gd[kk]) {
    00890                 id[kk] = id[kk] + 1;
    00891                 if (id[kk] == 1) npd[kk][1] = ii;
    00892                 npd[kk][2] = ii;
    00893             }
    00894         }
    00895         if (id[kk] != 0) {
    00896             kpts = npd[kk][2] - npd[kk][1] + 1;
    00897             sd[kk] = kpts*dt;
    00898         } else {
    00899             sd[kk] = 0.0;
    00900         }
    00901     }
    00902
    00903     sm->pgd = fabs(amin_dsp)>fabs(amax_dsp)? fabs(amin_dsp):fabs(amax_dsp);
    00904     sm->pgv = fabs(amin_vel)>fabs(amax_vel)? fabs(amin_vel):fabs(amax_vel);
    00905     sm->pga = fabs(amin_acc)>fabs(amax_acc)? fabs(amin_acc):fabs(amax_acc);
    00906
    00907     sm->tpgd = fabs(amin_dsp)>fabs(amax_dsp)? sm->t + dt*imin_dsp:sm->t + dt*imax_dsp;
    00908     sm->tpgv = fabs(amin_vel)>fabs(amax_vel)? sm->t + dt*imin_vel:sm->t + dt*imax_vel;
    00909     sm->tpga = fabs(amin_acc)>fabs(amax_acc)? sm->t + dt*imin_acc:sm->t + dt*imax_acc;
    00910
    00911     return 0;
    00912 }
    00913
    00914
    00915 /******************************************************************************
    00916  *  demean removes the mean from the n point series stored in array A.        *
    00917  *                                                                            *
    00918  ******************************************************************************/
    00919 void demean(float *A, int n)
    00920 {
    00921     int       i;
    00922     float    xm;
    00923
    00924     xm = 0.0;
    00925     for(i=0;i<n;i++) xm = xm + A[i];
    00926     xm = xm/n;
    00927     for(i=0;i<n;i++) A[i] = A[i] - xm;
    00928 }
    00929
    00930
    00931 /******************************************************************************
    00932  *  Butterworth locut filter order 2*nroll (nroll<=8)                         *
    00933  *   (see Kanasewich, Time Sequence Analysis in Geophysics,                   *
    00934  *   Third Edition, University of Alberta Press, 1981)                        *
    00935  *  written by W. B. Joyner 01/07/97                                          *
    00936  *                                                                            *
    00937  *  s[j] input = the time series to be filtered                               *
    00938  *      output = the filtered series                                          *
    00939  *      dimension of s[j] must be at least as large as                        *
    00940  *        nd+3.0*float(nroll)/(fcut*delt)                                     *
    00941  *  nd    = the number of points in the time series                           *
    00942  *  fcut  = the cutoff frequency                                              *
    00943  *  delt  = the timestep                                                      *
    00944  *  nroll = filter order                                                      *
    00945  *  causal if icaus.eq.1 - zero phase shift otherwise                         *
    00946  *                                                                            *
    00947  * The response is given by eq. 15.8-6 in Kanasewich:                         *
    00948  *  Y = sqrt((f/fcut)**(2*n)/(1+(f/fcut)**(2*n))),                            *
    00949  *                 where n = 2*nroll                                          *
    00950  *                                                                            *
    00951  * Dates: 01/07/97 - Written by Bill Joyner                                   *
    00952  *        12/17/99 - D. Boore added check for fcut = 0.0, in which case       *
    00953  *                   no filter is applied.  He also cleaned up the            *
    00954  *                   appearance of the code (indented statements in           *
    00955  *                   loops, etc.)                                             *
    00956  *        02/04/00 - Changed "n" to "nroll" to eliminate confusion with       *
    00957  *                   Kanesewich, who uses "n" as the order (=2*nroll)         *
    00958  *        03/01/00 - Ported to C by Jim Luetgert                              *
    00959  *                                                                            *
    00960  ******************************************************************************/
    00961 void locut(float *s, int nd, float fcut, float delt, int nroll, int icaus)
    00962 {
    00963     float    fact[8], b1[8], b2[8];
    00964     float    pi, w0, w1, w2, w3, w4, w5, xp, yp, x1, x2, y1, y2;
    00965     int       j, k, np2, npad;
    00966
    00967     if (fcut == 0.0) return;       /* Added by DMB  */
    00968
    00969     pi = 4.0*atan(1.0);
    00970     w0 = 2.0*pi*fcut;
    00971     w1 = 2.0*tan(w0*delt/2.0);
    00972     w2 = (w1/2.0)*(w1/2.0);
    00973     w3 = (w1*w1)/2.0 - 2.0;
    00974     w4 = 0.25*pi/nroll;
    00975
    00976     for(k=0;k<nroll;k++) {
    00977         w5 = w4*(2.0*k + 1.0);
    00978         fact[k] = 1.0/(1.0+sin(w5)*w1 + w2);
    00979         b1[k] = w3*fact[k];
    00980         b2[k] = (1.0-sin(w5)*w1 + w2)*fact[k];
    00981     }
    00982
    00983     np2 = nd;
    00984
    00985     if(icaus != 1) {
    00986         npad = 3.0*nroll/(fcut*delt);
    00987         np2 = nd+npad;
    00988         for(j=nd;j<np2;j++) s[j] = 0.0;
    00989     }
    00990
    00991     for(k=0;k<nroll;k++) {
    00992         x1 = x2 = y1 = y2 = 0.0;
    00993         for(j=0;j<np2;j++) {
    00994             xp = s[j];
    00995             yp = fact[k]*(xp-2.0*x1+x2) - b1[k]*y1 - b2[k]*y2;
    00996             s[j] = yp;
    00997             y2 = y1;
    00998             y1 = yp;
    00999             x2 = x1;
    01000             x1 = xp;
    01001         }
    01002     }
    01003
    01004     if(icaus != 1) {
    01005         for(k=0;k<nroll;k++) {
    01006             x1 = x2 = y1 = y2 = 0.0;
    01007             for(j=0;j<np2;j++) {
    01008                 xp = s[np2-j-1];
    01009                 yp = fact[k]*(xp-2.0*x1+x2) - b1[k]*y1 - b2[k]*y2;
    01010                 s[np2-j-1] = yp;
    01011                 y2 = y1;
    01012                 y1 = yp;
    01013                 x2 = x1;
    01014                 x1 = xp;
    01015             }
    01016         }
    01017     }
    01018
    01019     return;
    01020 }
    01021
    01022
    01023 /******************************************************************************
    01024  * rdrvaa                                                                     *
    01025  *                                                                            *
    01026  * This is a modified version of "Quake.For", originally                      *
    01027  * written by J.M. Roesset in 1971 and modified by                            *
    01028  * Stavros A. Anagnostopoulos, Oct. 1986.  The formulation is                 *
    01029  * that of Nigam and Jennings (BSSA, v. 59, 909-922, 1969).                   *
    01030  * Dates: 02/11/00 - Modified by David M. Boore, based on RD_CALC             *
    01031  *        03/01/00 - Ported to C by Jim Luetgert                              *
    01032  *                                                                            *
    01033  *   acc = acceleration time series                                           *
    01034  *    na = length of time series                                              *
    01035  * omega = 2*pi/per                                                           *
    01036  *  damp = fractional damping (e.g., 0.05)                                    *
    01037  *    dt = time spacing of input                                              *
    01038  *    rd = relative displacement of oscillator                                *
    01039  *    rv = relative velocity of oscillator                                    *
    01040  *    aa = absolute acceleration of oscillator                                *
    01041  ******************************************************************************/
    01042 void rdrvaa(float *acc, int na, float omega, float damp, float dt,
    01043             float *rd, float *rv, float *aa)
    01044 {
    01045     float    omt, d2, bom, d3, omd, om2, omdt, c1, c2, c3, c4, cc, ee;
    01046     float    s1, s2, s3, s4, s5, a11, a12, a21, a22, b11, b12, b21, b22;
    01047     float    y, ydot, y1, z, z1, z2, ra;
    01048     int       i;
    01049
    01050     omt  = omega*dt;
    01051     d2   = sqrt(1.0-damp*damp);
    01052     bom  = damp*omega;
    01053     d3   = 2.0*bom;
    01054     omd  = omega*d2;
    01055     om2  = omega*omega;
    01056     omdt = omd*dt;
    01057     c1 = 1.0/om2;
    01058     c2 = 2.0*damp/(om2*omt);
    01059     c3 = c1+c2;
    01060     c4 = 1.0/(omega*omt);
    01061     ee = exp(-damp*omt);
    01062     cc = cos(omdt)*ee;
    01063     s1 = sin(omdt)*ee/omd;
    01064     s2 = s1*bom;
    01065     s3 = s2 + cc;
    01066     s4 = c4*(1.0-s3);
    01067     s5 = s1*c4 + c2;
    01068
    01069     a11 =  s3;          a12 = s1;
    01070     a21 = -om2*s1;      a22 = cc - s2;
    01071
    01072     b11 =  s3*c3 - s5;  b12 = -c2*s3 + s5 - c1;
    01073     b21 = -s1+s4;       b22 = -s4;
    01074
    01075     y = ydot = *rd = *rv = *aa = 0.0;
    01076     for(i=0;i<na-1;i++) {
    01077         y1   = a11*y + a12*ydot + b11*acc[i] + b12*acc[i+1];
    01078         ydot = a21*y + a22*ydot + b21*acc[i] + b22*acc[i+1];
    01079         y = y1;    /* y is the oscillator output at time corresponding to index i   */
    01080         z = fabs(y);
    01081         if (z > *rd) *rd = z;
    01082         z1 = fabs(ydot);
    01083         if (z1 > *rv) *rv = z1;
    01084         ra = -d3*ydot -om2*y1;
    01085         z2 = fabs(ra);
    01086         if (z2 > *aa) *aa = z2;
    01087     }
    01088 }
    01089
    01090 /******************************************************************************
    01091  *   compute maximum amplitude and its associated period                      *
    01092  *                                                                            *
    01093  *   input:                                                                   *
    01094  *           npts   - number of points in timeseries                          *
    01095  *           dt     - sample spacing in sec                                   *
    01096  *           fc     - input timeseries                                        *
    01097  *   output:                                                                  *
    01098  *           amaxmm - raw maximum                                             *
    01099  *           pmax   - period of maximum                                       *
    01100  *           imax   - index of maxmimum point                                 *
    01101  *                                                                            *
    01102  ******************************************************************************/
    01103 void amaxper(int npts, float dt, float *fc, float *aminmm, float *amaxmm,
    01104                        float *pmax, int *imin, int *imax)
    01105 {
    01106     float    amin, amax, pp, pm, mean, frac;
    01107     int       i, j, jmin, jmax;
    01108
    01109     *imax = jmax = *imin = jmin = 0;
    01110     amax = amin = *amaxmm = *aminmm = fc[0];
    01111     *aminmm = *pmax = mean = 0.0;
    01112     for(i=0;i<npts;i++) {
    01113         mean = mean + fc[i]/npts;
    01114         if (fc[i] > amax) { jmax = i; amax = fc[i]; }
    01115         if (fc[i] < amin) { jmin = i; amin = fc[i]; }
    01116     }
    01117
    01118 /*     compute period of maximum    */
    01119
    01120     pp = pm = 0.0;
    01121     if (fc[jmax] > mean) {
    01122         j = jmax+1;
    01123         while(fc[j] > mean && j < npts) {
    01124             pp += dt;
    01125             j  += 1;
    01126         }
    01127         frac = dt*(mean-fc[j-1])/(fc[j]-fc[j-1]);
    01128         frac = 0.0;
    01129         pp = pp + frac;
    01130         j = jmax-1;
    01131         while(fc[j] > mean && j >= 0) {
    01132             pm += dt;
    01133             j  -= 1;
    01134         }
    01135         frac = dt*(mean-fc[j+1])/(fc[j]-fc[j+1]);
    01136         frac = 0.0;
    01137         pm = pm + frac;
    01138     } else {
    01139         j = jmax+1;
    01140         if(fc[j] < mean && j < npts) {
    01141             pp += dt;
    01142             j  += 1;
    01143         }
    01144         frac = dt*(mean-fc[j-1])/(fc[j]-fc[j-1]);
    01145         frac = 0.0;
    01146         pp = pp + frac;
    01147         j = jmax-1;
    01148         if(fc[j] < mean && j >= 0) {
    01149             pm += dt;
    01150             j  -= 1;
    01151         }
    01152         frac = dt*(mean-fc[j+1])/(fc[j]-fc[j+1]);
    01153         frac = 0.0;
    01154         pm = pm + frac;
    01155     }
    01156
    01157     *imin = jmin;
    01158     *imax = jmax;
    01159     *pmax = 2.0*(pm+pp);
    01160     *aminmm = amin;
    01161     *amaxmm = amax;
    01162
    01163     return;
    01164 }
    01165