00001 /* 00002 * THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE 00003 * CHECKED IT OUT USING THE COMMAND CHECKOUT. 00004 * 00005 * $Id: sac__2__ewevent_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.6 2002/03/22 20:02:38 lucky 00010 * Pre v6.1 checkin 00011 * 00012 * Revision 1.5 2001/07/27 20:41:58 lucky 00013 * *** empty log message *** 00014 * 00015 * Revision 1.4 2001/07/01 21:55:23 davidk 00016 * Cleanup of the Earthworm Database API and the applications that utilize it. 00017 * The "ewdb_api" was cleanup in preparation for Earthworm v6.0, which is 00018 * supposed to contain an API that will be backwards compatible in the 00019 * future. Functions were removed, parameters were changed, syntax was 00020 * rewritten, and other stuff was done to try to get the code to follow a 00021 * general format, so that it would be easier to read. 00022 * 00023 * Applications were modified to handle changed API calls and structures. 00024 * They were also modified to be compiler warning free on WinNT. 00025 * 00026 * Revision 1.3 2001/06/21 21:25:00 lucky 00027 * State of the code after LocalMag review was implemented and partially tested. 00028 * 00029 * Revision 1.2 2001/06/06 20:54:46 lucky 00030 * Changes made to support multitude magnitudes, as well as amplitude picks. This is w 00031 * in progress - checkin in for the sake of safety. 00032 * 00033 * Revision 1.1 2001/05/15 02:15:27 davidk 00034 * Initial revision 00035 * 00036 * Revision 1.5 2001/04/17 17:30:00 davidk 00037 * changed a bad subscript in a logging statement, due to a compiler warning on NT. 00038 * 00039 * Revision 1.4 2001/04/16 18:31:41 lucky 00040 * Added a check for a valid coda length before assining the SAC header 00041 * value to a db structure entry. 00042 * 00043 * Revision 1.3 2001/01/19 18:04:23 lucky 00044 * Added a check for mag label "DurCoda" as well as "Md" before looking 00045 * for a coda length 00046 * 00047 * Revision 1.2 2001/01/08 19:05:06 lucky 00048 * Fixed a bug in ProduceSAC_NextStationForEvent whereby only one 00049 * pick out of many would be written to the SAC header. 00050 * 00051 * Revision 1.1 2000/12/18 18:55:12 lucky 00052 * Initial revision 00053 * 00054 * Revision 1.12 2000/12/06 17:50:47 lucky 00055 * We now correctly keep track of the pick onset 00056 * 00057 * Revision 1.11 2000/11/08 18:46:33 lucky 00058 * Modified SAC_NextStationforEvent: let's stop pretending - we only 00059 * produce sac for the first Waveform of any SCN. 00060 * 00061 * Revision 1.10 2000/09/13 15:36:51 lucky 00062 * Fixed FirstMotion character handling; 00063 * Fixed handling of the coda length 00064 * 00065 * Revision 1.9 2000/09/07 21:17:45 lucky 00066 * Final version after the Review pages have been demonstrated. 00067 * 00068 * Revision 1.6 2000/08/31 16:10:24 lucky 00069 * Before mallocing snippet buffer, we need to calculate correct iByteLen field. 00070 * 00071 * Revision 1.5 2000/08/30 18:42:51 lucky 00072 * Added malloc of binSnippet before data is copied from the sacdata buffer 00073 * to the EWEvent structure. 00074 * 00075 * Revision 1.4 2000/08/30 17:41:57 lucky 00076 * InitEWEvent has been changed to include optional allocation of 00077 * pChanInfo space. This must be optional because GetEWEventInfo mallocs 00078 * space on the fly. 00079 * 00080 * Revision 1.3 2000/08/25 18:24:45 lucky 00081 * Fixed to work with S and P arrivals, as well as the new, explicitly 00082 * allocated pChanInfo. 00083 * 00084 * Revision 1.2 2000/08/19 00:10:12 lucky 00085 * Made SAC_reftime and SAC_Ktostr non-static -- needed in other modules 00086 * 00087 * Revision 1.1 2000/08/16 17:18:11 lucky 00088 * Initial revision 00089 * 00090 * Revision 1.13 2000/06/29 17:24:42 lucky 00091 * Added to Sac2EWEvent the translation of the coda duration from the sac header. 00092 * The coda duration is written to EW structs. 00093 * 00094 * Revision 1.12 2000/06/28 22:53:07 lucky 00095 * Increased the maximum number of Phases per event to 1000 to accomodate 00096 * extremely large Dewey-type events. 00097 * If BuildPhs fails, we record the failure, but do not exit. This way, 00098 * even if one arrival writing fails we still get the remaining ones 00099 * in the arc message. 00100 * 00101 * Revision 1.11 2000/06/28 22:22:34 lucky 00102 * Added a check in EWEvent2ArcMsg so that only the arrival-type channels 00103 * are written to Arc message. Recall:channels may be arrivals, station 00104 * magnitudes, or even just waveform data. 00105 * 00106 * Revision 1.10 2000/06/26 18:14:42 lucky 00107 * Fixed creation of hypoinverse arc files with NSN locations. 00108 * 00109 * Revision 1.9 2000/06/22 23:29:11 lucky 00110 * Fixed BuildPhs to work with USNSN messages' 00111 * 00112 * Revision 1.8 2000/06/20 21:42:51 lucky 00113 * *** empty log message *** 00114 * 00115 * Revision 1.7 2000/06/16 17:23:36 lucky 00116 * Set tObsPhase to be same as tCalcPhase when read from sac headers. 00117 * 00118 * Revision 1.6 2000/06/13 19:05:58 lucky 00119 * Fixed to work with the new hypo2ora, as well as the review stuff. 00120 * 00121 * Revision 1.5 2000/06/07 22:09:52 lucky 00122 * Added ArcMsg2EWEvent calls 00123 * 00124 * 00125 * Revision 1.1 2000/02/14 18:51:48 lucky 00126 * Initial revision 00127 * 00128 * 00129 */ 00130 00131 #include <earthworm.h> 00132 #include <sac_2_ewevent.h> 00133 00134 extern time_t timegm_ew (struct tm *); 00135 00136 double SAC_reftime (struct SAChead *); 00137 size_t SAC_Ktostr (char *, int, char *); 00138 00139 00140 /* Global Variable */ 00141 float WriteSAC_dGapThresh=1.5; 00142 00143 00144 /******************************************************** 00145 * 00146 * Sac2EWEvent: 00147 * Fills the EWEventInfoStruct structure pointed 00148 * to by pEvent, with the data from the SAC header 00149 * pointed to by sacheadp, with data for one 00150 * channel of trace data. 00151 * 00152 * Returns EW_SUCCESS if everything went OK, 00153 * EW_FAILURE otherwise. 00154 * 00155 * 00156 * Author: Lucky Vidmar 04/2000 00157 * 00158 ********************************************************/ 00159 int Sac2EWEvent (EWEventInfoStruct *pEvent, struct SAChead *sacheadp, 00160 char *sacdatap, int CopyData) 00161 { 00162 00163 double tref; 00164 int i, tmp, toalloc; 00165 int wt; 00166 EWDB_OriginStruct *pOrigin; 00167 EWDB_ArrivalStruct *pArrival; 00168 EWDB_StationMagStruct *pStaMag; 00169 EWDB_WaveformStruct *pWaveform; 00170 EWDB_StationStruct *pStation; 00171 EWChannelDataStruct *pChan; 00172 EWChannelDataStruct *pOldChan; 00173 00174 00175 00176 if ((pEvent == NULL) || (sacheadp == NULL)) 00177 { 00178 logit ("", "Invalid parameters passed in.\n"); 00179 return EW_FAILURE; 00180 } 00181 00182 00183 /* Check how many channels we processed */ 00184 if (pEvent->iNumChans >= pEvent->iNumAllocChans) 00185 { 00186 toalloc = (int) (pEvent->iNumAllocChans + CHAN_ALLOC_INCR); 00187 00188 logit ("", "Max Channels (%d) reached. Alloc to %d\n", 00189 pEvent->iNumAllocChans, toalloc); 00190 00191 pOldChan = pEvent->pChanInfo; 00192 if ((pEvent->pChanInfo = (EWChannelDataStruct *) malloc 00193 (toalloc * sizeof (EWChannelDataStruct))) == NULL) 00194 { 00195 logit ("", "Could not malloc %d pChanInfo entries.\n", toalloc); 00196 return (EW_FAILURE); 00197 } 00198 00199 memcpy (pEvent->pChanInfo, pOldChan, 00200 (pEvent->iNumAllocChans * sizeof (EWChannelDataStruct))); 00201 00202 for (tmp = pEvent->iNumChans; tmp < toalloc; tmp++) 00203 InitEWChan(&pEvent->pChanInfo[tmp]); 00204 00205 pEvent->iNumAllocChans = toalloc; 00206 00207 free (pOldChan); 00208 00209 } 00210 00211 00212 /* Origin structure 00213 *********************/ 00214 if ((pOrigin = &(pEvent->PrefOrigin)) == NULL) 00215 { 00216 logit ("", "Invalid pEvent -- can't get PrefOrigin.\n"); 00217 return EW_FAILURE; 00218 } 00219 00220 00221 /* Make sure that we have origin time 00222 **************************************/ 00223 if (sacheadp->o == (float) SACUNDEF) 00224 { 00225 logit ("", "No origin time in SAC header.\n"); 00226 return EW_FAILURE; 00227 } 00228 00229 /* Absolute event origin time 00230 ******************************/ 00231 if ((tref = SAC_reftime (sacheadp)) < 0) 00232 { 00233 logit ("", "Call to SAC_reftime failed.\n"); 00234 return EW_FAILURE; 00235 } 00236 00237 pOrigin->tOrigin = (double) sacheadp->o + tref; 00238 00239 /* Fill in the event location 00240 *****************************/ 00241 00242 /* Lattitude */ 00243 if (sacheadp->evla != (float) SACUNDEF) 00244 pOrigin->dLat = sacheadp->evla; 00245 00246 /* Longitude */ 00247 if (sacheadp->evlo != (float) SACUNDEF) 00248 pOrigin->dLon = sacheadp->evlo; 00249 00250 /* Depth */ 00251 if (sacheadp->evdp != (float) SACUNDEF) 00252 pOrigin->dDepth = sacheadp->evdp; 00253 00254 00255 00256 00257 /******* Channel Data *******/ 00258 00259 if ((pChan = &(pEvent->pChanInfo[pEvent->iNumChans])) == NULL) 00260 { 00261 logit ("", "Invalid pEvent -- can't get pChanInfo[%d].\n", pEvent->iNumChans); 00262 return EW_FAILURE; 00263 } 00264 00265 pChan->iNumArrivals = 0; 00266 pChan->iNumWaveforms = 1; 00267 00268 00269 /******* Station struct *******/ 00270 00271 if ((pStation = &(pChan->Station)) == NULL) 00272 { 00273 logit ("", "Invalid EWChannel -- can't get Station.\n"); 00274 return EW_FAILURE; 00275 } 00276 00277 SAC_Ktostr (pStation->Sta, 7, sacheadp->kstnm); 00278 SAC_Ktostr (pStation->Comp, 9, sacheadp->kcmpnm); 00279 SAC_Ktostr (pStation->Net, 9, sacheadp->knetwk); 00280 00281 /* Station Latitude */ 00282 if (sacheadp->stla != (float) SACUNDEF) 00283 pStation->Lat = sacheadp->stla; 00284 00285 /* Station Longitude */ 00286 if (sacheadp->stlo != (float) SACUNDEF) 00287 pStation->Lon = sacheadp->stlo; 00288 00289 /* Station Elevation */ 00290 if (sacheadp->stel != (float) SACUNDEF) 00291 pStation->Elev = sacheadp->stel; 00292 00293 00294 /* P-wave information */ 00295 if (sacheadp->a != (float) SACUNDEF) 00296 { 00297 i = pChan->iNumArrivals; 00298 00299 00300 /******* Arrival struct *******/ 00301 if ((pArrival = &(pChan->Arrivals[i])) == NULL) 00302 { 00303 logit ("", "Invalid EWChannel; can't get P-wave Arrivals.\n"); 00304 return EW_FAILURE; 00305 } 00306 00307 pArrival->pStation = pStation; 00308 00309 /* Arrival time */ 00310 if ((tref = SAC_reftime (sacheadp)) < 0) 00311 { 00312 logit ("", "Call to SAC_reftime failed.\n"); 00313 return EW_FAILURE; 00314 } 00315 00316 pArrival->tCalcPhase = (double) sacheadp->a + tref; 00317 pArrival->tObsPhase = (double) sacheadp->a + tref; 00318 pArrival->szObsPhase[0] = sacheadp->ka[0]; 00319 pArrival->szObsPhase[1] = sacheadp->ka[1]; 00320 pArrival->szObsPhase[2] = sacheadp->ka[2]; 00321 pArrival->szObsPhase[3] = '\0'; 00322 pArrival->cMotion = sacheadp->ka[1]; 00323 pArrival->cOnset = sacheadp->ka[3]; 00324 00325 /* Convert weight code into dSigma */ 00326 wt = sacheadp->ka[2] - '0'; 00327 00328 if (wt == 3) 00329 pArrival->dSigma = 0.08; 00330 else if (wt == 2) 00331 pArrival->dSigma = 0.05; 00332 else if (wt == 1) 00333 pArrival->dSigma = 0.03; 00334 else if (wt == 0) 00335 pArrival->dSigma = 0.02; 00336 else 00337 pArrival->dSigma = 99.99; 00338 00339 00340 pArrival->dDist = sacheadp->dist; 00341 pArrival->dAzm = sacheadp->az; 00342 00343 /******* StaMag struct *******/ 00344 00345 if ((pStaMag = &(pChan->Stamags[i])) == NULL) 00346 { 00347 logit ("", "Invalid EWChannel -- can't get StaMags.\n"); 00348 return EW_FAILURE; 00349 } 00350 00351 if (sacheadp->f != (float) SACUNDEF) 00352 { 00353 pStaMag->StaMagUnion.CodaDur.tCodaDurXtp = sacheadp->f - pArrival->tCalcPhase + tref; 00354 pStaMag->StaMagUnion.CodaDur.tCodaTermXtp = pArrival->tCalcPhase + 00355 pStaMag->StaMagUnion.CodaDur.tCodaDurXtp; 00356 } 00357 else 00358 { 00359 pStaMag->StaMagUnion.CodaDur.tCodaDurXtp = 0.0; 00360 pStaMag->StaMagUnion.CodaDur.tCodaTermXtp = pArrival->tCalcPhase; 00361 } 00362 00363 pStaMag->iMagType = MAGTYPE_DURATION; 00364 00365 pChan->iNumArrivals = pChan->iNumArrivals + 1; 00366 pChan->iNumStaMags = pChan->iNumStaMags + 1; 00367 } 00368 00369 00370 /* S-wave information */ 00371 if (sacheadp->t0 != (float) SACUNDEF) 00372 { 00373 i = pChan->iNumArrivals; 00374 00375 00376 /******* Arrival struct *******/ 00377 if ((pArrival = &(pChan->Arrivals[i])) == NULL) 00378 { 00379 logit ("", "Invalid EWChannel; can't get S-wave Arrivals.\n"); 00380 return EW_FAILURE; 00381 } 00382 00383 pArrival->pStation = pStation; 00384 00385 /* Arrival time */ 00386 if ((tref = SAC_reftime (sacheadp)) < 0) 00387 { 00388 logit ("", "Call to SAC_reftime failed.\n"); 00389 return EW_FAILURE; 00390 } 00391 00392 pArrival->tCalcPhase = (double) sacheadp->t0 + tref; 00393 pArrival->tObsPhase = (double) sacheadp->t0 + tref; 00394 pArrival->szObsPhase[0] = sacheadp->kt0[0]; 00395 pArrival->szObsPhase[1] = sacheadp->kt0[1]; 00396 pArrival->szObsPhase[2] = sacheadp->kt0[2]; 00397 pArrival->szObsPhase[3] = '\0'; 00398 pArrival->cMotion = sacheadp->kt0[1]; 00399 pArrival->cOnset = sacheadp->kt0[3]; 00400 00401 /* Convert weight code into dSigma */ 00402 wt = sacheadp->kt0[2] - '0'; 00403 00404 if (wt == 3) 00405 pArrival->dSigma = 0.08; 00406 else if (wt == 2) 00407 pArrival->dSigma = 0.05; 00408 else if (wt == 1) 00409 pArrival->dSigma = 0.03; 00410 else if (wt == 0) 00411 pArrival->dSigma = 0.02; 00412 else 00413 pArrival->dSigma = 99.99; 00414 00415 00416 pArrival->dDist = sacheadp->dist; 00417 pArrival->dAzm = sacheadp->az; 00418 00419 00420 /******* StaMag struct *******/ 00421 00422 if ((pStaMag = &(pChan->Stamags[i])) == NULL) 00423 { 00424 logit ("", "Invalid EWChannel -- can't get StaMags.\n"); 00425 return EW_FAILURE; 00426 } 00427 00428 if (sacheadp->f != (float) SACUNDEF) 00429 { 00430 pStaMag->StaMagUnion.CodaDur.tCodaDurXtp = sacheadp->f - pArrival->tCalcPhase + tref; 00431 pStaMag->StaMagUnion.CodaDur.tCodaTermXtp = pArrival->tCalcPhase + 00432 pStaMag->StaMagUnion.CodaDur.tCodaDurXtp; 00433 } 00434 else 00435 { 00436 pStaMag->StaMagUnion.CodaDur.tCodaDurXtp = 0.0; 00437 pStaMag->StaMagUnion.CodaDur.tCodaTermXtp = 0.0; 00438 } 00439 00440 pStaMag->iMagType = MAGTYPE_DURATION; 00441 00442 pChan->iNumArrivals = pChan->iNumArrivals + 1; 00443 pChan->iNumStaMags = pChan->iNumStaMags + 1; 00444 00445 } 00446 00447 00448 if (pChan->iNumArrivals == 0) 00449 { 00450 logit ("", "Warning: No picks found in SAC header.\n"); 00451 } 00452 00453 00454 /******* Waveform struct *******/ 00455 if (CopyData == 1) 00456 { 00457 if ((pWaveform = &(pChan->Waveforms[0])) == NULL) 00458 { 00459 logit ("", "Invalid EWChannel -- can't get Waveform.\n"); 00460 return EW_FAILURE; 00461 } 00462 00463 /* Allocate space of the snippet data */ 00464 pWaveform->iByteLen = sacheadp->npts * sizeof (float); 00465 if ((pWaveform->binSnippet = malloc (pWaveform->iByteLen)) == NULL) 00466 { 00467 logit ("", "Could not malloc binSnippet.\n"); 00468 return EW_FAILURE; 00469 } 00470 00471 memcpy (pWaveform->binSnippet, sacdatap, ((size_t) sacheadp->npts * sizeof (float))); 00472 } 00473 00474 00475 return EW_SUCCESS; 00476 } 00477 00478 00479 00480 00481 00482 00483 /******************************************************** 00484 * 00485 * SacHeaderInit: 00486 * Initialize the SAC header structure pointed to 00487 * by sacheadp -- be sure to malloc it before 00488 * calling this function. 00489 * 00490 * Also sets default (constant) values in the header 00491 * 00492 * Returns EW_SUCCESS if everything went OK, 00493 * EW_FAILURE otherwise. 00494 * 00495 * 00496 * Author: Lucky Vidmar 04/2000 00497 * 00498 ********************************************************/ 00499 int SacHeaderInit (struct SAChead *sacheadp) 00500 { 00501 00502 int i; 00503 struct SAChead2 *head2; /* use a simple structure here - we */ 00504 /* don't care what * the variables are - */ 00505 /* set them to 'undefined' */ 00506 if (sacheadp == NULL) 00507 { 00508 logit ("", "Invalid parameters passed in.\n"); 00509 return EW_FAILURE; 00510 } 00511 00512 /* change to a simpler format */ 00513 head2 = (struct SAChead2 *) sacheadp; 00514 00515 /* set all of the floats to 'undefined' */ 00516 for (i = 0; i < NUM_FLOAT; ++i) 00517 head2->SACfloat[i] = SACUNDEF; 00518 00519 /* set all of the ints to 'undefined' */ 00520 for (i = 0; i < MAXINT-5; ++i) 00521 head2->SACint[i] = SACUNDEF; 00522 00523 /* except for the logical integers - set them to 1 */ 00524 for ( ; i < MAXINT; ++i) 00525 head2->SACint[i] = 1; 00526 00527 /* set all of the strings to 'undefined' */ 00528 for (i = 0; i < MAXSTRING; ++i) 00529 (void) strncpy(head2->SACstring[i], SACSTRUNDEF,K_LEN); 00530 00531 /* SAC I.D. number */ 00532 head2->SACfloat[9] = SAC_I_D; 00533 00534 /* header version number */ 00535 head2->SACint[6] = SACVERSION; 00536 00537 00538 /* Set default values that remain constant */ 00539 sacheadp->b = 0; /* beginning time relative to reference time */ 00540 sacheadp->iftype = SAC_ITIME; /* File type is time series */ 00541 sacheadp->leven = 1; /* evenly spaced data */ 00542 sacheadp->idep = SAC_IUNKN; /* unknown independent data type */ 00543 sacheadp->iztype = SAC_IBEGINTIME; /* Reference time is Begin time */ 00544 sacheadp->ievtyp = SAC_IUNKN; /* event type */ 00545 00546 /* Label strings */ 00547 strncpy(sacheadp->ko, "Origin ", K_LEN); 00548 strncpy(sacheadp->ka, "Arrival ", K_LEN); 00549 strncpy(sacheadp->kt0, "S-Wave ", K_LEN); 00550 00551 return EW_SUCCESS; 00552 00553 } 00554 00555 00556 00557 /****************************************************** 00558 * SAC_reftime() reads the reference time in the SAC * 00559 * header and returns it in double seconds since * 00560 * January 1, 1970 00:00:00 GMT * 00561 * Returns: -1 on error * 00562 ******************************************************/ 00563 double SAC_reftime( struct SAChead *psac ) 00564 { 00565 struct tm stm; 00566 time_t tt; 00567 double tref; 00568 00569 /* Find seconds from 1970.01.01 to Jan 1 of the correct year 00570 ***********************************************************/ 00571 stm.tm_year = (int) psac->nzyear - 1900l; /* tm_year = yrs since 1900 */ 00572 stm.tm_mon = 0; /* January */ 00573 stm.tm_mday = 1; /* 1st */ 00574 stm.tm_hour = (int) psac->nzhour; 00575 stm.tm_min = (int) psac->nzmin; 00576 stm.tm_sec = (int) psac->nzsec; 00577 stm.tm_isdst = 0; 00578 stm.tm_wday = 0; /* not used by timegm */ 00579 stm.tm_yday = 0; /* not used by timegm */ 00580 00581 tt = timegm_ew (&stm); 00582 00583 if (tt < 0) 00584 { 00585 logit ("", "Call to timegm_ew failed.\n"); 00586 return (double) tt; 00587 } 00588 00589 /* Then add the number of seconds for the day-of-year plus milliseconds. 00590 ***********************************************************************/ 00591 tref = (double)tt + /* sec since 1970 to Jan 1 this year */ 00592 86400.*(psac->nzjday-1) + /* sec since Jan 1 in this year */ 00593 (double)psac->nzmsec/1000.; /* milliseconds */ 00594 00595 return( tref ); 00596 } 00597 00598 00599 /**************************************************************** 00600 * SAC_Ktostr() convert a SAC header K-field (which is filled * 00601 * with white-space) into a null-terminated character string * 00602 * Returns the length of the new string. * 00603 ****************************************************************/ 00604 size_t SAC_Ktostr( char *s1, int len1, char *s2 ) 00605 { 00606 int i; 00607 char tmp[K_LEN+1]; 00608 00609 /* NULL-fill the target 00610 **********************/ 00611 for( i=0; i<len1; i++ ) s1[i]='\0'; 00612 00613 /* Is the K-field NULL ? 00614 ***********************/ 00615 strncpy( tmp, s2, K_LEN ); 00616 tmp[K_LEN] = '\0'; 00617 if( strcmp( tmp, SACSTRUNDEF ) == 0 ) return( 0 ); 00618 00619 /* Null terminate after last non-space character 00620 ***********************************************/ 00621 for( i=K_LEN-1; i>=0; i-- ) /* start at the end of the string */ 00622 { 00623 if( tmp[i]!=' ' && tmp[i]!='\0' ) /* find last non-space char */ 00624 { 00625 tmp[i+1]='\0'; /* null-terminate after last non-space char */ 00626 break; 00627 } 00628 } 00629 if( i<0 ) return( 0 ); /* K-field was empty */ 00630 00631 strncpy( s1, tmp, len1-1 ); 00632 return( strlen( s1 ) ); 00633 } 00634 00635 00636 00637 /************************************************************************* 00638 * WriteSAC_init() * 00639 * Initializes the SAC writing environment. Uses the SACPABase * 00640 * routines for writing SAC data * 00641 *************************************************************************/ 00642 int WriteSAC_init(char * szOutDir,char * szOutputFormat, int bDebug) 00643 { 00644 /* Initialize SAC Putaway Base environment */ 00645 00646 if(bDebug) 00647 SACPABase_Debug(TRUE); 00648 return(SACPABase_init(MAX_SNIPPET_LEN+SACHEADERSIZE, szOutDir,TRUE,szOutputFormat)); 00649 } 00650 00651 00652 /************************************************************************* 00653 * WriteSAC_shutdown() * 00654 * Shuts down the SAC writing environment. Uses the SACPABase * 00655 * routines for writing SAC data * 00656 *************************************************************************/ 00657 int WriteSAC_shutdown(void) 00658 { 00659 /* Close SAC Putaway Base environment */ 00660 00661 return(SACPABase_close()); 00662 } 00663 00664 /************************************************************************* 00665 * WriteSAC_StartEvent () * 00666 * Initializes a SAC event, and records data that is related to * 00667 * the entire event. * 00668 * Uses the SACPABase routines for writing SAC data. * 00669 *************************************************************************/ 00670 int WriteSAC_Event_BAD (EWEventInfoStruct * pEvent, int WriteHtml) 00671 { 00672 int rc,i; 00673 00674 rc=WriteSAC_StartEvent(pEvent); 00675 if(rc != EW_SUCCESS) 00676 { 00677 logit("","WriteSAC_StartEvent() failed with rc=%d for idEvent %d\n", 00678 rc,pEvent->Event.idEvent); 00679 return(EW_FAILURE); 00680 } 00681 00682 for(i=0;i<pEvent->iNumChans;i++) 00683 { 00684 00685 if (WriteHtml == 1) 00686 { 00687 logit ("t", "STARTING CHANNEL %d of %d\n", 00688 i, pEvent->iNumChans); 00689 printf (" "); 00690 } 00691 00692 rc=ProduceSAC_NextStationForEvent(&(pEvent->pChanInfo[i])); 00693 if(rc != EW_SUCCESS) 00694 { 00695 logit("","ProduceSAC_NextStationForEvent() failed with rc=%d " 00696 "for chan record %d of %d.\n", 00697 rc,i,pEvent->iNumChans); 00698 } 00699 00700 00701 rc=WriteSAC_NextStationForEvent(&(pEvent->pChanInfo[i])); 00702 if(rc != EW_SUCCESS) 00703 { 00704 logit("","WriteSAC_NextStationForEvent() failed with rc=%d " 00705 "for chan record %d of %d.\n", 00706 rc,i,pEvent->iNumChans); 00707 } 00708 } /* end for i < iNumChans */ 00709 00710 return(WriteSAC_EndEvent()); 00711 } 00712 00713 /************************************************************************* 00714 * WriteSAC_StartEvent () * 00715 * Initializes a SAC event, and records data that is related to * 00716 * the entire event. * 00717 * Uses the SACPABase routines for writing SAC data. * 00718 *************************************************************************/ 00719 int WriteSAC_StartEvent(EWEventInfoStruct * pEvent) 00720 { 00721 char szEventID[20]; 00722 SAC_OriginStruct SAC_Origin; 00723 00724 sprintf(szEventID,"%d",pEvent->Event.idEvent); 00725 00726 SAC_Origin.dElev=pEvent->PrefOrigin.dDepth; 00727 SAC_Origin.dLat=pEvent->PrefOrigin.dLat; 00728 SAC_Origin.dLon=pEvent->PrefOrigin.dLon; 00729 SAC_Origin.tOrigin=pEvent->PrefOrigin.tOrigin; 00730 00731 return(SACPABase_next_ev(szEventID,SAC_Origin.tOrigin,&SAC_Origin)); 00732 } 00733 00734 /************************************************************************* 00735 * WriteSAC_EndEvent () * 00736 * Marks the end of a SAC event. * 00737 * Uses the SACPABase routines for writing SAC data. * 00738 *************************************************************************/ 00739 int WriteSAC_EndEvent(void) 00740 { 00741 return(SACPABase_end_ev()); 00742 } 00743 00744 /************************************************************************* 00745 * WriteSAC_StartEvent () * 00746 * Initializes a SAC event, and records data that is related to * 00747 * the entire event. * 00748 * Uses the SACPABase routines for writing SAC data. * 00749 *************************************************************************/ 00750 int ProduceSAC_NextStationForEvent(EWChannelDataStruct * pChannel) 00751 { 00752 int i,j,rc; 00753 TRACE_REQ TR; 00754 SAC_ArrivalStruct AS; 00755 SAC_StationStruct SS; 00756 SAC_ResponseStruct RS; 00757 EWDB_ArrivalStruct * pArrival; 00758 EWDB_WaveformStruct * pWaveform; 00759 /* SAC_AmpPickStruct AmpPick; 00760 EWDB_StationMagStruct * pStaMag; */ 00761 00762 00763 EWDB_StationStruct * pStation=&(pChannel->Station); 00764 00765 00766 if(pChannel->iNumWaveforms < 1) 00767 { 00768 logit("","ProduceSAC_NextStationForEvent() No waveforms for %s,%s,%s! Not much" 00769 " point in creating a SAC file for that station!\n", 00770 pStation->Sta,pStation->Comp,pStation->Net); 00771 return(EW_FAILURE); 00772 } 00773 00774 00775 pWaveform = &(pChannel->Waveforms[0]); 00776 if(pWaveform->iByteLen <= 0) 00777 { 00778 logit("","ProduceSAC_NextStationForEvent() Waveform for " 00779 "%s,%s,%s is 0 bytes long. \nNot much" 00780 " point in creating a SAC file for that station!\n", 00781 pStation->Sta,pStation->Comp,pStation->Net); 00782 return(EW_FAILURE); 00783 } 00784 00785 SACPABase_next_scn(pStation->Sta,pStation->Comp,pStation->Net); 00786 00787 strcpy(TR.sta,pStation->Sta); 00788 strcpy(TR.chan,pStation->Comp); 00789 strcpy(TR.net,pStation->Net); 00790 00791 TR.reqStarttime = TR.actStarttime = pWaveform->tStart; 00792 TR.reqEndtime = TR.actEndtime = pWaveform->tEnd; 00793 TR.pBuf=pWaveform->binSnippet; 00794 TR.bufLen = TR.actLen = pWaveform->iByteLen; 00795 00796 /* Note we are not setting the samprate parameter!!!!!! 00797 Hopefully this does not create a problem. We don't 00798 want to have to go through all of the header byte 00799 swapping here. davidk 000303 00800 Maybe we could just get it from pResponse 00801 *******************************************************/ 00802 00803 rc=SACPABase_write_trace(&TR,WriteSAC_dGapThresh); 00804 00805 if(rc != EW_SUCCESS) 00806 { 00807 logit("","SACPABase_write_trace() failed with rc=%d for %s,%s,%s\n", 00808 rc,pStation->Sta,pStation->Comp,pStation->Net); 00809 return(EW_FAILURE); 00810 } 00811 00812 /* Write arrival picks */ 00813 for (i = 0; i < pChannel->iNumArrivals; i++) 00814 { 00815 if ((pArrival = &pChannel->Arrivals[i]) == NULL) 00816 { 00817 logit ("", "Invalid EWChannel; Can't get P-wave Arrival.\n"); 00818 return (EW_FAILURE); 00819 } 00820 00821 if ((pArrival->szObsPhase[0] == 'P') || 00822 (pArrival->szObsPhase[0] == 'p')) 00823 { 00824 AS.tPhase=pArrival->tObsPhase; 00825 AS.cPhase='P'; 00826 AS.dDist=pArrival->dDist; 00827 AS.dAzm=pArrival->dAzm; 00828 AS.cFMotion=pArrival->cMotion; 00829 AS.cOnset=pArrival->cOnset; 00830 00831 if(pArrival->dSigma <=0.02) 00832 AS.iPhaseWt=0; 00833 else if(pArrival->dSigma <=0.03) 00834 AS.iPhaseWt=1; 00835 else if(pArrival->dSigma <=0.05) 00836 AS.iPhaseWt=2; 00837 else if(pArrival->dSigma <=0.08) 00838 AS.iPhaseWt=3; 00839 else 00840 AS.iPhaseWt=4; 00841 00842 00843 /* Find the coda length */ 00844 AS.dCodaLen=0; 00845 for (j = 0; j < pChannel->iNumStaMags && !AS.dCodaLen; j++) 00846 { 00847 if (pChannel->Stamags[j].iMagType == MAGTYPE_DURATION) 00848 { 00849 AS.dCodaLen = pChannel->Stamags[j].StaMagUnion.CodaDur.tCodaDurXtp; 00850 } 00851 } 00852 00853 rc=SACPABase_write_parametric(&AS, PWAVE); 00854 if(rc != EW_SUCCESS) 00855 { 00856 logit("","SACPABase_write_parametric (PWAVE) \n" 00857 "failed with rc=%d for %s,%s,%s\n", 00858 rc,pStation->Sta,pStation->Comp,pStation->Net); 00859 return(EW_FAILURE); 00860 } 00861 00862 } /* Processing Pwave */ 00863 00864 else if ((pArrival->szObsPhase[0] == 'S') || 00865 (pArrival->szObsPhase[0] == 's')) 00866 { 00867 AS.tPhase=pArrival->tObsPhase; 00868 AS.cPhase='S'; 00869 AS.dDist=pArrival->dDist; 00870 AS.dAzm=pArrival->dAzm; 00871 AS.cFMotion=pArrival->cMotion; 00872 AS.cOnset=pArrival->cOnset; 00873 00874 if(pArrival->dSigma <=0.02) 00875 AS.iPhaseWt=0; 00876 else if(pArrival->dSigma <=0.03) 00877 AS.iPhaseWt=1; 00878 else if(pArrival->dSigma <=0.05) 00879 AS.iPhaseWt=2; 00880 else if(pArrival->dSigma <=0.08) 00881 AS.iPhaseWt=3; 00882 else 00883 AS.iPhaseWt=4; 00884 00885 /* Find the coda length */ 00886 AS.dCodaLen=0; 00887 for (j = 0; j < pChannel->iNumStaMags && !AS.dCodaLen; j++) 00888 { 00889 if (pChannel->Stamags[j].iMagType == MAGTYPE_DURATION) 00890 { 00891 AS.dCodaLen = pChannel->Stamags[j].StaMagUnion.CodaDur.tCodaDurXtp; 00892 } 00893 } 00894 00895 rc=SACPABase_write_parametric(&AS, SWAVE); 00896 if(rc != EW_SUCCESS) 00897 { 00898 logit("","SACPABase_write_parametric (SWAVE) \n" 00899 "failed with rc=%d for %s,%s,%s\n", 00900 rc,pStation->Sta,pStation->Comp,pStation->Net); 00901 return(EW_FAILURE); 00902 } 00903 00904 } /* Processing Swave */ 00905 00906 else 00907 { 00908 logit ("", "Unknown arrival type ***%c***; continuing.\n", 00909 pArrival->szObsPhase[0]); 00910 } 00911 00912 } /* loop over arrivals */ 00913 00914 SS.dLat=pStation->Lat; 00915 SS.dLon=pStation->Lon; 00916 SS.dElev=pStation->Elev; 00917 SS.bResponseIsValid=pChannel->bResponseIsValid; 00918 00919 if(pChannel->bResponseIsValid) 00920 { 00921 00922 if ((pChannel->ResponseInfo.tfsFunc.iNumPoles > 0) || 00923 (pChannel->ResponseInfo.tfsFunc.iNumZeroes > 0)) 00924 { 00925 SS.pResponse=&RS; 00926 RS.dGain=pChannel->ResponseInfo.dGain; 00927 RS.iNumPoles=pChannel->ResponseInfo.tfsFunc.iNumPoles; 00928 RS.iNumZeroes=pChannel->ResponseInfo.tfsFunc.iNumZeroes; 00929 /* This assumes that SAC_PZNum and EWDBNum are of the 00930 same size and layout 00931 *****************************************/ 00932 memcpy(&(RS.Poles),&(pChannel->ResponseInfo.tfsFunc.Poles), 00933 RS.iNumPoles * sizeof(RS.Poles[0])); 00934 memcpy(&(RS.Zeroes),&(pChannel->ResponseInfo.tfsFunc.Zeroes), 00935 RS.iNumZeroes * sizeof(RS.Zeroes[0])); 00936 } 00937 } 00938 00939 rc=SACPABase_write_stainfo(&SS); 00940 if(rc != EW_SUCCESS) 00941 { 00942 logit("","ProduceSAC_NextStationForEvent():SACPABase_write_stainfo() \n" 00943 "failed with rc=%d for %s,%s,%s\n", 00944 rc,pStation->Sta,pStation->Comp,pStation->Net); 00945 return(EW_FAILURE); 00946 } 00947 00948 00949 return(EW_SUCCESS); 00950 } 00951 00952 /************************************************************************* 00953 * WriteSACEvent () * 00954 * Write the information from the sac header to the file * 00955 *************************************************************************/ 00956 int WriteSAC_NextStationForEvent(EWChannelDataStruct * pChannel) 00957 { 00958 00959 int rc; 00960 00961 rc=SACPABase_end_scn(); 00962 00963 if(rc != EW_SUCCESS) 00964 { 00965 logit("","ProduceSAC_NextStationForEvent():SACPABase_end_scn() \n" 00966 "failed with rc=%d for %s,%s,%s\n", 00967 rc, pChannel->Station.Sta, pChannel->Station.Comp, 00968 pChannel->Station.Net); 00969 return(EW_FAILURE); 00970 } 00971 00972 return EW_SUCCESS; 00973 } 00974 00975 00976 int WriteSAC_Set_dGapThresh(float IN_dGapThresh) 00977 { 00978 WriteSAC_dGapThresh=IN_dGapThresh; 00979 return(EW_SUCCESS); 00980 } 00981