source: trunk/src/seismic_processing/glass/src/modules/Glock/Glock.cpp @ 3212

Revision 3212, 16.4 KB checked in by paulf, 12 years ago (diff)

sync from hydra_proj circa Oct 2007

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/*
2 *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
3 *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
4 *
5 *    $Id$
6 *
7 *    Revision history:
8 *     $Log$
9 *     Revision 1.3  2007/12/27 18:05:23  paulf
10 *     sync from hydra_proj circa Oct 2007
11 *
12 *     Revision 1.6  2006/09/30 05:34:51  davidk
13 *     All origin-pick relation calculations were moved to OPCalc(Gap, etc).  Now Glock
14 *     calls OPCalc to do the calculation of Origin-Pick derived parameters, and
15 *     calls Solve to invert the pick data into a better origin solution.
16 *     So now it acts as more of a delegatory shell, instead of containing a bunch
17 *     of calculations.
18 *
19 *     Revision 1.5  2006/05/04 23:16:25  davidk
20 *     Moved pMod and parameter initialization from static tasks to constructor tasks,
21 *     the way they should be.
22 *
23 *     Revision 1.4  2006/01/30 23:16:22  davidk
24 *     Added int type-cast to get rid of compiler warning.
25 *
26 *     Revision 1.3  2005/10/17 06:34:31  davidk
27 *     Modified the way phase weighting is handled in the locator.  The person who
28 *     originally put it in, did not understand how the least square solver worked and
29 *     wrongly applied the weighting values.
30 *     Now weighting is done in integral increments ranging from 0 to 4.
31 *     Currently any phase with weight below (GLOCK_WEIGHT_GRANULARITY/2)
32 *     is not passed to the solver.
33 *
34 *     Revision 1.2  2005/10/13 22:21:58  davidk
35 *     Fixed a bug where the locator was performing multiple relocations with the same
36 *     starting location, instead of refining the location by using the results of the
37 *     previous location as the starting point for the next.
38 *     Made an aesthetic change to the calculation of the depth coefficient.
39 *
40 *     Revision 1.1.1.1  2005/06/22 19:30:48  michelle
41 *     new directory tree built from files in HYDRA_NEWDIR_2005-06-20 tagged hydra and earthworm projects
42 *
43 *     Revision 1.8  2005/02/28 23:42:11  davidk
44 *     Downgraded a debug message.
45 *
46 *     Revision 1.7  2005/02/15 20:02:46  davidk
47 *     Modified the gap calculation (for the origin) in the locator, so that it is limited
48 *     to P-class (within 110 degrees distance) and PKPdf phases.
49 *     This way some flukily associated S or ancillary phase doesn't affect the
50 *     quality of the gap metric.
51 *
52 *     Revision 1.6  2004/12/02 17:51:20  davidk
53 *     Incorporated phase-based LocWeight values, that cause different phases
54 *     to contribute different amounts to affecting the location. Utilized LocWeight
55 *     in weight values submittted to the Solver, and in determining RMS, so
56 *     that picks that aren't weighted in the solution to not adversly skew
57 *     the RMS value even if they have a high residual value. (This is particularly
58 *     important with phases like P'P', that have huge residual values, but should
59 *     be associated with the origin, yet not contribute to it's location.
60 *
61 *     Revision 1.5  2004/11/30 02:48:52  davidk
62 *     Found bug where input params for the dtdz vertical slowness parameter
63 *     were being input in degrees (same unit as dtdx).  Converted the output
64 *     from degrees to km, before using it to adjust the origin depth.
65 *
66 *     Inverted the dtdz coefficient so that it was negative instead, the same way
67 *     the other coefficients were.
68 *
69 *     Revision 1.4  2004/11/01 06:33:21  davidk
70 *     Modified to work with new hydra traveltime library.
71 *     Cleaned up struct ITravelTime references.
72 *
73 *     Revision 1.3  2004/09/16 01:09:38  davidk
74 *     Added a configable variable for the number of Iterations that the Locator-Solver runs
75 *     each time an origin is modified.  Initialized it to the previously hardcoded value of 3-net-1.
76 *
77 *     Cleaned up logging messages.
78 *     Downgraded some debugging messages.
79 *
80 *     Revision 1.2  2004/04/01 22:09:18  davidk
81 *     v1.11 update
82 *     Now utilizes OPCalc routines for calculating Origin/Pick association params.
83 *
84 *     Revision 1.3  2003/11/07 22:42:29  davidk
85 *     Added RCS Header.
86 *     Added function Locate_InitOrigin() that initializes a blank
87 *     origin Structure for use with the locator.
88 *     Added CompareOrigins() function for comparing the "quality"
89 *     of two origin estimates for the same quake.
90 *     Used function in code that formerly compared origins based on
91 *     RMS values.
92 *
93 *     Revision 1.3  2003/11/07 22:38:59  davidk
94 *     Added RCS Header.
95 *     Added two member functions: Locate_InitOrigin() and CompareOrigins().
96 *
97 *
98 **********************************************************/
99
100// glock.cpp
101
102#include <windows.h>
103#include <stdio.h>
104#include <math.h>
105#include <Debug.h>
106#include <ISolve.h>
107#include <opcalc.h>
108#include "glock.h"
109#include "GlockMod.h"
110//#include "date.h"
111//#include "str.h"
112extern "C" {
113#include "utility.h"
114}
115
116
117extern "C" {
118void __cdecl qsortd (double *base, unsigned num);
119}
120
121
122#define ITER_SUCCESS    0
123#define ITER_NODATA             1
124#define ITER_DEPTH              2
125
126#define RAD2DEG  57.29577951308
127#define DEG2RAD 0.01745329251994
128
129// Carl 120303
130static double pPickDistArray[OPCALC_MAX_PICKS];
131static double pPickAzmArray[OPCALC_MAX_PICKS];
132
133/*
134static int PickAzmComp(const void *az1, const void *az2) {
135        double azm1;
136        double azm2;
137        memmove(&azm1, az1, sizeof(double));
138        memmove(&azm2, az2, sizeof(double));
139        if(azm1 < azm2)
140                return -1;
141        if(azm1 > azm2)
142                return 1;
143        return 0;
144}
145
146// End Carl 120303
147
148static int PickDistComp(const void *elem1, const void *elem2) {
149        PICK *pck1;
150        PICK *pck2;
151        memmove(&pck1, elem1, sizeof(PICK *));
152        memmove(&pck2, elem2, sizeof(PICK *));
153        CDebug::Log(DEBUG_MINOR_INFO,"Comp() %.2f %.2f\n", pck1->dDelta, pck2->dDelta);
154        if(pck1->dDelta < pck2->dDelta)
155                return -1;
156        if(pck1->dDelta > pck2->dDelta)
157                return 1;
158        return 0;
159}
160***************************************************/
161
162//---------------------------------------------------------------------------------------CWave
163CGlock::CGlock(CMod * IN_pMod) 
164{
165        for(int i=0; i<4; i++)
166                bFree[i] = true;
167
168        this->pMod = IN_pMod;
169        CDebug::SetModuleName("CGlock");
170        pGlnt = pMod->pGlint;
171        pTT = pMod->pTT;
172        OPCalc_SetTravelTimePointer(pTT);
173        pSlv = pMod->pSolve;
174
175}
176
177//---------------------------------------------------------------------------------------~CWave
178CGlock::~CGlock() {
179}
180
181//---------------------------------------------------------------------------------------Locate
182// Locate earthquake with constraints. E.G. set tnez to "TZ" to allow free
183// origin time and depth, tnez set to "TNE" is a fixed depth solution.
184// T = Time, N = North, E = East, Z = Depth
185int CGlock::Locate(char *ent, char *mask) {
186        int i;
187        int res;
188
189  // Carl 120303  - reversed each one of the true/false definitions
190        for(i=0; i<4; i++)
191                bFree[i] = false;
192        for(i=0; i<(int)strlen(mask); i++) {
193                switch(mask[i]) {
194                case 't':
195                case 'T':
196                        bFree[0] = true;
197                        break;
198                case 'n':
199                case 'N':
200                        bFree[1] = true;
201                        break;
202                case 'e':
203                case 'E':
204                        bFree[2] = true;
205                        break;
206                case 'z':
207                case 'Z':
208                        bFree[3] = true;
209                        break;
210                default:
211                        return 7001;
212                }
213        }
214        res = Locate(ent);
215        for(i=0; i<4; i++)
216                bFree[i] = true;
217        return res;
218}
219
220//---------------------------------------------------------------------------------------Locate
221// Locate earthquake (retrieve data from glint, update event in glint)
222int CGlock::Locate(char *ent) {
223        PICK *pck;
224        double torg;
225        double dis;
226        double azm;
227  double dDistKM;
228        int res;
229        bool fix;
230        bool btemp;
231
232        pOrg = pGlnt->getOrigin(ent);
233        if(!pOrg)
234                return 1;
235  if(pOrg->nPh < 4)  // we need atleast 4 phases for hypocenter processing
236    return(1);
237
238        torg = pOrg->dT;
239  nPck = 0;
240  int iPickRef=0;
241        while(pck = pGlnt->getPicksFromOrigin(pOrg, &iPickRef)) 
242  {
243    // DK 011404  record the min and max pick id's used in this
244    //            location for debugging purposes
245    if(nPck == 0)
246    {
247      strcpy(pOrg->idMinPick, pck->idPick);
248      strcpy(pOrg->idMaxPick, pck->idPick);
249    }
250    else
251    {
252      if(strcmp(pOrg->idMinPick,pck->idPick) > 0)
253        strcpy(pOrg->idMinPick, pck->idPick);
254      else if(strcmp(pOrg->idMaxPick,pck->idPick) < 0)
255        strcpy(pOrg->idMaxPick, pck->idPick);
256    }
257
258                geo_to_km_deg(pOrg->dLat, pOrg->dLon, pck->dLat, pck->dLon, &dDistKM, &dis, &azm);
259//              CDebug::Log(DEBUG_MINOR_INFO,"Locate(): PCK:%s t:%.2f lat:%.4f lon:%.4f elev:%.4f dis:%.2f azm:%.2f\n",
260//                      pck->idPick, pck->dT-torg, pck->dLat, pck->dLon, pck->dElev, dis, azm);
261                pPick[nPck++] = pck;
262        }
263
264  // DK Added 081303
265  ORIGIN oBest, oCurrent, oNext;
266  memcpy(&oCurrent,pOrg,sizeof(ORIGIN));
267  Locate_InitOrigin(&oNext, pOrg->idOrigin);
268  Locate_InitOrigin(&oBest, pOrg->idOrigin);
269
270        fix = false;
271  /* Note we must loop through two times before we start looping
272     through locator iterations.  Something about the first time
273     through calcs residuals for the existing origin, and the
274     second time loads the solver's matrix.  So the 3rd Iterate()
275     call is the first one that results in an actual relocation.
276   **************************************************************/
277        for(nIter=0; nIter<(pMod->iNumLocatorIterations+2); nIter++) {
278                if(fix) {
279                        btemp = bFree[3];
280                        bFree[3] = false;
281    }  // end if(fix)
282    res = Iterate(true, &oCurrent, &oNext);
283          CDebug::Log(DEBUG_MINOR_INFO,"Locate(): Trying location: %.4f %.4f %.4f - %.2f (%.f/%.1f-%3d)\n", 
284                oCurrent.dLat, oCurrent.dLon, oCurrent.dZ, 
285                oCurrent.dT, oCurrent.dAff, oCurrent.dRms, oCurrent.nEq);
286                if(fix) {
287                        fix = false;
288                        bFree[3] = btemp;
289    }  // end if(fix)
290    if(CompareOrigins(&oBest, &oCurrent) > 0)
291      memcpy(&oBest,&oCurrent, sizeof(ORIGIN));
292
293                switch(res) {
294                case ITER_SUCCESS:
295                        break;
296                case ITER_NODATA:
297                        pOrg->nEq = 0;
298                        pOrg->dRms = 0;
299                        return 0;
300                        break;
301                case ITER_DEPTH:
302                        fix = true;
303                        break;
304                }
305    // make the next origin the current one and reprocess
306    memcpy(&oCurrent,&oNext, sizeof(ORIGIN));
307  }  // end for 3 iterations
308  /* process the last origin */
309  res = Iterate(false, &oCurrent, &oNext);
310          CDebug::Log(DEBUG_MINOR_INFO,"Locate(): Trying location: %.4f %.4f %.4f - %.2f (%.1f-%3d)\n", 
311                oCurrent.dLat, oCurrent.dLon, oCurrent.dZ, 
312                oCurrent.dT, oCurrent.dRms, oCurrent.nEq);
313    if(CompareOrigins(&oBest, &oCurrent) > 0)
314      memcpy(&oBest,&oCurrent, sizeof(ORIGIN));
315
316  /* reset the originpick information to the best. */
317  res = Iterate(false, &oBest, &oNext);
318
319  if(oBest.nPh > 4)
320    memcpy(pOrg,&oBest, sizeof(ORIGIN));
321  else
322  {
323          CDebug::Log(DEBUG_MAJOR_INFO,"Locate(): Unable to generate location using "
324                                 "atleast 4 phases for Origin %s nPh:%d nEq:%d\n",
325                oBest.idOrigin, oBest.nPh, oBest.nEq);
326    return(1);
327  }
328
329        CDebug::Log(DEBUG_MINOR_INFO,"Locate(): Chose location : %.4f %.4f %.4f - %.2f (%.1f-%3d)\n", 
330              oBest.dLat, oBest.dLon, oBest.dZ, 
331              oBest.dT, oBest.dRms, oBest.nEq);
332        return 0;
333} // end Locate()
334
335//---------------------------------------------------------------------------------------Iterate
336int CGlock::Iterate(bool solve, ORIGIN * poCurrent, ORIGIN * poNext) 
337{
338        /*
339  PICK *pck;
340        TTT *ttt;
341        double tres;
342        double coef[4];
343        int j;
344        double tobs;
345        double tcal;
346        double d;
347        double azm;
348        double toa;
349        double z;
350        double sum;
351        double torg;
352        double lat;
353        double lon;
354        double dep;
355  */
356
357  int iRetCode;
358
359        int i,j;
360        double *stp;
361        int nvar;
362        int ivar;
363       
364        nvar = 0;
365        for(i=0; i<4; i++)
366                if(bFree[i]) nvar++;
367        if(nvar < 1)
368                solve = false;
369        if(solve)
370                pSlv->Init(nvar);
371
372  iRetCode = ProcessLocation(solve, poCurrent, pPick, nPck);
373
374        nIter++;
375        if(poCurrent->nEq < 4)
376                return ITER_SUCCESS;
377        if(!solve)
378                return ITER_SUCCESS;
379
380        // Iterate solution
381        stp = pSlv->Solve();
382        j = 0;
383
384  memcpy(poNext,poCurrent,sizeof(ORIGIN));
385  poNext->nEq=0;
386  poNext->dRms=10000.0;
387
388        for(ivar=0; ivar<4; ivar++) {
389                if(bFree[ivar]) {
390                        switch(ivar) {
391                        case 0:
392                                poNext->dT      += stp[j++];
393                                break;
394                        case 1:
395                                poNext->dLat    += stp[j++];
396                                break;
397                        case 2:
398                                poNext->dLon    += stp[j++];
399                                break;
400                        case 3:
401                                poNext->dZ      += stp[j++] * DEG2KM;
402                                break;
403                        }
404                }
405        }
406        if(poNext->dZ < 0.0)
407                return ITER_DEPTH;
408        if(poNext->dZ > 800.0)
409                return ITER_DEPTH;
410
411  // DK 061903
412  while(poNext->dLon > 180.0)
413  {
414    CDebug::Log(DEBUG_MAJOR_INFO,"Locate(): adjusting origin(%.2f,%6.2f,%7.2f,%3.0f) original(%.2f,%6.2f,%7.2f,%3.0f)\n",
415                poNext->dT,poNext->dLat,poNext->dLon,poNext->dZ, poCurrent->dT,poCurrent->dLat,poCurrent->dLon,poCurrent->dZ);
416    poNext->dLon-=360.0;
417  }
418  while(poNext->dLon < -180.0)
419  {
420    CDebug::Log(DEBUG_MAJOR_INFO,"Locate(): adjusting origin(%.2f,%6.2f,%7.2f,%3.0f) original(%.2f,%6.2f,%7.2f,%3.0f)\n",
421                poNext->dT,poNext->dLat,poNext->dLon,poNext->dZ, poCurrent->dT,poCurrent->dLat,poCurrent->dLon,poCurrent->dZ);
422    poNext->dLon+=360.0;
423  }
424  while(poNext->dLat > 90.0)
425  {
426    CDebug::Log(DEBUG_MAJOR_INFO,"Locate(): adjusting origin(%.2f,%6.2f,%7.2f,%3.0f) original(%.2f,%6.2f,%7.2f,%3.0f)\n",
427                poNext->dT,poNext->dLat,poNext->dLon,poNext->dZ, poCurrent->dT,poCurrent->dLat,poCurrent->dLon,poCurrent->dZ);
428    poNext->dLat=180.0 - poNext->dLat;
429  }
430  while(poNext->dLat < -90.0)
431  {
432    CDebug::Log(DEBUG_MAJOR_INFO,"Locate(): adjusting origin(%.2f,%6.2f,%7.2f,%3.0f) original(%.2f,%6.2f,%7.2f,%3.0f)\n",
433                poNext->dT,poNext->dLat,poNext->dLon,poNext->dZ, poCurrent->dT,poCurrent->dLat,poCurrent->dLon,poCurrent->dZ);
434    poNext->dLat=-180.0 - poNext->dLat;
435  }
436
437        // Update solution
438
439  CDebug::Log(DEBUG_MINOR_INFO,"Locate(%s,%d): Solve moved origin(%.2f,%6.2f,%7.2f,%3.0f) original(%.2f,%6.2f,%7.2f,%3.0f)\n",
440              poCurrent->idOrigin,poCurrent->iOrigin, poNext->dT,poNext->dLat,poNext->dLon,poNext->dZ, poCurrent->dT,poCurrent->dLat,poCurrent->dLon,poCurrent->dZ);
441        return ITER_SUCCESS;
442
443}  // end Iterate()
444
445//---------------------------------------------------------------------------------------Affinity
446// Calculates affinity factors for each associated phase. For now it is just the ratio
447// of a stations distance to N times the modal distance, where N is to be determined.
448
449#define GLOCK_WEIGHT_GRANULARITY  0.25
450
451int CGlock::AddPickCoefficientsForSolve(PICK * pPick, const ORIGIN * pOrg)
452{
453  int j = 0;
454  double coef[4];
455  int i;
456
457  if(bFree[0]) 
458    coef[j++] = 1.0;
459  if(bFree[1]) 
460    coef[j++] = -cos(DEG2RAD*pPick->dAzm)*pPick->ttt.dtdx;
461  if(bFree[2]) 
462    coef[j++] = -sin(DEG2RAD*pPick->dAzm)*cos(DEG2RAD*pOrg->dLat) * pPick->ttt.dtdx;
463  if(bFree[3]) 
464    coef[j++] = -pPick->ttt.dtdz; 
465
466  i = (int)((pPick->ttt.dLocWeight + (GLOCK_WEIGHT_GRANULARITY/2)) / GLOCK_WEIGHT_GRANULARITY);
467  if(pSlv)
468  {
469    for(j=0; j < i; j++)
470      pSlv->Equation(coef, pPick->tRes);
471    return(0);
472  }
473  else
474  {
475    return(-1);
476  }
477}  // end AddPickCoefficientsForSolve()
478
479// Function called by Iterate to calculate Origin and Origin/Link
480// params, and to suggest better locations.
481int CGlock::ProcessLocation(bool solve, ORIGIN * pOrg, PICK ** pPick, int nPck) 
482{
483
484  // Initialize origin vars
485        int rc = 0;
486  int i;
487
488
489  rc = OPCalc_ProcessOrigin(pOrg, pPick, nPck);
490  if(rc)
491    return(rc);
492
493
494
495//  CDebug::Log(DEBUG_MAJOR_INFO,"Locate(I): Org(%6.3f,%8.3f,%5.1f - %.2f)\n",
496//              pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT);
497
498        for(i=0; i<nPck; i++) 
499  {
500    // calculate solve coefficients if desired
501    if(solve)
502      AddPickCoefficientsForSolve(pPick[i], pOrg);
503   
504  }  // end for each pick in the Pick Array
505
506  return(0);
507}  // End Process Location()
508
509
510// Function to initialize the comparison parameters for
511// relative solution quality.
512void  CGlock::Locate_InitOrigin(ORIGIN * pOrigin, const char * idOrigin)
513{
514  // set the comparison criteria to something really high/bad
515  memset(pOrigin,0,sizeof(ORIGIN));
516  pOrigin->dRms = 10000.0;
517  strncpy(pOrigin->idOrigin, idOrigin, sizeof(pOrigin->idOrigin));
518  pOrigin->idOrigin[sizeof(pOrigin->idOrigin)-1] = 0x00;
519}
520
521
522// Standard compare function.  returns -1 if Origin1 is "better"
523// than Origin2, 1 if vice-versa, and 0 if they are of equal
524// quality.
525int CGlock::CompareOrigins(ORIGIN * po1, ORIGIN * po2)
526{
527  CDebug::Log(DEBUG_MINOR_INFO,
528              "Comparing Origins:\n"
529              "1) %.2f/%.2f/%.2f - %.2f  Aff: %.2f  AfGp: %.2f AfNeq: %.2f  RMS: %.2f Gap: %.2f\n"
530              "2) %.2f/%.2f/%.2f - %.2f  Aff: %.2f  AfGp: %.2f AfNeq: %.2f  RMS: %.2f Gap: %.2f\n",
531              po1->dLat, po1->dLon, po1->dZ, po1->dT, po1->dAff, po1->dAffGap, po1->dAffNumArr, po1->dRms, po1->dGap,
532              po2->dLat, po2->dLon, po2->dZ, po2->dT, po2->dAff, po2->dAffGap, po2->dAffNumArr, po2->dRms, po2->dGap
533             );
534
535  if(po1->dAff > po2->dAff)
536    return(-1);
537  if(po1->dAff < po2->dAff)
538    return(1);
539  //else
540  return(0);
541}  // end CompareOrigins()
542
543
Note: See TracBrowser for help on using the repository browser.