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

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

minor changes in PPP
  • 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 string epoTimeStr = string(_epoTime);
423 const vector<t_pppParam*>& params = _parlist->params();
424
425 if (preProcessing) {
426 SLIP += 10.0;
427 }
428
429 for (unsigned ii = 0; ii < LCs.size(); ii++) {
430 const t_lc::type& tLC = LCs[ii];
431 if (t_lc::includesPhase(tLC)) {
432 for (unsigned iObs = 0; iObs < obsVector.size(); iObs++) {
433 const t_pppSatObs* obs = obsVector[iObs];
434
435 if (preProcessing && obs->prn() != refPrn) {continue;}
436
437 // Check set Slips and Jump Counters
438 // ---------------------------------
439 bool slip = false;
440 if (obs->slip()) {
441 LOG << epoTimeStr << "cycle slip set (obs) " << obs->prn().toString() << endl;
442 slip = true;
443 }
444
445 if (_slips[obs->prn()]._obsSlipCounter != -1 &&
446 _slips[obs->prn()]._obsSlipCounter > obs->slipCounter()) {
447 LOG << epoTimeStr << "cycle slip set (obsSlipCounter) " << obs->prn().toString() << endl;
448 slip = true;
449 }
450 if (!preProcessing) {
451 _slips[obs->prn()]._obsSlipCounter = obs->slipCounter();
452 }
453 if (_slips[obs->prn()]._biasJumpCounter != -1 &&
454 _slips[obs->prn()]._biasJumpCounter != obs->biasJumpCounter()) {
455 LOG << epoTimeStr << "cycle slip set (biasJumpCounter) " << obs->prn().toString() << endl;
456 slip = true;
457 }
458 if (!preProcessing) {
459 _slips[obs->prn()]._biasJumpCounter = obs->biasJumpCounter();
460 }
461 // Slip Set
462 // --------
463 if (slip) {
464 if (preProcessing) {
465 _obsPool->setRefSatChangeRequired(true);
466 }
467 else {
468 resetAmb(obs->prn(), obsVector);
469 }
470 }
471 // Check Pre-Fit Residuals
472 // -----------------------
473 else {
474 ColumnVector AA(params.size());
475 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
476 const t_pppParam* par = params[iPar];
477 AA[iPar] = par->partial(_epoTime, obs, tLC, refPrn);
478 }
479
480 double ll = obs->obsValue(tLC) - obs->cmpValue(tLC) - DotProduct(_x0, AA);
481 double vv = DotProduct(AA, _xFlt) - ll;
482
483 if (fabs(vv) > SLIP) {
484 if (preProcessing) {
485 _obsPool->setRefSatChangeRequired(true);
486 }
487 else {
488 LOG << epoTimeStr << " cycle slip detected " << t_lc::toString(tLC) << ' '
489 << obs->prn().toString() << ' ' << setw(8) << setprecision(4) << vv << endl;
490 resetAmb(obs->prn(), obsVector);
491 }
492 }
493 }
494 }
495 }
496 }
497 return success;
498}
499
500// Reset Ambiguity Parameter (cycle slip)
501////////////////////////////////////////////////////////////////////////////
502t_irc t_pppFilter::resetAmb(t_prn prn, const vector<t_pppSatObs*>& obsVector,
503 SymmetricMatrix* QSav, ColumnVector* xSav) {
504
505 t_irc irc = failure;
506 vector<t_pppParam*>& params = _parlist->params();
507 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
508 t_pppParam* par = params[iPar];
509 if (par->type() == t_pppParam::amb && par->prn() == prn) {
510 int ind = par->indexNew();
511 t_lc::type tLC = par->tLC();
512 LOG << string(_epoTime) << " RESET " << par->toString() << endl;
513 delete par; par = new t_pppParam(t_pppParam::amb, prn, tLC, &obsVector);
514 par->setIndex(ind);
515 params[iPar] = par;
516 for (unsigned ii = 1; ii <= params.size(); ii++) {
517 _QFlt(ii, ind+1) = 0.0;
518 if (QSav) {
519 (*QSav)(ii, ind+1) = 0.0;
520 }
521 }
522 _QFlt(ind+1,ind+1) = par->sigma0() * par->sigma0();
523 if (QSav) {
524 (*QSav)(ind+1,ind+1) = _QFlt(ind+1,ind+1);
525 }
526 _xFlt[ind] = 0.0;
527 if (xSav) {
528 (*xSav)[ind] = _xFlt[ind];
529 }
530 _x0[ind] = par->x0();
531 irc = success;
532 }
533 }
534
535 return irc;
536}
537
538// Compute various DOP Values
539////////////////////////////////////////////////////////////////////////////
540void t_pppFilter::cmpDOP(const vector<t_pppSatObs*>& obsVector) {
541
542 _dop.reset();
543
544 try {
545 const unsigned numPar = 4;
546 Matrix AA(obsVector.size(), numPar);
547 _numSat = 0;
548 for (unsigned ii = 0; ii < obsVector.size(); ii++) {
549 t_pppSatObs* obs = obsVector[ii];
550 char system = obs->prn().system();
551 t_prn refPrn = t_prn();
552 if (OPT->_refSatRequired) {
553 refPrn = _obsPool->getRefSatMapElement(system)->prn();
554 }
555 if (obs->isValid() && !obs->outlier()) {
556 ++_numSat;
557 for (unsigned iPar = 0; iPar < numPar; iPar++) {
558 const t_pppParam* par = _parlist->params()[iPar];
559 AA[_numSat-1][iPar] = par->partial(_epoTime, obs, t_lc::c1, refPrn);
560 }
561 }
562 }
563 if (_numSat < 4) {
564 return;
565 }
566 AA = AA.Rows(1, _numSat);
567 SymmetricMatrix NN; NN << AA.t() * AA;
568 SymmetricMatrix QQ = NN.i();
569
570 _dop.H = sqrt(QQ(1,1) + QQ(2,2));
571 _dop.V = sqrt(QQ(3,3));
572 _dop.P = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3));
573 _dop.T = sqrt(QQ(4,4));
574 _dop.G = sqrt(QQ(1,1) + QQ(2,2) + QQ(3,3) + QQ(4,4));
575 }
576 catch (...) {
577 }
578}
579
580// Compute various DOP Values
581////////////////////////////////////////////////////////////////////////////
582void t_pppFilter::predictCovCrdPart(const SymmetricMatrix& QFltOld) {
583
584 const vector<t_pppParam*>& params = _parlist->params();
585 if (params.size() < 3) {
586 return;
587 }
588
589 bool first = (params[0]->indexOld() < 0);
590
591 SymmetricMatrix Qneu(3); Qneu = 0.0;
592 for (unsigned ii = 0; ii < 3; ii++) {
593 const t_pppParam* par = params[ii];
594 if (first) {
595 Qneu[ii][ii] = par->sigma0() * par->sigma0();
596 }
597 else {
598 Qneu[ii][ii] = par->noise() * par->noise();
599 }
600 }
601
602 const t_pppStation* sta = PPP_CLIENT->staRover();
603 SymmetricMatrix Qxyz(3);
604 covariNEU_XYZ(Qneu, sta->ellApr().data(), Qxyz);
605
606 if (first) {
607 _QFlt.SymSubMatrix(1,3) = Qxyz;
608 }
609 else {
610 double dt = _epoTime - _firstEpoTime;
611 if (dt < OPT->_seedingTime) {
612 _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3);
613 }
614 else {
615 _QFlt.SymSubMatrix(1,3) = QFltOld.SymSubMatrix(1,3) + Qxyz;
616 }
617 }
618}
619
620// Compute datum transformation
621////////////////////////////////////////////////////////////////////////////
622void t_pppFilter::datumTransformation(const ColumnVector& xFltOld, const SymmetricMatrix& QFltOld) {
623 Matrix D21 = _datumTrafo->computeTrafoMatrix();
624 _QFlt << D21 * QFltOld * D21.t();
625 _xFlt = D21 * xFltOld;
626}
627
628
629// Reset Ambiguity Parameter (cycle slip)
630////////////////////////////////////////////////////////////////////////////
631t_irc t_pppFilter::addInfiniteNoise(t_pppParam::e_type para) {
632 t_irc irc = failure;
633 vector<t_pppParam*>& params = _parlist->params();
634 for (unsigned iPar = 0; iPar < params.size(); iPar++) {
635 t_pppParam* par = params[iPar];
636 if (par->type() == para) {
637 int ind = par->indexNew();
638 if (ind < 0) {
639 return irc;
640 }
641 _QFlt(ind+1,ind+1) += par->sigma0() * par->sigma0();
642 irc = success;
643 }
644 }
645 return irc;
646}
Note: See TracBrowser for help on using the repository browser.