- Timestamp:
- Aug 13, 2014, 1:55:04 PM (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/BNC/src/PPP/pppFilter.cpp
r5916 r5922 180 180 Matrix AA(maxObs, _parlist->nPar()); 181 181 ColumnVector ll(maxObs); 182 UpperTriangularMatrix Sl(maxObs); Sl= 0.0;182 DiagonalMatrix PP(maxObs); PP = 0.0; 183 183 184 184 int iObs = -1; 185 vector<t_pppSatObs*> 186 vector<t_lc::type> usedTypes;185 vector<t_pppSatObs*> usedObs; 186 vector<t_lc::type> usedTypes; 187 187 for (unsigned ii = 0; ii < obsVector.size(); ii++) { 188 188 t_pppSatObs* obs = obsVector[ii]; 189 189 if (!obs->outlier()) { 190 Matrix CC(LCs.size(), 4);191 190 for (unsigned jj = 0; jj < LCs.size(); jj++) { 192 191 const t_lc::type tLC = LCs[jj]; … … 198 197 AA[iObs][iPar] = par->partial(_epoTime, obs, tLC); 199 198 } 200 ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1)); 201 if (LCs.size() > 1) { 202 ColumnVector coeff(4); 203 obs->lc(tLC, 0.0, 0.0, 0.0, 0.0, &coeff); 204 CC[jj][0] = coeff[0] * obs->sigma(t_lc::l1); 205 CC[jj][1] = coeff[1] * obs->sigma(t_lc::l2); 206 CC[jj][2] = coeff[2] * obs->sigma(t_lc::c1); 207 CC[jj][3] = coeff[3] * obs->sigma(t_lc::c2); 208 } 209 else { 210 Sl[iObs][iObs] = obs->sigma(tLC); 211 } 212 } 213 if (LCs.size() > 1) { 214 SymmetricMatrix QQ; QQ << CC * CC.t(); 215 Sl.SymSubMatrix(iObs-LCs.size()+2, iObs+1) = Cholesky(QQ).t(); 199 ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1)); 200 PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC)); 216 201 } 217 202 } … … 225 210 AA = AA.Rows(1, iObs+1); 226 211 ll = ll.Rows(1, iObs+1); 227 Sl = Sl.SymSubMatrix(1, iObs+1);212 PP = PP.SymSubMatrix(1, iObs+1); 228 213 229 214 // Kalman update step 230 215 // ------------------ 231 DiagonalMatrix PP(iObs+1); PP = 0.0;232 for (int ii = 1; ii <= iObs+1; ii++) {233 PP(ii,ii) = 1.0 / (Sl(ii,ii) * Sl(ii,ii));234 }235 216 kalman(AA, ll, PP, _QFlt, _xFlt); 236 217
Note:
See TracChangeset
for help on using the changeset viewer.