source: trunk/src/seismic_processing/binder_ew/grid.c @ 7558

Revision 7558, 41.9 KB checked in by paulf, 10 months ago (diff)

cleaned up travel times in grid.c to be in seconds, not integer values converted back and forth

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/*
2 * grid.c : grid manager
3 */
4#include <math.h>
5#include <limits.h>
6#include <inttypes.h>
7#include <string.h>
8#include <stdio.h>
9#include <stdlib.h>
10#include <time.h>
11
12#include "chron3.h"
13#include "kom.h"
14#include "tlay.h"
15#include "site.h"
16#include "geo.h"
17#include "grid.h"
18#include "bind.h"
19#include "earthworm_simple_funcs.h"
20
21#define LAT(yy) ((yy) / grdFac.y + grdOrg.lat)
22#define LON(xx) ((xx) / grdFac.x + grdOrg.lon)
23#define X(LON) (grdFac.x * ((LON) - grdOrg.lon))
24#define Y(LAT) (grdFac.y * ((LAT) - grdOrg.lat))
25#define MAXDIS 100
26#define Panic(x) (logit("e","Panic point %d in grid.c\n",(x)),exit(-1))
27
28/* Function prototypes
29   *******************/
30int  GetInst( char *, unsigned char * );   /* getutil.c sys-independent  */
31void grid_init( void );
32int  grid_cmpsrt( const void *, const void * );
33int  grid_inc(int wt, unsigned char instid);
34void grid_dump(int ix, int iy, int iz);
35void grid_free( void );
36void bndr_link( long, long );
37void bndr_quake2k( long );
38
39/* Structure for sorting list of stacking picks on time
40 ******************************************************/
41typedef struct {
42     double t;       /* arrival time of pick                     */
43     short  ilist;   /* pre-sorted (original) index of this pick */
44} SRT;               /* structure added 960408:ldd */
45
46/* Output control
47   **************/
48static int logStack = 0;
49
50static int grid_debug = 0;
51/* the following two settings are OPTIONAL rejection parameters triggered by nearest_quake_dist command */
52static double nearest_quake_dist = 0.0;         /* epicentral km to nearest quake */
53#define NEAREST_QUAKE_TIME 2.0
54static double nearest_quake_time = NEAREST_QUAKE_TIME;  /* seconds to nearest quake */
55static int ignore_loc_code_in_same_station_decision = 0; /* in the is_same_SNL() function, if this is set to 1 ignore the location code as a tie breaker */
56
57/* Stacking parameters
58   *******************/
59#define MAXWT 10
60static int nWt = 0;
61static struct {
62        char          wt;       /* Weight code (0-4)                    */
63        char          inc;      /* Stack increment                      */
64        unsigned char instid;   /* Installation id; added 20020523:ldd  */
65} Wt[MAXWT];
66unsigned char InstWild = 0;     /* wildcard for installation id         */ /*020523:ldd*/
67char    InitInstWild   = 0;     /* initialization flag                  */ /*020523:ldd*/
68double facCluster      = 4.0;   /* Pick clustering cutoff               */
69double resGrid         = 1.0;   /* Tolerance for reassociating waifs    */
70double rmsGrid         = 1.0;   /* RMS threshold for new quake          */
71double rStack          = 100.0; /* Primary stacking dist. cutoff        */
72double GlitchNsec      = .035;  /* A glitch is defined as a group of at */ /*960408:ldd*/
73int    GlitchMinPk     = 4;     /*    least GlitchMinPk picks within    */ /*960408:ldd*/
74                                /*    GlitchNsec seconds                */ /*960408:ldd*/
75int    mThresh         = 6;     /* Association threshold                */
76long   mStack          = 20;    /* Maximum number of picks to stack     */
77unsigned int maskGrid  = 0;     /* Initialization mask                  */
78double  tStack          = 0.1;   /* Stack width (0.1 s)                  */
79int    stack_horizontals = 0;   /* set to 1 to allow horizontal components in stack - paulf */
80int    mPix;
81int    nPix;
82
83
84/* Grid definition parameters
85 ****************************/
86static int initGrid = 0;
87static int nX, nY, nZ;          /* Dimensions of stacking buffer        */
88static int nXS, nYS;            /* Dimensions of trav time template     */
89float dTime;            /* Time granularity (seconds)                   */
90float dSpace;           /* Space granularity (seconds)                  */
91int nThresh = 5;        /* Minimum stack value to associate             */
92int nFocus  = 5;        /* Maximum hits to be considered focused        */
93CART    grdTop;         /* Heigth, Width and Depth of grid              */
94CART    grdStp;         /* Grid spacing, km                             */
95CART    grdFac;         /* Geocentric to cartesian factor (km/degree)   */
96GEO     grdOrg;         /* Origin = (0,0,0) in cartesian coordinates    */
97GEO     grdMax;         /* Geocentric upper bound on grid               */
98
99/* Travel-time template
100 **********************/
101/* Trav[ix][iy][iz] : This array has  the travel time from the source
102                to the center of Stack[0][0][0] to the center of
103                Stack[ix][iy][iz]. Index is calculated as
104                iz + nZ*iy + nYS*nZ*ix (times are in .1 secs)           */
105double *lpTrav;
106int lpTravMax; /* the maximum index in lpTrav[]  */
107
108/* Stacking Grid
109 ***************/
110/* Grid[ix][iy][iz]:  The space time coordinate of a box in the
111                stack accumulator is taken to be the center of the box in
112                space-time. Index calculated as iz + nZ*iy + nY*nZ*ix.
113                The grid should be thought of as being
114                draped over the pick initiating the association         */
115short *lpGrid;
116int lpGridMax;
117static time_t sae=0;
118
119/* Functions to return index into 1-D array based on 3-D indices
120 ***************************************************************/
121#define IXGRID(IX, IY, IZ) ((IZ) + nZ * (IY) + nY * nZ * (IX))
122int IXTRAV(int ix, int iy, int iz)
123{
124        if(ix < 0)      ix = -ix;
125        if(iy < 0)      iy = -iy;
126        return (iz + (nZ*iy) + (nYS*nZ*ix));
127}
128
129/********************************************************************************
130 * grid_init() : Initialize grid calculations                                   *
131 *              (once per execution after commands are read)                    *
132 ********************************************************************************/
133void grid_init( void )
134{
135        double t, x, y, z;
136        double r;
137        double dtdr, dtdz;
138        long ix, iy, iz;
139        long ixyz;
140        long i;
141
142     /* Initialize the Stacking Grid! */
143        if(!initGrid) {
144                if ( maskGrid != 31 ) {
145                    logit( "e", "grid_init: Cannot initialize; "  );
146                    if( !(maskGrid&1)  ) logit( "e", "<grdlat> "  );
147                    if( !(maskGrid&2)  ) logit( "e", "<grdlon> "  );
148                    if( !(maskGrid&4)  ) logit( "e", "<grdz> "    );
149                    if( !(maskGrid&8)  ) logit( "e", "<dspace> "  );
150                    if( !(maskGrid&16) ) logit( "e", "<thresh> "  );
151                    logit( "e", "command(s) omitted; exiting!\n" );
152                    exit( -1 );
153                }
154                nThresh = mThresh;
155                grdFac.y = (float)(40000.0 / 360.0);
156                grdFac.x = grdFac.y
157                        * cos(0.5*(grdOrg.lat+grdMax.lat)*6.283185/360.0);
158                if(grdFac.x < (float)0.0)
159                        grdFac.x = -grdFac.x;
160                grdFac.z = (float)1.0;
161                initGrid = 1;
162                grid_cart(&grdMax, &grdTop);
163                nX = (int)(grdTop.x / dSpace + 0.5);
164                nY = (int)(grdTop.y / dSpace + 0.5);
165                nZ = (int)(grdTop.z / dSpace + 0.5);
166                if(nZ < 1)
167                        nZ = 1;
168                logit( "", "grid_init: entire stack grid dimensions: "
169                           "nx, ny, nz = %d %d %d\n", nX, nY, nZ); 
170                lpGridMax = nX * nY * nZ;
171                lpGrid = (short *)malloc( nX*nY*nZ*sizeof(short) );
172                for(i=0; i<nX*nY*nZ; i++) {
173                        lpGrid[i] = (short)0;
174                }
175                nXS = nYS = (int) rStack / dSpace;
176                logit( "", "grid_init: single station stack grid dimensions: "
177                           "nXS, nYS = %d %d\n", nXS, nYS);
178                lpTravMax = nXS * nYS * nZ;
179                lpTrav = (double *)malloc( nXS*nYS*nZ*sizeof(double) );
180                for(ix=0; ix<nXS; ix++) {
181                        x = dSpace * (double)ix;
182                        for(iy=0; iy<nYS; iy++) {
183                                y = dSpace * (double)iy;
184                                r = sqrt(x*x + y*y);
185                                for(iz=0; iz<nZ; iz++) {
186                                        z = dSpace * ((double)iz + 0.5);
187                                        t = t_lay(r, z,&dtdr, &dtdz);
188                                        ixyz = iz + nZ*iy + nYS*nZ*ix;
189                                        if ((ixyz<0) || (ixyz>=lpTravMax)) Panic( 11 );
190                                        lpTrav[ixyz] = t;
191                                }
192                        }
193                }
194                logit( "", "grid_init: stack grid initialized\n"); 
195
196             /* Do sanity checks on other configurable values */
197                if( mStack > mPick )  {
198                    logit( "","grid_init: stack depth (%ld) too large; "
199                           "resetting it to pick fifo length (%ld)\n", 
200                            mStack, mPick ); 
201                    mStack = mPick;
202                }
203        }
204        x = 0.0;
205}
206
207/**********************************************************************
208 * grid_com()  Process all recognized commands.                       *
209 **********************************************************************/
210int grid_com( void )
211{
212        if( !InitInstWild ) {
213           if( GetInst( "INST_WILDCARD", &InstWild ) != 0 ) {
214               logit("e","grid_com: INST_WILDCARD not valid; does earthworm_global.d"
215                        " file exist in EW_PARAMS?; exiting!\n" );
216               exit( 0 );
217           }
218           InitInstWild = 1;
219        }
220
221        if(k_its("log_stack")) {
222           logStack = 1;
223           return 1;
224        }
225        if(k_its("sae")) {
226           sae = k_int();
227           return 1;
228        }
229        if(k_its("ignore_loc_code_in_same_station_decision")) {
230           ignore_loc_code_in_same_station_decision = k_int();
231           return 1;
232        }
233        if(k_its("grid_debug")) {
234           grid_debug = k_int();
235           return 1;
236        }
237        if(k_its("nearest_quake_dist")) {
238           nearest_quake_dist = k_val();        /* required distance threshold */
239           nearest_quake_time = k_val();        /* optional time setting */
240           if (k_err()) nearest_quake_time = NEAREST_QUAKE_TIME; /* set to default */
241           return 1;
242        }
243
244        if(k_its("grdlat")) {
245           grdOrg.lat = (float)k_val();
246           grdMax.lat = (float)k_val();
247           maskGrid |= 1;
248           return 1;
249        }
250
251        if(k_its("grdlon")) {
252           grdOrg.lon = (float)k_val();
253           grdMax.lon = (float)k_val();
254           maskGrid |= 2;
255           return 1;
256        }
257
258        if(k_its("grdz")) {
259           grdOrg.z = (float)k_val();
260           grdMax.z = (float)k_val();
261           maskGrid |= 4;
262           return 1;
263        }
264
265        if(k_its("dspace")) {
266           dSpace = (float)k_val();
267           maskGrid |= 8;
268           return 1;
269        }
270        if(k_its("stack_horizontals")) {
271           stack_horizontals = 1;
272           return 1;
273        }
274
275        if(k_its("thresh")) {
276           mThresh = k_int();
277           maskGrid |= 16;
278           return 1;
279        }
280
281        if(k_its("rstack")) {
282           rStack = (float)k_val();
283           return 1;
284        }
285
286        if(k_its("tstack")) {
287           tStack = k_val(); 
288           return 1;
289        }
290
291        if(k_its("focus")) {
292           nFocus = k_int();
293           return 1;
294        }
295
296        if(k_its("rmsgrid")) {
297           rmsGrid = k_val();
298           return 1;
299        }
300
301        if(k_its("grid_wt")) {
302           if(nWt < MAXWT) {
303              Wt[nWt].wt     = k_int() + '0';
304              Wt[nWt].inc    = k_int();
305              Wt[nWt].instid = InstWild;
306              nWt++;
307              return 1;
308           }
309           logit( "e", "grid_com: Too many <grid_wt> commands! "
310                       "MAXWT = %d; exiting!\n", MAXWT );
311           exit( 0 );
312        }
313
314        if(k_its("grid_wt_instid")) {
315           char *str;
316           if( nWt < MAXWT ) {
317              Wt[nWt].wt  = k_int() + '0';
318              Wt[nWt].inc = k_int();
319              str = k_str();
320              if( str ) {
321                 if( GetInst( str, &Wt[nWt].instid ) != 0 ) {
322                    logit("e","grid_com: invalid installation id: %s in"
323                              "<grid_wt_instid> command; exiting!\n", str );
324                    exit( 0 );
325                 }
326              } else {
327                 logit("e","grid_com: no installation id in "
328                          "<grid_wt_instid> command; exiting!\n" );
329                 exit( 0 );
330              } 
331              logit("e","grid_wt_instid: %c %d %s:%u\n",
332                     Wt[nWt].wt, Wt[nWt].inc, str, Wt[nWt].instid );
333              nWt++;
334              return 1;
335           }
336           logit( "e", "grid_com: Too many <grid_wt_instid> commands! "
337                       "MAXWT = %d; exiting!\n", MAXWT );
338           exit( 0 );
339        }
340
341        if(k_its("stack")) {
342           mStack = k_long();
343           return 1;
344        }
345
346        if( k_its("define_glitch") ) {   /*new command; 960408:ldd*/
347           GlitchMinPk = k_int();
348           GlitchNsec  = k_val();
349           return 1;
350        }
351
352        return 0;
353}
354
355/********************************************************************************
356 * grid_cart(geo, cart)  Convert point from geocentric to cartesian coordinates *
357 ********************************************************************************/
358void grid_cart(GEO *geo, CART *cart)
359{
360        if(!initGrid)
361                grid_init();
362        cart->x = grdFac.x * (geo->lon - grdOrg.lon);
363        cart->y = grdFac.y * (geo->lat - grdOrg.lat);
364        cart->z = grdFac.z * (geo->z   - grdOrg.z);
365        cart->z = geo->z;
366}
367
368/********************************************************************************
369 * grid_geo(geo, cart) : Convert point from cartesian to geocentric coordinates *
370 ********************************************************************************/
371void grid_geo(GEO *geo, CART *cart)
372{
373        if(!initGrid)
374                grid_init();
375        geo->lon = cart->x / grdFac.x + grdOrg.lon;
376        geo->lat = cart->y / grdFac.y + grdOrg.lat;
377        geo->z   = cart->z;
378}
379
380/********************************************************************************
381 * grid_inc()    Return stacking increment based on pick quality and            *
382 *               installation id that produced the pick                         *
383 ********************************************************************************/
384int grid_inc(int wt, unsigned char instid)
385{
386        int inc;
387        int i;
388
389        inc = 0;
390        for(i=0; i<nWt; i++) {
391           if(Wt[i].wt == wt) {
392              if( Wt[i].instid == instid ) {   /* found rule for this installation */
393                 inc = Wt[i].inc;              /* set increment and return now!    */
394                 break;
395              }
396              if( Wt[i].instid == InstWild ) { /* wildcard: set increment but keep */
397                 inc = Wt[i].inc;              /* looking to see if there's a rule */
398              }                                /* for this specific installation!  */
399           }
400        }
401        return inc;
402}
403
404/********************************************************************************
405 * grid_dump()  Write a graphical representation of stack to log file           *
406 ********************************************************************************/
407void grid_dump(int ix, int iy, int iz)
408{
409        char txt[5000];
410        int ix1, ix2;
411        int iy1, iy2;
412        int jx, jy;
413        int is;
414        int i;
415
416        ix1 = 0;
417        ix2 = nX;
418        iy1 = 0;
419        iy2 = nY;
420        if(ix - nXS > ix1)      ix1 = ix - nXS;
421        if(ix + nXS < ix2)      ix2 = ix + nXS;
422        if(iy - nYS > iy1)      iy1 = iy - nYS;
423        if(iy + nYS < iy2)      iy2 = iy + nYS;
424        logit( "", "Grid counts, depth = %.1f km.\n", dSpace*(iz+0.5));
425        for(jy=iy2-1; jy>=iy1; jy--) {
426                i = 0;
427                for(jx=ix1; jx<ix2; jx++) {
428                        is = IXGRID(jx, jy, iz);
429                        if ((is<0) || (is>=lpGridMax)) Panic( 1 );
430                        if(lpGrid[is] > 9)
431                                txt[i] = 'a' + lpGrid[is] - 10;
432                        else
433                                txt[i] = '0' + lpGrid[is];
434                        if(jx == ix && jy == iy)
435                                txt[i] = '*';
436                        if(jx == 0) txt[i] = '|';
437                        if(jx == nX-1) txt[i] = '|';
438                        if(jy == 0) txt[i] = '-';
439                        if(jy == nY-1) txt[i] = '-';
440                        i++;
441                }
442                txt[i] = 0;
443                logit( "", "%s\n", txt);
444        }
445}
446
447/********************************************************************************
448 * grid_stack(lpick) : Stack a pick  (This is the guts!)                        *
449 ********************************************************************************/
450#define MAXLST 100
451int grid_stack(long lpick)
452{
453        CART    xyz;
454        GEO     geo;
455        double  tsum;
456        double  xsum;
457        double  ysum;
458        double  zsum;
459        int     ixpick;
460        int     site;           /* Stacking phase, site index           */
461        int     site_index;     /* Stackable pick, site index pointer   */
462        int     xd, yd;         /* Site dist. (temporary)               */
463        int     kdx, kdy;       /* Integer distance for focussing stack */
464        int     kdmin;          /* Current distance minimum             */
465        int     stkix1, stkix2; /* Stacking region for stacking pick    */
466        int     stkiy1, stkiy2;
467        int     ix1, ix2;       /* Stacking region, overlap of pTrav    */
468        int     iy1, iy2;
469        long    l1, l2;         /* Range of pick ids to stack           */
470        long    l, l_internal;  /* Iterating pick index                 */
471        int     ipick;          /* lpick % mPick                        */
472                                /* list changed to static by WMK        */
473        struct {                /* List of picks to stack, 0 is master  */
474                long  id;       /* Pick id                              */
475                double kt;      /* Offset in seconds                    */
476                short x;        /* Station location (grid coord)        */
477                short y;        /* Station location (grid coord)        */
478                char  use;      /* 1 if pk should be stacked; 0 if not  */ /*960408:ldd*/
479                char  horiz;    /* 1 horizontal component, 0 if vert */
480        } list[MAXLST];
481        SRT     srt[MAXLST];    /* sort structure for stacking list     */ /*960408:ldd*/
482        int     nlist;          /* Num. of picks in stacking list       */
483        double  tmax;           /* Time range of trav time template     */
484        double  ts;             /* Arrival time, stacking phase         */
485        double  t1, t2;         /* Stacking time overlap range          */
486        int     it;             /* Travel time iterator                 */
487        int     ix, iy, iz;     /* Spatial iterators                    */
488        double     kt;             /* Time accumator during stacking       */
489        int     i, i2;          /* Common iterator                      */
490        int     j;              /* Common iterator                      */
491        int     is;             /* Stack index                          */
492        int     inc;            /* Stack increment (wt dependent)       */
493        double     kt1, kt2, kt3;
494        int     idifx, idify;
495        int     reject;         /* used to reject a pick if higher weight one exists */
496        char    txt[80];
497
498#define MAXHIT 100
499        int mhit;               /* Maximum hits in stack                */
500        int nhit;               /* Number of hits in hit list           */
501        static struct {         /* hit changed to static by WMK         */
502                int x, y, z;
503        } hit[MAXHIT];          /* Hit list, stack indices              */
504
505        int npix;
506        static int ipix[MAXLST]; /* changed to static by WMK            */
507
508/* Initialization */
509        if(nWt < 1) {
510                nWt = 1;
511                Wt[0].wt = '0';
512                Wt[0].inc = 1;
513        }
514        if (sae != 0) {
515                if (time(NULL) > sae) {
516                    exit(2);
517                }
518        }
519
520/* Stacking phase, stack is constructed on the trav times of this phs   */
521        ipick = lpick % mPick;
522        site = pPick[ipick].site;
523        inc = grid_inc( pPick[ipick].wt, pPick[ipick].instid );
524        if ( inc < 1 ) return 0;
525        geo.lat = (float) Site[site].lat;
526        geo.lon = (float) Site[site].lon;
527        geo.z   = 0;
528        grid_cart(&geo, &xyz);
529        list[0].id  = lpick;
530        list[0].kt  = 0.0;
531        list[0].x   = (int)(xyz.x/dSpace);
532        list[0].y   = (int)(xyz.y/dSpace);
533        list[0].use = 1;                /*960408:ldd*/
534        list[0].horiz = 0;              /* starts out assuming vertical component */
535        srt[0].t     = pPick[ipick].t;  /*960408:ldd*/
536        srt[0].ilist = 0;               /*960408:ldd*/
537        if (stack_horizontals != 1 && is_component_horizontal(Site[site].comp[2],Site[site].net)) {
538                if (grid_debug)
539                        logit("t", "grid_stack: not using horizontal component to initiate a stack\n");
540                list[0].horiz = 1;
541                list[0].use = 0;
542                goto nada; /* changed in v1.0.5 to just not initiate on a horizontal */
543        }
544
545     /* Decide on range of picks to consider in stack */ 
546        l2 = lpick + mStack/2;
547        if(l2 > lPick)
548                l2 = lPick;
549        l1 = l2 - mStack;
550        if(l1 < 0)
551                l1 = 0;
552        it = IXTRAV(0, nYS-1, nZ-1);
553        if ((it<0) || (it>=lpTravMax))  Panic( 12 );
554        tmax = lpTrav[it];
555        ts = pPick[ipick].t;
556        t1 = ts - tmax;
557        t2 = ts + tmax;
558
559     /* Look thru picks in range; stash info for stackable picks */
560        nlist = 1;
561        for(l=l2-1; l>=l1; l--) {   /*go thru pick list backwards*/
562                i = l % mPick;
563                if(nlist >= MAXLST)             break;
564                if(l == lpick)                  continue;
565                if(pPick[i].quake)              continue;
566                if(pPick[i].t < t1)             continue;
567                if(pPick[i].t > t2)             continue;
568                site_index = pPick[i].site;
569                if (stack_horizontals != 1 && is_component_horizontal(Site[site_index].comp[2],Site[site_index].net)) {
570                     /* stacking was only designed for P arrivals, so don't use horizontals */
571                     if (grid_debug)
572                         logit("t", "grid_stack: not using pick on horizontal %s.%s.%s.%s from stashed pick info for stacking against\n",
573                                Site[site_index].name, Site[site_index].comp, Site[site_index].net, Site[site_index].loc);
574                     continue;
575                } 
576                /* the following are new checks to deal with multiple picks from the same station */
577                if (site == site_index) {
578                     /* don't add a pick from the same SCNL to the stackable list */
579                     if (grid_debug)
580                       logit("t", "grid_stack: not testing pick %s.%s.%s.%s wt(%c) for stacking against same SCNL %s.%s\n",
581                                Site[site_index].name, Site[site_index].comp, Site[site_index].net, Site[site_index].loc,
582                                pPick[i].wt, Site[site].net, Site[site].name);
583                     continue;
584                }
585                /* remember, sites are SCNL's so check S and N and L of pick against initiating phase,
586                   this is the case of multiple Z's on different components*/
587                if (is_same_SNL(site_index, site)) {
588                     if (grid_debug)
589                       logit("t", "grid_stack: not testing pick %s.%s.%s.%s wt(%c) for stacking against same station location %s.%s.%s\n",
590                                Site[site_index].name, Site[site_index].comp, Site[site_index].net, Site[site_index].loc,
591                                pPick[i].wt, Site[site].net, Site[site].name, Site[site].loc);
592                     continue;
593                }
594                /* if we reach here, we have a viable stackable pick (not a horizontal, comp not from same station as initiating pick */
595
596                /* now test the stackable pick against any others from this station, that are not horizontals! */
597                reject = 0;   /* set to 1 if this pick should be rejected because there is a better one from this sta in the list */
598                for (l_internal=l2-1; l_internal>=l1; l_internal--) {
599                     if (l_internal == l) continue; 
600                     j = l_internal % mPick;
601                     if(pPick[j].quake) continue;        /* already connected to another quake, don't stack against it */
602                     if(pPick[j].t < t1) continue;      /* pick outside time window */
603                     if(pPick[j].t > t2) continue;      /* pick outside time window */
604                     if (stack_horizontals != 1 && is_component_horizontal(Site[pPick[j].site].comp[2],Site[j].net)) continue;
605                     /* if it is a horizontal skip this check */
606                     /* okay, this is a viable pick for inclusion, but is it weighted higher then the current one being tested, "l"?? */
607                     /* test only if this is same SCNL but diff time pick <or> same SNL but diff time pick */
608                     if (pPick[j].site == site_index || is_same_SNL(pPick[j].site, site_index)) {
609                          /* remember wt here is 0,1,2,3 which is really quality, 0 better than 1 etc */
610                          /* is it case 1, higher weight? (note since wt=quality, < is used below) */
611                          if (pPick[j].wt < pPick[i].wt) reject = j;
612                          /* is it case 2, equal weight? */
613                          if (pPick[j].wt == pPick[i].wt && pPick[j].t > pPick[i].t) reject = j; /* choose earlier one */
614                          if (pPick[j].wt == pPick[i].wt && pPick[j].t == pPick[i].t) reject = j; /* just choose one they are equal */
615                          if (reject) break;
616                     }
617                }
618                if (reject) {
619                     /* weight was higher or equal to another one in the stack for same site */
620                     if (grid_debug)
621                       logit("t", "(DEBUG: grid_stack: not stacking against pick %s.%s.%s.%s with quality (%c) lower than a better pick (%c) %s.%s at same sta\n",
622                                Site[site_index].name, Site[site_index].comp, Site[site_index].net, Site[site_index].loc, 
623                                pPick[i].wt, pPick[reject].wt, Site[pPick[reject].site].comp, Site[pPick[reject].site].loc);
624                     continue;
625                }
626                geo.lat = (float) Site[site_index].lat;
627                geo.lon = (float) Site[site_index].lon;
628                geo.z   = 0;
629                grid_cart(&geo, &xyz);
630                list[nlist].id   = l;
631                list[nlist].kt   = pPick[i].t - ts;
632                list[nlist].x    = (int)(xyz.x/dSpace);
633                list[nlist].y    = (int)(xyz.y/dSpace);
634                list[nlist].use  = 1;                          /*960408:ldd*/
635                srt[nlist].t     = pPick[i].t;                 /*960408:ldd*/
636                srt[nlist].ilist = nlist;                      /*960408:ldd*/
637                xd = list[nlist].x - list[0].x;
638                if(xd < 0)      xd = -xd;
639                if(xd >= (2*nXS-1)) continue;
640                yd = list[nlist].y - list[0].y;
641                if(yd < 0)      yd = -yd;
642                if(yd >= (2*nYS-1)) continue;
643                nlist++;
644        }
645
646/* Look for glitches in list of picks to stack;
647   set glitch-pick use-flags to 0; block of code added 960408:ldd. */
648        if( GlitchMinPk != 0 )
649        {
650           qsort( srt, nlist, sizeof(SRT), grid_cmpsrt );
651           i  = 0;
652           i2 = GlitchMinPk - 1;
653           while( i2 < nlist ) {
654              if( (srt[i2].t-srt[i].t) <= GlitchNsec ) {  /* found a glitch! */
655                 for(j=i; j<=i2; j++) {           /* flag all picks in  */
656                    list[srt[j].ilist].use = 0;   /* glitch "don't use" */
657                 }
658              }
659              i++;
660              i2++;
661           }
662           if( list[0].use == 0 && list[0].horiz == 0) {    /* don't try to stack if master */
663              logit( "t",              /* pick belongs to a glitch     */
664                     "grid_stack, mhit = 0, glitch\n" );
665              goto nada;
666           } else if ( list[0].horiz == 1) {
667              logit( "t",              /* pick belongs to a horizontal     */
668                     "grid_stack, mhit = 0, horizontal\n" );
669              goto nada;
670           }
671        }
672 
673
674/* Sum the stack, scan for cell exceeding nThresh */
675        mhit = 0;
676        nhit = 0;
677
678/* Initialize grid elements to stacking-pick's increment value */
679/* The original loop over entire grid commented out; replaced by
680   loop over only the stack-pick's stacking grid.  950626:ldd. */
681/*      for(i=0; i<nX*nY*nZ; i++) */
682/*              lpGrid[i] = inc;  */
683        ix1 = -1;    /* fix: 20041110:LDD */
684        iy1 = -1;    /* fix: 20041110:LDD */
685        ix2 = nX;
686        iy2 = nY;
687        if(list[0].x - nXS > ix1)       ix1 = list[0].x - nXS;
688        if(list[0].x + nXS < ix2)       ix2 = list[0].x + nXS;
689        if(list[0].y - nYS > iy1)       iy1 = list[0].y - nYS;
690        if(list[0].y + nYS < iy2)       iy2 = list[0].y + nYS;
691        for(iz=0; iz<nZ; iz++) {
692            for(ix=ix1+1; ix<ix2; ix++) {         /* fix: 20041110:LDD */
693                for(iy=iy1+1; iy<iy2; iy++) {     /* fix: 20041110:LDD */
694                        is = IXGRID(ix, iy, iz);
695                        if ((is<0) || (is>=lpGridMax)) Panic( 2 );
696                        lpGrid[is] = inc;
697                }
698            }
699        }
700     /* Save boundaries of stack-pick's stacking grid */
701        stkix1 = ix1;  stkix2 = ix2;
702        stkiy1 = iy1;  stkiy2 = iy2;
703#ifdef DEBUG
704        { 
705        int ip = list[0].id % mPick;
706        int is = pPick[ip].site;
707        logit("","%s.%s.%s.%s inc: %d  x:%d ix: %d-%d  y:%d iy: %d-%d  iz: %d-%d\n",
708               Site[is].name,Site[is].comp,Site[is].net,Site[is].loc,
709               inc, list[0].x, ix1, ix2, list[0].y, iy1, iy2, 0, nZ ); */
710        }
711#endif /* DEBUG */
712/* end of loop over stack-pick's stacking grid.  950626:ldd. */
713
714/* loop thru other picks and increment bins in overlap of stacking-grids */
715        for(i=1; i<nlist; i++) {
716                int ip = list[i].id % mPick;
717#ifdef DEBUG
718                int s  = pPick[ip].site;
719#endif /* DEBUG */
720
721                if( list[i].use == 0 || list[i].horiz == 1 ) /* skip glitch-picks & horizontals */ /*960408:ldd*/
722                        continue;
723                inc = grid_inc(pPick[ip].wt, pPick[ip].instid);
724                if(inc < 1)
725                        continue;
726             /* find overlap between stack-pick's grid and this pick's grid */
727                ix1 = stkix1; 
728                iy1 = stkiy1;
729                ix2 = stkix2;
730                iy2 = stkiy2;
731                if(list[i].x - nXS > ix1)       ix1 = list[i].x - nXS; 
732                if(list[i].x + nXS < ix2)       ix2 = list[i].x + nXS;
733                if(list[i].y - nYS > iy1)       iy1 = list[i].y - nYS;
734                if(list[i].y + nYS < iy2)       iy2 = list[i].y + nYS;
735#ifdef DEBUG
736                { 
737                logit("","%s.%s.%s.%s inc: %d  x:%d ix: %d-%d  y:%d iy: %d-%d  iz: %d-%d\n",
738                      Site[s].name,Site[s].comp,Site[s].net,Site[s].loc,
739                      inc, list[i].x, ix1, ix2, list[i].y, iy1, iy2, 0, nZ );
740                } 
741#endif /* DEBUG */
742                for(iz=0; iz<nZ; iz++) {
743                    for(ix=ix1+1; ix<ix2; ix++) { 
744                       for(iy=iy1+1; iy<iy2; iy++) {
745                                kt1 = list[0].kt - list[i].kt;
746                                it = IXTRAV(list[0].x - ix,
747                                        list[0].y - iy, iz);
748                                if ((it<0) || (it>=lpTravMax))  Panic( 13 );
749                                kt2 = lpTrav[it];
750                                it = IXTRAV(list[i].x - ix,
751                                        list[i].y - iy, iz);
752                                if ((it<0) || (it>=lpTravMax))  Panic( 17 );
753                                kt3 = lpTrav[it];
754                                kt = kt1 - kt2 + kt3;
755                                is = IXGRID(ix, iy, iz);
756                                if ((is<0) || (is>=lpGridMax)) Panic( 3 );
757                                if(kt > -tStack && kt < tStack) {
758                                        lpGrid[is] += inc;
759#ifdef DEBUG
760                                        logit("","%s.%s.%s.%s hits lpGrid.%d.%d.%d %2d  resid:%.2f\n",
761                                              Site[s].name,Site[s].comp,                   
762                                              Site[s].net,Site[s].loc,                     
763                                              ix,iy,iz, lpGrid[is],(float)kt ); 
764#endif /* DEBUG */
765                                        if(lpGrid[is] < mhit)
766                                                continue;
767                                        if(lpGrid[is] > mhit) {
768                                                kdmin = INT_MAX;
769                                                mhit = lpGrid[is];
770                                                nhit = 0;
771                                        }
772                                        if(mhit < nThresh)
773                                                continue;
774                                        kdx = list[0].x - ix;
775                                        if(kdx < 0) kdx = - kdx;
776                                        if(kdx < kdmin) {
777                                                kdmin = kdx;
778                                                nhit = 0;
779                                        }
780                                        kdy = list[0].y - iy;
781                                        if(kdy < 0) kdy = - kdy;
782                                        if(kdy < kdmin) {
783                                                kdmin = kdy;
784                                                nhit = 0;
785                                        }
786                                        if(kdy > kdmin && kdx > kdmin)
787                                                continue;
788                                        if(nhit < MAXHIT) {
789                                                hit[nhit].x = ix;
790                                                hit[nhit].y = iy;
791                                                hit[nhit].z = iz;
792                                                nhit++;
793                                        }
794                                } 
795#ifdef DEBUG
796                                else { 
797                                   logit("","%s.%s.%s.%s miss lpGrid.%d.%d.%d %2d  resid:%.2f\n",
798                                      Site[s].name,Site[s].comp,                   
799                                      Site[s].net,Site[s].loc,                   
800                                      ix,iy,iz, lpGrid[is], kt );       
801                                } 
802#endif  /* END DEBUG loop */
803                        } /* end of loop over Y */
804                    } /* end of loop over X */
805                } /* end of loop over Z */
806        }
807        logit( "t", "grid_stack, mhit = %d\n", mhit);
808        if(mhit < nThresh)
809                goto nada;
810        if(nhit > nFocus) {
811                logit( "t", "grid_stack : null focus, nhit = %d\n", nhit);
812                goto nada;
813        }
814
815/* List phases contributing to stack maximum */
816        npix = 1;
817        ipix[0] = 0;    /* Including the master phase */
818        for(i=1; i<nlist; i++) {
819                if( list[i].use == 0 || list[i].horiz == 1) continue; /* skip glitch-picks or horizontal */ /*960408:ldd*/
820                for(j=0; j<nhit; j++) {
821                        ix = hit[j].x;
822                        iy = hit[j].y;
823                        iz = hit[j].z;
824
825                   /* Check distance from pick-cell to hit-cell; pick didn't  */
826                   /* contribute if it's beyond stacking limit. Added 2/13/96 */
827                        idifx = list[i].x - ix;
828                        idifx = (idifx < 0) ? -idifx : idifx;
829                        if ( idifx >= nXS ) continue;
830                        idify = list[i].y - iy;
831                        idify = (idify < 0) ? -idify : idify;
832                        if ( idify >= nYS ) continue;
833
834                   /* See if the travel-time residual falls in tStack range */
835                        kt1 = list[0].kt - list[i].kt;
836                        it = IXTRAV(list[0].x-ix, list[0].y-iy, iz);
837                        if ((it<0) || (it>=lpTravMax)) Panic( 14 );
838                        kt2 = lpTrav[it];
839                        it = IXTRAV(list[i].x-ix, list[i].y-iy, iz);
840                        if ((it<0) || (it>=lpTravMax))  Panic( 15 );
841                        kt3 = lpTrav[it];
842                        kt = kt1 - kt2 + kt3;
843                        if(kt > -tStack && kt < tStack) {
844                                ipix[npix++] = i; /* add the pick to the list that is used */
845                                break;
846                        }
847                }
848        }
849        if(npix < 4)
850                goto nada;
851
852/* Calculate initial location */
853        if(logStack)
854                grid_dump(list[0].x, list[0].y, hit[0].z);
855        tsum  = 0.0;
856        xsum  = 0.0;
857        ysum  = 0.0;
858        zsum  = 0.0;
859        for(i=0; i<nhit; i++) {
860                logit( "", "  hit[%d] = %d %d %d\n",
861                       i, hit[i].x, hit[i].y, hit[i].z);
862                it = IXTRAV(hit[i].x - list[0].x,
863                        hit[i].y - list[0].y, hit[i].z);
864                if ((it<0) || (it>=lpTravMax)) Panic( 16 );
865                tsum += ts - lpTrav[it];
866                xsum += hit[i].x;
867                ysum += hit[i].y;
868                zsum += hit[i].z;
869        }
870        i = lQuake % mQuake;
871        pQuake[i].t    = tsum/nhit;
872        pQuake[i].lat  = LAT(dSpace * (ysum/nhit + 0.5));
873        pQuake[i].lon  = LON(dSpace * (xsum/nhit + 0.5));
874        pQuake[i].z    = dSpace * (zsum/nhit + 0.5);
875        pQuake[i].rms  = 0.0;
876        pQuake[i].dmin = 0.0;
877        pQuake[i].ravg = 0.0;
878        pQuake[i].gap  = 0.0;
879        pQuake[i].npix = npix;
880        pQuake[i].nmod = 0;
881        pQuake[i].assessed   = 0;
882        pQuake[i].lpickfirst = 0;  /* filled in by bind_locate */
883        date20(pQuake[i].t, txt);
884        logit( "", "  starting location %s %.4f %.4f %.1f\n",
885                txt, pQuake[i].lat, pQuake[i].lon, pQuake[i].z);
886        logit( "", "grid_stack : npix = %d\n", npix);
887        if (is_quake_simultaneous(i)) {
888                pQuake[i].npix = 0; /* suggested fix by JinKooLee 5/29/2018 */
889                goto nada;      /* terminate nucleation, too close to other quake */
890        }
891        for(j=0; j<npix; j++) {
892                ixpick = ipix[j];
893                l = list[ixpick].id % mPick;
894                pPick[l].quake = lQuake;  /* assoc w/quake */
895                pPick[l].phase = 0;       /* label as P    */ /*bugfix 960408:ldd*/
896                bndr_link(lQuake, list[ixpick].id);
897        }
898        bndr_quake2k(lQuake);
899        lQuake++;
900        return 1;
901
902nada:
903        return 0;
904}
905
906/****************************************************************
907 *  grid_cmpsrt() compare 2 arrival times in a stacking list.   *
908 *                Function added 960408:ldd                     *
909 ****************************************************************/
910int grid_cmpsrt( const void *l1, const void *l2 )
911{
912        SRT *lst1;
913        SRT *lst2;
914
915        lst1 = (SRT *) l1;
916        lst2 = (SRT *) l2;
917        if(lst1->t < lst2->t)   return -1;
918        if(lst1->t > lst2->t)   return  1;
919        return 0;
920}
921
922/****************************************************************
923 * grid_free()  Free previously malloc'd memory                 *
924 ****************************************************************/
925void grid_free( void )
926{
927        free( (void *)lpGrid );
928        free( (void *)lpTrav );
929        return;
930}
931/****************************************************************/
932/* are the Site's Station, Network, and Location codes the same */
933/* Return 1 if true, 0 if false */
934/****************************************************************/
935int is_same_SNL(int i1, int i2) {
936
937int loc_is_the_same;
938         if (ignore_loc_code_in_same_station_decision) {
939           loc_is_the_same = 1; /* ignore the location code check */
940         } else {
941           loc_is_the_same = (strcmp(Site[i1].loc, Site[i2].loc) == 0); /* actually check loc code */
942         }   
943         if (strcmp(Site[i1].name, Site[i2].name) == 0 &&
944           strcmp(Site[i1].net, Site[i2].net) == 0 &&
945           loc_is_the_same) {
946           return 1;
947         }
948         return 0;
949}
950
951
952/* ****************
953  is_quake_simultaneous(i) - i is index of quake to test
954                returns n which is index of quake that meets test of nearest
955                or returns 0 - no nearest quake or test not activated
956***************** */
957int is_quake_simultaneous(int i) {
958int j;
959        if (nearest_quake_dist != 0.0) {
960                /* check to make sure no nearest neighbor quakes */
961                for(j=0;j<mQuake;j++) {
962                        /* first check if within time tolerance of other quakes */
963                        if (i != j && pQuake[j].npix !=0 && fabs(pQuake[j].t - pQuake[i].t)<nearest_quake_time) {
964                                /* then compute distance */
965                                double dist, azm;
966                                geo_to_km(pQuake[i].lat, pQuake[i].lon, pQuake[j].lat, pQuake[j].lon, &dist, &azm);
967                                /* debug */
968                                logit("t", "Debug: simultaneous event, time-diff %5.2fs and dist %6.2f km\n",
969                                                fabs(pQuake[j].t - pQuake[i].t), dist );
970                                if (dist < nearest_quake_dist) {
971                                        logit("t", "killing event, too close in time %5.2f < tolerance: %5.2fs and space to another by %6.2f km < tolerance %6.2f km\n",
972                                                fabs(pQuake[j].t - pQuake[i].t),nearest_quake_time, dist, nearest_quake_dist);
973                                         return(j);
974                                }
975                        }
976                }
977        }
978        return(0);
979}
Note: See TracBrowser for help on using the repository browser.