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