source: ntrip/trunk/BNC/src/PPP/pppFilter.cpp@ 9088

Last change on this file since 9088 was 9088, checked in by stuerze, 4 years ago

adjusted allocation of slip and LTI according to the respective RTCM version
  • Property svn:keywords set to Author Date Id Rev URL;svn:eol-style=native
  • Property svn:mime-type set to text/plain
File size: 20.1 KB
Line 
1/* -------------------------------------------------------------------------
2 * BKG NTRIP Client
3 * -------------------------------------------------------------------------
4 *
5 * Class: t_pppFilter
6 *
7 * Purpose: Filter Adjustment
8 *
9 * Author: L. Mervart
10 *
11 * Created: 29-Jul-2014
12 *
13 * Changes:
14 *
15 * -----------------------------------------------------------------------*/
16
17#include <iostream>
18#include <iomanip>
19#include <cmath>
20#include <newmat.h>
21#include <newmatio.h>
22#include <newmatap.h>
23
24#include "pppFilter.h"
25#include "bncutils.h"
26#include "pppParlist.h"
27#include "pppObsPool.h"
28#include "pppStation.h"
29#include "pppClient.h"
30
31using namespace BNC_PPP;
32using namespace std;
33
34// Constructor
35////////////////////////////////////////////////////////////////////////////
36t_pppFilter::t_pppFilter(t_pppObsPool* obsPool) {
37 _parlist = 0;
38 _numSat = 0;
39 _obsPool = obsPool;
40 _refPrn = t_prn();
41 _datumTrafo = new t_datumTrafo();
42}
43
44// Destructor
45////////////////////////////////////////////////////////////////////////////
46t_pppFilter::~t_pppFilter() {
47 delete _parlist;
48 delete _datumTrafo;
49}
50
51// Process Single Epoch
52////////////////////////////////////////////////////////////////////////////
53t_irc t_pppFilter::processEpoch(int num) {
54 _numSat = 0;
55 _numEpoProcessing = num;
56 const double maxSolGap = 60.0;
57
58 if (!_parlist) {
59 _parlist = new t_pppParlist();
60 }
61
62 // Vector of all Observations
63 // --------------------------
64 t_pppObsPool::t_epoch* epoch = _obsPool->lastEpoch();
65 if (!epoch) {
66 return failure;
67 }
68 vector<t_pppSatObs*>& allObs = epoch->obsVector();
69
70 // Time of the Epoch
71 // -----------------
72 _epoTime = epoch->epoTime();
73
74 if (!_firstEpoTime.valid() ||
75 !_lastEpoTimeOK.valid() ||
76 (maxSolGap > 0.0 && _epoTime - _lastEpoTimeOK > maxSolGap)) {
77 _firstEpoTime = _epoTime;
78 }
79
80 string epoTimeStr = string(_epoTime);
81
82 //--
83 // Set Parameters
84 // --------------
85 _parlist->set(_epoTime, allObs, _obsPool->getRefSatMap());
86 const vector<t_pppParam*>& params = _parlist->params();
87
88 // Status Vector, Variance-Covariance Matrix
89 // -----------------------------------------
90 ColumnVector xFltOld = _xFlt;
91 SymmetricMatrix QFltOld = _QFlt;
92
93 _QFlt.ReSize(_parlist->nPar()); _QFlt = 0.0;
94 _xFlt.ReSize(_parlist->nPar()); _xFlt = 0.0;
95 _x0.ReSize(_parlist->nPar()); _x0 = 0.0;
96 for (unsigned ii = 0; ii < params.size(); ii++) {
97 const t_pppParam* par1 = params[ii];
98 _x0[ii] = par1->x0();
99 int iOld = par1->indexOld();
100 if (iOld < 0) {
101 _QFlt[ii][ii] = par1->sigma0() * par1->sigma0(); // new parameter
102 }
103 else {
104 _QFlt[ii][ii] = QFltOld[iOld][iOld] + par1->noise() * par1->noise();
105 _xFlt[ii] = xFltOld[iOld];
106 for (unsigned jj = 0; jj < ii; jj++) {
107 const t_pppParam* par2 = params[jj];
108 int jOld = par2->indexOld();
109 if (jOld >= 0) {
110 _QFlt[ii][jj] = QFltOld(iOld+1,jOld+1);
111 }
112 }
113 }
114 }
115 predictCovCrdPart(QFltOld);
116
117 // Init Datum Trafo
118 // ----------------------------------------
119 if ((OPT->_obsModelType == OPT->DCMcodeBias ||
120 OPT->_obsModelType == OPT->DCMphaseBias) &&(_numEpoProcessing == 1)) {
121 _numAllUsedLCs = 0;
122 for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
123 char system = OPT->systems()[iSys];
124 _numAllUsedLCs += OPT->LCs(system).size();
125 if (OPT->_pseudoObsIono && epoch->pseudoObsIono() == false) {
126 _numAllUsedLCs -= 1; // GIM not used
127 if (OPT->_pseudoObsTropo) { // if no VTEC values /pseudo obs iono then no pseudo obs tropo
128 _numAllUsedLCs -= 1;
129 }
130 }
131 int modifyLCs = 0;
132 if (OPT->_pseudoObsTropo &&
133 OPT->_pseudoObsIono && epoch->pseudoObsIono() == true) {
134 modifyLCs = -1;
135 }
136 // max Obs
137 int maxObs = allObs.size() * (_numAllUsedLCs + modifyLCs);
138 if (OPT->_pseudoObsIono && epoch->pseudoObsIono() == true) {
139 maxObs -= 1; // pseudo obs iono with respect to refSat
140 }
141 if (OPT->_pseudoObsIono && epoch->pseudoObsIono() == true &&
142 OPT->_pseudoObsTropo) {
143 maxObs += 1; // pseudo obs tropo once per station
144 }
145 _datumTrafo->initAA(maxObs, _parlist->nPar());
146 }
147 }
148
149 // Pre-Process Satellite Systems separately
150 // ----------------------------------------
151 bool preProcessing = false;
152 if (OPT->_obsModelType == OPT->DCMcodeBias ||
153 OPT->_obsModelType == OPT->DCMphaseBias) {
154 preProcessing = true;
155 for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
156 (iSys) ? _datumTrafo->setFirstSystem(false) : _datumTrafo->setFirstSystem(true);
157 char system = OPT->systems()[iSys];
158 if (OPT->_refSatRequired) {
159 _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
160 if (_obsPool->hasHistoricalRefSat(_refPrn)) {
161 LOG << epoTimeStr << " Warning: prevent to process erroneous refSat again!";
162 return failure;
163 }
164 }
165 vector<t_pppSatObs*> obsVector;
166 for (unsigned jj = 0; jj < allObs.size(); jj++) {
167 if (allObs[jj]->prn().system() == system) {
168 obsVector.push_back(allObs[jj]);
169 }
170 }
171 if (processSystem(OPT->LCs(system), obsVector, _refPrn,
172 epoch->pseudoObsIono(), preProcessing) != success) {
173 return failure;
174 }
175 }
176 // refSat change required?
177 // -----------------------
178 if (_obsPool->refSatChangeRequired()) {
179 _xFlt = xFltOld;
180 _QFlt = QFltOld;
181 return success;
182 }
183 else if (!_obsPool->refSatChangeRequired() &&
184 _numEpoProcessing > 1) {
185 datumTransformation(xFltOld, QFltOld);
186 }
187 }
188
189 // Process Satellite Systems separately
190 // ------------------------------------
191 preProcessing = false;
192 for (unsigned iSys = 0; iSys < OPT->systems().size(); iSys++) {
193 char system = OPT->systems()[iSys];
194 if (OPT->_refSatRequired) {
195 _refPrn = (_obsPool->getRefSatMapElement(system))->prn();
196 }
197 unsigned int num = 0;
198 vector<t_pppSatObs*> obsVector;
199 for (unsigned jj = 0; jj < allObs.size(); jj++) {
200 if (allObs[jj]->prn().system() == system) {
201 obsVector.push_back(allObs[jj]);
202 num++;
203 }
204 }
205 LOG << epoTimeStr << " SATNUM " << system << ' ' << right << setw(2) << num << endl;
206 if (processSystem(OPT->LCs(system), obsVector, _refPrn,
207 epoch->pseudoObsIono(), preProcessing) != success) {
208 return failure;
209 }
210 }
211
212 // close epoch processing
213 // ----------------------
214 cmpDOP(allObs);
215 _parlist->printResult(_epoTime, _QFlt, _xFlt);
216 _lastEpoTimeOK = _epoTime; // remember time of last successful epoch processing
217
218 return success;
219}
220
221// Process Selected LCs
222////////////////////////////////////////////////////////////////////////////
223t_irc t_pppFilter::processSystem(const vector<t_lc::type>& LCs,
224 const vector<t_pppSatObs*>& obsVector,
225 const t_prn& refPrn,
226 bool pseudoObsIonoAvailable,
227 bool preProcessing) {
228 LOG.setf(ios::fixed);
229
230 // Detect Cycle Slips
231 // ------------------
232 if (detectCycleSlips(LCs, obsVector, refPrn, preProcessing) != success) {
233 return failure;
234 }
235
236 unsigned usedLCs = LCs.size();
237
238 ColumnVector xSav = _xFlt;
239 SymmetricMatrix QSav = _QFlt;
240 string epoTimeStr = string(_epoTime);
241 const vector<t_pppParam*>& params = _parlist->params();
242
243 if (OPT->_pseudoObsIono && !pseudoObsIonoAvailable) {
244 usedLCs -= 1; // GIM not used
245 if (OPT->_pseudoObsTropo) { // if no VTEC values /pseudo obs iono then no pseudo obs tropo
246 usedLCs -= 1;
247 }
248 }
249
250 int modifyLCs = 0;
251 if (OPT->_pseudoObsTropo &&
252 OPT->_pseudoObsIono && pseudoObsIonoAvailable) {
253 modifyLCs = -1;
254 }
255 // max Obs
256 unsigned maxObs = obsVector.size() * (usedLCs + modifyLCs);
257 if (OPT->_pseudoObsIono && pseudoObsIonoAvailable) {
258 maxObs -= 1; // pseudo obs iono with respect to refSat
259 }
260 if (OPT->_pseudoObsIono && pseudoObsIonoAvailable &&
261 OPT->_pseudoObsTropo) {
262 maxObs += 1; // pseudo obs tropo once per station
263 }
264
265 // Outlier Detection Loop
266 // ----------------------
267 for (unsigned iOutlier = 0; iOutlier < maxObs; iOutlier++) {
268
269 if (iOutlier > 0) {
270 _xFlt = xSav;
271 _QFlt = QSav;
272 }
273
274 // First-Design Matrix, Terms Observed-Computed, Weight Matrix
275 // -----------------------------------------------------------
276 Matrix AA(maxObs, _parlist->nPar());
277 ColumnVector ll(maxObs);
278 DiagonalMatrix PP(maxObs); PP = 0.0;
279
280 int iObs = -1;
281 vector<t_pppSatObs*> usedObs;
282 vector<t_lc::type> usedTypes;
283 for (unsigned ii = 0; ii < obsVector.size(); ii++) {
284 t_pppSatObs* obs = obsVector[ii];
285 unsigned hlp = ii+1;
286 if (!obs->outlier()) {
287 for (unsigned jj = 0; jj < usedLCs; jj++) {
288 const t_lc::type tLC = LCs[jj];
289 if (tLC == t_lc::GIM && obs->isReference()) {continue;}
290 if (tLC == t_lc::Tz0 && hlp != obsVector.size()) {continue;}
291 ++iObs;
292 usedObs.push_back(obs);
293 usedTypes.push_back(tLC);
294 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
295 const t_pppParam* par = params[iPar];
296 AA[iObs][iPar] = par->partial(_epoTime, obs, tLC, refPrn);
297 }
298 ll[iObs] = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA.Row(iObs+1));
299 PP[iObs] = 1.0 / (obs->sigma(tLC) * obs->sigma(tLC));
300 }
301 }
302 }
303
304 if ((preProcessing && ! iOutlier) &&
305 (OPT->_obsModelType == OPT->DCMcodeBias || OPT->_obsModelType == OPT->DCMphaseBias)) {
306 _datumTrafo->updateIndices(maxObs);
307 _datumTrafo->prepareAA(AA, _numEpoProcessing, _parlist->nPar());
308 if (_obsPool->refSatChangeRequired()) { // from detectCycleSlips()
309 return success;
310 }
311 }
312
313 // Check number of observations, truncate matrices
314 // -----------------------------------------------
315 if (iObs == -1) {
316 return failure;
317 }
318 AA = AA.Rows(1, iObs+1);
319 ll = ll.Rows(1, iObs+1);
320 PP = PP.SymSubMatrix(1, iObs+1);
321
322 // Kalman update step
323 // ------------------
324 kalman(AA, ll, PP, _QFlt, _xFlt);
325
326 // Check Residuals
327 // ---------------
328 ColumnVector vv = AA * _xFlt - ll;
329 double maxOutlier = 0.0;
330 int maxOutlierIndex = -1;
331 t_lc::type maxOutlierLC = t_lc::dummy;
332 for (unsigned ii = 0; ii < usedObs.size(); ii++) {
333 const t_lc::type tLC = usedTypes[ii];
334 double res = fabs(vv[ii]);
335 if (res > usedObs[ii]->maxRes(tLC)) {
336 if (res > fabs(maxOutlier)) {
337 maxOutlier = vv[ii];
338 maxOutlierIndex = ii;
339 maxOutlierLC = tLC;
340 }
341 }
342 }
343
344 // Mark outlier or break outlier detection loop
345 // --------------------------------------------
346 if (maxOutlierIndex > -1) {
347 t_pppSatObs* obs = usedObs[maxOutlierIndex];
348 t_pppParam* par = 0;
349 if (!preProcessing) {
350 LOG << epoTimeStr << " Outlier " << t_lc::toString(maxOutlierLC) << ' '
351 << obs->prn().toString() << ' '
352 << setw(8) << setprecision(4) << maxOutlier << endl;
353 }
354 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
355 t_pppParam* hlp = params[iPar];
356 if (hlp->type() == t_pppParam::amb &&
357 hlp->prn() == obs->prn() &&
358 hlp->tLC() == usedTypes[maxOutlierIndex]) {
359 par = hlp;
360 }
361 }
362 if (preProcessing) {
363 if (par && obs->prn() == refPrn) {
364 _obsPool->setRefSatChangeRequired(true);
365 break;
366 }
367 else if (par && obs->prn() != refPrn) {
368 par->setAmbResetCandidate();
369 }
370 }
371 else {
372 if (par) {
373 if (par->ambResetCandidate() || _obsPool->hasHistoricalRefSat(par->prn())) {
374 resetAmb(par->prn(), obsVector, &QSav, &xSav);
375 }
376 else {
377 par->setAmbResetCandidate();
378 obs->setOutlier();
379 }
380 }
381 else {
382 obs->setOutlier();
383 }
384 }
385 }
386 // Print Residuals
387 // ---------------
388 else {
389 if (!preProcessing) {
390 for (unsigned jj = 0; jj < LCs.size(); jj++) {
391 for (unsigned ii = 0; ii < usedObs.size(); ii++) {
392 const t_lc::type tLC = usedTypes[ii];
393 t_pppSatObs* obs = usedObs[ii];
394 if (tLC == LCs[jj]) {
395 obs->setRes(tLC, vv[ii]);
396 LOG << epoTimeStr << " RES "
397 << left << setw(3) << t_lc::toString(tLC) << right << ' '
398 << obs->prn().toString().substr(0,3) << ' '
399 << setw(8) << setprecision(4) << vv[ii] << endl;
400 }
401 }
402 }
403 }
404 break;
405 }
406 }
407
408 return success;
409}
410
411// Cycle-Slip Detection
412////////////////////////////////////////////////////////////////////////////
413t_irc t_pppFilter::detectCycleSlips(const vector<t_lc::type>& LCs,
414 const vector<t_pppSatObs*>& obsVector,
415 const t_prn& refPrn,
416 bool preProcessing) {
417 if (_obsPool->hasHistoricalRefSat(refPrn)) {
418 return success;
419 }
420
421 //double SLIP = 20.0; // slip threshold
422 double SLIP = 200.0; // slip threshold deactivated
423 string epoTimeStr = string(_epoTime);
424 const vector<t_pppParam*>& params = _parlist->params();
425
426 for (unsigned ii = 0; ii < LCs.size(); ii++) {
427 const t_lc::type& tLC = LCs[ii];
428 if (t_lc::includesPhase(tLC)) {
429 for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) {
430 const t_pppSatObs* obs = obsVector[iObs];
431
432 if (preProcessing && obs->prn() != refPrn) {continue;}
433
434 // Check set Slips and Jump Counters
435 // ---------------------------------
436 bool slip = false;
437 if (obs->slip()) {
438 LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString() << endl;
439 slip = true;
440 }
441
442 if (_slips[obs->prn()]._obsSlipCounter != -1 &&
443 _slips[obs->prn()]._obsSlipCounter != obs->slipCounter()) {
444 LOG << epoTimeStr << "cycle slip set (obsSlipCounter) " << obs->prn().toString() << endl;
445 slip = true;
446 }
447 if (!preProcessing) {
448 _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
449 }
450 if (_slips[obs->prn()]._biasJumpCounter != -1 &&
451 _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) {
452 LOG << epoTimeStr << "cycle slip set (biasJumpCounter) " << obs->prn().toString() << endl;
453 slip = true;
454 }
455 if (!preProcessing) {
456 _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
457 }
458 // Slip Set
459 // --------
460 if (slip) {
461 if (preProcessing) {
462 _obsPool->setRefSatChangeRequired(true);
463 }
464 else {
465 resetAmb(obs->prn(), obsVector);
466 }
467 }
468 // Check Pre-Fit Residuals
469 // -----------------------
470 else {
471 ColumnVector AA(params.size());
472 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
473 const t_pppParam* par = params[iPar];
474 AA[iPar] = par->partial(_epoTime, obs, tLC, refPrn);
475 }
476
477 double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
478 double vv = DotProduct(AA, _xFlt) - ll;
479
480 if (fabs(vv) > SLIP) {
481 if (preProcessing) {
482 _obsPool->setRefSatChangeRequired(true);
483 }
484 else {
485 LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
486 << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
487 resetAmb(obs->prn(), obsVector);
488 }
489 }
490 }
491 }
492 }
493 }
494 return success;
495}
496
497// Reset Ambiguity Parameter (cycle slip)
498////////////////////////////////////////////////////////////////////////////
499t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*>& obsVector,
500 SymmetricMatrix* QSav, ColumnVector* xSav) {
501
502 t_irc irc = failure;
503 vector<t_pppParam*>& params = _parlist->params();
504 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
505 t_pppParam* par = params[iPar];
506 if (par->type() == t_pppParam::amb && par->prn() == prn) {
507 int ind = par->indexNew();
508 t_lc::type tLC = par->tLC();
509 LOG << string(_epoTime) << " RESET " << par->toString() << endl;
510 delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector);
511 par->setIndex(ind);
512 params[iPar] = par;
513 for (unsigned ii = 1; ii <= params.size(); ii++) {
514 _QFlt(ii, ind+1) = 0.0;
515 if (QSav) {
516 (*QSav)(ii, ind+1) = 0.0;
517 }
518 }
519 _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0();
520 if (QSav) {
521 (*QSav)(ind+1,ind+1) = _QFlt(ind+1,ind+1);
522 }
523 _xFlt[ind] = 0.0;
524 if (xSav) {
525 (*xSav)[ind] = _xFlt[ind];
526 }
527 _x0[ind] = par->x0();
528 irc = success;
529 }
530 }
531
532 return irc;
533}
534
535// Compute various DOP Values
536////////////////////////////////////////////////////////////////////////////
537void t_pppFilter::cmpDOP(const vector<t_pppSatObs*>& obsVector) {
538
539 _dop.reset();
540
541 try {
542 const unsigned numPar = 4;
543 Matrix AA(obsVector.size(), numPar);
544 _numSat = 0;
545 for (unsigned ii = 0; ii < obsVector.size(); ii++) {
546 t_pppSatObs* obs = obsVector[ii];
547 char system = obs->prn().system();
548 t_prn refPrn = t_prn();
549 if (OPT->_refSatRequired) {
550 refPrn = _obsPool->getRefSatMapElement(system)->prn();
551 }
552 if (obs->isValid() && !obs->outlier()) {
553 ++_numSat;
554 for (unsigned iPar = 0; iPar < numPar; iPar++) {
555 const t_pppParam* par = _parlist->params()[iPar];
556 AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn);
557 }
558 }
559 }
560 if (_numSat < 4) {
561 return;
562 }
563 AA = AA.Rows(1, _numSat);
564 SymmetricMatrix NN; NN << AA.t() * AA;
565 SymmetricMatrix QQ = NN.i();
566
567 _dop.H = sqrt(QQ(1,1) + QQ(2,2));
568 _dop.V = sqrt(QQ(3,3));
569 _dop.P = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
570 _dop.T = sqrt(QQ(4,4));
571 _dop.G = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3) + QQ(4,4));
572 }
573 catch (...) {
574 }
575}
576
577// Compute various DOP Values
578////////////////////////////////////////////////////////////////////////////
579void t_pppFilter::predictCovCrdPart(const SymmetricMatrix& QFltOld) {
580
581 const vector<t_pppParam*>& params = _parlist->params();
582 if (params.size() < 3) {
583 return;
584 }
585
586 bool first = (params[0]->indexOld() < 0);
587
588 SymmetricMatrix Qneu(3); Qneu = 0.0;
589 for (unsigned ii = 0; ii < 3; ii++) {
590 const t_pppParam* par = params[ii];
591 if (first) {
592 Qneu[ii][ii] = par->sigma0() * par->sigma0();
593 }
594 else {
595 Qneu[ii][ii] = par->noise() * par->noise();
596 }
597 }
598
599 const t_pppStation* sta = PPP_CLIENT->staRover();
600 SymmetricMatrix Qxyz(3);
601 covariNEU_XYZ(Qneu, sta->ellApr().data(), Qxyz);
602
603 if (first) {
604 _QFlt.SymSubMatrix(1,3) = Qxyz;
605 }
606 else {
607 double dt = _epoTime - _firstEpoTime;
608 if (dt < OPT->_seedingTime) {
609 _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3);
610 }
611 else {
612 _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3) + Qxyz;
613 }
614 }
615}
616
617// Compute datum transformation
618////////////////////////////////////////////////////////////////////////////
619void t_pppFilter::datumTransformation(const ColumnVector& xFltOld, const SymmetricMatrix& QFltOld) {
620 Matrix D21 = _datumTrafo->computeTrafoMatrix();
621 _QFlt << D21 * QFltOld * D21.t();
622 _xFlt = D21 * xFltOld;
623}
624
625
626// Reset Ambiguity Parameter (cycle slip)
627////////////////////////////////////////////////////////////////////////////
628t_irc t_pppFilter::addInfiniteNoise(t_pppParam::e_type para) {
629 t_irc irc = failure;
630 vector<t_pppParam*>& params = _parlist->params();
631 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
632 t_pppParam* par = params[iPar];
633 if (par->type() == para) {
634 int ind = par->indexNew();
635 if (ind < 0) {
636 return irc;
637 }
638 _QFlt(ind+1,ind+1) += par->sigma0() * par->sigma0();
639 irc = success;
640 }
641 }
642 return irc;
643}
Note: See TracBrowser for help on using the repository browser.