[2176] | 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 | |
---|
[3212] | 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 | |
---|
[2176] | 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 | |
---|
[3212] | 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 | |
---|
[2176] | 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 |
---|
[3212] | 337 | CDebug::Log(DEBUG_MAJOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n", |
---|
[2176] | 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 |
---|
[3212] | 345 | pOrg->dAffGap = 4 * OPCalc_BellCurve(pOrg->dGap/360.0); |
---|
| 346 | if(pOrg->nEq < OPCALC_nMinStatPickCount && pOrg->dAffGap < 1.0) |
---|
[2176] | 347 | pOrg->dAffGap = 1.0; |
---|
[3212] | 348 | |
---|
[2176] | 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 |
---|
[3212] | 363 | pOrg->dAffNumArr = log10((double)(pOrg->nEq)); |
---|
[2176] | 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; |
---|
[3212] | 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 | |
---|
[2176] | 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 | |
---|