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 

Generated on Tue May 6 09:16:03 2003 for Earthworm Libs by doxygen1.3-rc3