source: SH_SHM/trunk/source/motif/seismics.c @ 334

Revision 334, 46.5 KB checked in by marcus, 12 years ago (diff)

r173 | klaus | 2011-02-23 12:24:31 +0100 (Mi, 23 Feb 2011) | 1 line

format of EMSC file has changed some time ago; adjusted reading program for identify_phase

Line 
1
2/* file seismics.c
3 *      ==========
4 *
5 * version 34, 16-Jun-2006
6 *
7 * seismic utilities
8 * K. Stammler, 25-Mar-93
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 <math.h>
39#include "basecnst.h"
40#ifdef BC_INC_STDLIB
41#include BC_INC_STDLIB
42#endif
43#include "sysbase.h"
44#include "erusrdef.h"
45#include "tcusrdef.h"
46#include "ptusrdef.h"
47#include "earthloc.h"
48#include "glusrdef.h"
49#include "seismics.h"
50#include "globalparams.h"
51
52
53
54/* for identify phases, two operation modes */
55#define cSiModeSearch 1
56#define cSiModeSelect 2
57
58/* maximum length of source string in char */
59#define cSiSrcLength 7
60
61/* actual tolerance values for identify_phase, values taken from parameter file */
62static float siv_tol_trav=(-1.0);
63static float siv_tol_slow=(-1.0);
64static float siv_tol_azim=(-1.0);
65/* optimum travel time fit for identify_phase */
66static float siv_trav_dt=1000.0;
67static char  siv_opt_phase[cBcShortStrLth+1];  /* best fitting phase */
68
69
70
71static char *siv_evlists[] = {
72        /*"neic-finger",*/ "sed-redpuma", "emsc", "isc", ""
73};
74
75#ifdef XXX
76static char *siv_infosrc[] = {
77        /*"neic-a",*/ "sed-rp", "emsc", ""
78};
79#endif
80
81/* phase list for phase identification */
82static char *siv_idphases[] = {
83        "P", "Pdiff", "PcP", "PP", "PPP", "PKPdf", "PKPab", "PKPbc", "SP",
84        "PKKPdf", "PKKPbc", "PKKPab", "SKPdf", "SKPab", "SKKPdf",
85        "S", "Sdiff", "ScS", "SS", "SKSdf", "SSS", "SKKSac", "SKKSdf", ""
86};
87
88/* ml distance correction (C.F.RICHTER (1958): Elementary
89 * Seismology. p 342, up to 600 km epicentral distance,
90 * G.SCHNEIDER (1975): Erdbeben. p.338, from 700 up to 1000 km
91 * epicentral distance)
92 */
93static float siv_ml_sigma[] = {
94        1.4,1.4,1.5,1.6,1.7,1.9,2.1,2.3,2.4,2.5,2.6,2.7,2.8,2.8,2.8,
95        2.8,2.9,2.9,3.0,3.0,3.0,3.1,3.1,3.2,3.2,3.3,3.3,3.4,3.4,3.5,
96        3.5,3.6,3.65,3.7,3.7,3.8,3.8,3.9,3.9,4.0,4.0,4.1,4.1,4.2,4.2,
97        4.3,4.3,4.3,4.4,4.4,4.5,4.5,4.5,4.6,4.6,4.6,4.6,4.7,4.7,4.7,
98        4.7,4.8,4.8,4.8,4.8,4.8,4.9,4.9,4.9,4.9,4.9,5.2,5.4,5.5,5.7
99};
100#define cSiMlSigmaLth 75
101
102
103/* prototypes of local routines */
104static void si_identify_phase_step( int mode, char onset[], char station[],
105        float slowness, float azimuth, char phase[], float *lat, float *lon,
106        float *dep, char origin[], char infosrc[], TSyStatus *status );
107static float si_get_maxtol( char phase[] );
108static void si_search_events( char type[], TSyBoolean newlist, char onset[],
109        float backsec, float trav[], float lat[], float lon[], float depth[],
110        char infsrc[][cSiSrcLength+1], int *evno, TSyStatus *status );
111
112
113/*---------------------------------------------------------------------------*/
114
115
116
117void si_fit_distance( char phase[], float slowness, float depth,
118        float tol, float *distance, STATUS *status )
119
120/* fits distance from slowness and depth
121 *
122 * parameters of routine
123 * char       phase[];      input; phase name
124 * float      slowness;     input; slowness in sec/deg
125 * float      depth;        input; depth in km
126 * float      tol;          input; tolerance in slowness
127 * float      *distance     output; distance in deg
128 */
129{
130        /* local variables */
131        float    min_dist, max_dist;   /* distance window to search */
132        float    curr_dist;            /* current distance */
133        float    curr_slow;            /* current slowness */
134        int      i;                    /* loop counter */
135
136        /* executable code */
137
138        min_dist = 2.0;
139        max_dist = 105.0;  /* for P only useful */
140
141        i = 0;
142        FOREVER  {
143                curr_dist = (min_dist + max_dist) / 2.0;
144                if  (GpGetInt(cGpI_debug_level) > 1)
145                        printf( "SHM-dbg2: curr dist %5.1f\n", curr_dist );
146                curr_slow = pt_slowness( phase, curr_dist, depth, status );
147                if  (Severe(status))  return;
148                if  (Abs(curr_slow-slowness) <= tol)  break;
149                if  (curr_slow > slowness)  {
150                        min_dist = curr_dist;
151                } else {
152                        max_dist = curr_dist;
153                } /*endif*/
154                if  (++i > 50)  {
155                        *status = SIE_NOCONVERG;
156                        return;
157                } /*endif*/
158        } /*endfor*/
159        *distance = curr_dist;
160
161} /* end of si_fit_distance */
162
163
164
165/*---------------------------------------------------------------------------*/
166
167
168
169#define MAGN_ML_MIN_DIST 16
170#define MAGN_ML_MAX_DIST 120
171
172
173
174float si_magn_mb( float ampl, float period, float distance, STATUS *status )
175
176/* Returns magnitude mb.  Passed amplitude and period must origin from a
177 * Z component.
178 *
179 * parameters of routine
180 * float      ampl;        input; amplitude in nm
181 * float      period;      input; amplitude in sec
182 * float      distance;    input; epicentral distance in deg
183 * STATUS     *status;     output; return status
184 *                         returns value of mb
185 */
186{
187
188        /* local data array */
189        static float mb_sigma[] = {
190                5.9,5.9,5.9,6.0,6.0,6.1,6.2,6.3,6.3,6.5,6.4,6.5,6.6,
191                6.6,6.6,6.7,6.7,6.7,6.7,6.7,6.6,6.5,6.5,6.4,6.4,6.5,6.5,6.5,6.5,
192                6.7,6.8,6.9,6.9,6.8,6.7,6.7,6.7,6.7,6.8,6.8,6.8,6.8,6.8,6.8,6.8,
193                6.6,7.0,6.9,7.0,7.0,7.0,7.0,7.0,7.0,6.9,6.9,6.9,6.9,6.8,6.8,6.9,
194                6.9,6.9,6.8,6.7,6.8,6.9,7.0,7.0,7.0,6.9,6.9,7.1,7.0,7.0,7.1,7.1,
195                7.2,7.1,7.2,7.3,7.4,7.5,7.5,7.4,7.3,7.4,7.5,7.6,7.7,7.8,7.9,7.9,
196                8.0,8.1,8.2,8.2,8.4,8.6,8.7,8.8,8.9,9.0,9.0,9.0
197        };
198
199        /* local variables */
200        int      idist;     /* rounded distance */
201        float    tmp;       /* scratch */
202
203        /* executable code */
204
205        idist = Nint( distance );
206        if  (idist < MAGN_ML_MIN_DIST || idist > MAGN_ML_MAX_DIST)  {
207                *status = SIE_MAGN_DIST_OOR;
208                return 0.0;
209        } /*endif*/
210
211        if  (ampl <= 0.0 || period <= 0.0)  {
212                *status = SIE_VALUE_OOR;
213                return 0.0;
214        } /*endif*/
215
216        /* rescale to microns */
217        ampl /= 1000.0;
218
219        tmp = log10( (double)(ampl/period) ) + mb_sigma[idist-MAGN_ML_MIN_DIST];
220
221        if  (GpGetInt(cGpI_debug_level) > 1)
222                printf( "SHM-dbg2: mb %3.1f, distance %5.1f, period %5.2f, ampl %f\n",
223                        tmp, distance, period, ampl );
224        return tmp;
225
226        /* this didn't work (in dbxtool):
227        return ( log10(ampl/period) + mb_sigma[idist-MAGN_ML_MIN_DIST] );
228        */
229
230} /* end of si_magn_mb */
231
232
233
234#undef MAGN_ML_MIN_DIST
235#undef MAGN_ML_MAX_DIST
236
237
238
239/*---------------------------------------------------------------------------*/
240
241
242
243float si_magn_ms_plain( float ampl, float period, float distance,
244        STATUS *status )
245
246/* returns MS magnitude, using Prague formula
247 *
248 * parameters of routine
249 * float      ampl;        input; amplitude in nm
250 * float      period;      input; amplitude in sec
251 * float      distance;    input; epicentral distance in deg
252 * STATUS     *status;     output; return status
253 *                         returns value of MS
254 */
255{
256        /* local data array */
257        float ms_min_period[] = {
258                4,4,5,5,5,6,6,7,7,9,10,12,13,14,15,16,16,16,17,17,18,18
259        };
260
261        /* local variables */
262        float    tmp;       /* scratch */
263        int      index;     /* index in ms_min_period */
264        int      idist;     /* rounded distance */
265        float    minper;    /* minimum period */
266
267        /* executable code */
268
269        /* check distance & period */
270        idist = Nint( distance );
271        if  (idist < 2)  {
272                *status = SIE_MAGN_DIST_OOR;
273                return 0.0;
274        } else if  (idist > 140)  {
275                minper = 18.0;
276        } else {
277                index = (idist > 10) ? Nint(distance/10.0)+7 : Nint(distance)-2;
278                minper = ms_min_period[index];
279        } /*endif*/
280        if  (period < minper)  {
281                *status = SIE_VALUE_OOR;
282                return 0.0;
283        } /*endif*/
284
285        /* check amplitude */
286        if  (ampl <= 0.0)  {
287                *status = SIE_VALUE_OOR;
288                return 0.0;
289        } /*endif*/
290
291        /* rescale to microns */
292        ampl /= 1000.0;
293
294        tmp = log10((double)(ampl/period)) + 1.66*log10((double)distance) + 3.3;
295
296        if  (GpGetInt(cGpI_debug_level) > 0)
297                printf( "SHM-dbg1: MS %3.1f, distance %5.1f, period %5.2f, ampl %f\n",
298                        tmp, distance, period, ampl );
299        return tmp;
300
301} /* end of si_magn_ms_plain */
302
303
304
305/*---------------------------------------------------------------------------*/
306
307
308
309#define MAGN_MS_MIN_DIST 1
310#define MAGN_MS_MAX_DIST 99
311#define MAGN_MS_MIN_PER 10
312#define MAGN_MS_MAX_PER 40
313
314
315
316float si_magn_ms_m_b( int type, float ampl, float period, float distance,
317        STATUS *status )
318
319/* returns MS magnitude, using Marshall-Basham formula
320 *
321 * parameters of routine
322 * int        type;        input; path type
323 * float      ampl;        input; amplitude in nm
324 * float      period;      input; amplitude in sec
325 * float      distance;    input; epicentral distance in deg
326 * STATUS     *status;     output; return status
327 *                         returns value of MS
328 */
329{
330        /* local data arrays */
331        /* Marshal-Basham MS: distance correction table */
332        float ms_delta[] = {
333                0.17,0.35,0.57,0.67,0.78,0.84,0.90,0.95,0.98,1.02,
334                1.05,1.08,1.11,1.13,1.15,1.17,1.19,1.21,1.23,1.25,
335                1.27,1.29,1.31,1.32,1.34,1.36,1.38,1.40,1.41,1.43,
336                1.44,1.45,1.47,1.48,1.50,1.52,1.54,1.55,1.56,1.57,
337                1.59,1.61,1.62,1.64,1.65,1.66,1.68,1.70,1.71,1.72,
338                1.74,1.75,1.76,1.77,1.78,1.80,1.82,1.83,1.84,1.85,
339                1.87,1.89,1.90,1.91,1.92,1.93,1.94,1.95,1.96,1.97,
340                1.98,1.99,2.01,2.02,2.03,2.04,2.05,2.06,2.07,2.08,
341                2.09,2.10,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,
342                2.18,2.19,2.20,2.21,2.22,2.23,2.24,2.25,2.26
343        };
344
345        /* Marshal Basham MS: path correction */
346        float ms_path[][31] = {
347                /* path: c-na */
348                {-0.75,-0.67,-0.61,-0.53,-0.46,-0.38,-0.30,-0.24,
349                 -0.16,-0.08, 0.00, 0.01, 0.03, 0.04, 0.05, 0.07,
350                  0.09, 0.11, 0.13, 0.14, 0.16, 0.17, 0.18, 0.20,
351                  0.21, 0.23, 0.25, 0.27, 0.28, 0.29, 0.31},
352                /* path: c-eu */
353                {-0.30,-0.27,-0.24, 0.21,-0.18,-0.15,-0.13,-0.10,
354                 -0.07,-0.04, 0.00, 0.03, 0.05, 0.07, 0.11, 0.14,
355                  0.18, 0.22, 0.24, 0.27, 0.30, 0.32, 0.33, 0.34,
356                  0.35, 0.36, 0.37, 0.38, 0.39, 0.40, 0.41},
357                /* path: c-o */
358                { 0.00, 0.01, 0.03, 0.04, 0.05, 0.07, 0.08, 0.09,
359                  0.10, 0.12, 0.13, 0.15, 0.16, 0.17, 0.18, 0.20,
360                  0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28,
361                  0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35},
362                /* path: 0 */
363                { 0.50, 0.45, 0.38, 0.33, 0.27, 0.20, 0.15, 0.09,
364                  0.03,-0.03,-0.09,-0.08,-0.07,-0.06,-0.05,-0.04,
365                 -0.03,-0.03,-0.02,-0.01, 0.00, 0.01, 0.02, 0.03,
366                  0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10}
367        };
368
369
370        /* local variables */
371        int      idist;        /* rounded distance */
372        int      iper;         /* rounded period */
373        int      path_idx;     /* path index */
374        float    tmp;          /* scratch */
375
376        /* executable code */
377
378        /* check distance */
379        idist = Nint( distance );
380        if  (idist < MAGN_MS_MIN_DIST || idist > MAGN_MS_MAX_DIST)  {
381                *status = SIE_MAGN_DIST_OOR;
382                return 0.0;
383        } /*endif*/
384
385        /* check period */
386        iper = Nint( period );
387        if  (iper < MAGN_MS_MIN_PER || iper > MAGN_MS_MAX_PER)  {
388                *status = SIE_VALUE_OOR;
389                return 0.0;
390        } /*endif*/
391
392        /* check amplitude */
393        if  (ampl <= 0.0)  {
394                *status = SIE_VALUE_OOR;
395                return 0.0;
396        } /*endif*/
397
398        /* get path index */
399        switch  (type)  {
400        case SIC_MAGN_MS_C_NA:  path_idx = 0;   break;
401        case SIC_MAGN_MS_C_EU:  path_idx = 1;   break;
402        case SIC_MAGN_MS_C_O:   path_idx = 2;   break;
403        case SIC_MAGN_MS_O:     path_idx = 3;   break;
404        default:  *status = SIE_ILL_MAGN_TYPE;  return 0.0;
405        } /*endswitch*/
406
407        /* rescale to microns */
408        ampl /= 1000.0;
409
410        tmp = log10((double)ampl) + ms_delta[idist-MAGN_MS_MIN_DIST]
411                + ms_path[path_idx][iper-MAGN_MS_MIN_PER] + 3.0;
412
413        if  (GpGetInt(cGpI_debug_level) > 0)
414                printf( "SHM-dbg1: MS %3.1f, distance %5.1f, period %5.2f, ampl %f\n",
415                        tmp, distance, period, ampl );
416        return tmp;
417
418} /* end of si_magn_ms_m_b */
419
420
421
422#undef MAGN_MS_MIN_DIST
423#undef MAGN_MS_MAX_DIST
424#undef MAGN_MS_MIN_PER
425#undef MAGN_MS_MAX_PER
426
427
428
429/*---------------------------------------------------------------------------*/
430
431
432
433#define MAGN_ML_MIN_DIST 1.0
434#define MAGN_ML_MAX_DIST 1000.0
435
436
437
438float si_magn_ml( float ampl, float distance, STATUS *status )
439
440/* Returns magnitude ml.  Amplitude and period must be read from a WOOD-ANDERSON
441 * simulated seismogram, horizontal component.
442 *
443 * parameters of routine
444 * float      ampl;        input; amplitude of WOOD-ANDERSON simulation
445 * float      distance;    input; epicentral distance in km
446 * STATUS     *status;     output; return status
447 *                         returns value of ml
448 */
449{
450        /* ml distance correction (C.F.RICHTER (1958): Elementary
451         * Seismology. p 342, up to 600 km epicentral distance,
452         * G.SCHNEIDER (1975): Erdbeben. p.338, from 700 up to 1000 km
453         * epicentral distance)
454         */
455
456        /* local variables */
457        int      idist;        /* rounded distance */
458        float    delta;        /* step size */
459        int      ioffset;      /* index offset */
460        float    grad;         /* gradient */
461        int      itmp;         /* scratch */
462        float    tmp;          /* scratch */
463
464        /* executable code */
465
466        /* check distance */
467        if  (distance < MAGN_ML_MIN_DIST || distance >= MAGN_ML_MAX_DIST)  {
468                *status = SIE_MAGN_DIST_OOR;
469                return 0.0;
470        } /*endif*/
471
472        /* get step size */
473        if  (distance > 600.0)  {
474                delta = 100.0;
475                ioffset = 64;
476        } else if  (distance > 100.0)  {
477                delta = 10.0;
478                ioffset = 10;
479        } else {
480                delta = 5.0;
481                ioffset = 0;
482        } /*endif*/
483
484        itmp = Nint( floor(distance/delta) );
485        idist = itmp + ioffset;
486        grad = ( siv_ml_sigma[idist+1] - siv_ml_sigma[idist] ) / delta;
487        tmp = distance - (float)itmp*delta;
488        tmp = siv_ml_sigma[idist] + tmp*grad;
489
490        tmp += log10((double)ampl);
491
492        if  (GpGetInt(cGpI_debug_level) > 0)
493                printf( "SHM-dbg1: ml %3.1f, distance %5.1f, ampl %f\n",
494                        tmp, distance, ampl );
495
496        return tmp;
497
498} /* end of si_magn_ml */
499
500
501
502#undef MAGN_ML_MIN_DIST
503#undef MAGN_ML_MAX_DIST
504
505
506
507/*---------------------------------------------------------------------------*/
508
509
510
511void si_ext_location( float lat, float lon, float depth, char refstation[],
512        float *distance, float *azimuth, float *p_slowness, STATUS *status )
513
514/* Returns distance, back-azimuth & P-slowness (if possible) of given
515 * hypocenter to refstation.
516 *
517 * parameters of routine
518 * float      lat, lon;       input; epicenter of event;
519 * float      depth;          input; event depth in km
520 * char       refstation[];   input; name of receiving station
521 * float      *distance;      output; distance in deg
522 * float      *azimuth;       output; back-azimuth in deg
523 * float      *p_slowness;    output; P-slowness or 0.0
524 * STATUS     *status;        output; return status
525 */
526{
527        /* local variables */
528        GLT_STATINF *statinf;    /* station info */
529        double      d_dist, d_azim, d_bazim;
530        float       p_slow;     /* local value of P-slowness */
531
532        /* executable code */
533
534        statinf = gl_store_station( refstation, TRUE, status );
535        if  (Severe(status))  return;
536        mb_locdiff( statinf->lat, statinf->lon, lat, lon,
537                &d_dist, &d_azim, &d_bazim );
538        *distance = d_dist;
539        *azimuth = d_bazim;
540        p_slow = pt_slowness( "P", *distance, depth, status );
541        if  (Severe(status))  {
542                *status = BC_NOERROR;
543                if  (GpGetInt(cGpI_debug_level) > 1)
544                        printf( "SHM-dbg2: P-slowness could not be computed\n" );
545        } else {
546                *p_slowness = p_slow;
547        } /*endif*/
548
549} /* end of si_ext_location */
550
551
552
553/*---------------------------------------------------------------------------*/
554
555
556
557void si_get_location( float distance, float azimuth, char refstation[],
558        float *lat, float *lon, int *ferid, char fername[], STATUS *status )
559
560/* computes location from given distance and azimuth relative to given
561 * reference station.  'fername' must have minimum length BC_LINELTH
562 *
563 * parameters of routine
564 * float     distance;       input; distance in deg
565 * float     azimuth;        input; back-azimuth in deg
566 * char      refstation[];   input; name of reference station
567 * float     *lat, *lon;     output; computed source location
568 * int       *ferid;         output; Flinn-Engdahl region number
569 * char      fername[];      output; name of Flinn-Engdahl region(>= BC_LINELTH)
570 * STATUS    *status;        output; return status
571 */
572{
573        /* local variables */
574        GLT_STATINF *statinf;    /* station info */
575
576        /* executable code */
577
578        statinf = gl_store_station( refstation, TRUE, status );
579        if  (Severe(status))  return;
580
581        mb_sphereadd( statinf->lat, statinf->lon, distance, azimuth, lat, lon );
582        mb_ferindex( *lat, *lon, ferid, status );
583        if  (Severe(status))  return;
584        mb_fername( *ferid, BC_LINELTH, fername, status );
585
586} /* end of si_get_location */
587
588
589
590/*---------------------------------------------------------------------------*/
591
592
593
594TSyBoolean si_is_first_onset( char mode[], char phase[] )
595
596/* Returns TRUE if 'phase' is a first onset phase.
597 *
598 * parameters of routine
599 * char       mode[];       input; 'tele' or not
600 * char       phase[];      input; name of phase
601 *                          returns TRUE if phase is first onset
602 */
603{
604        /* local variables */
605
606        /* executable code */
607
608        if  (strcmp(phase,"P") == 0)  return TRUE;
609        if  (strcmp(phase,"PKPab") == 0)  return TRUE;
610        if  (strcmp(phase,"PKPbc") == 0)  return TRUE;
611        if  (strcmp(phase,"PKPdf") == 0)  return TRUE;
612        if  (strcmp(phase,"PKP") == 0)  return TRUE;
613
614        if  (strcmp(mode,"tele") == 0)  return FALSE;
615
616        if  (strcmp(phase,"Pg") == 0)  return TRUE;
617        if  (strcmp(phase,"Pn") == 0)  return TRUE;
618
619        return FALSE;
620
621} /* end of si_is_first_onset */
622
623
624
625/*---------------------------------------------------------------------------*/
626
627
628/* length of search window in sec */
629#define SEARCHWDW 7200.0
630/* maximum number of event expected in search window */
631#define MAXEVNO 30
632
633
634void si_match_location( char type[], char phase[], char onset[], char station[],
635        float slowness, float azimuth, TSyBoolean newlist, float *lat, float *lon,
636        float *dep, char origin[], char infsrc[], TSyStatus *status )
637
638/* Finds matching locations in external event lists
639 *
640 * parameters of routine
641 * char       type[];      input; which event list
642 * char       phase[];     input; name of phase
643 * char       onset[];     input; onset time of phase
644 * char       station[];   input; station of phase reading
645 * float      slowness;    input; slowness of phase
646 * float      azimuth;     input; azimuth of phase
647 * TSyBoolean newlist;     input; retrieve new event list or use old one
648 * float      *lat;        output; latitude found
649 * float      *lon;        output; longitude found
650 * float      *dep;        output; depth found
651 * char       origin[];    output; origin time
652 * char       infsrc[];    output; information source
653 * TSyStatus  *status;     output; return status
654 */
655{
656        /* local variables */
657        float    trav[MAXEVNO];            /* travel time of phase */
658        float    evlat[MAXEVNO];           /* latitude of events */
659        float    evlon[MAXEVNO];           /* longitude of events */
660        float    depth[MAXEVNO];           /* depth of events */
661        char     infsrc_l[MAXEVNO][cSiSrcLength+1]; /* info source */
662        int      evno;                     /* number of events found */
663        int      i;                        /* counter */
664        GLT_STATINF *statinf;    /* station info */
665        double   d_dist, d_azim, d_bazim;  /* difference epi - station */
666        float    f_dist, f_azim;           /* difference epi - station */
667        float    ctrav, cslow;             /* current traveltime and slowness */
668        TSyStatus locstat;                 /* local status */
669        TSyBoolean found;                  /* event found */
670        float    diff_trav, diff_azim, diff_slow; /* differences */
671        float    obs_azimuth;              /* observed azimuth, +180 deg ? */
672
673        /* executable code */
674
675        /* initialise tolerance values if not yet done */
676        if  (siv_tol_trav < 0.0)  siv_tol_trav = GpGetFloat( cGpF_idphases_tol_trav );
677        if  (siv_tol_slow < 0.0)  siv_tol_slow = GpGetFloat( cGpF_idphases_tol_slow );
678        if  (siv_tol_azim < 0.0)  siv_tol_azim = GpGetFloat( cGpF_idphases_tol_azim );
679
680        /* find events in search window */
681        si_search_events( type, newlist, onset, SEARCHWDW, trav, evlat, evlon,
682                depth, infsrc_l, &evno, status );
683        if  (SySevere(status))  return;
684
685        statinf = gl_store_station( station, TRUE, status );
686        if  (SySevere(status))  return;
687
688        /* may change siv_tol_trav for surface waves, since LR is not searched
689    * in identify_phase */
690        if  (strcmp(phase,"LR") == 0)
691                siv_tol_trav = GpGetFloat( cGpF_idphases_tol_travsurf );
692
693        found = FALSE;
694        for  (i=0; i<evno; i++)  {
695
696                /* may have been changed in last loop */
697                obs_azimuth = azimuth;
698
699                /* compute traveltime, azimuth and slowness for each possible event */
700                mb_locdiff( statinf->lat, statinf->lon, evlat[i], evlon[i],
701                        &d_dist, &d_azim, &d_bazim );
702                f_dist = d_dist; f_azim = d_bazim;
703                locstat = cBcNoError;
704                if  (strcmp(phase,"PKP") == 0)  {
705                        ctrav = pt_travel( "PKPdf", f_dist, depth[i], &locstat );
706                        if  (locstat != cBcNoError)
707                                ctrav = pt_travel( "PKPbc", f_dist, depth[i], &locstat );
708                        if  (locstat != cBcNoError)
709                                ctrav = pt_travel( "PKPab", f_dist, depth[i], &locstat );
710                        if  (locstat != cBcNoError)  ctrav = -100.0;
711                        cslow = 0.0;
712                } else {
713                        ctrav = pt_travel( phase, f_dist, depth[i], &locstat );
714                        if  (locstat != cBcNoError)  ctrav = -100.0;
715                        cslow = pt_slowness( phase, f_dist, depth[i], &locstat );
716                        if  (locstat != cBcNoError)  cslow = -100.0;
717                        if  (cslow < 0.0 && cslow > -99.0)  {
718                                obs_azimuth = azimuth + 180.0;
719                                if  (obs_azimuth > 360.0)  obs_azimuth -= 360.0;
720                                cslow = -cslow;
721                        } /*endif*/
722                } /*endif*/
723
724                /* check for matching parameters */
725                /* difference in travel time */
726                if  (ctrav == -100.0)  {
727                        /* if no travel time could be computed, set to large value */
728                        diff_trav=9999.9;
729                } else {
730                        diff_trav = fabs( ctrav - trav[i] );
731                } /*endif*/
732                /* difference in slowness */
733                diff_slow = 0.0;
734                if  (strcmp(phase,"PKP") != 0 && strcmp(phase,"LR") != 0)
735                        if  (cslow == -100.0)  {
736                                /* if no slowness could be computed, set to large value */
737                                diff_slow = 99.9;
738                        } else {
739                                diff_slow = fabs( cslow - slowness );
740                        } /*endif*/
741                /* difference in azimuth */
742                if  (slowness < 2.5)  {
743                        /* there is no resolution in azimuth on very steep incidence angles */
744                        diff_azim = 0.0;
745                } else {
746                        diff_azim = fabs( (float)f_azim - obs_azimuth );
747                        if  (diff_azim > 180.0)  diff_azim = fabs( diff_azim - 360.0 );
748                } /*endif*/
749                if  (siv_tol_trav == 0.0)  {
750                        printf( "   %11s %5s: dt=%6.1f da=%5.1f ds=%4.1f",
751                                type, phase, diff_trav, diff_azim, diff_slow );
752                        printf( " ev%d: (%6.2f,%7.2f) %5.1f\n",
753                                i, evlat[i], evlon[i], depth[i] );
754                } /*endif*/
755
756#ifdef XXX
757                /* debug output */
758                printf( "%f %f %5.f %s: t:%f a:%f\n", evlat[i], evlon[i], depth[i],
759                        infsrc_l[i], diff_trav, diff_azim );
760#endif
761
762                /* store best diff_trav, if slowness and azimuth match */
763                if  (diff_slow < siv_tol_slow && diff_azim < siv_tol_azim)
764                        if  (diff_trav < siv_trav_dt)  {
765                                siv_trav_dt = diff_trav;
766                                strcpy( siv_opt_phase, phase );
767                        } /*endif*/
768
769                if  (diff_trav < siv_tol_trav && diff_slow < siv_tol_slow
770                        && diff_azim < siv_tol_azim)  {
771                        if  (found)  {
772                                *status = SIE_MORE_MATCHES;
773                                return;
774                        } /*endif*/
775                        found = TRUE;
776                        *lat = evlat[i];
777                        *lon = evlon[i];
778                        *dep = depth[i];
779                        strcpy( infsrc, infsrc_l[i] );
780                        tc_tadd( onset, -trav[i], origin, status );
781                        if  (SySevere(status))  return;
782                        printf( "   %s: accept phase %s on event %s (%5.2f,%6.2f) %5.1f from %s\n",
783                                type, phase, origin, *lat, *lon, *dep, infsrc );
784                        break;
785                } /*endif*/
786
787        } /*endfor*/
788
789        if  (!found)  *status = SIE_NO_EV_FOUND;
790
791} /* end of si_match_location */
792
793
794
795/*---------------------------------------------------------------------------*/
796
797
798#define TMPFILE "search_ev.dat"
799#define SEDFILE "sed_redpuma.dat"
800#define EMSCFILE "search_emsc.dat"
801#define ISCFILE "search_isc.dat"
802#define FINGERADDRESS "quake@ghtftp.cr.usgs.gov"
803#define SEDPAGE "http://www.seismo.ethz.ch/redpuma/redpuma_ami_list.html"
804/* #define EMSCPAGE "http://www.emsc-csem.org/cgi-bin/ALERT_all_messages.sh" */
805
806
807void si_search_events( char type[], TSyBoolean newlist, char onset[],
808        float backsec, float trav[], float lat[], float lon[], float depth[],
809        char infsrc[][cSiSrcLength+1], int *evno, TSyStatus *status )
810
811/* Searches for events in given event lists within a time window
812 *
813 * parameters of routine
814 * char       type[];         input; which list
815 * TSyBoolean newlist;        input; retrieve new list or use old one
816 * char       onset[];        input; end of search window (onset time)
817 * float      backsec;        input; time window in sec backward from onset
818 * float      trav[];         output; travel time (onset - origin)
819 * float      lat[];          output; latitudes
820 * float      lon[];          output; longitudes
821 * float      depth[];        output; depths
822 * char       infsrc[][cSiSrcLength+1]; output; information source
823 * int        *evno;          output; number of events
824 * TSyStatus  *status;        output; return status
825 */
826{
827        /* local variables */
828        char     cmd[cBcLongStrLth+1];     /* shell command */
829        char     tmpfile[cBcFileLth+1];    /* scratch file */
830        char     *env;                     /* pointer to environment */
831        FILE     *fp;                      /* pointer to file */
832        char     line[cBcLongStrLth+1];    /* current line of file */
833        char     *lptr;                    /* point to char in line */
834        TSyBoolean datalines;              /* line is a data line */
835        char     origtime[cBcTimeLth+1];   /* origin time string */
836        float    curtrav;                  /* travel time for current event */
837        char     tmpstr[cBcLineLth+1];     /* scratch string */
838        char     tmpstr2[cBcLineLth+1];    /* second scratch string */
839        char     excl_agency[cBcLineLth+1];/* exclusive agency */
840        char     tchar;                    /* scratch character */
841        float    curlat, curlon;           /* epicentre of current line */
842        float    curdep;                   /* current depth */
843        float    curmag;                   /* current magnitude */
844        int      i;                        /* counter */
845
846        /* executable code */
847
848        *evno = 0;
849
850        strcpy( excl_agency, GpGetString(cGpS_exclusive_agency) );
851
852        if  (strcmp(type,"neic-finger") == 0)  {
853
854                /* name of scratch file */
855                env = (char *)getenv( "SH_SCRATCH" );
856                sprintf( tmpfile, "%s/%s", env, TMPFILE );
857                if  (newlist)  {
858                        sprintf( cmd, "rm %s", tmpfile );
859                        if  (GpGetInt(cGpI_debug_level) > 2)
860                                printf( "SHM-dbg3: executing: %s\n", cmd );
861                        system( cmd );
862                        sprintf( cmd, "finger %s > %s", FINGERADDRESS, tmpfile );
863                        if  (GpGetInt(cGpI_debug_level) > 2)
864                                printf( "SHM-dbg3: executing: %s\n", cmd );
865                        system( cmd );
866                } /*endif*/
867                fp = fopen( tmpfile, "r" );
868                if  (fp == NULL)  {
869                        *status = SIE_OPEN_READ;
870                        err_setcontext( " ## file " );
871                        err_setcontext( tmpfile );
872                        return;
873                } /*endif*/
874                datalines = FALSE;
875                while  (fgets(line,cBcLongStrLth,fp) != NULL)  {
876                        if  (strncmp(line,"yy/mm/dd hh:mm:ss",17) == 0)  {
877                                datalines = TRUE;
878                                continue;
879                        } /*endif*/
880                        if  (!datalines)  continue;
881                        if  (*line == '\n' || *line == ' ')  continue;
882                        strncpy( origtime, line, 17 );
883                        origtime[17] = '\0';
884                        origtime[8] = '_';
885                        curtrav = tc_tdiff( onset, origtime, status );
886                        if  (SySevere(status))  {fclose(fp); return;}
887                        if  (curtrav >= 0.0 && curtrav <= backsec)  {
888                                if  (*evno == MAXEVNO)  {
889                                        *status = SIE_TOO_MANY_EV;
890                                        fclose( fp );
891                                        return;
892                                } /*endif*/
893                                trav[*evno] = curtrav;
894                                strncpy( tmpstr, line+19, 5 );
895                                tmpstr[5] = '\0';
896                                sscanf( tmpstr, "%f", lat+(*evno) );
897                                if  (line[24] == 'S')  lat[*evno] = -lat[*evno];
898                                strncpy( tmpstr, line+26, 6 );
899                                tmpstr[6] = '\0';
900                                sscanf( tmpstr, "%f", lon+(*evno) );
901                                if  (line[32] == 'W')  lon[*evno] = -lon[*evno];
902                                strncpy( tmpstr, line+34, 5 );
903                                tmpstr[5] = '\0';
904                                sscanf( tmpstr, "%f", depth+(*evno) );
905                                strcpy( infsrc[*evno], "neic-f" );
906                                (*evno)++;
907                        } /*endif*/
908                } /*endwhile*/
909                fclose( fp );
910
911        } else if  (strcmp(type,"sed-redpuma") == 0)  {
912
913                /* name of scratch file */
914                env = (char *)getenv( "SH_SCRATCH" );
915                sprintf( tmpfile, "%s/%s", env, SEDFILE );
916                if  (newlist)  {
917                        sprintf( cmd, "rm %s", tmpfile );
918                        if  (GpGetInt(cGpI_debug_level) > 2)
919                                printf( "SHM-dbg3: executing: %s\n", cmd );
920                        system( cmd );
921                        env = (char *)getenv( "SH_UTIL" );
922                        sprintf( cmd, "%s/get_html_text.csh %s %s", env, SEDPAGE, tmpfile );
923                        if  (GpGetInt(cGpI_debug_level) > 2)
924                                printf( "SHM-dbg3: executing: %s\n", cmd );
925                        system( cmd );
926                } /*endif*/
927                fp = fopen( tmpfile, "r" );
928                if  (fp == NULL)  {
929                        *status = SIE_OPEN_READ;
930                        err_setcontext( " ## file " );
931                        err_setcontext( tmpfile );
932                        return;
933                } /*endif*/
934                datalines = FALSE;
935                while  (fgets(line,cBcLongStrLth,fp) != NULL)  {
936                        if  (strncmp(line," D a t e  Time (UTC)",20) == 0)  {
937                                datalines = TRUE;
938                                continue;
939                        } /*endif*/
940                        if  (!datalines)  continue;
941                        if  (*line == '\n')  continue;
942                        if  (strlen(line) < 40)  continue;
943                        if  (line[0] < '0' || line[0] > '9')  continue;
944                        if  (line[9] != ' ')  continue;
945                        origtime[0] = line[0];
946                        origtime[1] = line[1];
947                        origtime[2] = '-';
948                        origtime[3] = line[2];
949                        origtime[4] = line[3];
950                        origtime[5] = line[4];
951                        origtime[6] = '-';
952                        strncpy( origtime+7, line+5, 4 );
953                        origtime[11] = '_';
954                        strncpy( origtime+12, line+10, 10 );
955                        origtime[22] = '\0';
956                        strncpy( tmpstr, line+44, 5);
957                        tmpstr[5] = '\0';
958                        if  (*excl_agency != '\0' && strcasecmp(excl_agency,tmpstr) != 0)  {
959                                /* if exclusive agency is set, accept only this */
960                                curtrav = -1.0;
961                        } else if  (strcmp(tmpstr,"A*NOR") == 0  || strcmp(tmpstr,"A*SED") == 0)  {
962                                curtrav = -1.0;  /* throw away */
963                        } else {
964                                curtrav = tc_tdiff( onset, origtime, status );
965                        } /*endif*/
966                        if  (SySevere(status))  {fclose(fp); return;}
967                        if  (curtrav >= 0.0 && curtrav <= backsec)  {
968                                if  (*evno == MAXEVNO)  {
969                                        *status = SIE_TOO_MANY_EV;
970                                        fclose( fp );
971                                        return;
972                                } /*endif*/
973                                /* reformat source string */
974                                tchar = *tmpstr;
975                                if  (strcmp(tmpstr+2,"EMS") == 0)  {
976                                        sprintf( tmpstr, "EMSC-%c", tchar );
977                                } else if  (strcmp(tmpstr+2,"GSR") == 0)  {
978                                        sprintf( tmpstr, "GSRC-%c", tchar );
979                                } else if  (strcmp(tmpstr+2,"NEI") == 0)  {
980                                        sprintf( tmpstr, "NEIC-%c", tchar );
981                                } else if  (strcmp(tmpstr+2,"SED") == 0)  {
982                                        sprintf( tmpstr, "SED-%c", tchar );
983                                } else if  (strcmp(tmpstr+2,"ODC") == 0)  {
984                                        sprintf( tmpstr, "ODC-%c", tchar );
985                                } /*endif*/
986                                strcpy( infsrc[*evno], tmpstr );
987                                trav[*evno] = curtrav;
988                                strncpy( tmpstr, line+21, 4 );
989                                tmpstr[4] = '\0';
990                                sscanf( tmpstr, "%f", lat+(*evno) );
991                                if  (line[25] == 'S')  lat[*evno] = -lat[*evno];
992                                strncpy( tmpstr, line+27, 5 );
993                                tmpstr[5] = '\0';
994                                sscanf( tmpstr, "%f", lon+(*evno) );
995                                if  (line[32] == 'W')  lon[*evno] = -lon[*evno];
996                                strncpy( tmpstr, line+33, 3 );
997                                tmpstr[3] = '\0';
998                                sscanf( tmpstr, "%f", depth+(*evno) );
999                                (*evno)++;
1000                        } /*endif*/
1001                } /*endwhile*/
1002                fclose( fp );
1003
1004        } else if  (strcmp(type,"emsc") == 0)  {
1005
1006                /* name of scratch file */
1007                env = (char *)getenv( "SH_SCRATCH" );
1008                sprintf( tmpfile, "%s/%s", env, EMSCFILE );
1009                if  (newlist)  {
1010                        sprintf( cmd, "rm %s", tmpfile );
1011                        if  (GpGetInt(cGpI_debug_level) > 2)
1012                                printf( "SHM-dbg3: executing: %s\n", cmd );
1013                        system( cmd );
1014                        env = (char *)getenv( "SH_UTIL" );
1015                        sprintf( cmd, "%s/quakeml/emsc.py -o %s",
1016                                env, tmpfile );
1017                        if  (GpGetInt(cGpI_debug_level) > 2)
1018                                printf( "SHM-dbg3: executing: %s\n", cmd );
1019                        system( cmd );
1020                } /*endif*/
1021                fp = fopen( tmpfile, "r" );
1022                if  (fp == NULL)  {
1023                        *status = SIE_OPEN_READ;
1024                        err_setcontext( " ## file " );
1025                        err_setcontext( tmpfile );
1026                        return;
1027                } /*endif*/
1028                /*datalines = FALSE; */
1029                datalines = TRUE;  /* only data lines now, after preprocessing of files */
1030                while  (fgets(line,cBcLongStrLth,fp) != NULL)  {
1031#ifdef XXX
1032                        if  (strncmp(line,"  DATE       TIME",17) == 0)  {
1033                                datalines = TRUE;
1034                                continue;
1035                        } /*endif*/
1036#endif
1037                        if  (!datalines)  continue;
1038                        lptr = line;
1039                        while  (*lptr == ' ')  lptr++;
1040                        if  (*lptr == '\n')  continue;
1041                        if  (strlen(lptr) < 56)  continue;
1042                        if  (lptr[0] < '0' || lptr[0] > '9')  continue;
1043                        if  (lptr[4] != '-')  continue;
1044                        if  (lptr[7] != '-')  continue;
1045                        if  (lptr[10] != ' ')  continue;
1046                        if  (lptr[11] != ' ')  continue;
1047                        if  (lptr[17] != ':')  continue;
1048                        /* origin time is in a fixed format */
1049                        origtime[0] = lptr[8];
1050                        origtime[1] = lptr[9];
1051                        origtime[2] = ',';
1052                        origtime[3] = lptr[5];
1053                        origtime[4] = lptr[6];
1054                        origtime[5] = ',';
1055                        origtime[6] = lptr[0];
1056                        origtime[7] = lptr[1];
1057                        origtime[8] = lptr[2];
1058                        origtime[9] = lptr[3];
1059                        origtime[10] = ',';
1060                        origtime[11] = lptr[12];
1061                        origtime[12] = lptr[13];
1062                        origtime[13] = ',';
1063                        origtime[14] = lptr[15];
1064                        origtime[15] = lptr[16];
1065                        origtime[16] = ',';
1066                        origtime[17] = lptr[18];
1067                        origtime[18] = lptr[19];
1068                        origtime[19] = ',';
1069                        origtime[20] = lptr[21];
1070                        origtime[21] = '0';
1071                        origtime[22] = '0';
1072                        origtime[23] = '\0';
1073                        lptr += 23;   /* now at latitude */
1074                        sscanf( lptr, "%f", &curlat );
1075                        while  (*lptr != ' ' && *lptr != '\n')  lptr++;  /* skip over latitude */
1076                        while  (*lptr == ' ')  lptr++;  /* now at latitude sign */
1077                        if  (*lptr == 'S' || *lptr == 's')  curlat *= -1.0;
1078                        lptr++;
1079                        while  (*lptr == ' ')  lptr++;  /* now at longitude */
1080                        sscanf( lptr, "%f", &curlon );
1081                        while  (*lptr != ' ' && *lptr != '\n')  lptr++;  /* skip over longitude */
1082                        while  (*lptr == ' ')  lptr++;  /* now at longitude sign */
1083                        if  (*lptr == 'W' || *lptr == 'w')  curlon *= -1.0;
1084                        lptr++;
1085                        while  (*lptr == ' ')  lptr++;  /* now see what we find here */
1086                        if  (*lptr >= '0' && *lptr <= '9')  {
1087                                /* here we found depth */
1088                                sscanf( lptr, "%f", &curdep );
1089                                while  (*lptr != ' ' && *lptr != '\n')  lptr++;  /* skip over depth */
1090                                while  (*lptr == ' ')  lptr++;
1091                                if  (*lptr == 'G' && lptr[1] == ' ')  lptr += 2;
1092                                if  (*lptr == 'f' && lptr[1] == ' ')  lptr += 2;
1093                        } else {
1094                                curdep = 33.0;
1095                        } /*endif*/
1096                        while  (*lptr == ' ')  lptr++;  /* now we should come to magnitude */
1097                        if  (*lptr >= '0' && *lptr <= '9')  {
1098                                /* here we found magnitude */
1099                                sscanf( lptr, "%f", &curmag );
1100                        } else {
1101                                curmag = -1.0;
1102                        } /*endif*/
1103                        while  (*lptr != ' ' && *lptr != '\n')  lptr++;  /* skip over mangitude */
1104                        while  (*lptr == ' ')  lptr++;  /* now we should come to magnitude type */
1105                        while  (*lptr != ' ' && *lptr != '\n')  lptr++;  /* skip over mangitude type */
1106                        tchar = *lptr;
1107                        if  (tchar == 'A')  {
1108                                /* && (
1109                                strcmp(tmpstr,"SED") == 0 ||
1110                                strcmp(tmpstr,"ODC") == 0 ||
1111                                strcmp(tmpstr,"EMSC") == 0 ||
1112                                strcmp(tmpstr,"GFZ") == 0 ))  {  */
1113                                curtrav = -1.0;  /* do not trust these automatic solutions */
1114                        } else {
1115                                curtrav = tc_tdiff( onset, origtime, status );
1116                        } /*endif*/
1117                        if  (SySevere(status))  {fclose(fp); return;}
1118                        if  (curtrav >= 0.0 && curtrav <= backsec)  {
1119                                if  (*evno == MAXEVNO)  {
1120                                        *status = SIE_TOO_MANY_EV;
1121                                        fclose( fp );
1122                                        return;
1123                                } /*endif*/
1124                                lptr = line + strlen(line) - 1;
1125                                while  (*lptr == ' ' || *lptr == '\n')  lptr--;
1126                                i = (int)(lptr - line);
1127                                while  (lptr > line && *lptr != ' ')  lptr--;
1128                                lptr++;
1129                                i -= (int)(lptr - line);
1130                                i++;
1131                                if  (i > cSiSrcLength-2)  i = cSiSrcLength-2;
1132                                strncpy( infsrc[*evno], lptr, i );
1133                                infsrc[*evno][i] = '\0';
1134                                /* if exclusive agency is set, accept only this */
1135                                if  (*excl_agency != '\0'
1136                                        && strcasecmp(excl_agency,infsrc[*evno]) != 0)  {
1137                                        if  (GpGetInt(cGpI_debug_level) > 3)
1138                                                printf( "SHM-dbg4: skipped agency %s\n", infsrc[*evno] );
1139                                        continue;
1140                                } /*endif*/
1141                                if  (tchar == 'M' || tchar == 'A')  {
1142                                        sprintf( tmpstr, "-%c", tchar );
1143                                        strcat( infsrc[*evno], tmpstr );
1144                                } /*endif*/
1145                                trav[*evno] = curtrav;
1146                                lat[*evno] = curlat;
1147                                lon[*evno] = curlon;
1148                                depth[*evno] = curdep;
1149                                (*evno)++;
1150                        } /*endif*/
1151                } /*endwhile*/
1152                fclose( fp );
1153
1154        } else if  (strcmp(type,"isc") == 0)  {
1155
1156                /* name of scratch file */
1157                env = (char *)getenv( "SH_SCRATCH" );
1158                sprintf( tmpfile, "%s/%s", env, ISCFILE );
1159                if  (newlist)  {
1160                        sprintf( cmd, "rm %s", tmpfile );
1161                        if  (GpGetInt(cGpI_debug_level) > 2)
1162                                printf( "SHM-dbg3: executing: %s\n", cmd );
1163                        system( cmd );
1164                        env = (char *)getenv( "SH_UTIL" );
1165                        tc_tadd( onset, -backsec, tmpstr2, status );
1166                        sprintf( cmd, "%s/get_isc_events.csh %s %s %s 4.0",
1167                                env, tmpstr2, onset, tmpfile );
1168                        if  (GpGetInt(cGpI_debug_level) > 2)
1169                                printf( "SHM-dbg3: executing: %s\n", cmd );
1170                        system( cmd );
1171                } /*endif*/
1172                fp = fopen( tmpfile, "r" );
1173                if  (fp == NULL)  {
1174                        *status = SIE_OPEN_READ;
1175                        err_setcontext( " ## file " );
1176                        err_setcontext( tmpfile );
1177                        return;
1178                } /*endif*/
1179                while  (fgets(line,cBcLongStrLth,fp) != NULL)  {
1180                        /* data lines start with the year, so either 19.. or 20.. */
1181                        if  (*line < '1' || *line > '2')  continue;
1182                        origtime[0] = line[8];
1183                        origtime[1] = line[9];
1184                        origtime[2] = ',';
1185                        origtime[3] = line[5];
1186                        origtime[4] = line[6];
1187                        origtime[5] = ',';
1188                        origtime[6] = line[0];
1189                        origtime[7] = line[1];
1190                        origtime[8] = line[2];
1191                        origtime[9] = line[3];
1192                        origtime[10] = ',';
1193                        origtime[11] = line[11];
1194                        origtime[12] = line[12];
1195                        origtime[13] = ',';
1196                        origtime[14] = line[14];
1197                        origtime[15] = line[15];
1198                        origtime[16] = ',';
1199                        origtime[17] = line[17];
1200                        origtime[18] = line[18];
1201                        origtime[19] = ',';
1202                        origtime[20] = (line[20] == ' ' ? '0' : line[20]);
1203                        origtime[21] = (line[21] == ' ' ? '0' : line[21]);
1204                        origtime[22] = '0';
1205                        origtime[23] = '\0';
1206                        if  (sscanf(line+36,"%f",&curlat) != 1)  continue;
1207                        if  (sscanf(line+45,"%f",&curlon) != 1)  continue;
1208                        strncpy( tmpstr2, line+71, 6 );
1209                        tmpstr2[6] = '\0';
1210                        if  (sscanf(tmpstr2,"%f",&curdep) != 1)  curdep = 33.0;
1211                        curtrav = tc_tdiff( onset, origtime, status );
1212                        if  (SySevere(status))  {fclose(fp); return;}
1213                        if  (curtrav >= 0.0 && curtrav <= backsec)  {
1214                                if  (*evno == MAXEVNO)  {
1215                                        *status = SIE_TOO_MANY_EV;
1216                                        fclose( fp );
1217                                        return;
1218                                } /*endif*/
1219                                lptr = line + 118;
1220                                while  (*lptr == ' ' || *lptr != '\n')  lptr++;
1221                                i = strlen( lptr );
1222                                if  (lptr[i] == '\n')  lptr[i--] = '\0';
1223                                if  (i > cSiSrcLength)  i = cSiSrcLength;
1224                                strncpy( infsrc[*evno], lptr, i );
1225                                infsrc[*evno][i] = '\0';
1226                                /* if exclusive agency is set, accept only this */
1227                                if  (*excl_agency != '\0'
1228                                        && strcasecmp(excl_agency,infsrc[*evno]) != 0)
1229                                        continue;
1230                                trav[*evno] = curtrav;
1231                                lat[*evno] = curlat;
1232                                lon[*evno] = curlon;
1233                                depth[*evno] = curdep;
1234                                (*evno)++;
1235                        } /*endif*/
1236                } /*endwhile*/
1237                fclose( fp );
1238        } else {
1239                *status = SIE_ILL_EVLIST;
1240                return;
1241        } /*endif*/
1242
1243} /* end of si_search_events */
1244
1245
1246
1247/*---------------------------------------------------------------------------*/
1248
1249
1250
1251void si_identify_phase( char onset[], char station[],
1252        float slowness, float azimuth, char phase[], float *lat, float *lon,
1253        float *dep, char origin[], char infosrc[], TSyStatus *status )
1254
1255/* identifies phase from onset time using external event lists
1256 *
1257 * parameters of routine
1258 * char       onset[];    input; onset time of phase
1259 * char       station[];  input; station name
1260 * float      slowness;   input; slowness of phase
1261 * float      azimuth;    input; back azimuth of phase
1262 * char       phase[];    output; phase name
1263 * float      *lat;       output; epicenter latitude
1264 * float      *lon;       output; epicenter longitude
1265 * float      *dep;       output; source depth
1266 * char       origin[];   output; origin time
1267 * char       infosrc[];  output; information source (which epi list)
1268 * TSyStatus  *status;    output; return status
1269 */
1270{
1271        /* local variables */
1272        TSyStatus   locstat;    /* local status */
1273
1274        /* executable code */
1275
1276        locstat = cBcNoError;
1277        si_identify_phase_step( cSiModeSearch, onset, station, slowness, azimuth,
1278                phase, lat, lon, dep, origin, infosrc, &locstat );
1279        si_identify_phase_step( cSiModeSelect, onset, station, slowness, azimuth,
1280                phase, lat, lon, dep, origin, infosrc, status );
1281
1282} /* end of si_identify_phase */
1283
1284
1285
1286/*---------------------------------------------------------------------------*/
1287
1288
1289
1290static void si_identify_phase_step( int mode, char onset[], char station[],
1291        float slowness, float azimuth, char phase[], float *lat, float *lon,
1292        float *dep, char origin[], char infosrc[], TSyStatus *status )
1293
1294/* Called by si_identify_phase, identifies phase from onset time using
1295 * external event lists
1296 *
1297 * parameters of routine
1298 * int        mode;       input; operation mode
1299 * char       onset[];    input; onset time of phase
1300 * char       station[];  input; station name
1301 * float      slowness;   input; slowness of phase
1302 * float      azimuth;    input; back azimuth of phase
1303 * char       phase[];    output; phase name
1304 * float      *lat;       output; epicenter latitude
1305 * float      *lon;       output; epicenter longitude
1306 * float      *dep;       output; source depth
1307 * char       origin[];   output; origin time
1308 * char       infosrc[];  output; information source (which epi list)
1309 * TSyStatus  *status;    output; return status
1310 */
1311{
1312        /* local variables */
1313        static TSyBoolean lists_created=FALSE;  /* lists are on scratch disk */
1314        int      i, l;         /* counters */
1315        TSyBoolean newlist;    /* create new list */
1316        TSyStatus  locstat;    /* local status */
1317        float      maxtol;     /* maximum tolerance time */
1318
1319        /* executable code */
1320
1321        if  (mode == cSiModeSearch)
1322                printf( "*SHM: identify onset %s on %s, slo %4.1f, baz %5.1f\n",
1323                        onset, station, slowness, azimuth );
1324
1325        /* set travel time tolerance according to operation mode */
1326        if  (mode == cSiModeSearch)  {
1327                /* if in search mode, set tolerance to zero, so nothing is accepted */
1328                siv_tol_trav = 0.0;
1329                siv_trav_dt = 1000.0;
1330                *siv_opt_phase = '\0';
1331        } else {
1332                /* if a reasonably fitting event was found, set the tolerance to the
1333                   best one found */
1334                maxtol = si_get_maxtol( siv_opt_phase );
1335                if  (siv_trav_dt < maxtol)  {
1336                        siv_tol_trav = siv_trav_dt + 0.1;
1337                } else {
1338                        siv_tol_trav = GpGetFloat( cGpF_idphases_tol_trav );
1339                } /*endif*/
1340        } /*endif*/
1341        if  (GpGetInt(cGpI_debug_level) > 2)
1342                printf( "SHM-dbg3: dt tol set to %4.1f\n", siv_tol_trav );
1343
1344        /* loop all available lists from different agencies */
1345        l = 0;
1346        while  (siv_evlists[l][0] != '\0')  {
1347
1348                /* check all phases */
1349                i = 0;
1350                if  (strcmp(siv_evlists[l],"isc") == 0)  {
1351                        /* ISC has a list over many years, each event may need a new request */
1352                        newlist = (mode == cSiModeSearch);
1353                } else {
1354                        /* from redpuma and EMSC only recent events are copied, one request
1355          * per SHM session is enough */
1356                        newlist = !lists_created;
1357                } /*endif*/
1358                while  (siv_idphases[i][0] != '\0')  {
1359                        locstat = cBcNoError;
1360                        si_match_location( siv_evlists[l], siv_idphases[i], onset, station,
1361                                slowness, azimuth, newlist, lat, lon, dep, origin, infosrc,
1362                                &locstat );
1363                        if  (locstat == cBcNoError)  {
1364                                /* phase found */
1365                                strcpy( phase, siv_idphases[i] );
1366                                /*strcpy( infosrc, siv_infosrc[l] );*/
1367                                /* return only unspecified PKP and SKS */
1368                                if  (strncmp(phase,"PKP",3) == 0)  phase[3] = '\0';
1369                                if  (strncmp(phase,"PKKP",4) == 0)  phase[4] = '\0';
1370                                if  (strncmp(phase,"SKKP",4) == 0)  phase[4] = '\0';
1371                                if  (strncmp(phase,"SKS",3) == 0)  phase[3] = '\0';
1372                                siv_tol_trav = GpGetFloat( cGpF_idphases_tol_trav );
1373                                return;
1374                        } /*endif*/
1375                        newlist = FALSE;
1376                        i++;
1377                } /*endwhile*/
1378
1379                l++;
1380
1381        } /*endwhile*/
1382
1383        lists_created = TRUE;
1384
1385        /* nothing found */
1386        *status = locstat;
1387        siv_tol_trav = GpGetFloat( cGpF_idphases_tol_trav );
1388
1389} /* end of si_identify_phase_step */
1390
1391
1392
1393/*---------------------------------------------------------------------------*/
1394
1395
1396
1397static float si_get_maxtol( char phase[] )
1398
1399/* Returns maximum tolerance time dependent on phase
1400 *
1401 * parameters of routine
1402 * char       phase[];   input; phase to give tolerance time
1403 */
1404{
1405        /* executable code */
1406
1407        if  (strcmp(phase,"P") == 0)  {
1408                return 10.0;
1409        } else if  (strcmp(phase,"Pdiff") == 0)  {
1410                return 20.0;
1411        } else if  (strcmp(phase,"PcP") == 0)  {
1412                return 10.0;
1413        } else if  (strcmp(phase,"PP") == 0)  {
1414                return 10.0;
1415        } else if  (strcmp(phase,"PPP") == 0)  {
1416                return 20.0;
1417        } else if  (strcmp(phase,"SP") == 0)  {
1418                return 30.0;
1419        } else if  (strncmp(phase,"PKP",3) == 0)  {
1420                return 15.0;
1421        } else if  (strcmp(phase,"S") == 0)  {
1422                return 30.0;
1423        } else if  (strcmp(phase,"Sdiff") == 0)  {
1424                return 40.0;
1425        } else if  (strcmp(phase,"ScS") == 0)  {
1426                return 40.0;
1427        } else if  (strcmp(phase,"SS") == 0)  {
1428                return 40.0;
1429        } else if  (strncmp(phase,"SKS",3) == 0)  {
1430                return 50.0;
1431        } else if  (strcmp(phase,"SSS") == 0)  {
1432                return 50.0;
1433        } else {
1434                return GpGetFloat( cGpF_idphases_tol_trav );
1435        } /*endif*/
1436
1437} /* end of si_get_maxtol */
1438
1439
1440
1441/*---------------------------------------------------------------------------*/
1442
1443
1444#define ISC_END_DATE "31-Dec-2004"
1445
1446
1447void si_lookup_agency( char ponset[], char agency[] )
1448
1449/* returns agency (either EMS or ISC) where a phase with onset time ponset
1450 * can be looked for.  Currently very simple solution: all phases before
1451 * ISC_END_DATE are searched at ISC more recent ones on EMSC
1452 *
1453 *
1454 * parameters of routine
1455 * char       ponset[];     input; phase onset time
1456 * char       agency[];     output; name of agency
1457 */
1458{
1459        /* local variables */
1460        TSyStatus locstat;   /* local status */
1461
1462        /* excutable code */
1463
1464        locstat = cBcNoError;
1465        if  (tc_tdiff(ponset,ISC_END_DATE,&locstat) < 0.0)  {
1466                strcpy( agency, "isc" );
1467        } else {
1468                strcpy( agency, "emsc" );
1469        } /*endif*/
1470
1471} /* end of si_lookup_agency */
1472
1473
1474
1475/*---------------------------------------------------------------------------*/
1476
1477
1478
1479void si_read_table( char fname[], char tablename[], TSyStatus *status )
1480
1481/* Read a table from a file, e.g. sigma-table for ml
1482 *
1483 * parameters of routine
1484 * char       fname[];     input; name of input file
1485 * char       tablename[]; input; name of table
1486 * TSyStatus  *status;     output; return status
1487 */
1488{
1489        /* local variables */
1490        float    *table;              /* pointer to table */
1491        int      length;              /* length of table */
1492        FILE     *fp;                 /* pointer to input file */
1493        char     line[cBcLineLth+1];  /* current line of file */
1494        int      i;                   /* line counter */
1495        float    xval;                /* xvalue, not checked */
1496        char     locfname[cBcFileLth+1]; /* local filename */
1497        char     *env;                /* pointer to environment variable */
1498
1499        /* executable code */
1500
1501        if  (strcmp(tablename,"ml-sigma") == 0)  {
1502                table = siv_ml_sigma;
1503                length = cSiMlSigmaLth;
1504        } else {
1505                *status = SIE_UNKNOWN_TABLE;
1506                err_setcontext( tablename );
1507                return;
1508        } /*endif*/
1509
1510        fp = fopen( fname, "r" );
1511        if  (fp == NULL)  {
1512                env = (char *)getenv( "SH_INPUTS" );
1513                if  (env == NULL)  {*status = SIE_OPEN_TABLE; fclose(fp); return;}
1514                if  (strlen(env)+strlen(fname) > cBcFileLth)  {
1515                        *status = SIE_OPEN_TABLE;
1516                        err_setcontext( fname );
1517                } /*endif*/
1518                strcpy( locfname, env );
1519                strcat( locfname, fname );
1520                fp = fopen( locfname, "r" );
1521                if  (fp == NULL)  {
1522                        *status = SIE_OPEN_TABLE;
1523                        err_setcontext( fname );
1524                        return;
1525                } /*endif*/
1526        } /*endif*/
1527
1528        i = 0;
1529        while  (i < length)  {
1530                if  (fgets(line,cBcLineLth,fp) == NULL)  {
1531                        *status = SIE_READ_TABLE;
1532                        fclose( fp );
1533                        return;
1534                } /*endif*/
1535                if  (*line == '!' || *line == '\n')  continue;
1536                if  (sscanf( line, "%f %f", &xval, table+i ) != 2)  {
1537                        *status = SIE_READ_TABLE;
1538                        fclose( fp );
1539                        return;
1540                } /*endif*/
1541                i++;
1542        } /*endwhile*/
1543
1544        fclose( fp );
1545
1546} /* end of si_read_table */
1547
1548
1549
1550/*---------------------------------------------------------------------------*/
Note: See TracBrowser for help on using the repository browser.