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 | extern "C" { |
---|
12 | void __cdecl qsortd (double *base, unsigned num); |
---|
13 | } |
---|
14 | |
---|
15 | |
---|
16 | int OPCalc_ProcessOriginSimple(ORIGIN * pOrg, PICK ** pPick, int nPck); |
---|
17 | |
---|
18 | |
---|
19 | // Lib Global variables - File Scope |
---|
20 | static ITravelTime * pTT = 0x00; // file-scope master traveltime pointer |
---|
21 | |
---|
22 | //---------------------------------------------------------------------------------------Zerg |
---|
23 | // Calculate the OPCalc_BellCurve function of x. |
---|
24 | // The height of the curve ranges from 0 to 1. |
---|
25 | // The curve is centered at 0, with a range from -1 to 1, |
---|
26 | // such that OPCalc_BellCurve(0) = 1 and OPCalc_BellCurve(-1) = OPCalc_BellCurve(1) = 0. |
---|
27 | // OPCalc_BellCurve is usually called as OPCalc_BellCurve(x/w), |
---|
28 | // where x is the deviation from the mean, and w is the half-width |
---|
29 | // of the distrubution desired. |
---|
30 | // OPCalc_BellCurve(-w) = OPCalc_BellCurve(w) = 0.0. |
---|
31 | // dBellCurve/dX(0) = dBellCurve/dX(1) = 0.0. |
---|
32 | // OPCalc_BellCurve is symetric about 0.0, and is roughly gaussian with a |
---|
33 | // half-width of 1. For x < -1 and x > 1, OPCalc_BellCurve is defined as 0. |
---|
34 | // In Carl terms, OPCalc_BellCurve is defined as Zerg. |
---|
35 | double OPCalc_BellCurve(double xx) |
---|
36 | { |
---|
37 | double x = xx; |
---|
38 | if(x < 0.0) |
---|
39 | x = -x; |
---|
40 | if(x > 1.0) |
---|
41 | return 0.0; |
---|
42 | return 2.0*x*x*x - 3.0*x*x + 1.0; |
---|
43 | } |
---|
44 | |
---|
45 | |
---|
46 | int __cdecl ComparePickAffinityAsc(const void * pvoid1, const void * pvoid2) |
---|
47 | { |
---|
48 | PICK * p1, *p2; |
---|
49 | |
---|
50 | p1 = *(PICK **) pvoid1; |
---|
51 | p2 = *(PICK **) pvoid2; |
---|
52 | |
---|
53 | if(p1->dAff < p2->dAff) |
---|
54 | return(-1); |
---|
55 | if(p1->dAff > p2->dAff) |
---|
56 | return(1); |
---|
57 | else |
---|
58 | return(0); |
---|
59 | } |
---|
60 | |
---|
61 | int __cdecl ComparePickResidDesc(const void * pvoid1, const void * pvoid2) |
---|
62 | { |
---|
63 | PICK * p1, *p2; |
---|
64 | |
---|
65 | p1 = *(PICK **) pvoid1; |
---|
66 | p2 = *(PICK **) pvoid2; |
---|
67 | |
---|
68 | if(fabs(p1->tRes) > fabs(p2->tRes)) |
---|
69 | return(-1); |
---|
70 | if(fabs(p1->tRes) < fabs(p2->tRes)) |
---|
71 | return(1); |
---|
72 | else |
---|
73 | return(0); |
---|
74 | } |
---|
75 | |
---|
76 | int __cdecl ComparePickiPickAsc(const void * pvoid1, const void * pvoid2) |
---|
77 | { |
---|
78 | PICK * p1, *p2; |
---|
79 | |
---|
80 | p1 = *(PICK **) pvoid1; |
---|
81 | p2 = *(PICK **) pvoid2; |
---|
82 | |
---|
83 | if(p1->iPick < p2->iPick) |
---|
84 | return(-1); |
---|
85 | if(p1->iPick > p2->iPick) |
---|
86 | return(1); |
---|
87 | else |
---|
88 | return(0); |
---|
89 | } |
---|
90 | |
---|
91 | |
---|
92 | int OPCalc_ProcessOrigin(ORIGIN * pOrg, PICK ** pPick, int nPck) |
---|
93 | { |
---|
94 | |
---|
95 | double dFullAff; |
---|
96 | int rc; |
---|
97 | int i; |
---|
98 | PICK * pck; |
---|
99 | |
---|
100 | /* process the origin with the full complement of picks |
---|
101 | *******************************************************/ |
---|
102 | if(rc=OPCalc_ProcessOriginSimple(pOrg, pPick, nPck)) |
---|
103 | return(rc); |
---|
104 | |
---|
105 | dFullAff = pOrg->dAff; |
---|
106 | |
---|
107 | // log the results |
---|
108 | CDebug::Log(DEBUG_MINOR_INFO,"OPC_PO: Full Origin: %.2f/%.2f/%.2f - %.2f Aff: %.2f AfGp: %.2f AfNeq: %.2f RMS: %.2f Gap: %.2f\n", |
---|
109 | pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT, pOrg->dAff, pOrg->dAffGap, pOrg->dAffNumArr, pOrg->dRms, pOrg->dGap, |
---|
110 | pOrg->dRms, pOrg->nEq); |
---|
111 | |
---|
112 | for(i=0; i<nPck; i++) |
---|
113 | { |
---|
114 | pck = pPick[i]; |
---|
115 | CDebug::Log(DEBUG_MINOR_INFO, |
---|
116 | "%3d) %-4s %-6s %5.1f %3.0f %3.0f %5.2f %6.1f %4.1f %3.1f/%3.1f/%3.1f %s\n", |
---|
117 | i, pck->sSite, pck->sPhase, |
---|
118 | pck->dDelta, pck->dAzm, pck->dToa * RAD2DEG, pck->tRes, pck->dT - pOrg->dT, |
---|
119 | pck->dAff, pck->dAffDis, pck->dAffRes, pck->dAffPPD, pck->sTag |
---|
120 | ); |
---|
121 | } |
---|
122 | |
---|
123 | // sort the picks by resid |
---|
124 | qsort(pPick, nPck, sizeof(PICK *), ComparePickResidDesc); |
---|
125 | |
---|
126 | /* when trimming picks, don't trim anything that is either: |
---|
127 | 1) outside of the worst 20% in terms of residual |
---|
128 | 2) inside 1 std. deviation (assuming normal distribution with stddev = 3.33 sec) |
---|
129 | We Make three BIG assumptions here: |
---|
130 | a) That the ORIGIN is an accurate approximation of the Event origin |
---|
131 | b) That the pick residual distribution is NORMAL. |
---|
132 | c) That the stddev for such a distribution is 3.33 seconds. |
---|
133 | **************************************************************/ |
---|
134 | #define OPCALC_MAX_RESID_TRIM_RATIO .20 |
---|
135 | #define OPCALC_PICK_RESID_DIST_STD_DEV_SEC 3.33 |
---|
136 | |
---|
137 | // trim out the worst picks (by res) |
---|
138 | for(i=0; i<(nPck*OPCALC_MAX_RESID_TRIM_RATIO); i++) |
---|
139 | { |
---|
140 | pck = pPick[i]; |
---|
141 | CDebug::Log(DEBUG_MINOR_INFO, "Pick filter: %d res: %.2f aff: %.2f (%.2f ~ %.2f)\n", |
---|
142 | i, pck->tRes, pck->dAff, |
---|
143 | fabs(pck->tRes), OPCALC_PICK_RESID_DIST_STD_DEV_SEC); |
---|
144 | if(fabs(pck->tRes) < OPCALC_PICK_RESID_DIST_STD_DEV_SEC) |
---|
145 | break; |
---|
146 | } |
---|
147 | |
---|
148 | /* process the origin with the full complement of picks |
---|
149 | *******************************************************/ |
---|
150 | if(i<nPck) |
---|
151 | { |
---|
152 | OPCalc_ProcessOriginSimple(pOrg, &pPick[i], nPck-i); |
---|
153 | |
---|
154 | // log the results |
---|
155 | CDebug::Log(DEBUG_MINOR_INFO,"OPC_PO: Trimmed Origin: %.2f/%.2f/%.2f - %.2f Aff: %.2f AfGp: %.2f AfNeq: %.2f RMS: %.2f Gap: %.2f\n", |
---|
156 | pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT, pOrg->dAff, pOrg->dAffGap, pOrg->dAffNumArr, pOrg->dRms, pOrg->dGap, |
---|
157 | pOrg->dRms, pOrg->nEq); |
---|
158 | |
---|
159 | for(; i<nPck; i++) |
---|
160 | { |
---|
161 | pck = pPick[i]; |
---|
162 | CDebug::Log(DEBUG_MINOR_INFO, |
---|
163 | "%3d) %-4s %-6s %5.1f %3.0f %3.0f %5.2f %6.1f %4.1f %3.1f/%3.1f/%3.1f %s\n", |
---|
164 | i, pck->sSite, pck->sPhase, |
---|
165 | pck->dDelta, pck->dAzm, pck->dToa * RAD2DEG, pck->tRes, pck->dT - pOrg->dT, |
---|
166 | pck->dAff, pck->dAffDis, pck->dAffRes, pck->dAffPPD, pck->sTag |
---|
167 | ); |
---|
168 | } |
---|
169 | |
---|
170 | // if the full was better, then reprocess the parameters using the FULL |
---|
171 | if(dFullAff > pOrg->dAff) |
---|
172 | { |
---|
173 | OPCalc_ProcessOriginSimple(pOrg, pPick, nPck); |
---|
174 | } |
---|
175 | } |
---|
176 | |
---|
177 | // sort the picks by resid |
---|
178 | qsort(pPick, nPck, sizeof(PICK *), ComparePickiPickAsc); |
---|
179 | |
---|
180 | return(0); |
---|
181 | |
---|
182 | } |
---|
183 | int OPCalc_ProcessOriginSimple(ORIGIN * pOrg, PICK ** pPick, int nPck) |
---|
184 | { |
---|
185 | |
---|
186 | double pPickDistArray[OPCALC_MAX_PICKS]; |
---|
187 | double pPickAzmArray[OPCALC_MAX_PICKS]; |
---|
188 | |
---|
189 | // Initialize origin vars |
---|
190 | int nEq = 0; |
---|
191 | int i; |
---|
192 | |
---|
193 | double sum = 0.0; |
---|
194 | double dWeight = 0.0; |
---|
195 | PICK * pck; |
---|
196 | int iRetCode; |
---|
197 | double dGap; |
---|
198 | |
---|
199 | PhaseType ptPhase; |
---|
200 | PhaseClass pcPhase; |
---|
201 | |
---|
202 | // Set default RMS to absurd level. Shouldn't this be an affinity setting. |
---|
203 | pOrg->dRms = 10000.0; |
---|
204 | |
---|
205 | |
---|
206 | // CDebug::Log(DEBUG_MAJOR_INFO,"Locate(I): Org(%6.3f,%8.3f,%5.1f - %.2f)\n", |
---|
207 | // pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT); |
---|
208 | |
---|
209 | for(i=0; i<nPck; i++) |
---|
210 | { |
---|
211 | pck = pPick[i]; |
---|
212 | iRetCode = OPCalc_CalculateOriginPickParams(pck,pOrg); |
---|
213 | if(iRetCode == -1) |
---|
214 | { |
---|
215 | CDebug::Log(DEBUG_MINOR_ERROR,"ProcessOrigin(): OPCalc_CalculateOriginPickParams failed with " |
---|
216 | "terminal error. See logfile. Returning!\n"); |
---|
217 | return(-1); |
---|
218 | } |
---|
219 | else if(iRetCode == 1) |
---|
220 | { |
---|
221 | // could not find TravelTimeTable entry for pOrg,pck |
---|
222 | CDebug::Log(DEBUG_MINOR_INFO,"ProcessOrigin(): OPCalc_CalculateOriginPickParams() could not " |
---|
223 | "find TTT entry for Org,Pick!\n"); |
---|
224 | continue; |
---|
225 | } |
---|
226 | |
---|
227 | // accumulate statistics |
---|
228 | sum += pck->tRes*pck->tRes*pck->ttt.dLocWeight; |
---|
229 | |
---|
230 | dWeight += pck->ttt.dLocWeight; |
---|
231 | |
---|
232 | // copy the delta/azimuth of the pick to proper arrays for |
---|
233 | // further calculations |
---|
234 | pPickDistArray[nEq] = pck->dDelta; |
---|
235 | |
---|
236 | // DK 021105 Change gap calculation to only use P-Class phases |
---|
237 | // that way if we pull in random S, PPP, and SKiKP phases, |
---|
238 | // they won't affect the gap "quality" of the solution. |
---|
239 | ptPhase = GetPhaseType(pck->ttt.szPhase); |
---|
240 | pcPhase = GetPhaseClass(ptPhase); |
---|
241 | if(( (pcPhase == PHASECLASS_P) && (pck->dDelta < OPCALC_MAX_BELIEVABLE_DISTANCE_FOR_P_PHASE )) || |
---|
242 | (ptPhase == PHASE_PKPdf) |
---|
243 | ) |
---|
244 | pPickAzmArray[nEq] = pck->dAzm; |
---|
245 | else |
---|
246 | pPickAzmArray[nEq] = -1.0; |
---|
247 | |
---|
248 | |
---|
249 | // increment the number of picks used to locate the origin |
---|
250 | nEq++; |
---|
251 | |
---|
252 | // keep the arrays from being overrun |
---|
253 | if(nEq > OPCALC_MAX_PICKS) |
---|
254 | break; |
---|
255 | } // end for each pick in the Pick Array |
---|
256 | |
---|
257 | // if we didn't use atleast one pick. Give up and go home. |
---|
258 | if(nEq < 1) |
---|
259 | return(1); |
---|
260 | |
---|
261 | |
---|
262 | // Record the number of equations used to calc the location |
---|
263 | pOrg->nEq = nEq; |
---|
264 | |
---|
265 | // Record the residual standard deviation |
---|
266 | pOrg->dRms = sqrt(sum/dWeight); |
---|
267 | |
---|
268 | /* Calculate Pick Azimuth statistics for the Origin (Gap) |
---|
269 | ****************************************************/ |
---|
270 | // Sort the pick-azimuth array |
---|
271 | qsortd(pPickAzmArray, nEq); |
---|
272 | |
---|
273 | // If the Pick wasn't a P-class Pick, then we gave it a negative azm, |
---|
274 | // so igonre all of the negative azm picks (which are all at the beginning |
---|
275 | // since the array has been sorted |
---|
276 | for(i=0; i<nEq-2; i++) |
---|
277 | { |
---|
278 | if(pPickAzmArray[i] >= 0.0) |
---|
279 | break; |
---|
280 | } |
---|
281 | |
---|
282 | CDebug::Log(DEBUG_MAJOR_INFO, |
---|
283 | "ProcessOrigin(): Calculating Gap for Org %d(%5.1f, %6.1f, %3.0f) based on %d valid picks.\n ", |
---|
284 | pOrg->iOrigin, pOrg->dLat, pOrg->dLon, pOrg->dZ, nEq - i); |
---|
285 | // Calculate the gap for the last-to-first wrap |
---|
286 | dGap = pPickAzmArray[i] + 360.0 - pPickAzmArray[nEq-1]; |
---|
287 | for(; i<nEq-1; i++) |
---|
288 | { |
---|
289 | // for each gap between two picks, record only if largest so far |
---|
290 | if(pPickAzmArray[i+1] - pPickAzmArray[i] > dGap) |
---|
291 | dGap = pPickAzmArray[i+1] - pPickAzmArray[i]; |
---|
292 | } |
---|
293 | pOrg->dGap = dGap; |
---|
294 | |
---|
295 | /* Calculate Pick Distance statistics for the Origin |
---|
296 | ****************************************************/ |
---|
297 | // Sort the pick-distance array |
---|
298 | qsortd(pPickDistArray, nEq); |
---|
299 | |
---|
300 | // Calculate the minimum station distance for the origin |
---|
301 | pOrg->dDelMin = pPickDistArray[0]; |
---|
302 | |
---|
303 | // Calculate the maximum station distance for the origin |
---|
304 | pOrg->dDelMax = pPickDistArray[nEq-1]; |
---|
305 | |
---|
306 | // Calculate the median station distance for the origin |
---|
307 | pOrg->dDelMod = 0.5*(pPickDistArray[nEq/2] + pPickDistArray[(nEq-1)/2]); |
---|
308 | |
---|
309 | |
---|
310 | /* Calculate the Affinity for the Origin (and Picks) |
---|
311 | ****************************************************/ |
---|
312 | iRetCode = OPCalc_CalcOriginAffinity(pOrg, pPick, nPck); |
---|
313 | if(iRetCode) |
---|
314 | { |
---|
315 | CDebug::Log(DEBUG_MINOR_ERROR, |
---|
316 | "Error: ProcessOrigin(): OPCalc_CalcOriginAffinity() return (%d)\n", |
---|
317 | iRetCode); |
---|
318 | return(iRetCode); |
---|
319 | } |
---|
320 | |
---|
321 | return(0); |
---|
322 | } // End Process Location() |
---|
323 | |
---|
324 | |
---|
325 | int OPCalc_CalcOriginAffinity(ORIGIN * pOrg, PICK ** ppPick, int nPck) |
---|
326 | { |
---|
327 | // PICK *pck; |
---|
328 | int i; |
---|
329 | int iRetCode; |
---|
330 | double dAffinitySum = 0.0; // Sum of usable-pick affinities |
---|
331 | int iAffinityNum = 0; // Number of usable-picks |
---|
332 | |
---|
333 | if(nPck < 1) |
---|
334 | return(-1); |
---|
335 | |
---|
336 | // Log the calculated values |
---|
337 | CDebug::Log(DEBUG_MAJOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n", |
---|
338 | pOrg->dDelMin, pOrg->dDelMod, pOrg->dDelMax); |
---|
339 | |
---|
340 | // Gap affinity |
---|
341 | // DK CHANGE 011504. Changing the range of the |
---|
342 | // dAffGap statistic |
---|
343 | // pOrg->dAffGap = 1.0+OPCalc_BellCurve(pOrg->dGap/360.0); |
---|
344 | // dAffGap range is 0.0 - 2.0 |
---|
345 | pOrg->dAffGap = 4 * OPCalc_BellCurve(pOrg->dGap/360.0); |
---|
346 | if(pOrg->nEq < OPCALC_nMinStatPickCount && pOrg->dAffGap < 1.0) |
---|
347 | pOrg->dAffGap = 1.0; |
---|
348 | |
---|
349 | // pOrg->dAffGap = 3 * OPCalc_BellCurve(pOrg->dGap/360.0); |
---|
350 | |
---|
351 | if(pOrg->dAffGap > 2.0) |
---|
352 | pOrg->dAffGap = 2.0; |
---|
353 | else if(pOrg->dAffGap < 0.5 && pOrg->dAffGap > 0.1) |
---|
354 | pOrg->dAffGap = 0.5; // don't let AffGap drop below 0.5 unless |
---|
355 | //the Gap is AWFUL!!!! |
---|
356 | |
---|
357 | |
---|
358 | |
---|
359 | // nEq affinity (macho statistic) |
---|
360 | if(pOrg->nEq < OPCALC_nMinStatPickCount) |
---|
361 | pOrg->dAffNumArr = 1.0; |
---|
362 | else |
---|
363 | pOrg->dAffNumArr = log10((double)(pOrg->nEq)); |
---|
364 | |
---|
365 | for(i=0; i<nPck; i++) |
---|
366 | { |
---|
367 | iRetCode = OPCalc_CalcPickAffinity(pOrg, ppPick[i]); |
---|
368 | CDebug::Log(DEBUG_MINOR_INFO,"Affinity(): Pck[%d] Delta:%.2f Aff:%.2f\n", i, ppPick[i]->dDelta, ppPick[i]->dAff); |
---|
369 | |
---|
370 | if(iRetCode < 0) |
---|
371 | { |
---|
372 | CDebug::Log(DEBUG_MINOR_ERROR, |
---|
373 | "OPCalc_CalcOriginAffinity(): Fatal Error(%d) calculating" |
---|
374 | " affinity for pick[%d]\n", |
---|
375 | iRetCode,i); |
---|
376 | return(-1); |
---|
377 | } |
---|
378 | if(iRetCode > 0) |
---|
379 | continue; |
---|
380 | |
---|
381 | dAffinitySum+=ppPick[i]->dAff; |
---|
382 | iAffinityNum++; |
---|
383 | } |
---|
384 | |
---|
385 | if(iAffinityNum > 0) |
---|
386 | { |
---|
387 | pOrg->dAff = dAffinitySum / iAffinityNum; |
---|
388 | |
---|
389 | /* now lets do something more advanced. |
---|
390 | We are trying to calculate the "quality" of this origin. |
---|
391 | We're not trying to recalculate where the origin should be, |
---|
392 | only trying to caculate a "quality" number. |
---|
393 | Under the current basic calculations, we can run into problems |
---|
394 | where a couple of random picks can have residuals close enough in |
---|
395 | to associate, but far enough out to lower the quality of the origin. |
---|
396 | A problem is, that as soon as you remove a pick from the origin |
---|
397 | calculation, you have to do a lot of quality recalculation. Removing |
---|
398 | one pick potentially affects both AffGap affinity component as well |
---|
399 | as the AffNumArr. |
---|
400 | What we'd like is to find some sort of peak quality value for the Origin, |
---|
401 | (a reverse saddle). We'll try to get there by taking the full result |
---|
402 | (all-picks), and remove the picks with the lowest affinity until we hit |
---|
403 | a point where the origin affinity starts to decrease. |
---|
404 | Unfortunately this can't all be done in the OPCalc routines, because the |
---|
405 | Gap is currently calculated in the Glock module ONLY. |
---|
406 | DK 0920'06 |
---|
407 | ************************************************************************/ |
---|
408 | |
---|
409 | for(i=0; i<nPck; i++) |
---|
410 | { |
---|
411 | } |
---|
412 | |
---|
413 | |
---|
414 | |
---|
415 | } |
---|
416 | else |
---|
417 | { |
---|
418 | pOrg->dAff = 0; |
---|
419 | } |
---|
420 | return(0); |
---|
421 | } // end OPCalc_CalcOriginAffinity() |
---|
422 | |
---|
423 | |
---|
424 | int OPCalc_CalcPickAffinity(ORIGIN * pOrg, PICK * pPick) |
---|
425 | { |
---|
426 | |
---|
427 | double dDistRatio; // ratio between distance of current pick and median distance |
---|
428 | double dMedianDistEst; |
---|
429 | |
---|
430 | int iRetCode; |
---|
431 | |
---|
432 | // Do full calculations for picks not already associated with the |
---|
433 | // origin. (distance, azimuth, traveltime, TakeOffAngle) |
---|
434 | if(pPick->iState == GLINT_STATE_WAIF) |
---|
435 | { |
---|
436 | iRetCode = OPCalc_CalculateOriginPickParams(pPick,pOrg); |
---|
437 | if(iRetCode) |
---|
438 | return(iRetCode); |
---|
439 | } |
---|
440 | |
---|
441 | // Residual affinity |
---|
442 | // Range 0.0 - 2.0 where 0.0 is no residual and 2.0 is |
---|
443 | // any absolute residual >= 10 seconds (OPCALC_dResidualCutOff) |
---|
444 | pPick->dAffRes = 2.0 * OPCalc_BellCurve(pPick->tRes/(pPick->ttt.dResidWidth)); |
---|
445 | |
---|
446 | // Distance affinity |
---|
447 | if(OPCALC_fAffinityStatistics & AFFINITY_USE_DISTANCE) |
---|
448 | { |
---|
449 | |
---|
450 | // hack in a calculation that increases the value used for the median distance if the number of |
---|
451 | // phases is large. |
---|
452 | if(pOrg->nPh > (pOrg->dDelMod*2.0)) |
---|
453 | dMedianDistEst = (double) (pOrg->nPh/2); |
---|
454 | else |
---|
455 | dMedianDistEst = pOrg->dDelMod; |
---|
456 | |
---|
457 | // Calculate the distance of this pick relative |
---|
458 | // to the distance of median pick used in the solution |
---|
459 | if(pOrg->nEq < OPCALC_nMinStatPickCount) |
---|
460 | { |
---|
461 | dDistRatio = pPick->dDelta / (OPCALC_dCutoffMultiple*dMedianDistEst * 1.5); |
---|
462 | } |
---|
463 | else |
---|
464 | { |
---|
465 | dDistRatio = pPick->dDelta / (OPCALC_dCutoffMultiple*dMedianDistEst); |
---|
466 | } |
---|
467 | // Distance Affinity |
---|
468 | // Penalize picks (based upon distance) when |
---|
469 | // they are more than dCutoffMultiple times the |
---|
470 | // median pick distance from the Origin. |
---|
471 | // Range 0.0 - 2.0 where 0.0 is beyond the Modal distance cutoff and 2.0 is |
---|
472 | // a pick at the hypocenter. |
---|
473 | pPick->dAffDis = 2 * OPCalc_BellCurve(dDistRatio); |
---|
474 | |
---|
475 | if(pOrg->nEq < 4) |
---|
476 | pPick->dAffDis = 1.0; |
---|
477 | } |
---|
478 | else |
---|
479 | { |
---|
480 | pPick->dAffDis = 1.0; |
---|
481 | } |
---|
482 | |
---|
483 | if(OPCALC_fAffinityStatistics & AFFINITY_USE_PPD) |
---|
484 | { |
---|
485 | // Pick Probability Affinity |
---|
486 | // Calculate a Pick Probability Affinity statistic. |
---|
487 | // Range 0.5 - 5.0 (highest observed value is 3.0) |
---|
488 | // where the value is ((log 2 (X+1)) + 1) / 2 |
---|
489 | // where X is the number of times that a phase was picked |
---|
490 | // within the Pick's travel time entry per 10000 picks. |
---|
491 | pPick->dAffPPD = pPick->ttt.dAssocStrength / 2.0; |
---|
492 | if(pPick->dAffPPD > 2.0) |
---|
493 | pPick->dAffPPD = 2.0; |
---|
494 | } |
---|
495 | else |
---|
496 | { |
---|
497 | pPick->dAffPPD = 1.0; |
---|
498 | } |
---|
499 | |
---|
500 | // Composite affinity statistic |
---|
501 | pPick->dAff = pOrg->dAffGap * pOrg->dAffNumArr * pPick->dAffRes * pPick->dAffDis * pPick->dAffPPD; |
---|
502 | |
---|
503 | return(0); |
---|
504 | } // end CalcPickAffinity() |
---|
505 | |
---|
506 | |
---|
507 | void OPCalc_SetTravelTimePointer(ITravelTime * pTravTime) |
---|
508 | { |
---|
509 | pTT = pTravTime; |
---|
510 | } |
---|
511 | |
---|
512 | |
---|
513 | int OPCalc_CalculateOriginPickParams(PICK * pPick, ORIGIN * pOrg) |
---|
514 | { |
---|
515 | TTEntry ttt; |
---|
516 | TTEntry * pttt; |
---|
517 | double dDistKM; |
---|
518 | |
---|
519 | if(!(pTT && pPick && pOrg)) |
---|
520 | { |
---|
521 | CDebug::Log(DEBUG_MINOR_ERROR,"OPCalc_CalculateOriginPickParams(): Null Pointer: " |
---|
522 | "pPick(%u) pOrg(%u) pTT(%u)\n", |
---|
523 | pPick, pOrg, pTT); |
---|
524 | return(-1); |
---|
525 | } |
---|
526 | |
---|
527 | if(pPick->dLat < -90.0 || pPick->dLat > 90.0 || pPick->dLon < -180.0 || pPick->dLon > 180.0) |
---|
528 | return(-1); |
---|
529 | |
---|
530 | geo_to_km_deg(pOrg->dLat, pOrg->dLon, pPick->dLat, pPick->dLon, |
---|
531 | &dDistKM, &pPick->dDelta, &pPick->dAzm); |
---|
532 | pttt = pTT->TBest(pOrg->dZ, pPick->dDelta, pPick->dT - pOrg->dT, &ttt); |
---|
533 | if(!pttt) |
---|
534 | { |
---|
535 | pPick->bTTTIsValid = false; |
---|
536 | pPick->tRes = 9999.9; |
---|
537 | pPick->ttt.szPhase = Phases[PHASE_Unknown].szName; |
---|
538 | |
---|
539 | return(1); |
---|
540 | } |
---|
541 | |
---|
542 | pPick->dTrav = pttt->dTPhase; |
---|
543 | pPick->dToa = pttt->dTOA; |
---|
544 | pPick->tRes = (pPick->dT - pOrg->dT) - pttt->dTPhase; |
---|
545 | strcpy(pPick->sPhase, pttt->szPhase); |
---|
546 | memcpy(&pPick->ttt, pttt, sizeof(TTEntry)); |
---|
547 | pPick->bTTTIsValid = true; |
---|
548 | CDebug::Log(DEBUG_MINOR_INFO,"Iterate(): Eq[%d] d:%.2f z:%.2f tobs:%.2f azm:%.2f ttt:%d\n", |
---|
549 | pOrg->nEq, pPick->dDelta, pOrg->dZ, pPick->dT - pOrg->dT, pPick->dAzm, pttt); |
---|
550 | |
---|
551 | return(0); |
---|
552 | } // end OPCalc_CalculateOriginPickParams() |
---|
553 | |
---|
554 | |
---|