12 const std::vector<float> &obsVal,
13 const std::vector<float> &obsErr,
14 const std::vector<float> &bkgVal,
15 const std::vector<float> &bkgErr,
16 const std::vector<float> &PdBad,
17 const bool ModelLevels,
18 std::vector<int> &flags,
19 std::vector<float> &PGE,
20 std::vector<float> &PGEBd,
22 const std::vector<float> *obsVal2,
23 const std::vector<float> *bkgVal2)
25 const float missingValueFloat = util::missingValue(1.0f);
27 const double PGEMult = 1000.0;
29 const double PGEMDI = 1.111;
39 const double SDiffCrit = obsVal2 && bkgVal2 ?
43 const size_t numLevelsToCheck = obsVal.size();
46 PGEBd.assign(obsVal.size(), PGEMDI);
57 const bool obsErrEmpty = obsErr.empty();
58 const bool bkgErrEmpty = bkgErr.empty();
60 for (
size_t jlev = 0; jlev < numLevelsToCheck; ++jlev) {
62 if (!obsErrEmpty && !bkgErrEmpty &&
63 obsErr[jlev] >= 0 && bkgErr[jlev] >= 0) {
64 ErrVar = std::pow(ObErrMult * obsErr[jlev], 2) +
65 std::pow(BkgErrMult * bkgErr[jlev], 2);
67 ErrVar = missingValueFloat;
70 if (ErrVarMax > 0.0) {
71 ErrVar = std::min(ErrVar, ErrVarMax);
75 if (obsVal[jlev] != missingValueFloat &&
76 bkgVal[jlev] != missingValueFloat &&
77 ErrVar != missingValueFloat) {
78 if (obsVal2 && bkgVal2 &&
79 (*obsVal2)[jlev] != missingValueFloat &&
80 (*bkgVal2)[jlev] != missingValueFloat) {
81 SDiff = (std::pow(obsVal[jlev] - bkgVal[jlev], 2) +
82 std::pow((*obsVal2)[jlev] - (*bkgVal2)[jlev], 2)) /
static_cast<double>(ErrVar);
84 PdGood = std::exp(-0.5 * std::min(SDiff, 2.0 * ExpArgMax)) / (2.0 * M_PI * ErrVar);
86 SDiff = std::pow(obsVal[jlev] - bkgVal[jlev], 2) / ErrVar;
88 PdGood = std::exp(-0.5 * std::min(SDiff, 2.0 * ExpArgMax)) /
89 std::sqrt(2.0 * M_PI * ErrVar);
93 PGEBk = (PdBad[jlev] * PGE[jlev]) /
94 (PdBad[jlev] * PGE[jlev] + PdGood * (1.0 - PGE[jlev]));
98 if (PGEBk >= PGECrit) {
108 PGE[jlev] = trunc(PGEBk * PGEMult) + PGE[jlev];
112 (SDiff >= SDiffCrit ||