Context Navigation

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 */
62	static float siv_tol_trav=(-1.0);
63	static float siv_tol_slow=(-1.0);
64	static float siv_tol_azim=(-1.0);
65	/* optimum travel time fit for identify_phase */
66	static float siv_trav_dt=1000.0;
67	static char siv_opt_phase[cBcShortStrLth+1]; /* best fitting phase */
68
69
70
71	static char *siv_evlists[] = {
72	/"neic-finger",/ "sed-redpuma", "emsc", "isc", ""
73	};
74
75	#ifdef XXX
76	static char *siv_infosrc[] = {
77	/"neic-a",/ "sed-rp", "emsc", ""
78	};
79	#endif
80
81	/* phase list for phase identification */
82	static 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	*/
93	static 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 */
104	static 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 );
107	static float si_get_maxtol( char phase[] );
108	static 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
117	void 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
174	float 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
243	float 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
316	float 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
438	float 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
511	void 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
557	void 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
594	TSyBoolean 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
634	void 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
807	void 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,"ANOR") == 0 \|\| strcmp(tmpstr,"ASED") == 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
1251	void 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
1290	static 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
1397	static 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
1447	void 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
1479	void 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.

Download in other formats: