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 

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