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

Revision 4964, 42.0 KB checked in by paulf, 8 years ago (diff)

cleaned up windows NT compile

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