source: trunk/src/seismic_processing/glass/src/modules/OPCalc/opcalc.cpp @ 2176

Revision 2176, 7.0 KB checked in by paulf, 13 years ago (diff)

added from hydra_proj, new version

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#include <math.h>
2#include <Debug.h>
3#include <string.h>
4
5
6# define _OPCALC_SO_EXPORT  __declspec( dllexport)
7
8#include <opcalc.h>
9#include <opcalc_const.h>
10
11// Lib Global variables - File Scope
12static ITravelTime * pTT = 0x00;      // file-scope master traveltime pointer
13
14//---------------------------------------------------------------------------------------Zerg
15// Calculate the OPCalc_BellCurve function of x.
16// The height of the curve ranges from 0 to 1.
17// The curve is centered at 0, with a range from -1 to 1,
18// such that OPCalc_BellCurve(0) = 1 and OPCalc_BellCurve(-1) = OPCalc_BellCurve(1) = 0.
19// OPCalc_BellCurve is usually called as OPCalc_BellCurve(x/w),
20// where x is the deviation from the mean, and w is the half-width
21// of the distrubution desired.
22// OPCalc_BellCurve(-w) = OPCalc_BellCurve(w) = 0.0.
23// dBellCurve/dX(0) = dBellCurve/dX(1) = 0.0.
24// OPCalc_BellCurve is symetric about 0.0, and is roughly gaussian with a
25// half-width of 1. For x < -1 and x > 1, OPCalc_BellCurve is defined as 0.
26// In Carl terms, OPCalc_BellCurve is defined as Zerg.
27double OPCalc_BellCurve(double xx) 
28{
29        double x = xx;
30        if(x < 0.0)
31                x = -x;
32        if(x > 1.0)
33                return 0.0;
34        return 2.0*x*x*x - 3.0*x*x + 1.0;
35}
36
37
38int OPCalc_CalcOriginAffinity(ORIGIN * pOrg, PICK ** ppPick, int nPck) 
39{
40//      PICK *pck;
41        int i;
42  int iRetCode;
43  double dAffinitySum = 0.0;  // Sum of usable-pick affinities
44  int    iAffinityNum = 0;    // Number of usable-picks
45
46  if(nPck < 1)
47    return(-1);
48
49  // Log the calculated values
50        CDebug::Log(DEBUG_MINOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n",
51                pOrg->dDelMin, pOrg->dDelMod, pOrg->dDelMax);
52
53        // Gap affinity
54                // DK CHANGE 011504.  Changing the range of the
55    //                    dAffGap statistic
56    //         pOrg->dAffGap = 1.0+OPCalc_BellCurve(pOrg->dGap/360.0);
57    // dAffGap range is 0.0 - 2.0
58        if(pOrg->nEq < OPCALC_nMinStatPickCount)
59                pOrg->dAffGap = 1.0;
60  else
61    pOrg->dAffGap = 4 * OPCalc_BellCurve(pOrg->dGap/360.0); 
62//    pOrg->dAffGap = 3 * OPCalc_BellCurve(pOrg->dGap/360.0);
63
64  if(pOrg->dAffGap > 2.0)
65    pOrg->dAffGap = 2.0;
66  else if(pOrg->dAffGap < 0.5 && pOrg->dAffGap > 0.1)
67    pOrg->dAffGap = 0.5;  // don't let AffGap drop below 0.5 unless
68                          //the Gap is AWFUL!!!!
69
70
71
72        // nEq affinity (macho statistic)
73        if(pOrg->nEq < OPCALC_nMinStatPickCount)
74                pOrg->dAffNumArr = 1.0;
75        else
76                pOrg->dAffNumArr = log10((float)(pOrg->nEq));
77
78        for(i=0; i<nPck; i++) 
79  {
80    iRetCode = OPCalc_CalcPickAffinity(pOrg, ppPick[i]);
81          CDebug::Log(DEBUG_MINOR_INFO,"Affinity(): Pck[%d] Delta:%.2f Aff:%.2f\n", i, ppPick[i]->dDelta, ppPick[i]->dAff);
82
83    if(iRetCode < 0)
84    {
85          CDebug::Log(DEBUG_MINOR_ERROR,
86                  "OPCalc_CalcOriginAffinity(): Fatal Error(%d) calculating"
87                  " affinity for pick[%d]\n",
88                  iRetCode,i);
89      return(-1);
90    }
91    if(iRetCode > 0)
92      continue;
93
94    dAffinitySum+=ppPick[i]->dAff;
95    iAffinityNum++;
96  }
97
98  if(iAffinityNum > 0)
99  {
100    pOrg->dAff = dAffinitySum / iAffinityNum;
101  }
102  else
103  {
104    pOrg->dAff = 0;
105  }
106  return(0);
107}  // end OPCalc_CalcOriginAffinity()
108
109
110int OPCalc_CalcPickAffinity(ORIGIN * pOrg, PICK * pPick) 
111{
112
113  double dDistRatio;  // ratio between distance of current pick and median distance
114  double dMedianDistEst;
115
116  int iRetCode;
117
118        // Do full calculations for picks not already associated with the
119  // origin.  (distance, azimuth, traveltime, TakeOffAngle)
120        if(pPick->iState == GLINT_STATE_WAIF) 
121  {
122    iRetCode = OPCalc_CalculateOriginPickParams(pPick,pOrg);
123    if(iRetCode)
124      return(iRetCode);
125        }
126
127        // Residual affinity
128  // Range   0.0 - 2.0  where 0.0 is no residual and 2.0 is
129  //         any absolute residual >= 10 seconds (OPCALC_dResidualCutOff)
130        pPick->dAffRes = 2.0 * OPCalc_BellCurve(pPick->tRes/(pPick->ttt.dResidWidth));
131
132        // Distance affinity
133  if(OPCALC_fAffinityStatistics & AFFINITY_USE_DISTANCE)
134  {
135
136    // hack in a calculation that increases the value used for the median distance if the number of
137    // phases is large.
138    if(pOrg->nPh > (pOrg->dDelMod*2.0))
139      dMedianDistEst = (double) (pOrg->nPh/2);
140    else
141      dMedianDistEst = pOrg->dDelMod;
142   
143    // Calculate the distance of this pick relative
144    // to the distance of median pick used in the solution
145    if(pOrg->nEq < OPCALC_nMinStatPickCount) 
146    {
147        dDistRatio = pPick->dDelta / (OPCALC_dCutoffMultiple*dMedianDistEst * 1.5);
148    } 
149    else 
150    {
151        dDistRatio = pPick->dDelta / (OPCALC_dCutoffMultiple*dMedianDistEst);
152    }
153    // Distance Affinity
154    //       Penalize picks (based upon distance) when
155    //       they are more than dCutoffMultiple times the
156    //       median pick distance from the Origin.
157    // Range   0.0 - 2.0  where 0.0 is beyond the Modal distance cutoff and 2.0 is
158    //         a pick at the hypocenter.
159    pPick->dAffDis = 2 * OPCalc_BellCurve(dDistRatio);
160
161    if(pOrg->nEq < 4)
162      pPick->dAffDis = 1.0;
163  }
164  else
165  {
166    pPick->dAffDis = 1.0;
167  }
168
169  if(OPCALC_fAffinityStatistics & AFFINITY_USE_PPD)
170  {
171    // Pick Probability Affinity
172    // Calculate a Pick Probability Affinity statistic.
173    // Range   0.5 - 5.0 (highest observed value is 3.0)
174    // where the value is ((log 2 (X+1)) + 1) / 2
175    // where X is the number of times that a phase was picked
176    // within the Pick's travel time entry per 10000 picks.
177    pPick->dAffPPD = pPick->ttt.dAssocStrength / 2.0;
178    if(pPick->dAffPPD > 2.0)
179      pPick->dAffPPD = 2.0;
180  }
181  else
182  {
183    pPick->dAffPPD = 1.0;
184  }
185
186        // Composite affinity statistic
187        pPick->dAff = pOrg->dAffGap * pOrg->dAffNumArr * pPick->dAffRes * pPick->dAffDis * pPick->dAffPPD;
188
189  return(0);
190}  // end CalcPickAffinity()
191
192
193void OPCalc_SetTravelTimePointer(ITravelTime * pTravTime)
194{
195  pTT = pTravTime;
196}
197
198
199int OPCalc_CalculateOriginPickParams(PICK * pPick, ORIGIN * pOrg)
200{
201  TTEntry   ttt;
202  TTEntry * pttt;
203  double    dDistKM;
204
205  if(!(pTT && pPick && pOrg))
206  {
207                CDebug::Log(DEBUG_MINOR_ERROR,"OPCalc_CalculateOriginPickParams(): Null Pointer: "
208                                  "pPick(%u) pOrg(%u) pTT(%u)\n", 
209                pPick, pOrg, pTT);
210    return(-1);
211  }
212
213  if(pPick->dLat < -90.0 || pPick->dLat > 90.0 || pPick->dLon < -180.0 || pPick->dLon > 180.0)
214    return(-1);
215
216        geo_to_km_deg(pOrg->dLat, pOrg->dLon, pPick->dLat, pPick->dLon, 
217                &dDistKM, &pPick->dDelta, &pPick->dAzm);
218        pttt = pTT->TBest(pOrg->dZ, pPick->dDelta, pPick->dT - pOrg->dT, &ttt);
219        if(!pttt) 
220  {
221    pPick->bTTTIsValid = false;
222    pPick->tRes = 9999.9;
223    pPick->ttt.szPhase = Phases[PHASE_Unknown].szName;
224
225    return(1);
226  }
227
228        pPick->dTrav = pttt->dTPhase;
229        pPick->dToa = pttt->dTOA;
230        pPick->tRes = (pPick->dT - pOrg->dT) - pttt->dTPhase;
231  strcpy(pPick->sPhase, pttt->szPhase);
232  memcpy(&pPick->ttt, pttt, sizeof(TTEntry));
233  pPick->bTTTIsValid = true;
234  CDebug::Log(DEBUG_MINOR_INFO,"Iterate(): Eq[%d] d:%.2f z:%.2f tobs:%.2f azm:%.2f ttt:%d\n", 
235              pOrg->nEq, pPick->dDelta, pOrg->dZ, pPick->dT - pOrg->dT, pPick->dAzm, pttt);
236
237  return(0);
238}  // end OPCalc_CalculateOriginPickParams()
239
240
Note: See TracBrowser for help on using the repository browser.