Ignore:
Timestamp:
12/27/07 13:05:23 (12 years ago)
Author:
paulf
Message:

sync from hydra_proj circa Oct 2007

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/seismic_processing/glass/src/modules/OPCalc/opcalc.cpp

    r2176 r3212  
    88#include <opcalc.h> 
    99#include <opcalc_const.h> 
     10 
     11extern "C" { 
     12void __cdecl qsortd (double *base, unsigned num); 
     13} 
     14 
     15 
     16int OPCalc_ProcessOriginSimple(ORIGIN * pOrg, PICK ** pPick, int nPck); 
     17 
    1018 
    1119// Lib Global variables - File Scope 
     
    3644 
    3745 
     46int __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 
     61int __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   
     76int __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 
     92int 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} 
     183int 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 
    38325int OPCalc_CalcOriginAffinity(ORIGIN * pOrg, PICK ** ppPick, int nPck)  
    39326{ 
     
    48335 
    49336  // Log the calculated values 
    50         CDebug::Log(DEBUG_MINOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n", 
     337        CDebug::Log(DEBUG_MAJOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n", 
    51338                pOrg->dDelMin, pOrg->dDelMod, pOrg->dDelMax); 
    52339 
     
    56343    //         pOrg->dAffGap = 1.0+OPCalc_BellCurve(pOrg->dGap/360.0); 
    57344    // dAffGap range is 0.0 - 2.0 
    58         if(pOrg->nEq < OPCALC_nMinStatPickCount) 
     345  pOrg->dAffGap = 4 * OPCalc_BellCurve(pOrg->dGap/360.0);  
     346        if(pOrg->nEq < OPCALC_nMinStatPickCount && pOrg->dAffGap < 1.0) 
    59347                pOrg->dAffGap = 1.0; 
    60   else 
    61     pOrg->dAffGap = 4 * OPCalc_BellCurve(pOrg->dGap/360.0);  
     348 
    62349//    pOrg->dAffGap = 3 * OPCalc_BellCurve(pOrg->dGap/360.0);  
    63350 
     
    74361                pOrg->dAffNumArr = 1.0; 
    75362        else 
    76                 pOrg->dAffNumArr = log10((float)(pOrg->nEq)); 
     363                pOrg->dAffNumArr = log10((double)(pOrg->nEq)); 
    77364 
    78365        for(i=0; i<nPck; i++)  
     
    99386  { 
    100387    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 
    101415  } 
    102416  else 
Note: See TracChangeset for help on using the changeset viewer.