Source
1406
1406
diff = ppa - pa[k];
1407
1407
t = diff >= 0 ? 0.5 : -0.5;
1408
1408
Int npi = Int(diff/C::pi + t);
1409
1409
fitpa[k] = pa[k] + npi*C::pi;
1410
1410
}
1411
1411
fitpa[0] = pa[0];
1412
1412
// Do least squares fit
1413
1413
if (!_rmLsqFit (pars, wsq, fitpa, paerr)) {
1414
1414
return false;
1415
1415
}
1416
-
if (pars(4) < chiSq) {
1417
-
chiSq = pars(4);
1416
+
if (pars[4] < chiSq) {
1417
+
chiSq = pars[4];
1418
1418
nTurns = h; // Number of turns
1419
-
rmFitted = pars(0); // Fitted RM
1420
-
rmErrFitted = pars(1); // Error in RM
1421
-
pa0Fitted = pars(2); // Fitted intrinsic angle
1422
-
pa0ErrFitted = pars(3); // Error in angle
1423
-
rChiSqFitted = pars(4); // Recued chi squared
1419
+
rmFitted = pars[0]; // Fitted RM
1420
+
rmErrFitted = pars[1]; // Error in RM
1421
+
pa0Fitted = pars[2]; // Fitted intrinsic angle
1422
+
pa0ErrFitted = pars[3]; // Error in angle
1423
+
rChiSqFitted = pars[4]; // Recued chi squared
1424
1424
if (n > 2) {
1425
1425
rChiSqFitted /= Float(n - 2);
1426
1426
}
1427
1427
}
1428
1428
}
1429
1429
return true;
1430
1430
}
1431
1431
1432
-
Bool ImagePolarimetry::_rmSupplementaryFit(Float& nTurns, Float& rmFitted, Float& rmErrFitted,
1433
-
Float& pa0Fitted, Float& pa0ErrFitted,
1434
-
Float& rChiSqFitted, const Vector<Float>& wsq,
1435
-
const Vector<Float>& pa, const Vector<Float>& paerr)
1436
-
{
1437
-
1438
-
// For supplementary points find lowest residual RM
1439
-
1440
-
const uInt n = wsq.nelements();
1441
-
//
1442
-
Float absRM = 1e30;
1443
-
Vector<Float> fitpa(pa.copy());
1444
-
Vector<Float> pars;
1445
-
for (Int i=-2; i<3; i++) {
1446
-
fitpa(n-1) = pa(n-1) + C::pi*i;
1447
-
1448
-
// Do least squares fit
1449
-
1450
-
if (! _rmLsqFit (pars, wsq, fitpa, paerr)) return false;
1451
-
1452
-
// Save solution with lowest absolute RM
1453
-
1454
-
if (abs(pars(0)) < absRM) {
1455
-
absRM = abs(pars(0));
1456
-
//
1457
-
nTurns = i; // nTurns
1458
-
rmFitted = pars(0); // Fitted RM
1459
-
rmErrFitted = pars(1); // Error in RM
1460
-
pa0Fitted = pars(2); // Fitted intrinsic angle
1461
-
pa0ErrFitted = pars(3); // Error in angle
1462
-
rChiSqFitted = pars(4); // Reduced chi squared
1463
-
if (n > 2) rChiSqFitted /= Float(n - 2);
1464
-
}
1465
-
1466
-
}
1467
-
//
1468
-
return true;
1432
+
Bool ImagePolarimetry::_rmSupplementaryFit(
1433
+
Float& nTurns, Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted,
1434
+
Float& pa0ErrFitted, Float& rChiSqFitted, const Vector<Float>& wsq,
1435
+
const Vector<Float>& pa, const Vector<Float>& paerr
1436
+
) {
1437
+
// For supplementary points find lowest residual RM
1438
+
const uInt n = wsq.size();
1439
+
Float absRM = 1e30;
1440
+
Vector<Float> fitpa(pa.copy());
1441
+
Vector<Float> pars;
1442
+
for (Int i=-2; i<3; ++i) {
1443
+
fitpa[n-1] = pa[n-1] + C::pi*i;
1444
+
// Do least squares fit
1445
+
if (! _rmLsqFit (pars, wsq, fitpa, paerr)) {
1446
+
return false;
1447
+
}
1448
+
// Save solution with lowest absolute RM
1449
+
if (abs(pars[0]) < absRM) {
1450
+
absRM = abs(pars[0]);
1451
+
nTurns = i; // nTurns
1452
+
rmFitted = pars(0); // Fitted RM
1453
+
rmErrFitted = pars[1]; // Error in RM
1454
+
pa0Fitted = pars[2]; // Fitted intrinsic angle
1455
+
pa0ErrFitted = pars[3]; // Error in angle
1456
+
rChiSqFitted = pars[4]; // Reduced chi squared
1457
+
if (n > 2) {
1458
+
rChiSqFitted /= Float(n - 2);
1459
+
}
1460
+
}
1461
+
}
1462
+
return true;
1469
1463
}
1470
1464
1471
-
1472
-
1473
-
1474
1465
Bool ImagePolarimetry::_rmLsqFit (Vector<Float>& pars, const Vector<Float>& wsq,
1475
1466
const Vector<Float> pa, const Vector<Float>& paerr) const
1476
1467
{
1477
1468
1478
1469
// Perform fit on unmasked data
1479
1470
1480
1471
static Vector<Float> solution;
1481
1472
try {
1482
1473
solution = _fitter->fit(wsq, pa, paerr);
1483
1474
} catch (AipsError x) {