source: SH_SHM/trunk/source/residual.c @ 16

Revision 16, 11.6 KB checked in by marcus, 15 years ago (diff)

r1 | svn | 2007-12-13 11:10:29 +0100 (Do, 13 Dez 2007) | 2 lines

Initial import
Line 
1
2/* file residual.c
3 *      ==========
4 *
5 * version 4, 22-May-2006
6 *
7 * Reading interpolated residual times from residual tables.
8 * K. Stammler, 3-Feb-2001
9 */
10
11
12/*
13 *
14 *  SeismicHandler, seismic analysis software
15 *  Copyright (C) 1992,  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#include <stdio.h>
36#include <string.h>
37#include "basecnst.h"
38#ifdef BC_INC_STDLIB
39#include BC_INC_STDLIB
40#endif
41#include "sysbase.h"
42#include "erusrdef.h"
43#include "utusrdef.h"
44#include "residual.h"
45
46
47
48/* global variables */
49TRsTable    *vrs_roottable=NULL;   /* pointer to root table */
50
51
52/* prototypes of local routines */
53static void RsEnlistTable( TRsTable *tab );
54static void RsDumpTable( TRsTable *tab );
55
56
57
58/*----------------------------------------------------------------------------*/
59
60
61
62void RsReadTables( char fname[], TSyStatus *status )
63
64/* Read in residual tables from file
65 *
66 * parameters of routine
67 * char       file[];       input; filename
68 * STATUS     *status;      output; return status
69 */
70{
71        /* local variables */
72        char     locname[cBcFileLth+1]; /* local file name */
73        char     *env;                /* pointer to environment variable */
74        FILE     *fp;                 /* pointer to input file */
75        char     line[cBcLineLth+1];  /* current line of file */
76        char     tmpstr[cBcLineLth+1];/* scratch string */
77        TRsTable *curtab;             /* pointer to current table */
78        int      i;                   /* counter */
79
80        /* executable code */
81
82        /* generate filename if not given */
83        if  (*fname == '\0' || strcmp(fname,"default") == 0)  {
84                env = (char *)getenv( "SH_INPUTS" );
85                if  ((strlen(env)+strlen(cRsDefaultName)) > cBcFileLth)  {
86                        *status = eRsLongFilename;
87                        return;
88                } /*endif*/
89                strcpy( locname, env );
90                strcat( locname, cRsDefaultName );
91        } else {
92                if  (strlen(fname) > cBcFileLth)  {
93                        *status = eRsLongFilename;
94                        return;
95                } /*endif*/
96                strcpy( locname, fname );
97        } /*endif*/
98
99        /* open file */
100        fp = fopen( locname, "r" );
101        if  (fp == NULL)  {
102                *status = eRsFileNotFound;
103                return;
104        } /*endif*/
105
106        /* read through file */
107        curtab = NULL;
108        while  (fgets(line,cBcLineLth,fp) != NULL)  {
109
110                /* ignore comment and empty lines */
111                if  (*line == '!')  continue;
112                if  (*line == '\n')  continue;
113
114                /* create new table on station */
115                if  (strncmp(line,"station:",8) == 0)  {
116                        if  (curtab != NULL)  {
117                                *status = eRsTableNotClosed;
118                                sy_deallocmem( curtab );
119                                return;
120                        } /*endif*/
121                        /* create new table */
122                        curtab = (TRsTable *)sy_allocmem( 1, sizeof(TRsTable), status );
123                        if  (SySevere(status))  return;
124                        RsInitTable( curtab );
125                        sscanf( line, "station: %s\n", tmpstr );
126                        if  (strlen(tmpstr) > cRsStationLth)  {
127                                *status = eRsLongName;
128                                return;
129                        } /*endif*/
130                        ut_cap( tmpstr );
131                        strcpy( curtab->station, tmpstr );
132                } /*endif*/
133
134                /* read phase list */
135                if  (strncmp(line,"phases:",7) == 0)  {
136                        if  (curtab == NULL)  {
137                                *status = eRsSyntaxError;
138                                return;
139                        } /*endif*/
140                        sscanf( line, "phases: %s\n", tmpstr );
141                        if  (strlen(tmpstr) > cRsPhaseListLth-2)  {
142                                *status = eRsLongPhaseList;
143                                return;
144                        } /*endif*/
145                        curtab->phases[0] = ',';
146                        strcpy( curtab->phases+1, tmpstr );
147                        strcat( curtab->phases, "," );
148                } /*endif*/
149
150                /* read slowness range */
151                if  (strncmp(line,"slowness:",9) == 0)  {
152                        if  (curtab == NULL)  {
153                                *status = eRsSyntaxError;
154                                return;
155                        } /*endif*/
156                        if  (sscanf(line,"slowness:%f\n-%f\n",
157                                &(curtab->minslow),&(curtab->maxslow)) != 2)  {
158                                *status = eRsReadSlowness;
159                                return;
160                        } /*endif*/
161                } /*endif*/
162
163                /* read length, azim & resid */
164                if  (strncmp(line,"entries:",8) == 0)  {
165                        if  (curtab == NULL)  {
166                                *status = eRsSyntaxError;
167                                return;
168                        } /*endif*/
169                        if  (sscanf(line,"entries:%d\n",&(curtab->tablth)) != 1)  {
170                                *status = eRsReadTablth;
171                                return;
172                        } /*endif*/
173                        fscanf( fp, "%s", tmpstr );
174                        if  (strcmp(tmpstr,"azim:") != 0)  {
175                                *status = eRsSyntaxError;
176                                err_setcontext( " ## 'entries' must be followed by 'azim' " );
177                        } /*endif*/
178                        curtab->azim = (float *)sy_allocmem(
179                                curtab->tablth, (int)sizeof(float), status );
180                        if  (SySevere(status))  return;
181                        for  (i=0; i<(curtab->tablth); i++)  {
182                                if  (fscanf(fp,"%f\n",(curtab->azim)+i) != 1)  {
183                                        *status = eRsSyntaxError;
184                                        err_setcontext( " ## reading azim list " );
185                                        return;
186                                } /*endif*/
187                        } /*endfor*/
188                        fscanf( fp, "%s", tmpstr );
189                        if  (strcmp(tmpstr,"resid:") != 0)  {
190                                *status = eRsSyntaxError;
191                                err_setcontext( " ## 'azim' must be followed by 'resid' " );
192                        } /*endif*/
193                        curtab->resid = (float *)sy_allocmem(
194                                curtab->tablth, (int)sizeof(float), status );
195                        if  (SySevere(status))  return;
196                        for  (i=0; i<(curtab->tablth); i++)  {
197                                if  (fscanf(fp,"%f\n",(curtab->resid)+i) != 1)  {
198                                        *status = eRsSyntaxError;
199                                        err_setcontext( " ## reading resid list " );
200                                        return;
201                                } /*endif*/
202                        } /*endfor*/
203                } /*endif*/
204
205                /* insert table when finished */
206                if  (strncmp(line,"end",3) == 0)  {
207                        /* check for completeness */
208                        if  (curtab == NULL)  {*status = eRsIncompleteTable; return;}
209                        if  (curtab->station[0] == '\0')
210                                {*status = eRsIncompleteTable; err_setcontext("station"); return;}
211                        if  (curtab->phases[0] == '\0')
212                                {*status = eRsIncompleteTable; err_setcontext("phases"); return;}
213                        if  (curtab->minslow < 0.0)
214                                {*status = eRsIncompleteTable; err_setcontext("minslow"); return;}
215                        if  (curtab->maxslow < 0.0)
216                                {*status = eRsIncompleteTable; err_setcontext("maxslow"); return;}
217                        if  (curtab->tablth <= 0)
218                                {*status = eRsIncompleteTable; err_setcontext("tablth"); return;}
219                        if  (curtab->azim == NULL)
220                                {*status = eRsIncompleteTable; err_setcontext("azim"); return;}
221                        if  (curtab->resid == NULL)
222                                {*status = eRsIncompleteTable; err_setcontext("resid"); return;}
223                        RsEnlistTable( curtab );
224                        curtab = NULL;
225                } /*endif*/
226
227        } /*endwhile*/
228
229        if  (curtab != NULL)  {
230                *status = eRsIncompleteTable;
231                err_setcontext( curtab->station );
232                sy_deallocmem( curtab );
233                return;
234        } /*endif*/
235
236        fclose( fp );
237
238} /* end of RsReadTables */
239
240
241
242/*----------------------------------------------------------------------------*/
243
244
245
246void RsInitTable( TRsTable *tab )
247
248/* initializes table
249 *
250 * parameters of routine
251 * TRsTable   *tab;      output; table to initialize
252 */
253{
254        /* executable code */
255
256        tab->station[0] = '\0';
257        tab->minslow = cRsEmptySlowness;
258        tab->maxslow = cRsEmptySlowness;
259        tab->tablth = 0;
260        tab->azim = NULL;
261        tab->resid = NULL;
262        tab->next = NULL;
263
264} /* end of RsInitTable */
265
266
267
268/*----------------------------------------------------------------------------*/
269
270
271
272void RsDumpTables( void )
273
274/* Dumps all tables to stdout
275 *
276 * no parameters
277 */
278{
279        /* local variables */
280        TRsTable *t;     /* moving pointer */
281
282        /* executable code */
283
284        printf( "\nResidual Table List Dump\n\n" );
285
286        t = vrs_roottable;
287        while  (t != NULL)  {
288                RsDumpTable( t );
289                t = t->next;
290        } /*endwhile*/
291
292        printf( "\n\n" );
293
294} /* end of RsDumpTables */
295
296
297
298/*----------------------------------------------------------------------------*/
299
300
301
302float RsResidual( char station[], char phase[], float slowness, float azimuth,
303        TSyStatus *status )
304
305/* Returns residual time in s for given station, phase, slowness and
306 * back-azimuth
307 *
308 * parameters of routine
309 * char       station[];      input; station name
310 * char       phase[];        input; phase name
311 * float      slowness;       input; slowness in s/deg
312 * float      azimuth;        input; back-azimuth n deg
313 * TSyStatus  status;         output; return status
314 */
315{
316        /* local variables */
317        char     locstat[cRsStationLth+1]; /* local station name */
318        TRsTable *t;           /* matching table */
319        TSyBoolean found;      /* found right table */
320        int      i;            /* index */
321        float    resid;        /* interpolated residual */
322        char     pstr[cRsPhaseListLth+1]; /* search string */
323
324        /* executable code */
325
326        /* create phase search string with commas */
327        *pstr = ',';
328        strcpy( pstr+1, phase );
329        strcat( pstr, "," );
330
331        /* create uppercase station name */
332        if  (strlen(station) > cRsStationLth)  {
333                *status = eRsLongName;
334                return;
335        } /*endif*/
336        strcpy( locstat, station );
337        ut_cap( locstat );
338
339        /* find matching table */
340        t = vrs_roottable;
341        while  (t != NULL)  {
342                found = TRUE;
343                if  (strcmp(t->station,locstat) != 0)  found = FALSE;
344                if  (strstr(t->phases,pstr) == NULL)  found = FALSE;
345                if  (slowness < t->minslow)  found = FALSE;
346                if  (slowness >= t->maxslow)  found = FALSE;
347                if  (found)  break;
348                t = t->next;
349        } /*endwhile*/
350
351        if  (t == NULL)  {
352                *status = eRsNoTableFound;
353                err_setcontext( " ## station: " );
354                err_setcontext( locstat );
355                err_setcontext( "  phase: " );
356                err_setcontext( phase );
357                err_setcontext( "  slowness: " );
358                err_setcontext_r( slowness );
359                return 0.0;
360        } /*endif*/
361
362        /*
363        printf( "\nfound table:\n" );
364        RsDumpTable( t );
365        */
366
367        if  (azimuth < t->azim[0] || azimuth > t->azim[(t->tablth)-1])  {
368                *status = eRsAzimOutOfRange;
369                err_setcontext( " ## station: " );
370                err_setcontext( locstat );
371                err_setcontext( "  slowness: " );
372                err_setcontext_r( slowness );
373                err_setcontext( "  azimuth: " );
374                err_setcontext_r( azimuth );
375                return 0.0;
376        } /*endif*/
377
378        i = 1;
379        while  (t->azim[i] < azimuth)  {
380                i++;
381                if  (i >= t->tablth)  {
382                        *status = eRsProgramBug;
383                        err_setcontext( " ## find azimuth index" );
384                        return;
385                } /*endif*/
386        } /*endwhile*/
387
388        /*
389        printf( "azimuth %5.1f is between %5.1f and %5.1f\n",
390                azimuth, t->azim[i-1], t->azim[i] );
391        */
392
393        resid = t->resid[i-1] + (azimuth - t->azim[i-1])
394                * ((t->resid[i] - t->resid[i-1]) / (t->azim[i] - t->azim[i-1]));
395
396        return resid;
397
398} /* end of RsResidual */
399
400
401
402/*----------------------------------------------------------------------------*/
403
404
405
406static void RsEnlistTable( TRsTable *tab )
407
408/* Appends table to the table list
409 *
410 * parameters of routine
411 * TRsTable   *tab;       input; table to append
412 */
413{
414        /* local variables */
415        TRsTable *t;              /* moving pointer */
416
417        /* executable code */
418
419        if  (vrs_roottable == NULL)  {
420                vrs_roottable = tab;
421                return;
422        } /*endif*/
423
424        t = vrs_roottable;
425        while  (t->next != NULL)
426                t = t->next;
427
428        t->next = tab;
429
430} /* end of RsEnlistTable */
431
432
433
434/*----------------------------------------------------------------------------*/
435
436
437
438static void RsDumpTable( TRsTable *tab )
439
440/* Dumps residual table to stdout
441 *
442 * parameters of routine
443 * TRsTable   *tab;        input; table to dump
444 */
445{
446        /* local varaibles */
447        int      i;                /* counter */
448
449        /* executable code */
450
451        if  (tab == NULL)  {
452                printf( "table == NULL\n" );
453                return;
454        } /*endif*/
455
456        printf( "station:   >%s<\n", tab->station );
457        printf( "phases:    %s\n", tab->phases );
458        printf( "slowness:  %4.2f - %4.2f\n", tab->minslow, tab->maxslow );
459        printf( "entries:   %d\n", tab->tablth );
460        printf( "azim:      " );
461        for  (i=0; i<(tab->tablth); i++)
462                printf( "%5.1f ", tab->azim[i] );
463        printf( "\n" );
464        printf( "resid:     " );
465        for  (i=0; i<(tab->tablth); i++)
466                printf( "%5.2f ", tab->resid[i] );
467        printf( "\n\n" );
468
469} /* RsDumpTable */
470
471
472
473/*----------------------------------------------------------------------------*/
474
Note: See TracBrowser for help on using the repository browser.