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

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

r71 | svn | 2008-12-11 09:20:53 +0100 (Do, 11 Dez 2008) | 1 line

fixed bug in check_evt (missing braces at if)
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        } /*endif*/
117        efv_showinitials = FALSE;
118        if  (argc > 3)
119                efv_showinitials = (strcmp(argv[3],"showinitials") == 0);
120
121        /* efv_infile is used for messages, directory will be removed */
122        i = strlen( infile ) - 1;
123        while  (infile[i] != '/' && i > 0)
124                i--;
125        strcpy( efv_infile, infile+i );
126
127        status = EveNOERROR;
128        EvReformatEventfile( infile, outfile, EvGetEvent,
129                EfPutCheck, &status );
130
131        return 0;
132
133} /* end of main */
134
135
136
137/*---------------------------------------------------------------------------*/
138
139
140
141void EfPutCheck( FILE *out, EvEventT *event, EvStatusT *status )
142
143/* Prints error messages of event to file 'out'.
144 *
145 * parameters of routine
146 * FILE       *out;          input; output file
147 * EvEventT   *event;        input; event information
148 * EvStatusT  *status;       output; return status
149 */
150{
151        /* local variables */
152        static BOOLEAN locqual_checked=FALSE; /* location quality checked ? */
153        BOOLEAN  amplitude_measured;      /* amplitude measured */
154        BOOLEAN  is_l_phase;              /* phase is L (or (L)) */
155        BOOLEAN  is_pkp_phase;            /* phase is PKP */
156        BOOLEAN  is_p_phase;              /* phase is P or (P) */
157        BOOLEAN  loc_given;               /* location given in this record */
158        NTIME    ntime;                   /* numeric origin time */
159
160        /* executable code */
161
162        /* final check if event is NULL */
163        if  (event == NULL)  {
164                EfFinalCheck( out );
165                locqual_checked = FALSE;
166                return;
167        } /*endif*/
168
169        /* fprintf( out, "--> check phase %s\n", event->phase ); */
170
171        /* is a location given here ? */
172        loc_given = (event->latitude > 0.0 && event->longitude > 0.0);
173        if  (loc_given)  {
174                efv_latitude = event->latitude;
175                efv_longitude = event->longitude;
176        } /*endif*/
177
178        /* is this an L phase ? */
179        is_l_phase = (strcmp(event->phase,"L") == 0
180                || strcmp(event->phase,"(L)") == 0);
181
182        /* is this a P phase ? */
183        is_p_phase = (strcmp(event->phase,"P") == 0
184                || strcmp(event->phase,"(P)") == 0);
185
186        /* is this a PKP phase ? */
187        is_pkp_phase = (strstr(event->phase,"PKP") != NULL);
188        if  (is_pkp_phase)
189                efv_foundTele = TRUE;
190
191        /* a depth phase can found only in teleseismic events */
192        if  (event->phase[0] == 'p' || event->phase[0] == 's')
193                efv_foundTele = TRUE;
194
195        /* check for local and regional phases */
196        if  (strcmp(event->phase,"Pg") == 0)
197                efv_foundPg = TRUE;
198        if  (strcmp(event->phase,"Pn") == 0)
199                efv_foundPnSn = TRUE;
200        if  (strcmp(event->phase,"Sn") == 0)
201                efv_foundPnSn = TRUE;
202
203        if  (event->event_type != EvcEventTypeUndefined)
204                efv_evtype = event->event_type;
205
206        /* setup global variables used in final check */
207        /* ------------------------------------------ */
208
209        /* is it a GRF station ? */
210        if  (strncmp(event->station,"GR",2) != 0
211                || event->station[2] < 'A' || event->station[2] > 'C')
212                efv_only_grf = FALSE;
213        if  (strcmp(event->station,"GRA1") == 0)
214                efv_gra1_picked = TRUE;
215
216        /* is a first arrival phase there ? */
217        if  (EfPhaseIsFirstArrival(event->phase))
218                efv_first_arrival = TRUE;
219
220        /* is there a location done with uncorrected beam ? */
221        if  ((loc_given && event->b_slowness != EvEMPTY_SLOWNESS &&
222                event->l_slowness == EvEMPTY_SLOWNESS
223                && event->loc_quality == EvcLocQualReliable
224                && strcmp(event->source,"SZGRF") == 0)
225                || (event->loc_method == EvcLocMethUncorrBeam))  {
226                efv_uncorr_loc = TRUE;
227                strcpy( efv_uncorrphase, event->phase );
228        } /*endif*/
229
230        /* is there an origin time from 1970? */
231        if  (loc_given)  {
232                tc_t2n( event->origin_time, &ntime, status );
233                if  (SySevere(status))  {
234                        efv_illegal_origin = TRUE;
235                } else if  (ntime.year == 1970)  {
236                        efv_illegal_origin = TRUE;
237                } /*endif*/
238        } /*endif*/
239
240        /* look for telex phase */
241        if  (strchr(event->phase_flags,EvFLAG_TELEX) != NULL)
242                efv_telex_phase = TRUE;
243
244        /* look for depth */
245        if  (efv_depth < -99.0 && event->depth_type != EvcDepthUndefined)
246                efv_depth = event->depth;
247
248        /* find plain mistakes */
249        /* ------------------- */
250
251        /* special for L phases */
252        if  (is_l_phase)  {
253                /* is there an impulsive L phase ? */
254                if  (event->onset_type == EvcOnsetImpulsive)
255                        fprintf( out, "F001 %s: phase L with impulsive onset\n", efv_infile );
256                /* is a legal filter used ? */
257                if  (strcmp(event->filter,"G_WWSSN_SP") == 0
258                        || strcmp(event->filter,"STANDARD_BP") == 0
259                        || event->filter[0] == '\0')
260                        fprintf( out, "F002 %s: L phase reading with illegal filter %s\n",
261                                efv_infile, event->filter );
262                /* is there slowness and azimuth ? if yes that's real nonsense */
263                if  (event->b_slowness != EvEMPTY_SLOWNESS
264                        || event->b_azimuth != EvEMPTY_AZIMUTH
265                        || event->l_slowness != EvEMPTY_SLOWNESS
266                        || event->l_azimuth != EvEMPTY_AZIMUTH)
267                        fprintf( out, "F003 %s: slowness and/or azimuth for phase %s\n",
268                                efv_infile, event->phase );
269                /* is amplitude & period determined ? */
270                if  (event->amplitude == EvEMPTY_AMPLITUDE
271                        || event->period == EvEMPTY_PERIOD)
272                        fprintf( out, "F004 %s: no amplitude / period for phase %s\n",
273                                efv_infile, event->phase );
274        } /*endif*/
275
276        if  (is_pkp_phase)  {
277                /* check slowness range */
278                if  (event->b_slowness != EvEMPTY_SLOWNESS
279                        && event->b_slowness > MAX_PKP_SLOWNESS)
280                        fprintf( out, "F005 %s: illegal beam slowness %3.1f for phase %s\n",
281                                efv_infile, event->b_slowness, event->phase );
282                if  (event->l_slowness != EvEMPTY_SLOWNESS
283                        && event->l_slowness > MAX_PKP_SLOWNESS)
284                        fprintf( out, "F006 %s: illegal corr. slowness %3.1f for phase %s\n",
285                                efv_infile, event->l_slowness, event->phase );
286                efv_subtype = EVTYPE_PKP;
287                efv_pkp_count++;
288        } /*endif*/
289
290        if  (is_p_phase)  {
291                /* check slowness range */
292                if  (event->b_slowness != EvEMPTY_SLOWNESS
293                        && event->b_slowness < MIN_P_SLOWNESS)
294                        if  (event->l_slowness == EvEMPTY_SLOWNESS)
295                                fprintf( out, "F007 %s: illegal beam slowness %3.1f for phase %s\n",
296                                        efv_infile, event->b_slowness, event->phase );
297                if  (event->l_slowness != EvEMPTY_SLOWNESS
298                        && event->l_slowness < MIN_P_SLOWNESS)
299                        fprintf( out, "F008 %s: illegal corr. slowness %3.1f for phase %s\n",
300                                efv_infile, event->l_slowness, event->phase );
301                efv_subtype = EVTYPE_P;
302        } /*endif*/
303
304        if  (EfIsLocalPhase(event->phase))  efv_subtype = EVTYPE_LOCAL;
305
306        /* amplitude measurement must be complete */
307        amplitude_measured = (event->period != EvEMPTY_PERIOD
308                || event->amplitude != EvEMPTY_AMPLITUDE
309                || event->amplitude_time != EvEMPTY_AMPLITUDE_TIME
310                || event->amplitude_vel != EvEMPTY_AMPLITUDE);
311        if  (amplitude_measured)
312                if  (event->period == EvEMPTY_PERIOD
313                        || event->amplitude == EvEMPTY_AMPLITUDE
314                        || (event->amplitude_time == EvEMPTY_AMPLITUDE_TIME && !is_l_phase)
315                        || event->amplitude_vel == EvEMPTY_AMPLITUDE
316                        || event->ap_source == EvcApSourceUndefined)
317                        fprintf( out,
318                                "F009 %s: incomplete amplitude determination at phase %s, station %s\n",
319                                efv_infile, event->phase, event->station );
320
321        /* check for magnitudes determined */
322        if  (event->mag[EvcMagMb] != EvEMPTY_MAGNITUDE)  {
323                efv_mb_found = TRUE;
324                if  (event->mag[EvcMagMb] > efv_magn)  efv_magn = event->mag[EvcMagMb];
325        } /*endif*/
326        if  (event->mag[EvcMagMs] != EvEMPTY_MAGNITUDE)  {
327                efv_ms_found = TRUE;
328                if  (event->mag[EvcMagMs] > efv_magn)  efv_magn = event->mag[EvcMagMs];
329        } /*endif*/
330        if  (event->mag[EvcMagMbb] != EvEMPTY_MAGNITUDE) efv_mbb_found = TRUE;
331
332        /* magnitudes mb,MS are valid only with amplitude */
333        if  (event->mag[EvcMagMs] != EvEMPTY_MAGNITUDE
334                || event->mag[EvcMagMb] != EvEMPTY_MAGNITUDE)
335                if  (!amplitude_measured)
336                        fprintf( out, "F010 %s: magnitude without amplitude on phase %s\n",
337                                efv_infile, event->phase );
338
339        /* a reliable location with PKP is not possible */
340        if  (event->loc_quality == EvcLocQualReliable && is_pkp_phase
341                && strcmp(event->source,"SZGRF") == 0 && loc_given)
342                fprintf( out, "F011 %s: reliable location with phase %s not possible\n",
343                        efv_infile, event->phase );
344
345        /* event type tele and local phases Pn, Pg, Sn, Sg */
346        if  (event->event_type == EvcEventTypeTeleQuake
347                && EfIsLocalPhase(event->phase))
348                fprintf( out, "F012 %s: incompatible event type TELE with phase %s\n",
349                        efv_infile, event->phase );
350
351        /* is there a location quality specified in the first phase ? */
352        if  (loc_given && !locqual_checked && event->phase[0] != '\0')  {
353                if  (event->loc_quality == EvcLocQualUndefined)
354                        fprintf( out, "F013 %s: no location quality specified\n", efv_infile );
355                locqual_checked = TRUE;
356        } /*endif*/
357
358        if  (event->analyst[0] != '\0')
359                strcpy( efv_analyst, event->analyst );
360
361} /* end of EfPutCheck */
362
363
364
365/*---------------------------------------------------------------------------*/
366
367
368
369static BOOLEAN EfPhaseIsFirstArrival( char name[] )
370
371/* Checks whether phase is first arrival
372 *
373 * parameters of routine
374 * char       name[];        input; name of phase
375 *                           returns TRUE if phase is possible first arrival
376 */
377{
378        /* executable code */
379
380        if  (strcmp(name,"P") == 0)  return TRUE;
381        if  (strcmp(name,"(P)") == 0)  return TRUE;
382        if  (strncmp(name,"PKP",3) == 0)  return TRUE;
383        if  (strncmp(name,"(PKP",4) == 0)  return TRUE;
384        if  (strncmp(name,"Pn",2) == 0)  return TRUE;
385        if  (strncmp(name,"Pg",2) == 0)  return TRUE;
386        if  (strcmp(name,"(Pn)") == 0)  return TRUE;
387        if  (strcmp(name,"(Pg)") == 0)  return TRUE;
388        if  (strcmp(name,"Pdiff") == 0)  return TRUE;
389        if  (strcmp(name,"(Pdiff)") == 0)  return TRUE;
390        if  (strcmp(name,"Pdif") == 0)  return TRUE;
391        if  (strcmp(name,"(Pdif)") == 0)  return TRUE;
392
393        return FALSE;
394
395} /* end of EfPhaseIsFirstArrival */
396
397
398
399/*---------------------------------------------------------------------------*/
400
401
402
403static void EfFinalCheck( FILE *out )
404
405/* Prints result of final check to output file 'out'
406 *
407 * parameters of routine
408 * FILE       *out;         input; pointer to output file
409 */
410{
411        /* local variables */
412        char     expl[cBcLongStrLth+1];   /* explanation string */
413
414        /* executable code */
415
416        /* fprintf( out, "--> final check\n" ); */
417        if  (efv_showinitials)
418                sprintf( expl, "%s (%s)", efv_infile, efv_analyst );
419        else
420                strcpy( expl, efv_infile );
421
422        /* a first arrival must be specified */
423        if  (!efv_first_arrival)
424                fprintf( out, "F014 %s: no first arrival phase selected\n", expl );
425
426        /* beam with uncorrected slowness at GRF is at least suspicious */
427        if  (efv_only_grf && efv_uncorr_loc)
428                fprintf( out, "F015 %s: GRF location with uncorrected beam at phase %s\n",
429                        expl, efv_uncorrphase );
430
431        /* only GRF picked but no GRA1? */
432        if  (efv_only_grf && !efv_gra1_picked)
433                fprintf( out, "F035 %s: only GRF picks but none at GRA1\n", expl );
434
435        /* no telex phase */
436        /* disabled 12-Nov-2008 K.S.
437        if  (!efv_telex_phase)
438                fprintf( out, "F016 %s: no telex phase found\n", expl );
439        */
440
441        /* mb determined? */
442        if  (efv_subtype == EVTYPE_P && !efv_mb_found)
443                fprintf( out, "F017 %s: no Magnitude mb determined for P phase\n", expl );
444
445        /* Ms determined? */
446        if  (efv_subtype >= EVTYPE_P && efv_depth <= 50.0 && efv_magn > 5.6 &&
447                !efv_ms_found)
448                fprintf( out, "F018 %s: no Magnitude Ms determined for magn %3.1f event\n",
449                        expl, efv_magn );
450        if  (efv_subtype == EVTYPE_PKP && efv_depth <= 50.0 && efv_pkp_count > 3 &&
451                !efv_ms_found)
452                fprintf( out, "F019 %s: no Magnitude Ms determined for PKP event\n",
453                        expl );
454        if  (efv_ms_found && efv_depth > 50.0)
455                fprintf( out, "F020 %s: Magnitude Ms determined for depth %5.1f event\n",
456                        expl, efv_depth );
457
458        /* Mbb determined? */
459        if  (efv_mb_found && efv_magn > 6.2 && !efv_mbb_found)
460                fprintf( out, "F021 %s: no Broadband Magnitude determined for magn %3.1f event\n",
461                        expl, efv_magn );
462
463        /* depth ok? */
464        if  (efv_evtype == EvcEventTypeMining && Abs(efv_depth-1.0) > 0.1)
465                fprintf( out, "F022 %s: depth %5.1f given for mining event\n",
466                        expl, efv_depth );
467        if  (efv_evtype == EvcEventTypeBlast && Abs(efv_depth) > 0.1)
468                fprintf( out, "F023 %s: depth %5.1f given for blast\n",
469                        expl, efv_depth );
470        if  (efv_evtype == EvcEventTypeTeleQuake && (efv_depth < 10.0))
471                fprintf( out, "F026 %s: depth %5.1f for teleseismic event\n",
472                        expl, efv_depth );
473        if  ((efv_evtype == EvcEventTypeLocalQuake
474                || efv_evtype == EvcEventTypeRegioQuake) && (efv_depth < 2.0))
475                fprintf( out, "F027 %s: depth %5.1f for local/regional event\n",
476                        expl, efv_depth );
477        if  (efv_evtype == EvcEventTypeLocalQuake && efv_depth > 50.0)
478                fprintf( out, "F036 %s: depth %5.1f given for local event\n",
479                        expl, efv_depth );
480        if  (efv_evtype == EvcEventTypeUndefined)
481                fprintf( out, "F024 %s: no event type specified\n", expl );
482        if  (efv_depth < -99.0)
483                fprintf( out, "F025 %s: no depth specified\n", expl );
484
485        /* proper event type? */
486        if  (efv_foundPg && efv_evtype == EvcEventTypeTeleQuake)
487                fprintf( out, "F028 %s: improper teleseismic event type: found Pg\n",
488                        expl );
489        if  (efv_foundPnSn && efv_evtype == EvcEventTypeTeleQuake)
490                fprintf( out, "F029 %s: improper teleseismic event type: found Pn/Sn\n",
491                        expl );
492        if  (efv_foundTele && efv_evtype == EvcEventTypeLocalQuake)
493                fprintf( out, "F030 %s: improper local event type: found tele phase\n",
494                        expl );
495        if  (efv_foundTele && efv_evtype == EvcEventTypeRegioQuake)
496                fprintf( out, "F031 %s: improper regional event type: found tele phase\n",
497                        expl );
498
499        /* check distance and event type */
500        if  (efv_latitude != 0.0 || efv_longitude != 0.0)  {
501                double dist, azim, bazim;   /* result of distance computation */
502                float fdist;    /* single precision distance */
503                mb_locdiff( (double)efv_reflat, (double)efv_reflon, (double)efv_latitude,
504                        (double)efv_longitude, &dist, &azim, &bazim );
505                fdist = dist;
506                if  (dist > 15.0 && efv_evtype == EvcEventTypeLocalQuake)
507                        fprintf( out, "F032 %s: improper local event type (distance %f deg)\n",
508                                expl, dist );
509                if  (dist > 30.0 && efv_evtype == EvcEventTypeRegioQuake)
510                        fprintf( out, "F033 %s: improper regional event type (distance %f deg)\n",
511                                expl, dist );
512                if  (dist < 25.0 && efv_evtype == EvcEventTypeTeleQuake)
513                        fprintf( out, "F034 %s: improper teleseismic event type (distance %f deg)\n",
514                                expl, dist );
515        } /*endif*/
516
517} /* end of EfFInalCheck */
518
519
520
521/*---------------------------------------------------------------------------*/
522
523
524
525static BOOLEAN EfIsLocalPhase( char phase[] )
526
527/* checks whether 'phase' is a local phase
528 *
529 * parameters of routine
530 * char       phase[];         input; phase name
531 */
532{
533        /* local variables */
534
535        /* executable code */
536
537        if  (strcmp(phase,"Pg") == 0)  return TRUE;
538        if  (strcmp(phase,"Pn") == 0)  return TRUE;
539        if  (strcmp(phase,"Sg") == 0)  return TRUE;
540        if  (strcmp(phase,"Sn") == 0)  return TRUE;
541        if  (strcmp(phase,"(Pg)") == 0)  return TRUE;
542        if  (strcmp(phase,"(Pn)") == 0)  return TRUE;
543        if  (strcmp(phase,"(Sg)") == 0)  return TRUE;
544        if  (strcmp(phase,"(Sn)") == 0)  return TRUE;
545
546        return FALSE;
547
548} /* end of EfIsLocalPhase */
549
550
551
552/*---------------------------------------------------------------------------*/
553
Note: See TracBrowser for help on using the repository browser.