source: SH_SHM/trunk/util/check_evt.c @ 87

Revision 87, 18.2 KB checked in by marcus, 14 years ago (diff)

r69 | svn | 2008-12-07 15:30:25 +0100 (So, 07 Dez 2008) | 1 line

more tests in check_evt
Line 
1
2/* file check_evt.c
3 *      ===========
4 *
5 * version 10, 16-Jan-2007
6 *
7 * Checks consistency of created output file of SHM
8 * K. Stammler, 21-Jul-94
9 */
10
11
12/*
13 *
14 *  SeismicHandler, seismic analysis software
15 *  Copyright (C) 1996,  Klaus Stammler, Federal Institute for Geosciences
16 *                                       and Natural Resources (BGR), Germany
17 *
18 *  This program is free software; you can redistribute it and/or modify
19 *  it under the terms of the GNU General Public License as published by
20 *  the Free Software Foundation; either version 2 of the License, or
21 *  (at your option) any later version.
22 *
23 *  This program is distributed in the hope that it will be useful,
24 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
25 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
26 *  GNU General Public License for more details.
27 *
28 *  You should have received a copy of the GNU General Public License
29 *  along with this program; if not, write to the Free Software
30 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
31 *
32 */
33
34
35
36#include <stdio.h>
37#include <string.h>
38#include BASECNST
39#include BC_SYSBASE
40#include "eventdsc.h"
41#include "tcusrdef.h"
42#include "earthloc.h"
43
44
45#define MAX_PKP_SLOWNESS 4.5
46#define MIN_P_SLOWNESS 4.3
47
48#define EVTYPE_UNDEF 0
49#define EVTYPE_LOCAL 1
50#define EVTYPE_P     2
51#define EVTYPE_PKP   3
52
53/* next error number currently 031 */
54
55
56
57/* global variables */
58static float   efv_reflat=49.0;           /* reference latitude */
59static float   efv_reflon=11.0;           /* reference longitude */
60static char    efv_infile[BC_FILELTH+1];  /* name of input file */
61static char    efv_uncorrphase[EvPHASELTH+1]; /* name of phase on uncorr beam */
62static BOOLEAN efv_first_arrival=FALSE;   /* first arrival phase selected */
63static BOOLEAN efv_only_grf=TRUE;         /* only GRF stations used */
64static BOOLEAN efv_gra1_picked=FALSE;     /* pick on GRA1 found */
65static BOOLEAN efv_uncorr_loc=FALSE;      /* uncorrected beam for location */
66static BOOLEAN efv_telex_phase=FALSE;     /* telex phase found */
67static BOOLEAN efv_mb_found=FALSE;        /* mb determination done */
68static BOOLEAN efv_ms_found=FALSE;        /* ms determination done */
69static BOOLEAN efv_mbb_found=FALSE;       /* mbb determination done */
70static int     efv_subtype=EVTYPE_UNDEF;  /* event type */
71static EvEventTypeT efv_evtype=EvcEventTypeUndefined; /* event type */
72static float   efv_depth=(-100.0);        /* depth of event */
73static float   efv_magn=0.0;              /* max. magnitude found */
74static int     efv_pkp_count=0;           /* number of PKP phases */
75static BOOLEAN efv_illegal_origin;        /* illegal origin time */
76static BOOLEAN efv_foundPg=FALSE;         /* Pg phase found */
77static BOOLEAN efv_foundPnSn=FALSE;       /* Pn/Sn found */
78static BOOLEAN efv_foundTele=FALSE;       /* teleseismic phase found */
79static float   efv_latitude=0.0;          /* epicentre latitude */
80static float   efv_longitude=0.0;         /* epicentre longitude */
81static char    efv_analyst[EvANALYSTLTH+1]=""; /* analysts initials */
82static BOOLEAN efv_showinitials=FALSE;    /* show initials */
83
84
85/* prototypes of local routines */
86void EfPutCheck( FILE *out, EvEventT *event, EvStatusT *status );
87static BOOLEAN EfPhaseIsFirstArrival( char name[] );
88static void EfFinalCheck( FILE *out );
89static BOOLEAN EfIsLocalPhase( char phase[] );
90
91
92
93int main( int argc, char *argv[] )
94{
95        /* local variables */
96        char      infile[BC_FILELTH+1];    /* input file */
97        char      outfile[BC_FILELTH+1];   /* output file */
98        EvStatusT status;                  /* return status */
99        int       i;                       /* counter */
100
101        /* executable code */
102
103        if  (argc < 2)  {
104                fprintf( stderr, "Usage: %s <input> [<output>] [<showinitials>]\n",
105                        argv[0] );
106                return 1;
107        } /*endif*/
108
109        strcpy( infile, argv[1] );
110        if  (argc > 2)
111                strcpy( outfile, argv[2] );
112                if  (*outfile == '\0')
113                        strcpy( outfile, "TT" );
114        else
115                strcpy( outfile, "TT" );
116        efv_showinitials = FALSE;
117        if  (argc > 3)
118                efv_showinitials = (strcmp(argv[3],"showinitials") == 0);
119
120        /* efv_infile is used for messages, directory will be removed */
121        i = strlen( infile ) - 1;
122        while  (infile[i] != '/' && i > 0)
123                i--;
124        strcpy( efv_infile, infile+i );
125
126        status = EveNOERROR;
127        EvReformatEventfile( infile, outfile, EvGetEvent,
128                EfPutCheck, &status );
129
130        return 0;
131
132} /* end of main */
133
134
135
136/*---------------------------------------------------------------------------*/
137
138
139
140void EfPutCheck( FILE *out, EvEventT *event, EvStatusT *status )
141
142/* Prints error messages of event to file 'out'.
143 *
144 * parameters of routine
145 * FILE       *out;          input; output file
146 * EvEventT   *event;        input; event information
147 * EvStatusT  *status;       output; return status
148 */
149{
150        /* local variables */
151        static BOOLEAN locqual_checked=FALSE; /* location quality checked ? */
152        BOOLEAN  amplitude_measured;      /* amplitude measured */
153        BOOLEAN  is_l_phase;              /* phase is L (or (L)) */
154        BOOLEAN  is_pkp_phase;            /* phase is PKP */
155        BOOLEAN  is_p_phase;              /* phase is P or (P) */
156        BOOLEAN  loc_given;               /* location given in this record */
157        NTIME    ntime;                   /* numeric origin time */
158
159        /* executable code */
160
161        /* final check if event is NULL */
162        if  (event == NULL)  {
163                EfFinalCheck( out );
164                locqual_checked = FALSE;
165                return;
166        } /*endif*/
167
168        /* fprintf( out, "--> check phase %s\n", event->phase ); */
169
170        /* is a location given here ? */
171        loc_given = (event->latitude > 0.0 && event->longitude > 0.0);
172        if  (loc_given)  {
173                efv_latitude = event->latitude;
174                efv_longitude = event->longitude;
175        } /*endif*/
176
177        /* is this an L phase ? */
178        is_l_phase = (strcmp(event->phase,"L") == 0
179                || strcmp(event->phase,"(L)") == 0);
180
181        /* is this a P phase ? */
182        is_p_phase = (strcmp(event->phase,"P") == 0
183                || strcmp(event->phase,"(P)") == 0);
184
185        /* is this a PKP phase ? */
186        is_pkp_phase = (strstr(event->phase,"PKP") != NULL);
187        if  (is_pkp_phase)
188                efv_foundTele = TRUE;
189
190        /* a depth phase can found only in teleseismic events */
191        if  (event->phase[0] == 'p' || event->phase[0] == 's')
192                efv_foundTele = TRUE;
193
194        /* check for local and regional phases */
195        if  (strcmp(event->phase,"Pg") == 0)
196                efv_foundPg = TRUE;
197        if  (strcmp(event->phase,"Pn") == 0)
198                efv_foundPnSn = TRUE;
199        if  (strcmp(event->phase,"Sn") == 0)
200                efv_foundPnSn = TRUE;
201
202        if  (event->event_type != EvcEventTypeUndefined)
203                efv_evtype = event->event_type;
204
205        /* setup global variables used in final check */
206        /* ------------------------------------------ */
207
208        /* is it a GRF station ? */
209        if  (strncmp(event->station,"GR",2) != 0
210                || event->station[2] < 'A' || event->station[2] > 'C')
211                efv_only_grf = FALSE;
212        if  (strcmp(event->station,"GRA1") == 0)
213                efv_gra1_picked = TRUE;
214
215        /* is a first arrival phase there ? */
216        if  (EfPhaseIsFirstArrival(event->phase))
217                efv_first_arrival = TRUE;
218
219        /* is there a location done with uncorrected beam ? */
220        if  ((loc_given && event->b_slowness != EvEMPTY_SLOWNESS &&
221                event->l_slowness == EvEMPTY_SLOWNESS
222                && event->loc_quality == EvcLocQualReliable
223                && strcmp(event->source,"SZGRF") == 0)
224                || (event->loc_method == EvcLocMethUncorrBeam))  {
225                efv_uncorr_loc = TRUE;
226                strcpy( efv_uncorrphase, event->phase );
227        } /*endif*/
228
229        /* is there an origin time from 1970? */
230        if  (loc_given)  {
231                tc_t2n( event->origin_time, &ntime, status );
232                if  (SySevere(status))  {
233                        efv_illegal_origin = TRUE;
234                } else if  (ntime.year == 1970)  {
235                        efv_illegal_origin = TRUE;
236                } /*endif*/
237        } /*endif*/
238
239        /* look for telex phase */
240        if  (strchr(event->phase_flags,EvFLAG_TELEX) != NULL)
241                efv_telex_phase = TRUE;
242
243        /* look for depth */
244        if  (efv_depth < -99.0 && event->depth_type != EvcDepthUndefined)
245                efv_depth = event->depth;
246
247        /* find plain mistakes */
248        /* ------------------- */
249
250        /* special for L phases */
251        if  (is_l_phase)  {
252                /* is there an impulsive L phase ? */
253                if  (event->onset_type == EvcOnsetImpulsive)
254                        fprintf( out, "F001 %s: phase L with impulsive onset\n", efv_infile );
255                /* is a legal filter used ? */
256                if  (strcmp(event->filter,"G_WWSSN_SP") == 0
257                        || strcmp(event->filter,"STANDARD_BP") == 0
258                        || event->filter[0] == '\0')
259                        fprintf( out, "F002 %s: L phase reading with illegal filter %s\n",
260                                efv_infile, event->filter );
261                /* is there slowness and azimuth ? if yes that's real nonsense */
262                if  (event->b_slowness != EvEMPTY_SLOWNESS
263                        || event->b_azimuth != EvEMPTY_AZIMUTH
264                        || event->l_slowness != EvEMPTY_SLOWNESS
265                        || event->l_azimuth != EvEMPTY_AZIMUTH)
266                        fprintf( out, "F003 %s: slowness and/or azimuth for phase %s\n",
267                                efv_infile, event->phase );
268                /* is amplitude & period determined ? */
269                if  (event->amplitude == EvEMPTY_AMPLITUDE
270                        || event->period == EvEMPTY_PERIOD)
271                        fprintf( out, "F004 %s: no amplitude / period for phase %s\n",
272                                efv_infile, event->phase );
273        } /*endif*/
274
275        if  (is_pkp_phase)  {
276                /* check slowness range */
277                if  (event->b_slowness != EvEMPTY_SLOWNESS
278                        && event->b_slowness > MAX_PKP_SLOWNESS)
279                        fprintf( out, "F005 %s: illegal beam slowness %3.1f for phase %s\n",
280                                efv_infile, event->b_slowness, event->phase );
281                if  (event->l_slowness != EvEMPTY_SLOWNESS
282                        && event->l_slowness > MAX_PKP_SLOWNESS)
283                        fprintf( out, "F006 %s: illegal corr. slowness %3.1f for phase %s\n",
284                                efv_infile, event->l_slowness, event->phase );
285                efv_subtype = EVTYPE_PKP;
286                efv_pkp_count++;
287        } /*endif*/
288
289        if  (is_p_phase)  {
290                /* check slowness range */
291                if  (event->b_slowness != EvEMPTY_SLOWNESS
292                        && event->b_slowness < MIN_P_SLOWNESS)
293                        if  (event->l_slowness == EvEMPTY_SLOWNESS)
294                                fprintf( out, "F007 %s: illegal beam slowness %3.1f for phase %s\n",
295                                        efv_infile, event->b_slowness, event->phase );
296                if  (event->l_slowness != EvEMPTY_SLOWNESS
297                        && event->l_slowness < MIN_P_SLOWNESS)
298                        fprintf( out, "F008 %s: illegal corr. slowness %3.1f for phase %s\n",
299                                efv_infile, event->l_slowness, event->phase );
300                efv_subtype = EVTYPE_P;
301        } /*endif*/
302
303        if  (EfIsLocalPhase(event->phase))  efv_subtype = EVTYPE_LOCAL;
304
305        /* amplitude measurement must be complete */
306        amplitude_measured = (event->period != EvEMPTY_PERIOD
307                || event->amplitude != EvEMPTY_AMPLITUDE
308                || event->amplitude_time != EvEMPTY_AMPLITUDE_TIME
309                || event->amplitude_vel != EvEMPTY_AMPLITUDE);
310        if  (amplitude_measured)
311                if  (event->period == EvEMPTY_PERIOD
312                        || event->amplitude == EvEMPTY_AMPLITUDE
313                        || (event->amplitude_time == EvEMPTY_AMPLITUDE_TIME && !is_l_phase)
314                        || event->amplitude_vel == EvEMPTY_AMPLITUDE
315                        || event->ap_source == EvcApSourceUndefined)
316                        fprintf( out,
317                                "F009 %s: incomplete amplitude determination at phase %s, station %s\n",
318                                efv_infile, event->phase, event->station );
319
320        /* check for magnitudes determined */
321        if  (event->mag[EvcMagMb] != EvEMPTY_MAGNITUDE)  {
322                efv_mb_found = TRUE;
323                if  (event->mag[EvcMagMb] > efv_magn)  efv_magn = event->mag[EvcMagMb];
324        } /*endif*/
325        if  (event->mag[EvcMagMs] != EvEMPTY_MAGNITUDE)  {
326                efv_ms_found = TRUE;
327                if  (event->mag[EvcMagMs] > efv_magn)  efv_magn = event->mag[EvcMagMs];
328        } /*endif*/
329        if  (event->mag[EvcMagMbb] != EvEMPTY_MAGNITUDE) efv_mbb_found = TRUE;
330
331        /* magnitudes mb,MS are valid only with amplitude */
332        if  (event->mag[EvcMagMs] != EvEMPTY_MAGNITUDE
333                || event->mag[EvcMagMb] != EvEMPTY_MAGNITUDE)
334                if  (!amplitude_measured)
335                        fprintf( out, "F010 %s: magnitude without amplitude on phase %s\n",
336                                efv_infile, event->phase );
337
338        /* a reliable location with PKP is not possible */
339        if  (event->loc_quality == EvcLocQualReliable && is_pkp_phase
340                && strcmp(event->source,"SZGRF") == 0 && loc_given)
341                fprintf( out, "F011 %s: reliable location with phase %s not possible\n",
342                        efv_infile, event->phase );
343
344        /* event type tele and local phases Pn, Pg, Sn, Sg */
345        if  (event->event_type == EvcEventTypeTeleQuake
346                && EfIsLocalPhase(event->phase))
347                fprintf( out, "F012 %s: incompatible event type TELE with phase %s\n",
348                        efv_infile, event->phase );
349
350        /* is there a location quality specified in the first phase ? */
351        if  (loc_given && !locqual_checked && event->phase[0] != '\0')  {
352                if  (event->loc_quality == EvcLocQualUndefined)
353                        fprintf( out, "F013 %s: no location quality specified\n", efv_infile );
354                locqual_checked = TRUE;
355        } /*endif*/
356
357        if  (event->analyst[0] != '\0')
358                strcpy( efv_analyst, event->analyst );
359
360} /* end of EfPutCheck */
361
362
363
364/*---------------------------------------------------------------------------*/
365
366
367
368static BOOLEAN EfPhaseIsFirstArrival( char name[] )
369
370/* Checks whether phase is first arrival
371 *
372 * parameters of routine
373 * char       name[];        input; name of phase
374 *                           returns TRUE if phase is possible first arrival
375 */
376{
377        /* executable code */
378
379        if  (strcmp(name,"P") == 0)  return TRUE;
380        if  (strcmp(name,"(P)") == 0)  return TRUE;
381        if  (strncmp(name,"PKP",3) == 0)  return TRUE;
382        if  (strncmp(name,"(PKP",4) == 0)  return TRUE;
383        if  (strncmp(name,"Pn",2) == 0)  return TRUE;
384        if  (strncmp(name,"Pg",2) == 0)  return TRUE;
385        if  (strcmp(name,"(Pn)") == 0)  return TRUE;
386        if  (strcmp(name,"(Pg)") == 0)  return TRUE;
387        if  (strcmp(name,"Pdiff") == 0)  return TRUE;
388        if  (strcmp(name,"(Pdiff)") == 0)  return TRUE;
389        if  (strcmp(name,"Pdif") == 0)  return TRUE;
390        if  (strcmp(name,"(Pdif)") == 0)  return TRUE;
391
392        return FALSE;
393
394} /* end of EfPhaseIsFirstArrival */
395
396
397
398/*---------------------------------------------------------------------------*/
399
400
401
402static void EfFinalCheck( FILE *out )
403
404/* Prints result of final check to output file 'out'
405 *
406 * parameters of routine
407 * FILE       *out;         input; pointer to output file
408 */
409{
410        /* local variables */
411        char     expl[cBcLongStrLth+1];   /* explanation string */
412
413        /* executable code */
414
415        /* fprintf( out, "--> final check\n" ); */
416        if  (efv_showinitials)
417                sprintf( expl, "%s (%s)", efv_infile, efv_analyst );
418        else
419                strcpy( expl, efv_infile );
420
421        /* a first arrival must be specified */
422        if  (!efv_first_arrival)
423                fprintf( out, "F014 %s: no first arrival phase selected\n", expl );
424
425        /* beam with uncorrected slowness at GRF is at least suspicious */
426        if  (efv_only_grf && efv_uncorr_loc)
427                fprintf( out, "F015 %s: GRF location with uncorrected beam at phase %s\n",
428                        expl, efv_uncorrphase );
429
430        /* only GRF picked but no GRA1? */
431        if  (efv_only_grf && !efv_gra1_picked)
432                fprintf( out, "F035 %s: only GRF picks but none at GRA1\n", expl );
433
434        /* no telex phase */
435        /* disabled 12-Nov-2008 K.S.
436        if  (!efv_telex_phase)
437                fprintf( out, "F016 %s: no telex phase found\n", expl );
438        */
439
440        /* mb determined? */
441        if  (efv_subtype == EVTYPE_P && !efv_mb_found)
442                fprintf( out, "F017 %s: no Magnitude mb determined for P phase\n", expl );
443
444        /* Ms determined? */
445        if  (efv_subtype >= EVTYPE_P && efv_depth <= 50.0 && efv_magn > 5.6 &&
446                !efv_ms_found)
447                fprintf( out, "F018 %s: no Magnitude Ms determined for magn %3.1f event\n",
448                        expl, efv_magn );
449        if  (efv_subtype == EVTYPE_PKP && efv_depth <= 50.0 && efv_pkp_count > 3 &&
450                !efv_ms_found)
451                fprintf( out, "F019 %s: no Magnitude Ms determined for PKP event\n",
452                        expl );
453        if  (efv_ms_found && efv_depth > 50.0)
454                fprintf( out, "F020 %s: Magnitude Ms determined for depth %5.1f event\n",
455                        expl, efv_depth );
456
457        /* Mbb determined? */
458        if  (efv_mb_found && efv_magn > 6.2 && !efv_mbb_found)
459                fprintf( out, "F021 %s: no Broadband Magnitude determined for magn %3.1f event\n",
460                        expl, efv_magn );
461
462        /* depth ok? */
463        if  (efv_evtype == EvcEventTypeMining && Abs(efv_depth-1.0) > 0.1)
464                fprintf( out, "F022 %s: depth %5.1f given for mining event\n",
465                        expl, efv_depth );
466        if  (efv_evtype == EvcEventTypeBlast && Abs(efv_depth) > 0.1)
467                fprintf( out, "F023 %s: depth %5.1f given for blast\n",
468                        expl, efv_depth );
469        if  (efv_evtype == EvcEventTypeTeleQuake && (efv_depth < 10.0))
470                fprintf( out, "F026 %s: depth %5.1f for teleseismic event\n",
471                        expl, efv_depth );
472        if  ((efv_evtype == EvcEventTypeLocalQuake
473                || efv_evtype == EvcEventTypeRegioQuake) && (efv_depth < 2.0))
474                fprintf( out, "F027 %s: depth %5.1f for local/regional event\n",
475                        expl, efv_depth );
476        if  (efv_evtype == EvcEventTypeLocalQuake && efv_depth > 50.0)
477                fprintf( out, "F036 %s: depth %5.1f given for local event\n",
478                        expl, efv_depth );
479        if  (efv_evtype == EvcEventTypeUndefined)
480                fprintf( out, "F024 %s: no event type specified\n", expl );
481        if  (efv_depth < -99.0)
482                fprintf( out, "F025 %s: no depth specified\n", expl );
483
484        /* proper event type? */
485        if  (efv_foundPg && efv_evtype == EvcEventTypeTeleQuake)
486                fprintf( out, "F028 %s: improper teleseismic event type: found Pg\n",
487                        expl );
488        if  (efv_foundPnSn && efv_evtype == EvcEventTypeTeleQuake)
489                fprintf( out, "F029 %s: improper teleseismic event type: found Pn/Sn\n",
490                        expl );
491        if  (efv_foundTele && efv_evtype == EvcEventTypeLocalQuake)
492                fprintf( out, "F030 %s: improper local event type: found tele phase\n",
493                        expl );
494        if  (efv_foundTele && efv_evtype == EvcEventTypeRegioQuake)
495                fprintf( out, "F031 %s: improper regional event type: found tele phase\n",
496                        expl );
497
498        /* check distance and event type */
499        if  (efv_latitude != 0.0 || efv_longitude != 0.0)  {
500                double dist, azim, bazim;   /* result of distance computation */
501                float fdist;    /* single precision distance */
502                mb_locdiff( (double)efv_reflat, (double)efv_reflon, (double)efv_latitude,
503                        (double)efv_longitude, &dist, &azim, &bazim );
504                fdist = dist;
505                if  (dist > 15.0 && efv_evtype == EvcEventTypeLocalQuake)
506                        fprintf( out, "F032 %s: improper local event type (distance %f deg)\n",
507                                expl, dist );
508                if  (dist > 30.0 && efv_evtype == EvcEventTypeRegioQuake)
509                        fprintf( out, "F033 %s: improper regional event type (distance %f deg)\n",
510                                expl, dist );
511                if  (dist < 25.0 && efv_evtype == EvcEventTypeTeleQuake)
512                        fprintf( out, "F034 %s: improper teleseismic event type (distance %f deg)\n",
513                                expl, dist );
514        } /*endif*/
515
516} /* end of EfFInalCheck */
517
518
519
520/*---------------------------------------------------------------------------*/
521
522
523
524static BOOLEAN EfIsLocalPhase( char phase[] )
525
526/* checks whether 'phase' is a local phase
527 *
528 * parameters of routine
529 * char       phase[];         input; phase name
530 */
531{
532        /* local variables */
533
534        /* executable code */
535
536        if  (strcmp(phase,"Pg") == 0)  return TRUE;
537        if  (strcmp(phase,"Pn") == 0)  return TRUE;
538        if  (strcmp(phase,"Sg") == 0)  return TRUE;
539        if  (strcmp(phase,"Sn") == 0)  return TRUE;
540        if  (strcmp(phase,"(Pg)") == 0)  return TRUE;
541        if  (strcmp(phase,"(Pn)") == 0)  return TRUE;
542        if  (strcmp(phase,"(Sg)") == 0)  return TRUE;
543        if  (strcmp(phase,"(Sn)") == 0)  return TRUE;
544
545        return FALSE;
546
547} /* end of EfIsLocalPhase */
548
549
550
551/*---------------------------------------------------------------------------*/
552
Note: See TracBrowser for help on using the repository browser.