source: trunk/src/seismic_processing/gmew/gm_util.c @ 8019

Revision 8019, 35.4 KB checked in by kevin, 3 months ago (diff)

Added Arias Intensity and ColorPGA option

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/*
2 *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
3 *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
4 *
5 *    $Id$
6 *
7 *    Revision history:
8 *     $Log$
9 *     Revision 1.10  2006/03/15 14:21:54  paulf
10 *     SCNL version of gmew v0.2.0
11 *
12 *     Revision 1.9  2002/09/26 19:29:37  dhanych
13 *     fixed 'BadBitMap' to 'BadBitmap'
14 *
15 *     Revision 1.8  2002/09/25 17:33:20  dhanych
16 *     added state checking at start of makeGM() to reject traces for error
17 *     conditions found by prepTrace(); e.g. gaps, signal-to-noise.
18 *
19 *     Revision 1.7  2002/02/28 17:03:36  lucky
20 *     Moved gma.c and gma.h to libsrc and main include
21 *
22 *     Revision 1.6  2001/07/18 19:41:36  lombard
23 *     Changed XMLDir, TempDir and MappingFile in GMPARAMS struct from string
24 *     arrays to string pointers. Changed gm_config.c and gm_util.c to support thes
25 *     changes. This solved a problem where the GMPARAMS structure was getting
26 *     corrupted when a pointer to it was passed into getGMFromTrace().
27 *     It's not clear why this was necessary; purify didn't complain.
28 *
29 *     Revision 1.5  2001/07/18 19:18:25  lucky
30 *     *** empty log message ***
31 *
32 *     Revision 1.4  2001/06/10 21:29:07  lombard
33 *     fixed memory leak in endEvent.
34 *
35 *     Revision 1.3  2001/04/21 19:18:32  lombard
36 *     removed prohibition of using Z components
37 *
38 *     Revision 1.2  2001/04/11 21:12:20  lombard
39 *     ../localmag/site.h ../localmag/lm_site.h
40 *
41 *     Revision 1.1  2001/03/30 19:14:25  lombard
42 *     Initial revision
43 *
44 *
45 *
46 */
47/*
48 * gm_util.c: A bunch of utility functions used for gmew.
49 * 
50 *    Pete Lombard; January, 2001
51 */
52
53
54#include <stdio.h>
55#include <stdlib.h>
56#include <string.h>
57#include <math.h>
58#include <errno.h>
59#include <earthworm.h>
60#include <chron3.h>
61#include <fft_prep.h>
62#include <read_arc.h>
63#include <rw_strongmotionII.h>
64#include <tlay.h>
65#include <transfer.h>
66#include <gma.h>
67#include "../localmag/lm_misc.h"
68#include "gm.h"
69#include "gm_sac.h"
70#include "gm_util.h"
71#include "gm_ws.h"
72#include "gm_xml.h"
73#include "../localmag/lm_site.h"
74
75#define GRAVITY 980.665 /* G cm/s/s  */
76
77/* Standard Spectral Response periods and damping value */
78/*#define NSP 3*/
79#define NSD 1
80static int NSP = 0;
81static double spectPer[SM_MAX_RSA];
82static double spectDamp[NSD] = {0.05};
83
84/* global variables for gm_util.c: */
85static DATABUF gTrace;    /* the trace buffer for raw and processed data  */
86static double *gWork;     /* the work array for convertWave               */
87static double *gWorkFFT;  /* the work array for fft99                     */
88static ResponseStruct gScnPZ;  /* Poles, zeros, and gain for the SCNL      */
89static SM_INFO gSM;       /* Strong-motion infor structure                */
90
91/* Internal Function Prototypes */
92static int getRespPZ(char *, char *, char *, char *, GMPARAMS *, EVENT *);
93
94/*
95 * getGMFromTrace: selector function to call the desired procedure
96 *           for getting ground-motion from a trace source such as wave_server,
97 *           SAC file, etc.
98 *    Returns: 0 on success
99 *            -1 on fatal errors
100 */
101int getGMFromTrace( EVENT *pEv, GMPARAMS *pgmParams)
102{
103  int rc;
104 
105  /* Initialize the SAC trace saver. We assume that this is the   *
106   * only type of trace saver; if more are added, then we'll need *
107   * an initSave() routine that will select the desired saver.    */
108  if (pEv->eventId[0] && pgmParams->saveTrace == GM_ST_SAC)
109    initSACSave(pEv, pgmParams);
110
111  /* Initialize the XML writer */
112  if (pEv->eventId[0] && pgmParams->XMLDir && strlen(pgmParams->XMLDir) > 0)
113    (void) Start_XML_file( pgmParams, pEv );
114
115  switch(pgmParams->traceSource)
116  {
117  case GM_TS_WS:
118    rc = getGMFromWS( pEv, pgmParams, &gTrace, &gSM);
119    break;
120#ifdef EWDB
121  case GM_TS_EWDB:
122    rc = getGMFromEWDB( pEv, pgmParams, &gTrace, &gSM);
123    break;
124#endif
125  default:
126    logit("et", "gmew getGMFromTrace: unknown source: %d\n",
127          pgmParams->traceSource);
128    rc = -1;
129  }
130
131  if (pEv->eventId[0] && pgmParams->saveTrace == GM_ST_SAC)
132    termSACSave(pEv, pgmParams);
133 
134  /* Terminate the XML writer */
135  if (pEv->eventId[0] && pgmParams->XMLDir && strlen(pgmParams->XMLDir) > 0)
136  {
137    if (Close_XML_file( )== 0 && pgmParams->sendActivate == 1) 
138    {
139       Send_XML_activate(pEv, pgmParams);
140    }
141  }
142
143  return rc;
144}
145
146/*
147 * addCompToEvent: Add a new component to the EVENT structure.
148 *     Checks that SCN meets the selection criteria of the lists add and pDel,
149 *     the station is within the distance limits, and that the limit
150 *     on number of stations is not exceeded.
151 *   Returns: 0 on success
152 *           +1 if comp is a duplicate of one already in the list
153 *           +2 if comp is not a selected component
154 *           +3 if station number limit is reached
155 *           +4 if component direction is not understood
156 *           +5 if station is outside distance limit
157 *           +6 if station location cannot be determined
158 *           -1 on failure (error allocating memory)
159 */
160int addCompToEvent( char *sta, char *comp, char *net, char *loc, EVENT *pEvt, 
161                    GMPARAMS *pgmParams, STA **ppSta, COMP3 **ppComp)
162{
163  COMP1 *this, *last;
164  PSCNLSEL thisSel;
165  SCNLPAR keySCNL;
166  double dist, lat, lon;
167  int i, dir;
168  int foundIt = 0, new = 0;
169  int selected = 0;
170  int mSta, mComp, mNet, mLoc;   /* `match' flags */
171  char orientation;
172 
173  /* Is this SCNL on the selection "Add" list?                         *
174   * If Add list is empty, assume that all stations are wanted.       */
175  if ( (thisSel = pgmParams->pAdd) == (PSCNLSEL)NULL )
176    selected = 1;
177  else
178  {
179    while (thisSel != (PSCNLSEL)NULL)
180    {
181      mSta = mComp = mNet = mLoc = 0;
182      if (thisSel->sta[0] == '*' || strcmp(thisSel->sta, sta) == 0) mSta = 1;
183      if (thisSel->net[0] == '*' || strcmp(thisSel->net, net) == 0) mNet = 1;
184      if (thisSel->comp[0] == '*' || 
185          memcmp(thisSel->comp, comp, strlen(thisSel->comp)) == 0) mComp = 1;
186      if (thisSel->loc[0] == '*' || 
187          memcmp(thisSel->loc, loc, strlen(thisSel->loc)) == 0) mLoc = 1;
188      if ( (selected = mSta && mNet && mComp && mLoc) == 1)
189        break;
190     
191      thisSel = thisSel->next;
192    }
193  }
194  if (!selected)
195  {
196    if (pgmParams->debug & GM_DBG_SEL)
197      logit("et", "addCompToEvent: <%s.%s.%s.%s> not in select list: (%d.%d.%d.%d)\n",
198            sta, comp, net, loc, mSta, mComp, mNet, mLoc);
199    return +2;
200  }
201
202  /* Is this SCNL on the selection "Del" list? */
203  thisSel = pgmParams->pDel;
204  while (thisSel != (PSCNLSEL)NULL)
205  {
206    mSta = mComp = mNet =  mLoc =0;
207    if (thisSel->sta[0] == '*' || strcmp(thisSel->sta, sta) == 0) mSta = 1;
208    if (thisSel->net[0] == '*' || strcmp(thisSel->net, net) == 0) mNet = 1;
209    if (thisSel->comp[0] == '*' || 
210        memcmp(thisSel->comp, comp, strlen(thisSel->comp)) == 0) mComp = 1;
211    if (thisSel->loc[0] == '*' || 
212        memcmp(thisSel->loc, loc, strlen(thisSel->loc)) == 0) mLoc = 1;
213    if ( (mSta && mNet && mComp && mLoc) == 1)
214    {
215      if (pgmParams->debug & GM_DBG_SEL)
216        logit("et", "addCompToEvent: <%s.%s.%s.%s> in delete list: (%d.%d.%d.%d)\n",
217              sta, comp, net, loc, mSta, mComp, mNet, mLoc);
218      return +2;
219    }
220    thisSel = thisSel->next;
221  }
222
223  orientation = comp[2]; /* the 3rd char in the seed channel name holds orientation info ZNE, 123 etc */
224
225  /* try map a number to a channel string */
226  if (strlen(pgmParams->ChannelNumberMap) > 1) 
227  {
228    int num;
229    num = atoi(&comp[2]);
230    if (num >= 1 && num <= 3) 
231    {
232      orientation = pgmParams->ChannelNumberMap[num];
233      logit("", "addCompToEvent: mapping orientation number %d to orientation %c SCNL now <%s.%s.%s.%s>\n",
234          num, orientation, sta, comp, net, loc);
235    }
236  }
237 
238  /* What direction is this component? We do this here, since it could
239   * possibly fail. We don't want to allocate a new COMP1 structure (below)
240   * and then not use it because this test failed. */
241  switch(orientation)
242  {
243  case 'E':
244    dir = GM_E;
245    break;
246  case 'N':
247    dir = GM_N;
248    break;
249  case 'Z':
250    dir = GM_Z;
251    break;
252  default:
253    logit("et", "addCompToEvent: unknown component direction for <%s.%s.%s.%s>\n",
254          sta, comp, net, loc);
255    return +4;
256  }
257
258  /* Search the station list for a match of sta, net */
259  for (i = 0; i < pEvt->numSta; i++)
260  {
261    if (strncmp(pEvt->Sta[i].sta, sta, 6) == 0 &&
262        strncmp(pEvt->Sta[i].net, net, 3) == 0 )
263    {
264        foundIt = 1;
265        break; 
266    }
267  }
268  /* Our desired sta/net is at index `i' regardless of foundIt value */
269  if (foundIt == 0)  /* No match; add to the list if it isn't full */
270  {
271    if ( pEvt->eventId[0]!=0 && pEvt->eventId[0]!='-' ) {
272                /* Find its epicentral distance; make sure its in range */
273                dist = getStaDist(sta, comp, net, loc, &lat, &lon, pEvt, pgmParams);
274                if ( dist < 0.0 )
275                  return +6;
276                else if (dist > pgmParams->maxDist)
277                {
278                  if (pgmParams->debug & GM_DBG_SEL)
279                        logit("et", "<%s.%s.%s.%s> distance (%e) exceeds limit\n", sta, comp, net, loc,
280                                  dist);
281                  return +5;
282                }
283        }
284    if (pEvt->numSta < pgmParams->maxSta)
285    {
286      strncpy(pEvt->Sta[i].sta, sta, TRACE_STA_LEN);
287      strncpy(pEvt->Sta[i].net, net, TRACE_NET_LEN);
288      if ( pEvt->eventId[0] ) {
289                  pEvt->Sta[i].lat = lat;
290                  pEvt->Sta[i].lon = lon;
291                  pEvt->Sta[i].dist = dist;
292          }
293          pEvt->Sta[i].comp = NULL;
294      pEvt->numSta++;
295    }
296    else
297    {
298      logit("et", 
299            "addCompToEvent: station limit reached; skipping <%s.%s.%s.%s>\n",
300            sta, comp, net, loc );
301      pEvt->Sta[i].dist = -1.0;
302      return +3;
303    }
304  }
305  *ppSta = &pEvt->Sta[i];
306 
307  if ( pEvt->Sta[i].comp == (COMP1 *)NULL)
308  { 
309      /* No components for this station, so create a new one */
310    if ( (pEvt->Sta[i].comp = (COMP1 *)calloc(sizeof(COMP1), 1)) == 
311         (COMP1 *)NULL)
312    {
313      logit("et", "addCompToEvent: out of memory for COMP1\n");
314      return -1;
315    }
316    this = pEvt->Sta[i].comp;
317    new = 1;
318  }
319  else
320  {
321    this = pEvt->Sta[i].comp;
322    foundIt = 0;
323    /* Walk the COMP list */
324    do
325    {    /* Search for match of the first 2 letters in component name */
326      if (memcmp(this->n2, comp, 2) == 0)
327      {
328        foundIt = 1;
329        break;
330      }
331      last = this;
332    } 
333    while ( (this = this->next) != (COMP1 *)NULL);
334   
335    if (foundIt == 0)
336    {           /* No match, so add a new COMP1 structure to the list */
337      if  ( (last->next = (COMP1 *)calloc(sizeof(COMP1), 1)) == 
338            (COMP1 *)NULL)
339        return -1;
340      new = 1;
341      this = last->next;
342    }
343  }
344  /* Now `this' is pointing at our COMP1 structure */
345
346  if (new == 1)
347  {      /* Fill in the new COMP1 name */
348    this->n2[0] = comp[0]; 
349    this->n2[1] = comp[1];
350  }
351  else if (this->c3[dir].name[0] != 0 && pgmParams->allowDuplicates == 0)
352  {  /* This component/direction is a duplicate */
353    if ( (pgmParams->debug & GM_DBG_SEL) > 0 )
354      logit("et", "addCompToEvent: <%s.%s.%s.%s> rejected as a duplicate\n",
355            sta, comp, net, loc);
356    return +1;
357  }
358
359  /* This component/direction is new for this event */
360  strncpy(this->c3[dir].name, comp, 3);
361  strncpy(this->c3[dir].loc, loc, 3);
362  this->c3[dir].BadBitmap = 0; /* this may just be a different location code, zero out values from prior comps */
363  *ppComp = &this->c3[dir];
364
365  /* Search for an SCNLPAR structure for this SCNL */
366  if (pgmParams->numSCNLPar > 0)
367  {
368    memset(&keySCNL, 0, sizeof(SCNLPAR));
369    strcpy(keySCNL.sta, sta);
370    strcpy(keySCNL.comp, comp);
371    strcpy(keySCNL.net, net);
372    strcpy(keySCNL.loc, loc);
373    this->c3[dir].pSCNLPar = (SCNLPAR *)bsearch(&keySCNL, pgmParams->pSCNLPar,
374                                              pgmParams->numSCNLPar,
375                                              sizeof(SCNLPAR), CompareSCNLPARs);
376  }
377  else
378    this->c3[dir].pSCNLPar = (SCNLPAR *)NULL;
379 
380  return 0;
381}
382
383
384/*
385 * getStaDist: Determine the epicentral distance for a station, given
386 *          the station name "sta", the network name "net", the component
387 *          name "comp", and the EVENT structure "pEvt".
388 *          Currently, staDist gets station locations from the
389 *          site table initialized with the site_load() calls from
390 *          site.c. This initialization must be done before calling staDist().
391 *   Returns: the epicentral distance if the station could be found
392 *           -1.0 if the station location could not be determined.
393 */
394double getStaDist( char *sta, char *comp, char *net, char *loc, double *pLat, 
395                   double *pLon, EVENT *pEvt, GMPARAMS *pgmParams)
396{
397  double r = -1.0;           /* Epicentral distance */
398  SITE *pSite;               /* Site table entry */
399 
400  switch(pgmParams->staLoc)
401  {
402  case GM_SL_HYP:
403    if ( (pSite = find_site( sta, comp, net, loc )) == 
404         (SITE *)NULL)
405    {
406      logit("et", "getStaDist: <%s.%s.%s.%s> - Not in station list.\n", 
407            sta, comp, net, loc);
408      return r;
409    }
410    *pLat = pSite->lat;
411    *pLon = pSite->lon;
412    if ( pEvt->eventId[0] )
413                r = utmcal(pEvt->lat, pEvt->lon, pSite->lat, pSite->lon);
414        else
415                r = 0.0;
416    break;
417#ifdef EWDB
418  case GM_SL_EWDB:
419    /* Get station location from EWDB */
420    break;
421#endif
422  }
423
424
425  return r;
426}
427
428
429
430/*
431 * EstPhaseArrivals: estimate the P and S arrival times for this
432 *    station; set the values for all components of this station.
433 *    Uses the "tlay" routines for travel-time calculations using
434 *    a single multi-layered regional velocity model.
435 */ 
436void EstPhaseArrivals(STA *pSta, EVENT *pEvt, int debug)
437{
438  TPHASE treg[4];
439  int nreg;
440 
441  nreg = t_region(pSta->dist, pEvt->depth, treg);
442  if (nreg == 2)
443  {   
444    pSta->p_est = treg[0].t;
445    pSta->s_est = treg[1].t;
446  }
447  else
448  {
449    pSta->p_est = (treg[0].t < treg[2].t) ? treg[0].t : treg[2].t;
450    pSta->s_est = (treg[1].t < treg[3].t) ? treg[1].t : treg[3].t;
451  }
452 
453  if (debug)
454  {
455    logit("et","Est P: %10.4f", treg[0].t);
456    if (nreg == 2)
457      logit("e", "\n");
458    else
459      logit("e", "  (%10.4f)\n", treg[2].t);
460    logit("et", "Est S: %10.4f", treg[1].t);
461    if (nreg == 2)
462      logit("e", "\n");
463    else
464      logit("e", "  (%10.4f)\n", treg[3].t);
465  }
466
467  pSta->p_est += pEvt->origin_time;
468  pSta->s_est += pEvt->origin_time;
469 
470  return;
471}
472
473
474/*
475 * Return the desired start time for traces from pSta.
476 */
477double traceStartTime(STA *pSta, GMPARAMS *pgmParams)
478{
479  return pSta->p_est - pgmParams->traceStart;
480}
481
482
483/*
484 * Return the desired end time for traces from pSta.
485 */
486double traceEndTime(STA *pSta, GMPARAMS *pgmParams)
487{
488  return pSta->s_est + pgmParams->traceEnd;
489}
490
491/*
492 * Return the desired peak-search start time for traces from pSta.
493 */
494double peakSearchStart(STA *pSta, GMPARAMS *pgmParams)
495{
496  double start;
497 
498  start = pgmParams->peakSearchStart * (pSta->s_est - pSta->p_est);
499  start = (start > pgmParams->peakSearchStartMin) ? start :
500    pgmParams->peakSearchStartMin;
501  return pSta->s_est - start;
502}
503
504
505/*
506 * Return the desired peak-search end time for traces from pSta.
507 */
508double peakSearchEnd(STA *pSta, GMPARAMS *pgmParams)
509{
510  double end;
511 
512  end = pgmParams->peakSearchEnd * (pSta->s_est - pSta->p_est);
513  end = (end > pgmParams->peakSearchEndMin) ? end :
514    pgmParams->peakSearchEndMin;
515  return pSta->s_est + end;
516}
517
518
519
520/*
521 * MakeGM: Transfer a raw trace into synthetic GM traces.
522 *         The trace data is in the DATABUF trace, for the SCN identified
523 *         by STA and COMP. Instrument response (pole/zero/gain) is
524 *         obtained from the specified response source.
525 *   Returns: 0 on success
526 *           +1 on non-fatal error (no response data for this SCN)
527 *           -1 on fatal errors 
528 */
529int makeGM(STA *pSta, COMP3 *pComp, GMPARAMS *pgmParams, EVENT *pEvt, 
530            DATABUF *pTrace)
531{
532  int i, rc;
533  double tTaper, taperFreqs[4] = {1.0, 1.0, 10.0, 10.0};
534  long padlen, nfft;
535 
536  /* reject bad traces (gaps, S/N ratio, etc.)  */
537  if ( pComp->BadBitmap )
538     return 1;
539 
540  /* Get the instrument response */
541  if ( (rc = getRespPZ(pSta->sta, pComp->name, pSta->net, pComp->loc, pgmParams, 
542                       pEvt)) != 0)
543  {
544    logit("et", "makeGM: no response data for <%s.%s.%s.%s>; skipping\n",
545          pSta->sta, pComp->name, pSta->net, pComp->loc);
546    return rc;
547  }
548
549  if (pComp->pSCNLPar == (SCNLPAR *)NULL)
550  {  /* Use the defaults if not set externally */
551    /* Set the low-frequency taper band to 0.05 and 0.1 hz.*/
552    taperFreqs[0] = 0.05;
553    taperFreqs[1] = 0.1;
554
555    /* Set the high-frequency taper band to 90% and 100% of Nyquist */
556    taperFreqs[3] = 0.5/pTrace->delta;
557    taperFreqs[2] = 0.45/pTrace->delta;
558
559    /* Set acceleration time-series taper length to default 0 */
560    tTaper = 0.0;
561  }
562  else
563  {
564    for (i = 0; i < 4; i++)
565      taperFreqs[i] = pComp->pSCNLPar->fTaper[i];
566    tTaper = pComp->pSCNLPar->taperTime;
567  }
568
569  if (pgmParams->debug & (GM_DBG_PZG | GM_DBG_TRS | GM_DBG_ARS) )
570    printf("\nResponse for <%s.%s.%s.%s> delta %e\nftaper %e %e %e %e\n", 
571           pSta->sta, pComp->name, pSta->net, pComp->loc, pTrace->delta, taperFreqs[0],
572           taperFreqs[1], taperFreqs[2], taperFreqs[3]);
573
574  padlen = -1;  /* let convertWave figure the padding */
575  rc = gma(pTrace->rawData, pTrace->nRaw, pTrace->delta, &gScnPZ, taperFreqs, 
576           tTaper, &padlen, &nfft, spectPer, NSP, spectDamp, NSD, 
577           pTrace->procData, pTrace->lenProc, gWork, gWorkFFT);
578  if (rc < 0)
579  {
580    switch(rc)
581    {
582    case -1:
583      logit("et", "gma failed: out of memory\n");
584      return -1;
585      break;
586    case -3:
587      logit("et", "gma failed: invalid arguments\n");
588      return -1;
589      break;
590    case -4:
591      logit("et", "gma failed: FFT error; nfft: %ld\n", nfft);
592      return -1;
593      break;
594    default:
595      logit("et", "gma failed: unknown error %d\n", rc);
596      return -1;
597    }
598  }
599  /* Do we need to adjust the end of the processed trace? */
600  if (nfft - padlen < pTrace->nRaw)
601  {    /* gma had to chop some of the end */
602    pTrace->endtime -= (pTrace->nRaw - (nfft - padlen)) * pTrace->delta;
603    pTrace->nProc = nfft - padlen;
604  }
605  else
606    pTrace->nProc = pTrace->nRaw;
607
608  pTrace->padLen = padlen;
609 
610  return 0;
611}
612
613/*
614 * prepTrace: prepare trace data for processing. Preps include: check for
615 *    gaps in peak-search window, compute and remove mean, fill gaps,
616 *    and check for clipping.
617 *    Also fills in the peak-search start and end times.
618 */
619void prepTrace(DATABUF *pTrace, STA *pSta, COMP3 *pComp, GMPARAMS *pgmParams)
620{
621  long i, iStart, iEnd, npoints;
622  double mean, dMax, dMin, gapEnd, clipLimit, preP, postP;
623  GAP *pGap;
624
625  if (pComp->pSCNLPar != (SCNLPAR *)NULL)
626    clipLimit = pComp->pSCNLPar->clipLimit;
627  else
628    clipLimit = 7.55e6; /* 90% of 2^23; clip for 24-bit digitizer */
629 
630  /* Fill in the peak-search window limits */
631  pComp->peakWinStart = peakSearchStart(pSta, pgmParams);
632  pComp->peakWinEnd = peakSearchEnd(pSta, pgmParams);
633
634  /* Find mean value of non-gap data */
635  pGap = pTrace->gapList;
636  i = 0L;
637  npoints = 0L;
638  mean = 0.0;
639  dMax = -4.0e6;
640  dMin = 4.0e6;
641  /*
642   * Loop over all the data, skipping any gaps. Note that a `gap' will not
643   * be declared at the end of the data, so the counter `i' will always
644   * get to pTrace->nRaw.
645   */
646  do
647  {
648    if (pGap == (GAP *)NULL)
649    {
650      iEnd = pTrace->nRaw;
651    }
652    else
653    {
654      iEnd = pGap->firstSamp - 1;
655
656      /* Test for gap within peak-search window */
657      gapEnd = pGap->starttime + (pGap->lastSamp - pGap->firstSamp + 1) *
658        pTrace->delta;
659      if ( (pGap->starttime >= pComp->peakWinStart && 
660            pGap->starttime <= pComp->peakWinEnd) ||
661           (gapEnd >= pComp->peakWinStart && gapEnd <= pComp->peakWinEnd) ||
662           (pGap->starttime < pComp->peakWinStart && 
663            gapEnd > pComp->peakWinEnd) )
664      {
665        logit("et", "trace from <%s.%s.%s.%s> has gap in peak-search window\n",
666              pSta->sta, pComp->name, pSta->net, pComp->loc);
667        pComp->BadBitmap |= GM_BAD_GAP;
668      }
669    }
670    for (; i < iEnd; i++)
671    {
672      mean += pTrace->rawData[i];
673      dMax = (dMax > pTrace->rawData[i]) ? dMax : pTrace->rawData[i];
674      dMin = (dMin < pTrace->rawData[i]) ? dMin : pTrace->rawData[i];
675      npoints++;
676    }
677    if (pGap != (GAP *)NULL)     /* Move the counter over this gap */   
678    {
679      i = pGap->lastSamp + 1;
680      pGap = pGap->next;
681    }
682  } while (i < pTrace->nRaw );
683 
684  mean /= (double)npoints;
685  if (dMax > clipLimit || dMin < -clipLimit)
686  {
687    logit("et", "trace from <%s.%s.%s.%s> may be clipped\n", pSta->sta,
688          pComp->name, pSta->net, pComp->loc);
689    pComp->BadBitmap |= GM_BAD_CLIP;
690  }
691 
692  /* Now remove the mean, and set points inside gaps to zero */
693  i = 0;
694  do
695  {
696    if (pGap == (GAP *)NULL)
697      iEnd = pTrace->nRaw;
698    else
699      iEnd = pGap->firstSamp - 1;
700
701    for (; i < iEnd; i++)
702      pTrace->rawData[i] -= mean;
703
704    if (pGap != (GAP *)NULL)     /* Fill in the gap with zeros */   
705    {
706      for ( ;i < pGap->lastSamp + 1; i++)
707        pTrace->rawData[i] = 0.0;
708      pGap = pGap->next;
709    }
710  } while (i < pTrace->nRaw );
711
712  /* See how much noise is in the pre-P-arrival data.       *
713   * We hope and pray it isn't the coda of a previous event */
714  preP = postP = 0.0;
715  iEnd = (long)( 0.5 + (pSta->p_est - pTrace->starttime) / 
716                 pTrace->delta );
717  for (i = 0; i < iEnd; i++)
718    preP += pTrace->rawData[i] * pTrace->rawData[i];
719
720  if (iEnd > 0)
721    preP /= (double)iEnd;
722 
723  iStart = (long)( 0.5 + (pComp->peakWinStart - pTrace->starttime) / 
724                   pTrace->delta );
725  if (iStart < 0 )
726    iStart = 0;
727  iEnd = (long)( 0.5 + (pComp->peakWinEnd - pTrace->starttime ) / 
728                 pTrace->delta);
729  if (iEnd > pTrace->nRaw)
730    iEnd = pTrace->nRaw;
731  for (i = iStart; i < iEnd; i++)
732    postP += pTrace->rawData[i] * pTrace->rawData[i];
733
734  if (iEnd > iStart)
735    postP /= (double)(iEnd - iStart);
736 
737  if (postP < preP * pgmParams->snrThresh)
738  {
739    logit("et", "prepTrace: <%s.%s.%s.%s> event signal (%5g) too small for pre-event threshold (%5g)\n",
740          pSta->sta, pComp->name, pSta->net, pComp->loc, postP, 
741          pgmParams->snrThresh * preP);
742    pComp->BadBitmap |= GM_LOW_SNR;
743  }
744
745  return;
746}
747
748/* calcAriasIntensity: Calculate the AI (Arias Intensity) from acceleration data.
749*/
750double calcAriasIntensity(double delta, double *acc, int acclen)
751{
752  int i;
753  double acc_value, ai, sum = 0.;
754  // integrate acceleration squared
755  for (i = 0; i < acclen; i++)
756  {
757    acc_value = acc[i];
758    // square the acceleration
759    acc_value *= acc_value;
760    if (i == 0 || i == acclen - 1)
761      sum += acc_value;
762    else
763      sum += 2. * acc_value;
764  }
765  sum *= (delta/2.);
766  ai = sum;
767  ai *= (M_PI / (2. * GRAVITY));
768  return ai;
769}
770
771/* getPeakGM: Find the maximum absolute values of acceleration, velocity,
772 *            displacement and Spectra Response from the processed data.
773 *            This function looks at the data within the "peak search"
774 *            window, defined by peakWinStart and peakWinEnd in the COMP3
775 *            structure, for the peak acceleration, velocity and displacement.
776 *            It looks from peakWinStart to the end of the padded data for
777 *            the peak Spectral Response values.
778 *            The peak values are used to fill in a SM_INFO structure
779 *            which is returned to be used by the caller.
780 */
781void getPeakGM(DATABUF *pTrace, COMP3 *pComp, STA *pSta, GMPARAMS *pgmParams,
782               EVENT *pEvt, SM_INFO *pSms)
783{
784  double *acc, *vel, *disp, *psa;
785  double mAcc, mVel, mDisp, mPsa;
786  long iAcc, iVel, iDisp, iPsa;
787  long i, iStart, iEnd;  /* start and end of peak-search window */
788  int isp;
789 
790  memset(pSms, 0, sizeof(SM_INFO));
791  strcpy(pSms->sta, pSta->sta);
792  strcpy(pSms->comp, pComp->name);
793  strcpy(pSms->net, pSta->net);
794  strcpy(pSms->loc, pComp->loc);
795  /* pSms->loc[0] = '\0'; */
796  strcpy(pSms->qid, pEvt->eventId);
797  strcpy(pSms->qauthor, pEvt->authorId);
798 
799  pSms->talt = 0.0;
800  pSms->altcode = SM_ALTCODE_NONE;
801  pSms->tload = -1.0;
802 
803  /* Search for the extrema in acc, vel, disp */
804  if ( pEvt->eventId[0]==0 ) {
805        pComp->peakWinStart = pTrace->starttime;
806        pComp->peakWinEnd = pComp->peakWinStart + pgmParams->alarmDuration;
807  }
808  iStart = (long)( 0.5 + (pComp->peakWinStart - pTrace->starttime) / 
809                                   pTrace->delta );
810  if (iStart < 0 )
811        iStart = 0;
812  iEnd = (long)( 0.5 + (pComp->peakWinEnd - pTrace->starttime ) / 
813                                 pTrace->delta);
814  if (iEnd > pTrace->nProc)
815        iEnd = pTrace->nProc;
816 
817  if (pgmParams->debug & GM_DBG_TIME)
818  {
819    logit("e", "trace start: %10.4lf end: %10.4lf (%ld)\n", 
820          pTrace->starttime - pEvt->origin_time,
821          pTrace->endtime - pEvt->origin_time, pTrace->nProc);
822    logit("e", "search start: %10.4lf (%ld) end: %10.4lf (%ld)\n", 
823          pComp->peakWinStart - pEvt->origin_time, iStart, 
824          pComp->peakWinEnd - pEvt->origin_time, iEnd);
825  }
826
827  acc = pTrace->procData;
828  vel = &pTrace->procData[pTrace->lenProc];
829  disp = &pTrace->procData[pTrace->lenProc * 2];
830 
831  mAcc = mVel = mDisp = 0.0;
832  iAcc = iVel = iDisp = -1;
833  for (i = iStart; i < iEnd; i++)
834  {
835    if (fabs(acc[i]) > mAcc)
836    {
837      mAcc = fabs(acc[i]);
838      iAcc = i;
839    }
840    if (fabs(vel[i]) > mVel)
841    {
842      mVel = fabs(vel[i]);
843      iVel = i;
844    }
845    if (fabs(disp[i]) > mDisp)
846    {
847      mDisp = fabs(disp[i]);
848      iDisp = i;
849    }
850  }
851  pSms->pga = mAcc;
852  pSms->tpga = pTrace->starttime + pTrace->delta * iAcc;
853  pSms->t = pSms->tpga;
854 
855  pSms->pgv = mVel;
856  pSms->tpgv = pTrace->starttime + pTrace->delta * iVel;
857  if (pSms->tpgv < pSms->t) pSms->t = pSms->tpgv;
858 
859  pSms->pgd = mDisp;
860  pSms->tpgd = pTrace->starttime + pTrace->delta * iDisp;
861  if (pSms->tpgd < pSms->t) pSms->t = pSms->tpgd;
862 
863  /* Search for peak Spectral Response values.             *
864   * This section is hard-coded to use the first           *
865   * spectDamp value, which is 0.05 as needed by ShakeMap. */ 
866  iStart = (long)( 0.5 + (pComp->peakWinStart - pTrace->starttime) / 
867                   pTrace->delta );
868  if (iStart < 0 )
869    iStart = 0;
870  /* Search all the way out through the padded data; *
871   * This is needed since the peak Spectral Response *
872   * may occur  after the forcing is completed.      */
873  iEnd = pTrace->nProc + pTrace->padLen;
874 
875  for (isp = 0; isp < NSP; isp++)
876  {
877    mPsa = 0.0;
878    psa = &pTrace->procData[pTrace->lenProc * (3 + isp)];
879    for (i = iStart; i < iEnd; i++)
880    {
881      if (mPsa < fabs(psa[i])) 
882      {
883        mPsa = fabs(psa[i]);
884        iPsa = i;
885      }
886    }
887    pSms->rsa[isp] = mPsa;
888    pSms->pdrsa[isp] = spectPer[isp];
889    /* Record the time in case we're saving results in SAC files */
890    pComp->RSAPeakTime[isp] = pTrace->starttime + pTrace->delta * iPsa;
891  }
892  pSms->nrsa = NSP;
893 
894  if (pgmParams->ariasIntensity)
895  {
896     pSms->ai = calcAriasIntensity(pTrace->delta, acc, pTrace->lenProc);
897  }
898
899  /* Write some XML */
900  if (pgmParams->XMLDir && strlen(pgmParams->XMLDir) > 0)
901    next_XML_SCNL( pgmParams, pSta, pSms );
902
903  return;
904}
905
906
907/*
908 * cleanTrace: initialize the global trace buffer, freeing old GAP structures.
909 *    Returns: nothing
910 */
911void cleanTrace( DATABUF *pTrace )
912{
913  GAP *pGap;
914 
915  pTrace->nRaw = 0L;
916  pTrace->nProc = 0L;
917  pTrace->delta = 0.0;
918  pTrace->starttime = 0.0;
919  pTrace->endtime = 0.0;
920  pTrace->nGaps = 0;
921 
922  /* Clear out the gap list */
923  pTrace->nGaps = 0;
924  while ( (pGap = pTrace->gapList) != (GAP *)NULL)
925  {
926    pTrace->gapList = pGap->next;
927    free(pGap);
928  }
929  return;
930}
931
932/*
933 * initBufs: allocate the three arrays needed for handling trace data
934 *           The argument reqLen is the reqeusted length; this value
935 *           is used for the raw data array. That value or larger,
936 *           if needed to find a multiple of the FFT factors, is used
937 *           for the size of the processed data array. The work array
938 *           is sized as needed by the gma routine.
939 *           This routine is intended to be called at startup instead
940 *           of later in the process life, to minimize memory growth
941 *           during process lifetime.
942 *   Returns: 0 on success
943 *           -1 when out of memory
944 */
945int initBufs( long reqLen )
946{
947  long nfft;
948  FACT *pF;
949 
950  /* Array for input data */
951  if ( (gTrace.rawData = (double *)malloc(reqLen * sizeof(double))) ==
952       (double*)NULL)
953  {
954    logit("et", "initBuffs: out of memory for rawData\n");
955    return -1;
956  }
957  gTrace.lenRaw = reqLen;
958 
959  if ( (nfft = prepFFT( reqLen, &pF)) < 0)
960  {
961    logit("et", "initBuffs: out of memory for FFT factors\n");
962    return -1;
963  }
964 
965  /* Number of processed data arrays: 3 for acc, vel and disp, plus *
966   * the ones needed for Spectral Response */
967  gTrace.numProc = 3 + NSD * NSP;
968  if ( (gTrace.procData = 
969        (double *)malloc(gTrace.numProc * (nfft + FFT_EXTRA) * 
970                         sizeof(double))) == (double*)NULL)
971  {
972    logit("et", "initBuffs: out of memory for procData\n");
973    return -1;
974  }
975  gTrace.lenProc = nfft + FFT_EXTRA;
976 
977  /* Work array for frequency response functions in gma:  *
978   * 3 for acc, vel, disp, plus one for Spectral Response */
979  if ( (gWork = (double *)malloc(4 * (nfft + 2) * 
980                                 sizeof(double))) == (double*)NULL)
981  {
982    logit("et", "initBuffs: out of memory for gWork\n");
983    return -1;
984  }
985
986  /* Work array for Temperton FFTs */
987  if ( (gWorkFFT = 
988        (double *)malloc(gTrace.numProc * (nfft + FFT_EXTRA) * 
989                         sizeof(double))) == (double*)NULL)
990  {
991    logit("et", "initBuffs: out of memory for gWorkFFT\n");
992    return -1;
993  }
994
995  return 0;
996}
997
998
999/*************************************************************
1000 *                   CompareSCNLPARs()                       *
1001 *                                                           *
1002 *  This function is passed to qsort() and bsearch().        *
1003 *************************************************************/
1004int CompareSCNLPARs( const void *s1, const void *s2 )
1005{
1006   int rc;
1007   SCNLPAR *t1 = (SCNLPAR *) s1;
1008   SCNLPAR *t2 = (SCNLPAR *) s2;
1009
1010   rc = strcmp( t1->sta, t2->sta );
1011   if ( rc != 0 ) return(rc);
1012   rc = strcmp( t1->comp, t2->comp);
1013   if ( rc != 0 ) return(rc);
1014   rc = strcmp( t1->net,  t2->net );
1015   if ( rc != 0 ) return(rc);
1016   rc = strcmp( t1->loc,  t2->loc );
1017   return(rc);
1018}
1019
1020/*************************************************************
1021 *                   CompareStaDist()                        *
1022 *                                                           *
1023 *  This function can be passed to qsort() and bsearch().    *
1024 *  Compare STA distances (from epicenter); closest first.   *
1025 *************************************************************/
1026int CompareStaDist( const void *s1, const void *s2 )
1027{
1028   STA *t1 = (STA *) s1;
1029   STA *t2 = (STA *) s2;
1030
1031   if (t1->dist < t2->dist)
1032     return -1;
1033   else if (t1->dist > t2->dist)
1034     return 1;
1035
1036   return 0;
1037}
1038
1039
1040/****************** Locally Used Functions Below Here *********************/
1041
1042/*************************************************************
1043 *                   CompareDoubles()                        *
1044 *                                                           *
1045 *  This function can be passed to qsort() and bsearch().    *
1046 *  Compare two doubles for ordering in increasing order     *
1047 *************************************************************/
1048int CompareDoubles( const void *s1, const void *s2 )
1049{
1050  double *d1 = (double *)s1;
1051  double *d2 = (double *)s2;
1052
1053  if (*d1 < *d2)
1054     return -1;
1055  else if (*d1 > *d2)
1056    return 1;
1057 
1058  return 0;
1059}
1060
1061/****************** Locally Used Functions Below Here *********************/
1062
1063/*
1064 * GetRespPZ: get response function (as poles and zeros) for the
1065 *   instrument sta, comp, net from the location specified in GMPARAMS.
1066 *   Response is filed in static responseStruct gScnPZ for later use
1067 *   by makeGM().
1068 *  Resturns: 0 on success
1069 *           -1 on fatal errors
1070 *           +1 on non-fatal errors
1071 */
1072static int getRespPZ(char *sta, char *comp, char *net, char *loc, GMPARAMS *pgmParams, 
1073              EVENT *pEvt)
1074{
1075  char respfile[PATH_MAX];
1076  int len, rc;
1077 
1078  /* Make sure the respone structure doesn't have any old memory allocated */
1079  cleanPZ(&gScnPZ);
1080  memset(respfile, 0, PATH_MAX);
1081 
1082  switch(pgmParams->respSource)
1083  {
1084#ifdef EWDB
1085  case GM_RS_EWDB:
1086    rc = getRespEWDB(sta, comp, net, &gScnPZ);
1087    return rc;
1088    break;
1089#endif
1090  case GM_RS_FILE:  /* fill in directory; rest of action is after switch() */
1091    strcpy(respfile, pgmParams->respDir);
1092    break;
1093  default:
1094    logit("et", "getRespPZ: unknown response source <%d>\n", 
1095          pgmParams->respSource);
1096    return -1;
1097  }
1098 
1099  /* Continue preparing respfile for PZ access: */
1100  len = strlen(respfile);
1101  if (respfile[len - 1] != '/' || respfile[len - 1] != '\\')
1102  {
1103    strcat(respfile, "/");
1104    len++;
1105  }
1106  if (fmtFilename(sta, comp, net, loc, respfile, PATH_MAX - len, 
1107                  pgmParams->respNameFormat) < 0)
1108  {
1109    logit("et", "getRespPZ: error formating respfile <%s>\n", respfile);
1110    return -1;
1111  }
1112  if ( (rc = readPZ(respfile, &gScnPZ)) < 0)
1113  {
1114    switch(rc)
1115    {
1116    case -1:
1117      logit("et", "getRespPZ: out of memory building pole-zero struct\n");
1118      break;
1119    case -2:
1120      logit("et", "getRespPZ: error parsing pole-zero file <%s> for <%s.%s.%s.%s>\n", 
1121            respfile, sta, comp, net, loc);
1122      return +1;  /* Maybe the other resp files will work */
1123      break;
1124    case -3:
1125      logit("et", "getRespPZ: invalid arguments passed to readPZ\n");
1126      break;
1127    case -4:
1128      logit("et", "getRespPZ: error opening pole-zero file <%s> for <%s.%s.%s.%s>\n", 
1129            respfile, sta, comp, net, loc);
1130      return +1;  /* Maybe the other resp files will work */
1131      break;
1132    default:
1133      logit("et", "getRespPZ: unknown error <%d> from readPZ(%s)\n", rc,
1134            respfile);
1135    }
1136    return -1;
1137  }
1138  /* all done! */
1139  return 0;
1140}
1141
1142/*
1143 * endEvent: Terminate the event and clean up. Calls cleanup routines
1144 *           for those data sinks that need it, then cleans up the
1145 *           EVENT structure in preparation for the next event.
1146 */
1147void endEvent(EVENT *pEvt, GMPARAMS *pgmParams)
1148{
1149  int i;
1150  COMP1 *pC1;
1151
1152  /* Call `endEvent' for those data sinks that may need it: EWDB? */
1153 
1154  /* Clean up the EVENT structure in preparation for next event */
1155  for (i = 0; i < pgmParams->maxSta; i++)
1156  {
1157    while ( (pC1 = pEvt->Sta[i].comp) != (COMP1 *)NULL)
1158    {
1159      pEvt->Sta[i].comp = pC1->next;
1160      free(pC1);
1161    }
1162  }
1163  pEvt->numSta = 0;
1164  memset(pEvt->Sta, 0, pgmParams->maxSta * sizeof(STA));
1165
1166  return;
1167}
1168
1169/*
1170 * procArc: parse the summary line from a hyp2000 archive message into
1171 *          the EVENT structure
1172 *   Returns: 0 on success
1173 *           -1 on parse errors
1174 */
1175int procArc( char *msg, EVENT *pEvt, MSG_LOGO iLogo, MSG_LOGO fLogo)
1176{
1177  struct Hsum Sum;          /* Hyp2000 summary data    */
1178  char shdw = '\0';   /* to store a fake shadow card   */
1179
1180  /* Process the hypocenter card (1st line of msg) */
1181  if ( read_hyp( msg, &shdw, &Sum ) < 0 )
1182  {  /* Should we send an error message? */
1183    logit( "t", "procArc: error parsing HYPO summary message\n" );
1184    return( -1 );
1185  }
1186  sprintf(pEvt->eventId, "%lu", Sum.qid);
1187  sprintf(pEvt->authorId, "%03u%03u%03u:%03u%03u%03u", iLogo.type, iLogo.mod,
1188          iLogo.instid, fLogo.type, fLogo.mod, fLogo.instid);
1189  pEvt->lon = (double)Sum.lon;
1190  pEvt->lat = (double)Sum.lat;
1191  pEvt->depth = (double)Sum.z;
1192  pEvt->origin_time = Sum.ot - GSEC1970; /* 1600 to 1970, in seconds */
1193
1194  return 0;
1195}
1196
1197/*
1198 * addSRaddSR2list: add a new spectra response period to the list of those to be
1199 *           processed
1200 *   Returns: 0 on success
1201 *            1 if maximum size of list exceeded
1202 */
1203int addSR2list( double newSR ) {
1204        if ( NSP == SM_MAX_RSA ) {
1205                logit( "w", "Maximum of %d spectral periods allowed; ignoring '%lf'\n", 
1206                        SM_MAX_RSA, newSR );
1207                return 1;
1208        }
1209        spectPer[NSP++] = newSR;
1210        return 0;
1211}
1212
1213
1214/*
1215 * checkSRlist: if list of spectra responses to be processed is empty,
1216 *           fill w/ default list
1217 *   Returns: 0 on success
1218 *            1 if maximum size of list exceeded
1219 */
1220void checkSRlist() {
1221        if ( NSP == 0 ) {
1222                spectPer[NSP++] = 0.3;
1223                spectPer[NSP++] = 1.0;
1224                spectPer[NSP++] = 3.0;
1225        }
1226}
Note: See TracBrowser for help on using the repository browser.