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

# sac_2_ewevent.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: 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