source: trunk/src/seismic_processing/pick_wcatwc/pick_wcatwc.c @ 4963

Revision 4963, 41.9 KB checked in by paulf, 8 years ago (diff)

patched pick_wcatwc from latest SVN release of EB

Line 
1      /*****************************************************************
2       *                         pick_wcatwc.c                         *
3       *                                                               *
4       *  This program is a conversion of the P-picker used at         *
5       *  at the West Coast/Alaska Tsunami Warning Center since the    *
6       *  early 1980s.  The basic algorithm for the P-picker was       *
7       *  written by Veith in 1978 and is described in Technical Note  *
8       *  1/78, "Seismic Signal Detection Algorithm" by Karl F. Veith, *
9       *  Teledyne GeoTech.  This technique is an adaptation of the    *
10       *  original Rex Allen picker.  Several variations to the        *
11       *  Veith method have been added over the years, such as         *
12       *  magnitude determinations and auto-calibrations.              *
13       *                                                               *
14       *  In general, the signal must go through                       *
15       *  three processing stages before a pick is declared.  The      *
16       *  first stage is a simple test which looks for higher than     *
17       *  normal signal amplitudes.  Actually, it is not the amplitude *
18       *  that is tested, but the cummulative difference between       *
19       *  samples (the MDF).  This is compared to the background MDF   *
20       *  every sample.                                                *
21       *  The background MDF is computed with a moving averages        *
22       *  technique.  Every LTASeconds the average MDF is computed and *
23       *  this is averaged with the existing MDF to produce a new MDF. *
24       *  The second phase of signal processing consists of two tests. *
25       *  Test 1 specifies that the MDF must exceed a trigger          *
26       *  threshold for a period of time after Phase 1 is passed twice *
27       *  as often as it is under the threshold.  Test 2 states        *
28       *  that the true amplitude must exceed the signal-to-noise      *
29       *  ratio at some time while Test 1 is being conducted.  Phase 3 *
30       *  consists of three tests.  Test 1 requires MDF greater than   *
31       *  the trigger threshold at least 6 times in opposing directions*
32       *  (i.e. it must see three full cycles of signal).  Test 2      *
33       *  requires the frequency of the above oscillations to be       *
34       *  greater than MinFreq.  Test 3 requires the MDF to be above   *
35       *  the trigger for at least half the time it takes to pass the  *
36       *  above two tests.  If these three Phases are passed, a P-pick *
37       *  is declared.                                                 *
38       *                                                               *
39       *  At this point, the original MDF start time is saved as       *
40       *  the event's P-time.  The picker continues to evaluate the    *
41       *  signal for magnitude data.  Both Mb and Ml's are computed    *
42       *  here along with sine-wave calibration processing.  (Sine-wave*
43       *  cals are automatically fed into WC/ATWC network analog       *
44       *  stations twice a day.  These cals are analyzed to update     *
45       *  station magnification at each cal.)  After the magnitude     *
46       *  processing has finished, evaluation on this channel ceases   *
47       *  and all variables are reinitialized.                         *
48       *                                                               *
49       *  The picker evaluates data at the incoming sample rate.  No   *
50       *  re-sampling is necessary.  However, the picker does not work *
51       *  well on anything but short-period data.  So, broadband data  *
52       *  must be filtered before P-picking.  Filtering is provided by *
53       *  this program.   Response for the default filter has been     *
54       *  computed and is used to convert to ground motion.  If other  *
55       *  filter parameters are used, the response must be calculated  *
56       *  and used in the ground motion calculations.  The raw         *
57       *  broadband data is still used, though, for Mwp computations.  *
58       *                                                               *
59       *  In this re-write of the P-picker, first-motion determination *
60       *  is attempted.  Coda magnitude information is not computed    *
61       *  here.  Gap interpolation is performed as in pick_ew.         *
62       *                                                               *
63       *  One other feature was added in this implementaion of the     *
64       *  Veith picker.  That is an alarm feature.  This feature looks *
65       *  for potentially large quakes recorded at any site.  When     *
66       *  the alarm criteria is passed, an alarm messages is fed       *
67       *  into the alarm ring.                                         *
68       *                                                               *
69       *  A no data alarm was also added to pick_wcatwc in 8/01.  An   *
70       *  alarm time is specified in the .d file.  If this time passes *
71       *  with data data arriving in the ring, an alarm message will   *
72       *  be issued to the alarm ring.                                 *
73       *                                                               *
74       *  The structure of this program was modeled after pick_ew.     *
75       *  Many of the functions which are applicable to any P-picker   *
76       *  were also taken from pick_ew.                                *
77       *                                                               *
78       *  May, 2009: Added Matlab neural net initialization for P-pick *
79       *  November, 2007: Expanded no data or bad time alarms          *
80       *  December, 2004: Added multi-station alarm feature and        *
81       *                  PStrength computations.                      *
82       *  October, 2005: Added options to run this module in           *
83       *                 conjunction with ATPlayer.                    *
84       *                                                               *
85       *  Written by Paul Whitmore, December, 2000                     *
86       *                                                               *
87       ****************************************************************/
88           
89#include <stdio.h>
90#include <stdlib.h>
91#include <string.h>
92#include <sys/types.h>
93#include <time.h>
94#ifdef _WINNT
95#include <malloc.h> // jmc
96#include <mat.h>
97#endif
98
99#include <earthworm.h>
100#include <transport.h>
101#include <trace_buf.h>
102#include <swap.h>
103#include <kom.h>
104#include "pick_wcatwc.h"
105#ifdef _USE_NEURALNET_
106#include <libNeuralNet.h>
107#endif
108
109      /***********************************************************
110       *              The main program starts here.              *
111       *                                                         *
112       *  Argument:                                              *
113       *     argv[1] = Name of picker configuration file         *
114       ***********************************************************/
115
116/* Global Variables
117   ****************/
118GPARM                   Gparm;           /* Configuration file parameters */
119STATION                 *StaArray;       /* Station array                 */
120char                    *WaveBuf;        /* Pointer to waveform buffer    */
121ALARMSTRUCT             *AS;             /* Pointer to Alarm buffer       */
122int                      Nsta;           /* Number of stations in list    */
123
124int main( int argc, char **argv )
125{
126   char          *configfile;     /* Pointer to name of config file */
127   static double dLastEndTime;    /* End time of last packet */
128   EWH           Ewh;             /* Parameters from earthworm.h */
129   MSG_LOGO      getlogo;         /* Logo of requested waveforms */
130   MSG_LOGO      hrtlogo;         /* Logo of outgoing heartbeats */
131   int           i, iTemp;        /* Loop counter */
132   long          InBufl;          /* Maximum message size in bytes */
133   static int    iNoDataAlarmIssued=0; /* 1, 2... after "no data" alarm issued*/
134   static int    iNumRegions;     /* Number of alarm regions */
135   char          line[40];        /* Heartbeat message */
136   int           lineLen;         /* Length of heartbeat message */
137   time_t        lLastData;       /* 1/1/70 time of last data provided */
138   MSG_LOGO      logo;            /* Logo of retrieved msg */
139   long          MsgLen;          /* Size of retrieved message */
140   pid_t         myPid;           /* Process id of this process */
141   long          RawBufl;         /* Raw data buffer size (for Mwp) */
142   time_t        then;            /* Previous heartbeat time */
143   time_t        diffTime;
144   TRACE2_HEADER *WaveHead;       /* Pointer to waveform header */
145   long          *WaveLong;       /* Long pointer to waveform data */
146   long          WaveRaw[MAX_TRACEBUF_SIZ];/* Unfiltered waveform data */
147   short         *WaveShort;      /* Short pointer to waveform data */
148   static int    MclInitFlag = 0;
149   char *paramdir;
150   char  FullTablePath[512];
151   
152/* Call the MCR and library initialization functions (Added 5/09 dln) */
153#ifdef _USE_NEURALNET_
154   if (MclInitFlag == 0 )
155   {
156      MclInitFlag = 1;
157      if ( !mclInitializeApplication( NULL, 0 ) )
158      {
159         fprintf( stderr, "Could not initialize the matlab NN application.\n" );
160         exit(1);               
161      }
162      if ( !libNeuralNetInitialize() )
163      {
164         fprintf( stderr, "Could not initialize the matlab nn library.\n" );
165         logit( "", "Could not initialize the matlab nn library.\n" );
166         logit( "" , " errno: %d  \n", errno );
167      }
168   } 
169#endif
170
171/* Check command line arguments
172   ****************************/
173   if ( argc != 2 )
174   {
175      fprintf( stderr, "Usage: pick_wcatwc <configfile>\n" );
176      return -1;
177   }
178   configfile = argv[1];
179
180/* Environmental Parameter to Read Stations
181   ****************************************/
182   FullTablePath[0] = 0;
183#ifdef _WINNT
184   paramdir = getenv( "EW_PARAMS" ); 
185   if ( paramdir == (char *)NULL ) return(-1);
186   strcpy( FullTablePath, paramdir  );
187   strcat( FullTablePath, "\\");
188#endif
189   strcat( FullTablePath, configfile );  // for unix just use relative dir to EW_PARAMS since that is where all is run
190
191/* Get parameters from the configuration files
192   *******************************************/
193   if ( GetConfig( FullTablePath, &Gparm ) == -1 )
194   {
195      fprintf( stderr, "pick_wcatwc: GetConfig() failed. Exiting.\n" );
196      return -1;
197   }
198
199/* Look up info in the earthworm.h tables
200   **************************************/
201   if ( GetEwh( &Ewh ) < 0 )
202   {
203      fprintf( stderr, "pick_wcatwc: GetEwh() failed. Exiting.\n" );
204      return -1;
205   }
206
207#ifdef _WINNT
208   if( SetConsoleCtrlHandler( (PHANDLER_ROUTINE) CtrlHandler, TRUE ) ) 
209      fprintf( stderr, "The Control Handler is installed. \n" ); 
210   else
211      fprintf( stderr, "ERROR: Could not set control handler \n" ); 
212#endif
213
214/* Specify logos of incoming waveforms and outgoing heartbeats
215   ***********************************************************/
216   getlogo.instid = Ewh.GetThisInstId;
217   getlogo.mod    = Ewh.GetThisModId;
218   getlogo.type   = Ewh.TypeWaveform;
219   hrtlogo.instid = Ewh.MyInstId;
220   hrtlogo.mod    = Gparm.MyModId;
221   hrtlogo.type   = Ewh.TypeHeartBeat;
222
223/* Initialize name of log-file & open it
224   *************************************/
225   logit_init( configfile, Gparm.MyModId, 256, 1 );
226
227/* Get our own pid for restart purposes
228   ************************************/
229   myPid = getpid();
230   if ( myPid == -1 )
231   {
232      logit( "e", "pick_wcatwc: Can't get my pid. Exiting.\n" );
233      return -1;
234   }
235
236 
237/* Log the configuration parameters
238   ********************************/
239   LogConfig( &Gparm );
240   
241/* Read in Alarm parameters if turned on
242   *************************************/
243   if ( Gparm.TwoStnAlarmOn == 1 )
244      if ( ReadAlarmParams( &AS, &iNumRegions, Gparm.TwoStnAlarmFile ) == 0 )
245      {
246         logit( "e", "pick_wcatwc: can't read alarm file %s; Exiting.\n",
247                Gparm.TwoStnAlarmFile );
248         return -1;
249      }             
250         
251/* Allocate the waveform buffer
252   ****************************/
253   InBufl = MAX_TRACEBUF_SIZ*2 + sizeof(long)*(Gparm.MaxGap-1);
254   WaveBuf = (char *) malloc( (size_t) InBufl );
255   if ( WaveBuf == NULL )
256   {
257      logit( "et", "pick_wcatwc: Cannot allocate waveform buffer\n" );
258      free( AS );
259      return -1;
260   }
261
262/* Point to header and data portions of waveform message
263   *****************************************************/
264   WaveHead  = (TRACE2_HEADER *) WaveBuf;
265   WaveLong  = (long *) (WaveBuf + sizeof(TRACE2_HEADER));
266   WaveShort = (short *) (WaveBuf + sizeof(TRACE2_HEADER));
267
268/* Read the station list and return the number of stations found.
269   Allocate the station list array.
270   **************************************************************/
271   if ( ReadStationList( &StaArray, &Nsta, Gparm.StaFile, Gparm.StaDataFile,
272                         Gparm.ResponseFile, MAX_STATIONS, 1 ) == -1 )
273   {
274      logit( "", "pick_wcatwc: ReadStationList() failed. Exiting.\n" );
275      free( WaveBuf );
276      free( AS );
277      return -1;
278   }
279   if ( Nsta == 0 )
280   {
281      logit( "et", "pick_wcatwc: Empty station list. Exiting." );
282      free( AS );
283      free( WaveBuf );
284      free( StaArray );
285      return -1;
286   }
287   logit( "t", "pick_wcatwc: Picking %d stations", Nsta );
288                 
289/* Convert iAlarmStatus to iAlarmPage and iAlarmSpeak
290   **************************************************/
291   for ( i=0; i<Nsta; i++ )
292   {
293      StaArray[i].iAlarmSpeak = 0;
294      StaArray[i].iAlarmPage = 0;
295      if ( StaArray[i].iAlarmStatus == 1 )
296      {
297         StaArray[i].iAlarmPage = 1;       
298         StaArray[i].iAlarmSpeak = 0;     
299      }
300      if ( StaArray[i].iAlarmStatus == 2 )
301      {
302         StaArray[i].iAlarmPage = 0;       
303         StaArray[i].iAlarmSpeak = 1;     
304         StaArray[i].iAlarmStatus = 1;
305      }
306      if ( StaArray[i].iAlarmStatus == 3 )
307      {
308         StaArray[i].iAlarmPage = 1;       
309         StaArray[i].iAlarmSpeak = 1;     
310         StaArray[i].iAlarmStatus = 1;
311      }
312   }
313
314/* Log the station list
315   ********************/
316   LogStaListP( StaArray, Nsta );
317
318/* Attach to existing transport rings
319   **********************************/
320   if ( Gparm.OutKey != Gparm.InKey )
321   {
322      tport_attach( &Gparm.InRegion,  Gparm.InKey );
323      tport_attach( &Gparm.OutRegion, Gparm.OutKey );
324      tport_attach( &Gparm.AlarmRegion, Gparm.AlarmKey );
325   }
326   else
327   {
328      tport_attach( &Gparm.InRegion, Gparm.InKey );
329      Gparm.OutRegion = Gparm.InRegion;
330   }
331
332/* Flush the input ring
333   ********************/
334   while ( tport_getmsg( &Gparm.InRegion, &getlogo, 1, &logo, &MsgLen,
335                          WaveBuf, MAX_TRACEBUF_SIZ) != GET_NONE );
336
337/* Send 1st heartbeat to the transport ring
338   ****************************************/
339   time( &then );
340   sprintf( line, "%ld %d\n", (long) then, myPid );
341   lineLen = strlen( line );
342   if ( tport_putmsg( &Gparm.InRegion, &hrtlogo, lineLen, line ) !=
343        PUT_OK )
344   {
345      logit( "et", "pick_wcatwc: Error sending 1st heartbeat. Exiting." );
346      if ( Gparm.OutKey != Gparm.InKey )
347      {
348         tport_detach( &Gparm.InRegion );
349         tport_detach( &Gparm.OutRegion );
350         tport_detach( &Gparm.AlarmRegion );
351      }
352      else
353         tport_detach( &Gparm.InRegion );
354      free( AS );
355      free( WaveBuf );
356      free( StaArray );
357      return 0;
358   }                                                                 
359   lLastData = then;
360
361/* Loop to read waveform messages and invoke the picker
362   ****************************************************/
363   dLastEndTime = 0.;   
364
365   while ( tport_getflag( &Gparm.InRegion ) != TERMINATE )
366   {
367      char    type[3];
368      static  STATION *Sta;     /* Pointer to the station being processed */
369      int     rc;               /* Return code from tport_getmsg() */
370      static  time_t  now;      /* Current time */
371      double  GapSizeD;         /* Number of missing samples (double) */
372      long    GapSize;          /* Number of missing samples (integer) */
373
374/* Send a heartbeat to the transport ring
375   **************************************/
376      time( &now );
377      if ( (now - then) >= Gparm.HeartbeatInt )
378      { 
379/* Have we received data through ring without AlarmTime seconds elapsing?
380   **********************************************************************/
381         then = now;
382         sprintf( line, "%ld %d\n", (long) now, myPid );
383         lineLen = strlen( line );
384         if ( tport_putmsg( &Gparm.InRegion, &hrtlogo, lineLen, line ) !=
385              PUT_OK )
386         {
387            logit( "et", "pick_wcatwc: Error sending heartbeat. Exiting." );
388            break;
389         }                           
390      }                     
391/* Send No Data Alarm if the clock was set way back in time, or no data
392   has been received in a long time.
393   ********************************************************************/
394      if ( strlen( Gparm.ATPLineupFileBB ) < 3 ) 
395          {
396                  // JFP -- Needed to do this to get of Warning.
397             diffTime = now - lLastData;
398                 if( diffTime < 0 ) diffTime = diffTime * -1;
399         if ( diffTime > Gparm.AlarmTime )
400                 // ----------------------------------------------
401         {
402            if ( iNoDataAlarmIssued == 0 )     
403            {
404               logit( "t", "No Data Alarm Activated\n" );
405               ReportAlarm( StaArray, Gparm.MyModId, Gparm.AlarmRegion,
406                Ewh.TypeAlarm, Ewh.MyInstId, 5, "", 0 );
407            }
408            iNoDataAlarmIssued++;
409         }
410          }
411/* Get a waveform buffer from transport region
412   *******************************************/
413      rc = tport_getmsg( &Gparm.InRegion, &getlogo, 1, &logo, &MsgLen,
414                         WaveBuf, MAX_TRACEBUF_SIZ);
415
416      if ( rc == GET_NONE )
417      {
418         sleep_ew( 200 );
419         continue;
420      }
421
422      if ( rc == GET_NOTRACK )
423         logit( "et", "pick_wcatwc: Tracking error.\n");
424
425      if ( rc == GET_MISS_LAPPED )
426         logit( "et", "pick_wcatwc: Got lapped on the ring.\n");
427
428      if ( rc == GET_MISS_SEQGAP )
429         logit( "et", "pick_wcatwc: Gap in sequence numbers.\n");
430
431      if ( rc == GET_MISS )
432         logit( "et", "pick_wcatwc: Missed messages.\n");
433
434      if ( rc == GET_TOOBIG )
435      {
436         logit( "et", "pick_wcatwc: Retrieved message too big (%d) for msg.\n",
437                MsgLen );
438         continue;
439      }
440
441/* If necessary, swap bytes in the message
442   ***************************************/
443// JFP
444#ifdef _ADD_TRACE2_HEADER_
445      if ( WaveMsg2MakeLocal( WaveHead ) < 0 )
446#else
447      if ( WaveMsgMakeLocal( WaveHead ) < 0 )
448#endif
449// ---
450      {
451         logit( "et", "pick_wcatwc: Unknown waveform type.\n" );
452         continue;
453      }
454
455#ifdef _WINNT
456/* If we are in the ATPlayer version of pick_, see if we should re-init
457   ********************************************************************/
458      if ( strlen( Gparm.ATPLineupFileBB ) > 2 )     /* Then we are */
459         if ( fabs( WaveHead->starttime-(int) dLastEndTime ) > 1800 ) /* Big gap */
460         {
461            free( StaArray );                                 
462            Nsta = ReadLineupFile( Gparm.ATPLineupFileBB, &StaArray, 
463                                   MAX_STATIONS );
464            logit( "", "reset StaArray Nsta = %d\n", Nsta );     
465            if ( Nsta < 2 )
466            {
467               logit( "", "Bad Lineup File read %s\n", Gparm.ATPLineupFileBB );
468               continue;
469            }       
470            for ( i=0; i<Nsta; i++ )
471            {
472               free( StaArray[i].plRawData );
473               free( StaArray[i].pdRawDispData );
474               free( StaArray[i].pdRawIDispData );
475               free( StaArray[i].plPickBuf );
476               free( StaArray[i].plPickBufRaw );
477            }
478         }       
479#endif
480
481/* Look up SCN number in the station list
482   **************************************/
483      Sta = NULL;                                                                 
484      for ( i=0; i<Nsta; i++ )
485         if ( !strcmp( WaveHead->sta,  StaArray[i].szStation ) &&
486              !strcmp( WaveHead->chan, StaArray[i].szChannel ) &&
487              !strcmp( WaveHead->loc,  StaArray[i].szLocation ) &&
488              !strcmp( WaveHead->net,  StaArray[i].szNetID ) )
489         {
490            Sta = (STATION *) &StaArray[i];
491            break;
492         }
493
494      if ( Sta == NULL )      /* SCN not found */
495         continue;
496                 
497/* Check if the time stamp is reasonable.  If it is ahead of the present
498   1/1/70 time, it is not reasonable. (+1. to account for int).
499   *********************************************************************/
500      if ( WaveHead->endtime > (double) now+1. )
501      {
502         logit( "t", "%s %s %s %s endtime (%lf) ahead of present (%ld)\n",
503                Sta->szStation, Sta->szChannel, Sta->szNetID, Sta->szLocation, WaveHead->endtime, (long) now );
504         continue;
505      }
506
507/* Do this the first time we get a message with this SCN
508   *****************************************************/
509      if ( Sta->iFirst == 1 )
510      {
511         Sta->lPickIndex = 1;
512         Sta->dSampRate = WaveHead->samprate;
513         Sta->dEndTime = WaveHead->endtime;
514         Sta->dStartTime = WaveHead->starttime;
515         Sta->dDataEndTime = WaveHead->endtime;
516         Sta->iFirst = 0;
517         ResetFilter( Sta );
518                 
519/* Allocate memory for raw data buffer (assume samprate will not change for
520   each station and that samprate is greater than 1)
521   *************************************************/
522         RawBufl = sizeof (long) * (long) (WaveHead->samprate+0.001) *
523                           Gparm.MwpSeconds;
524         if ( RawBufl/sizeof (long) > MAXMWPARRAY )
525         {
526            logit( "et", "Mwp array too large for %s\n", Sta->szStation );
527            RawBufl = (MAXMWPARRAY-1) * sizeof (long);
528         }
529                 // JFP -- OVERRUN
530                 // make it MAXMWPARRAY for now.
531         RawBufl = (MAXMWPARRAY-1) * sizeof (long);
532                 // ---
533         Sta->plRawData = (long *) malloc( (size_t) RawBufl );
534         if ( Sta->plRawData == NULL )
535         {
536            logit( "et", "pick_wcatwc: Can't allocate raw data buffer for %s\n",
537                               Sta->szStation );
538            free( AS );
539            free( WaveBuf );              /* Free-up allocated memory */
540            free( StaArray );
541            return -1;
542         }
543
544         RawBufl = (MAX_NN_SAMPS) * sizeof (long);
545         Sta->plPickBuf = (long *) malloc( (size_t) RawBufl );
546         if ( Sta->plPickBuf == NULL )
547         {
548            logit( "et", "pick_wcatwc: Can't allocate pick data buff for %s\n",
549                               Sta->szStation );
550            free( AS );
551            free( WaveBuf );              /* Free-up allocated memory */
552            for ( i=0; i<Nsta; i++ ) 
553               free( StaArray[i].plRawData );
554            free( StaArray );
555            return -1;
556         }
557
558         RawBufl = (MAX_NN_SAMPS) * sizeof (long);
559         Sta->plPickBufRaw = (long *) malloc( (size_t) RawBufl );
560         if ( Sta->plPickBufRaw == NULL )
561         {
562            logit( "et", "pick_wcatwc: Can't allocate pick data buff for %s\n",
563                               Sta->szStation );
564            free( AS );
565            free( WaveBuf );              /* Free-up allocated memory */
566            for ( i=0; i<Nsta; i++ ) free( StaArray[i].plRawData );
567            for ( i=0; i<Nsta; i++ ) free( StaArray[i].plPickBuf );
568            free( StaArray );
569            return -1;
570         }
571                 
572/* Allocate memory for Mwp displacement data buffer
573   ************************************************/
574         RawBufl = sizeof (double) * (long) (WaveHead->samprate+0.001) *
575                           Gparm.MwpSeconds;
576         if ( RawBufl/sizeof (double) > MAXMWPARRAY )
577         {
578            logit( "et", "Mwp array too large for %s\n", Sta->szStation );
579            RawBufl = (MAXMWPARRAY-1) * sizeof (double);
580         }
581                 // JFP -- OVERRUN
582                 // make it MAXMWPARRAY for now.
583         RawBufl = (MAXMWPARRAY-1) * sizeof (long);
584                 // ---
585         Sta->pdRawDispData = (double *) malloc( (size_t) RawBufl );
586         if ( Sta->pdRawDispData == NULL )
587         {
588            logit( "et", "pick_wcatwc: Can't allocate disp. data buff for %s\n",
589                               Sta->szStation );
590            free( AS );
591            free( WaveBuf );              /* Free-up allocated memory */
592            for ( i=0; i<Nsta; i++ ) 
593            {
594               free( StaArray[i].plRawData );
595               free( StaArray[i].plPickBuf );
596               free( StaArray[i].plPickBufRaw );
597            }
598            free( StaArray );
599            return -1;
600         }
601                 
602/* Allocate memory for Mwp integrated displacement data buffer
603   ***********************************************************/
604         Sta->pdRawIDispData = (double *) malloc( (size_t) RawBufl );
605         if ( Sta->pdRawIDispData == NULL )
606         {
607            logit( "et", "pick_wcatwc: Cannot allocate int. disp. data buffer for %s\n",
608                               Sta->szStation );
609            free( AS );
610            free( WaveBuf );              /* Free-up allocated memory */
611            for ( i=0; i<Nsta; i++ ) 
612            {
613               free( StaArray[i].plRawData );
614               free( StaArray[i].pdRawDispData );
615               free( StaArray[i].plPickBuf );
616               free( StaArray[i].plPickBufRaw );
617            }
618            free( StaArray );
619            return -1;
620         }
621                                 
622                 // JFP -- Allocate for other pointers as well.
623         Sta->plFiltCircBuff = (long *) malloc( 100 * sizeof( long) );
624         Sta->plRawCircBuff = (long *) malloc( 100 * sizeof( long) );
625                 // -------
626
627         dLastEndTime = WaveHead->endtime;
628         continue;
629      }
630     
631/* If data is not in order, throw it out
632   *************************************/
633      if ( Sta->dEndTime > WaveHead->starttime+0.001 ) continue;
634      lLastData = now;         /* We got some data, so update last alarm time */
635      iNoDataAlarmIssued = 0;
636
637/* If the samples are shorts, make them longs
638   ******************************************/
639      strcpy( type, WaveHead->datatype );
640      if ( (strcmp( type,"i2" ) == 0) || (strcmp( type,"s2" ) == 0) )
641         for ( i = WaveHead->nsamp-1; i>-1; i-- )
642            WaveLong[i] = (long) WaveShort[i];
643
644/* Save long data as raw, unfiltered data
645   **************************************/
646      for ( i = WaveHead->nsamp - 1; i > -1; i-- ) WaveRaw[i] = WaveLong[i];     
647         
648/* Compute the number of samples since the end of the previous message.
649   If (GapSize == 1), no data has been lost between messages.
650   If (1 < GapSize <= Gparm.MaxGap), data will be interpolated.
651   If (GapSize > Gparm.MaxGap), the picker will go into restart mode.
652   *******************************************************************/
653      GapSizeD = WaveHead->samprate * (WaveHead->starttime - Sta->dEndTime);
654
655      if ( GapSizeD < 0. )          /* Invalid. Time going backwards. */
656         GapSize = 0;
657      else
658         GapSize  = (long) (GapSizeD + 0.5);
659
660/* Interpolate missing samples and prepend them to the current message
661   *******************************************************************/
662      if ( (GapSize > 1) && (GapSize <= Gparm.MaxGap) )
663         Interpolate( Sta, WaveBuf, GapSize );
664
665/* Announce large sample gaps
666   **************************/
667      if ( GapSize > Gparm.MaxGap )
668      {
669         int      lineLen;
670         time_t   errTime;
671         char     errmsg[80];
672         MSG_LOGO logo;
673
674         time( &errTime );
675         sprintf( errmsg,
676               "%ld %d Found %4ld sample gap. Restarting station %-5s%-2s%-3s %s\n",
677               (long)errTime, PK_RESTART, GapSize, Sta->szStation, Sta->szNetID,
678               Sta->szChannel, Sta->szLocation );
679         lineLen = strlen( errmsg );
680         logo.type   = Ewh.TypeError;
681         logo.mod    = Gparm.MyModId;
682         logo.instid = Ewh.MyInstId;
683         tport_putmsg( &Gparm.InRegion, &logo, lineLen, errmsg );
684         if ( Gparm.Debug )
685            logit( "t", "pick_wcatwc: Restarting %-5s%-2s %-3s %s. GapSize = %d\n",
686                     Sta->szStation, Sta->szNetID, Sta->szChannel, Sta->szLocation, GapSize ); 
687      }
688
689/* For big gaps, enter restart mode. In restart mode, calculate
690   moving averages without picking.  Start picking again after a
691   specified number of samples has been processed.
692   *************************************************************/
693      if ( GapSize > Gparm.MaxGap )
694      {
695         Sta->iPickStatus = 1;
696         ResetFilter( Sta );
697         InitVar( Sta );
698         if ( Sta->iAlarmStatus >= 1 ) Sta->iAlarmStatus = 1;
699      }
700         
701/* Short-period bandpass filter the data if specified in station file -
702   Here filter the data without removing DC offset.  This offset is taken out
703   later.  This may occasionally cause pick at data start or after gaps if
704   a station has a large DC offset (filter response to offset).
705   (This is usually necessary for broadband data, not so for SP data.
706   Taper the signal for restarts, but not for continuous signal.).
707   ******************************************************************/
708      if ( Sta->iFiltStatus == 1 )
709         FilterPacket( WaveLong, Sta, WaveHead->nsamp, WaveHead->samprate,
710                       Gparm.LowCutFilter, Gparm.HighCutFilter,
711                       3.*(1./Gparm.LowCutFilter) );
712
713/* Save time and amplitude of the end of the current message
714   *********************************************************/
715      Sta->lEndData = WaveLong[WaveHead->nsamp - 1];
716      Sta->dEndTime = WaveHead->endtime;
717      Sta->dDataEndTime = WaveHead->endtime;
718      Sta->dDataStartTime = WaveHead->starttime;
719      dLastEndTime = WaveHead->endtime;
720         
721/* Run data through MovingAvg function and picker
722   **********************************************/       
723      PickV( Sta, WaveHead->starttime, Gparm.AlarmOn, Gparm.TwoStnAlarmOn,
724             Gparm.AlarmTimeout, Gparm.MinFreq, Gparm.LTASeconds,
725             Gparm.MwpSeconds, Gparm.MwpSigNoise, Gparm.LGSeconds,
726             Gparm.MbCycles, Gparm.dSNLocal, Gparm.dMinFLoc, Gparm.MyModId,
727             Gparm.AlarmRegion, Gparm.OutRegion, Ewh.TypeAlarm, Ewh.TypePickTWC,
728             Ewh.MyInstId, WaveRaw, WaveLong, AS, iNumRegions, 1, &iTemp,
729             Gparm.NeuralNet );
730   }
731
732/* Detach from the ring buffers
733   ****************************/
734   if ( Gparm.OutKey != Gparm.InKey )
735   {
736      tport_detach( &Gparm.InRegion );
737      tport_detach( &Gparm.OutRegion );
738      tport_detach( &Gparm.AlarmRegion );
739   }
740   else
741      tport_detach( &Gparm.InRegion );
742   free( AS );
743   free( WaveBuf );                     /* Free-up allocated memory */
744   for ( i=0; i<Nsta; i++ ) 
745   {
746      free( StaArray[i].plRawData );
747      free( StaArray[i].pdRawDispData );
748      free( StaArray[i].pdRawIDispData );
749      free( StaArray[i].plPickBuf );
750      free( StaArray[i].plPickBufRaw );
751   }
752   
753#ifdef _USE_NEURALNET_
754   libNeuralNetTerminate();    //added 5/5/2009 by dln
755   mclTerminateApplication();  //added 5/5/2009 by dln
756#endif
757
758   free( StaArray );
759   logit( "t", "Termination requested. Exiting.\n" );
760   return 0;
761}
762
763#ifdef _WINNT
764// JFP --
765// Added a new function to control when the program is closed.
766// ------
767BOOL CtrlHandler( DWORD fdwCtrlType ) 
768{ 
769   switch( fdwCtrlType ) 
770   {   
771// Handle the CTRL-C signal.
772      case CTRL_C_EVENT:       // printf( "Ctrl-C event\n\n" );
773// CTRL-CLOSE: confirm that the user wants to exit.
774      case CTRL_CLOSE_EVENT:   // printf( "Ctrl-Close event\n\n" );
775// Pass other signals to the next handler.
776      case CTRL_BREAK_EVENT:   // printf( "Ctrl-Break event\n\n" );
777#if 0
778/* Detach from the ring buffers */
779/* **************************   */
780          if ( Gparm.OutKey != Gparm.InKey )
781          {
782             tport_detach( &Gparm.InRegion );
783             tport_detach( &Gparm.OutRegion );
784             tport_detach( &Gparm.AlarmRegion );
785          }
786          else
787             tport_detach( &Gparm.InRegion );
788             free( WaveBuf );                     /* Free-up allocated memory */
789             free( AS );
790             for ( i=0; i<Nsta; i++ )
791             {
792                free( StaArray[i].plRawData );
793                free( StaArray[i].pdRawDispData );
794                free( StaArray[i].pdRawIDispData );
795                free( StaArray[i].plPickBuf );
796                free( StaArray[i].plPickBufRaw );
797             }
798   
799#ifdef _USE_NEURALNET_
800          libNeuralNetTerminate();    //added 5/5/2009 by dln
801          mclTerminateApplication();  //added 5/5/2009 by dln
802#endif
803          free( StaArray );
804#endif
805          logit( "t", "Termination requested. Exiting.\n" );
806          exit( 0 );
807          return( TRUE ); 
808   
809      case CTRL_LOGOFF_EVENT:
810         Beep( 1000, 200 ); 
811         printf( "Ctrl-Logoff event\n\n" );
812         return FALSE; 
813 
814      case CTRL_SHUTDOWN_EVENT:
815         Beep( 750, 500 ); 
816         printf( "Ctrl-Shutdown event\n\n" );
817         return FALSE; 
818
819      default: 
820         return FALSE; 
821   } 
822} 
823#endif
824
825      /*******************************************************
826       *                      GetEwh()                       *
827       *                                                     *
828       *      Get parameters from the earthworm.d file.      *
829       *******************************************************/
830
831int GetEwh( EWH *Ewh )
832{
833   if ( GetLocalInst( &Ewh->MyInstId ) != 0 )
834   {
835      fprintf( stderr, "pick_wcatwc: Error getting MyInstId.\n" );
836      return -1;
837   }
838
839   if ( GetInst( "INST_WILDCARD", &Ewh->GetThisInstId ) != 0 )
840   {
841      fprintf( stderr, "pick_wcatwc: Error getting GetThisInstId.\n" );
842      return -2;
843   }
844   if ( GetModId( "MOD_WILDCARD", &Ewh->GetThisModId ) != 0 )
845   {
846      fprintf( stderr, "pick_wcatwc: Error getting GetThisModId.\n" );
847      return -3;
848   }
849   if ( GetType( "TYPE_HEARTBEAT", &Ewh->TypeHeartBeat ) != 0 )
850   {
851      fprintf( stderr, "pick_wcatwc: Error getting TypeHeartbeat.\n" );
852      return -4;
853   }
854   if ( GetType( "TYPE_ERROR", &Ewh->TypeError ) != 0 )
855   {
856      fprintf( stderr, "pick_wcatwc: Error getting TypeError.\n" );
857      return -5;
858   }
859   if ( GetType( "TYPE_ALARM", &Ewh->TypeAlarm ) != 0 )
860   {
861      fprintf( stderr, "pick_wcatwc: Error getting TypeAlarm.\n" );
862      return -6;
863   }
864   if ( GetType( "TYPE_PICKTWC", &Ewh->TypePickTWC ) != 0 )
865   {
866      fprintf( stderr, "pick_wcatwc: Error getting TypePickTWC.\n" );
867      return -7;
868   }
869   // if ( GetType( "TYPE_TRACEBUF", &Ewh->TypeWaveform ) != 0 )
870   if ( GetType( "TYPE_TRACEBUF2", &Ewh->TypeWaveform ) != 0 )
871   {
872      fprintf( stderr, "pick_wcatwc: Error getting TYPE_TRACEBUF.\n" );
873      return -8;
874   }
875   if ( GetType( "TYPE_PICK_GLOBAL", &Ewh->TypePickGlobal ) != 0 )
876   {
877      fprintf( stderr, "pick_wcatwc: Error getting TypePickGlobal.\n" );
878      return -9;
879   }
880   return 0;
881}
882
883   /*************************************************************
884    *                       Interpolate()                       *
885    *                                                           *
886    *  Interpolate samples and insert them at the beginning of  *
887    *  the waveform.                                            *
888    *************************************************************/
889
890void Interpolate( STATION *Sta, char *WaveBuf, int GapSize )
891{
892   int      i;
893   int      j = 0;
894   int      nInterp = GapSize - 1;
895// JFP
896#ifdef _ADD_TRACE2_HEADER_
897   TRACE2_HEADER *WaveHead  = (TRACE2_HEADER *) WaveBuf;
898   long         *WaveLong  = (long *) (WaveBuf + sizeof(TRACE2_HEADER));
899#else
900   TRACE_HEADER *WaveHead  = (TRACE_HEADER *) WaveBuf;
901   long         *WaveLong  = (long *) (WaveBuf + sizeof(TRACE_HEADER));
902#endif
903// ---
904   double   SampleInterval = 1. / WaveHead->samprate;
905   double   delta = (double)(WaveLong[0] - Sta->lEndData) / GapSize;
906
907   for ( i = WaveHead->nsamp - 1; i >= 0; i-- )
908      WaveLong[i + nInterp] = WaveLong[i];
909
910   for ( i = 0; i < nInterp; i++ )
911      WaveLong[i] = (long) (Sta->lEndData + (++j * delta) + 0.5);
912
913   WaveHead->nsamp += nInterp;
914   WaveHead->starttime = Sta->dEndTime + SampleInterval;
915}
916
917 /***********************************************************************
918  *                             LogStaListP()                           *
919  *                                                                     *
920  *                         Log the station list                        *
921  ***********************************************************************/
922
923void LogStaListP( STATION *Sta, int Nsta )
924{
925   int i;
926
927   logit( "", "\nStation List:\n" );
928   for ( i = 0; i < Nsta; i++ )
929   {
930      logit( "", "%4s",     Sta[i].szStation );
931      logit( "", " %3s",    Sta[i].szChannel );
932      logit( "", " %2s",    Sta[i].szNetID );
933      logit( "", " %2s",    Sta[i].szLocation );
934      logit( "", " %1d",    Sta[i].iPickStatus );
935      logit( "", " %1d",    Sta[i].iFiltStatus );
936      logit( "", " %2d",    Sta[i].iSignalToNoise );
937      logit( "", " %1d",    Sta[i].iAlarmStatus );
938      logit( "", " %10.8lf",Sta[i].dAlarmAmp );
939      logit( "", " %6.1lf", Sta[i].dAlarmDur );
940      logit( "", " %6.3lf", Sta[i].dAlarmMinFreq );
941      logit( "", " %1d",    Sta[i].iComputeMwp );
942      logit( "", "\n" );
943   }
944   logit( "", "\n" );
945}
946
947   /*************************************************************
948    *                     ReadAlarmParams()                     *
949    *                                                           *
950    *  Fill up AlarmStructur with parameters from control file. *
951    *                                                           *
952    *  Variables:                                               *
953    *   pAS             Alarm Structure                         *
954    *   piNumRegs       Number of alarm regions                 *
955    *   pszFile         Alarm parameter file                    *
956    *                                                           *
957    *  Return:                                                  *
958    *   1=ok, 0=problem                                         *
959    *                                                           *
960    *************************************************************/
961
962int ReadAlarmParams( ALARMSTRUCT **pAS, int *piNumRegs, char *pszFile )
963{
964   int     i, ii, j, k, nfiles;
965   long    InBufl;                    /* buffer size */
966   static  ALARMSTRUCT *ptAS;         /* temp pointer */
967   
968/* Open file and read if open ok */   
969   nfiles = k_open( pszFile );
970   if ( nfiles == 0 )
971   {
972      logit( "", "Read problem in ReadAlarmParams (%s)\n", pszFile );
973      return( 0 );
974   }
975   
976   while ( nfiles > 0 )          /* While there are config files open */
977   {
978      while ( k_rd() )           /* Read next line */
979      {
980         char *com;
981         char *str;
982
983         com = k_str();          /* Get the first token from line */
984
985         if ( !com ) continue;             /* Ignore blank lines */
986         if ( com[0] == '#' ) continue;    /* Ignore comments */
987
988/* Read configuration parameters
989   *****************************/
990         if ( k_its( "iNumRegions" ) ) 
991         {
992            *piNumRegs = k_int();
993/* Allocate and init the Alarm buffer
994   **********************************/
995            InBufl = sizeof( ALARMSTRUCT ) * *piNumRegs; 
996            ptAS = (ALARMSTRUCT *) calloc( *piNumRegs, sizeof( ALARMSTRUCT ) );
997            if ( ptAS == NULL )
998            {
999               logit( "et", "pick_wcatwc: Cannot allocate alarm buffer\n");
1000               nfiles = k_close();
1001               return( 0 );
1002            }
1003
1004            for ( i=0; i<*piNumRegs; i++ )
1005            {
1006               for ( ii=0; ii<MAX_ALARM_STN; ii++ )
1007                  strcpy( ptAS[i].szStnAlarm[ii], "\0" );
1008               ptAS[i].dLastTime = 0.;
1009               ptAS[i].iNumPicksCnt = 0;
1010               for ( j=0; j<5; j++ )
1011               {
1012                  k_rd();
1013                  com = k_str();                    /* Read more */
1014                  if ( k_its( "Region" ) )
1015                  { if ( str = k_str() ) strcpy( ptAS[i].szRegionName, str ); }
1016                  else if ( k_its( "AlarmThresh" ) )
1017                  {  ptAS[i].iAlarmThresh = k_int(); }
1018                  else if ( k_its( "PStrength" ) )
1019                  {  ptAS[i].dThresh = k_val(); }
1020                  else if ( k_its( "MaxTime" ) )
1021                  {  ptAS[i].dMaxTime = k_val(); }
1022                  else if ( k_its( "NumStnInReg" ) )
1023                  {
1024                     ptAS[i].iNumStnInReg = k_int();
1025                     if ( ptAS[i].iNumStnInReg >= MAX_ALARM_STN )
1026                     {
1027                        logit( "et", "Too many alarm stns - %ld\n",
1028                                ptAS[i].iNumStnInReg );
1029                        nfiles = k_close();
1030                        free( ptAS );
1031                        return( 0 );
1032                     }
1033                     for ( k=0; k<ptAS[i].iNumStnInReg; k++ )
1034                     {
1035                        k_rd();
1036                        com = k_str();                    /* Read more */
1037                        if ( k_its( "Station" ) )
1038                        {  if ( str = k_str() ) strcpy( ptAS[i].szStation[k], str ); }
1039                     }
1040                  }
1041               }                                 
1042            }
1043         }
1044      }             
1045      nfiles = k_close();
1046   }
1047   
1048/* Log input data
1049   **************/   
1050   logit( "", "NumAlarmRegions = %ld\n", *piNumRegs );
1051   for ( i=0; i<*piNumRegs; i++ )
1052   {
1053      logit( "", "%s\n", ptAS[i].szRegionName );
1054      logit( "", "Alarm Threshold %ld\n", ptAS[i].iAlarmThresh );
1055      logit( "", "PStrength Threshold %lf\n", ptAS[i].dThresh );
1056      logit( "", "MaxTime %lf\n", ptAS[i].dMaxTime );
1057      logit( "", "NumStnInReg %ld\n", ptAS[i].iNumStnInReg );
1058      for ( k=0; k<ptAS[i].iNumStnInReg; k++ )
1059         logit( "", "%s\n", ptAS[i].szStation[k] );
1060   }
1061   *pAS = ptAS;
1062     
1063   return( 1 );
1064}
Note: See TracBrowser for help on using the repository browser.