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

Revision 7593, 34.8 KB checked in by paulf, 11 months ago (diff)

added ChannelNumberMap? to gmew to map numbers to components....before this 2 and 3 were default mapped to N and E respectively

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