Source
1166
1166
}
1167
1167
1168
1168
void ImagePolarimetry::_fiddleStokesCoordinate(
1169
1169
ImageInterface<Complex>& ie, Stokes::StokesTypes type
1170
1170
) const {
1171
1171
CoordinateSystem cSys = ie.coordinates();
1172
1172
_fiddleStokesCoordinate(cSys, type);
1173
1173
ie.setCoordinateInfo(cSys);
1174
1174
}
1175
1175
1176
-
void ImagePolarimetry::_fiddleTimeCoordinate(ImageInterface<Complex>& ie, const Quantum<Double>& f,
1177
-
Int coord) const
1178
-
{
1179
-
LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1180
-
//
1181
-
CoordinateSystem cSys = ie.coordinates();
1182
-
Coordinate* pC = cSys.coordinate(coord).clone();
1183
-
AlwaysAssert(pC->nPixelAxes()==1,AipsError);
1184
-
AlwaysAssert(pC->type()==Coordinate::LINEAR,AipsError);
1185
-
//
1186
-
Vector<String> axisUnits = pC->worldAxisUnits();
1187
-
axisUnits = String("s");
1188
-
if (!pC->setWorldAxisUnits(axisUnits)) {
1189
-
os << "Failed to set TimeCoordinate units to seconds because " << pC->errorMessage() << LogIO::EXCEPTION;
1190
-
}
1191
-
1192
-
// Find factor to convert from time (s) to rad/m/m
1193
-
1194
-
Vector<Double> inc = pC->increment();
1195
-
Double ff = f.getValue(Unit("Hz"));
1196
-
Double lambda = QC::c( ).getValue(Unit("m/s")) / ff;
1197
-
Double fac = -C::pi * ff / 2.0 / lambda / lambda;
1198
-
inc *= fac;
1199
-
//
1200
-
Vector<String> axisNames(1);
1201
-
axisNames = String("RotationMeasure");
1202
-
axisUnits = String("rad/m/m");
1203
-
Vector<Double> refVal(1,0.0);
1204
-
//
1205
-
LinearCoordinate lC(axisNames, axisUnits, refVal, inc,
1206
-
pC->linearTransform().copy(), pC->referencePixel().copy());
1207
-
//
1208
-
cSys.replaceCoordinate(lC, coord);
1209
-
ie.setCoordinateInfo(cSys);
1210
-
delete pC;
1176
+
void ImagePolarimetry::_fiddleTimeCoordinate(
1177
+
ImageInterface<Complex>& ie, const Quantum<Double>& f, Int coord
1178
+
) const {
1179
+
LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1180
+
CoordinateSystem cSys = ie.coordinates();
1181
+
auto_ptr<Coordinate> pC(cSys.coordinate(coord).clone());
1182
+
AlwaysAssert(pC->nPixelAxes()==1,AipsError);
1183
+
AlwaysAssert(pC->type()==Coordinate::LINEAR,AipsError);
1184
+
auto axisUnits = pC->worldAxisUnits();
1185
+
axisUnits = String("s");
1186
+
ThrowIf(
1187
+
! pC->setWorldAxisUnits(axisUnits),
1188
+
"Failed to set TimeCoordinate units to seconds because "
1189
+
+ pC->errorMessage()
1190
+
);
1191
+
// Find factor to convert from time (s) to rad/m/m
1192
+
auto inc = pC->increment();
1193
+
const auto ff = f.getValue(Unit("Hz"));
1194
+
const auto lambda = QC::c( ).getValue(Unit("m/s")) / ff;
1195
+
const auto fac = -C::pi * ff / 2.0 / lambda / lambda;
1196
+
inc *= fac;
1197
+
Vector<String> axisNames(1);
1198
+
axisNames = String("RotationMeasure");
1199
+
axisUnits = String("rad/m/m");
1200
+
Vector<Double> refVal(1,0.0);
1201
+
LinearCoordinate lC(
1202
+
axisNames, axisUnits, refVal, inc,
1203
+
pC->linearTransform().copy(), pC->referencePixel().copy()
1204
+
);
1205
+
cSys.replaceCoordinate(lC, coord);
1206
+
ie.setCoordinateInfo(cSys);
1211
1207
}
1212
1208
1213
-
1214
-
Quantum<Double> ImagePolarimetry::_findCentralFrequency(const Coordinate& coord, Int shape) const
1215
-
{
1216
-
AlwaysAssert(coord.nPixelAxes()==1,AipsError);
1217
-
//
1218
-
Vector<Double> pixel(1);
1219
-
Vector<Double> world;
1220
-
pixel(0) = Double(shape - 1) / 2.0;
1221
-
if (!coord.toWorld(world, pixel)) {
1222
-
LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1223
-
os << "Failed to convert pixel to world for SpectralCoordinate because "
1224
-
<< coord.errorMessage() << LogIO::EXCEPTION;
1225
-
}
1226
-
Vector<String> units = coord.worldAxisUnits();
1227
-
return Quantum<Double>(world(0), units(0));
1209
+
Quantum<Double> ImagePolarimetry::_findCentralFrequency(
1210
+
const Coordinate& coord, Int shape
1211
+
) const {
1212
+
AlwaysAssert(coord.nPixelAxes()==1,AipsError);
1213
+
Vector<Double> pixel(1);
1214
+
Vector<Double> world;
1215
+
pixel(0) = Double(shape - 1) / 2.0;
1216
+
ThrowIf(
1217
+
! coord.toWorld(world, pixel),
1218
+
"Failed to convert pixel to world for SpectralCoordinate because "
1219
+
+ coord.errorMessage()
1220
+
);
1221
+
const auto units = coord.worldAxisUnits();
1222
+
return Quantum<Double>(world(0), units(0));
1228
1223
}
1229
1224
1230
-
1231
-
Int ImagePolarimetry::_findSpectralCoordinate(const CoordinateSystem& cSys, LogIO& os,
1232
-
Bool fail) const
1233
-
{
1234
-
Int afterCoord = -1;
1235
-
Int coord = cSys.findCoordinate(Coordinate::SPECTRAL, afterCoord);
1236
-
if (coord<0) {
1237
-
if (fail) os << "No spectral coordinate in this image" << LogIO::EXCEPTION;
1238
-
}
1239
-
if (afterCoord>0) {
1240
-
os << LogIO::WARN << "This image has more than one spectral coordinate; only first used"
1241
-
<< LogIO::POST;
1242
-
}
1243
-
return coord;
1244
-
}
1245
-
1246
-
Bool ImagePolarimetry::_findRotationMeasure (Float& rmFitted, Float& rmErrFitted,
1247
-
Float& pa0Fitted, Float& pa0ErrFitted,
1248
-
Float& rChiSqFitted, Float& nTurns,
1249
-
const Vector<uInt>& sortidx,
1250
-
const Vector<Float>& wsq2, const Vector<Float>& pa2,
1251
-
const Array<Bool>& paMask2,
1252
-
const Array<Float>& paerr2,
1253
-
Float rmFg, Float rmMax, Float maxPaErr,
1254
-
const String& posString)
1255
-
//
1256
-
// wsq is lambda squared in m**2 in increasing wavelength order
1257
-
// pa is position angle in radians
1258
-
// paerr is pa error in radians
1259
-
// maxPaErr is maximum tolerated error in position angle
1260
-
// rmfg is a user specified foreground RM rad/m/m
1261
-
// rmmax is a user specified maximum RM
1262
-
//
1263
-
{
1264
-
static Vector<Float> paerr;
1265
-
static Vector<Float> pa;
1266
-
static Vector<Float> wsq;
1267
-
1268
-
// Abandon if less than 2 points
1269
-
1270
-
uInt n = sortidx.nelements();
1271
-
rmFitted = rmErrFitted = pa0Fitted = pa0ErrFitted = rChiSqFitted = 0.0;
1272
-
if (n<2) return false;
1273
-
1274
-
// Sort into decreasing frequency order and correct for foreground rotation
1275
-
// Remember wsq already sorted. Discard points that are too noisy or masked
1276
-
1277
-
const Vector<Float>& paerr1(paerr2.nonDegenerate(0));
1278
-
const Vector<Bool>& paMask1(paMask2.nonDegenerate(0));
1279
-
paerr.resize(n);
1280
-
pa.resize(n);
1281
-
wsq.resize(n);
1282
-
//
1283
-
uInt j = 0;
1284
-
for (uInt i=0; i<n; i++) {
1285
-
if (abs(paerr1(sortidx(i)))<maxPaErr && paMask1(sortidx(i))) {
1286
-
pa(j) = pa2(sortidx(i)) - rmFg*wsq2(i);
1287
-
paerr(j) = paerr1(sortidx(i));
1288
-
wsq(j) = wsq2(i);
1289
-
j++;
1290
-
}
1291
-
}
1292
-
n = j;
1293
-
if (n<=1) return false;
1294
-
//
1295
-
pa.resize(n,true);
1296
-
paerr.resize(n,true);
1297
-
wsq.resize(n, true);
1298
-
1299
-
// Treat supplementary and primary points separately
1300
-
1301
-
Bool ok = false;
1302
-
if (n==2) {
1303
-
ok = _rmSupplementaryFit(nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted,
1304
-
rChiSqFitted, wsq, pa, paerr);
1305
-
} else {
1306
-
ok = _rmPrimaryFit(nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted,
1307
-
rChiSqFitted, wsq, pa, paerr, rmMax, /*plotter,*/ posString);
1308
-
}
1309
-
1310
-
// Put position angle into the range 0->pi
1311
-
1312
-
static MVAngle tmpMVA1;
1313
-
if (ok) {
1314
-
MVAngle tmpMVA0(pa0Fitted);
1315
-
tmpMVA1 = tmpMVA0.binorm(0.0);
1316
-
pa0Fitted = tmpMVA1.radian();
1317
-
1318
-
// Add foreground back on
1319
-
1320
-
rmFitted += rmFg;
1321
-
}
1322
-
return ok;
1225
+
Int ImagePolarimetry::_findSpectralCoordinate(
1226
+
const CoordinateSystem& cSys, LogIO& os, Bool fail
1227
+
) const {
1228
+
Int afterCoord = -1;
1229
+
Int coord = cSys.findCoordinate(Coordinate::SPECTRAL, afterCoord);
1230
+
ThrowIf(coord < 0 && fail, "No spectral coordinate in this image");
1231
+
if (afterCoord>0) {
1232
+
os << LogIO::WARN << "This image has more than one spectral "
1233
+
<< "coordinate; only first used" << LogIO::POST;
1234
+
}
1235
+
return coord;
1323
1236
}
1324
1237
1325
-
void ImagePolarimetry::_hasQU () const
1326
-
{
1327
-
Bool has = _stokes[ImagePolarimetry::Q]!=0 &&
1328
-
_stokes[ImagePolarimetry::U]!=0;
1329
-
if (!has) {
1330
-
LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1331
-
os << "This image does not have Stokes Q and U which are required for this function" << LogIO::EXCEPTION;
1332
-
}
1238
+
Bool ImagePolarimetry::_findRotationMeasure(
1239
+
Float& rmFitted, Float& rmErrFitted, Float& pa0Fitted, Float& pa0ErrFitted,
1240
+
Float& rChiSqFitted, Float& nTurns, const Vector<uInt>& sortidx,
1241
+
const Vector<Float>& wsq2, const Vector<Float>& pa2,
1242
+
const Array<Bool>& paMask2, const Array<Float>& paerr2, Float rmFg,
1243
+
Float rmMax, Float maxPaErr, const String& posString
1244
+
) {
1245
+
// wsq is lambda squared in m**2 in increasing wavelength order
1246
+
// pa is position angle in radians
1247
+
// paerr is pa error in radians
1248
+
// maxPaErr is maximum tolerated error in position angle
1249
+
// rmfg is a user specified foreground RM rad/m/m
1250
+
// rmmax is a user specified maximum RM
1251
+
static Vector<Float> paerr;
1252
+
static Vector<Float> pa;
1253
+
static Vector<Float> wsq;
1254
+
// Abandon if less than 2 points
1255
+
uInt n = sortidx.size();
1256
+
if (n<2) {
1257
+
return false;
1258
+
}
1259
+
rmFitted = rmErrFitted = pa0Fitted = pa0ErrFitted = rChiSqFitted = 0.0;
1260
+
// Sort into decreasing frequency order and correct for foreground rotation
1261
+
// Remember wsq already sorted. Discard points that are too noisy or masked
1262
+
const Vector<Float>& paerr1(paerr2.nonDegenerate(0));
1263
+
const Vector<Bool>& paMask1(paMask2.nonDegenerate(0));
1264
+
paerr.resize(n);
1265
+
pa.resize(n);
1266
+
wsq.resize(n);
1267
+
uInt j = 0;
1268
+
for (uInt i=0; i<n; ++i) {
1269
+
if (abs(paerr1(sortidx(i))) < maxPaErr && paMask1(sortidx(i))) {
1270
+
pa(j) = pa2(sortidx(i)) - rmFg*wsq2(i);
1271
+
paerr(j) = paerr1(sortidx(i));
1272
+
wsq(j) = wsq2(i);
1273
+
++j;
1274
+
}
1275
+
}
1276
+
n = j;
1277
+
if (n<=1) {
1278
+
return false;
1279
+
}
1280
+
pa.resize(n, true);
1281
+
paerr.resize(n, true);
1282
+
wsq.resize(n, true);
1283
+
// Treat supplementary and primary points separately
1284
+
Bool ok = n == 2
1285
+
? _rmSupplementaryFit(
1286
+
nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted,
1287
+
rChiSqFitted, wsq, pa, paerr
1288
+
)
1289
+
: _rmPrimaryFit(
1290
+
nTurns, rmFitted, rmErrFitted, pa0Fitted, pa0ErrFitted,
1291
+
rChiSqFitted, wsq, pa, paerr, rmMax, /*plotter,*/ posString
1292
+
);
1293
+
// Put position angle into the range 0->pi
1294
+
static MVAngle tmpMVA1;
1295
+
if (ok) {
1296
+
MVAngle tmpMVA0(pa0Fitted);
1297
+
tmpMVA1 = tmpMVA0.binorm(0.0);
1298
+
pa0Fitted = tmpMVA1.radian();
1299
+
// Add foreground back on
1300
+
rmFitted += rmFg;
1301
+
}
1302
+
return ok;
1333
1303
}
1334
1304
1335
-
1336
-
ImageExpr<Float> ImagePolarimetry::_makeStokesExpr(ImageInterface<Float>* imPtr,
1337
-
Stokes::StokesTypes type, const String& name) const
1338
-
{
1339
-
LogIO os(LogOrigin("ImagePolarimetry", __func__, WHERE));
1340
-
if (imPtr==0) {
1341
-
os << "This image does not have Stokes " << Stokes::name(type) << LogIO::EXCEPTION;
1342
-
}
1343
-
1344
-
// Make node.
1345
-
1346
-
LatticeExprNode node(*imPtr);
1347
-
1348
-
// Make expression
1349
-
1350
-
LatticeExpr<Float> le(node);
1351
-
ImageExpr<Float> ie(le, name);
1352
-
ie.setUnits(_image->units());
1353
-
_fiddleStokesCoordinate(ie, type);
1354
-
1355
-
return ie;
1305
+
void ImagePolarimetry::_hasQU () const {
1306
+
ThrowIf(
1307
+
! (_stokes[Q] && _stokes[U]),
1308
+
"This image does not have Stokes Q and U which are "
1309
+
"required for this function"
1310
+
);
1356
1311
}
1357
1312
1313
+
ImageExpr<Float> ImagePolarimetry::_makeStokesExpr(
1314
+
ImageInterface<Float>* imPtr, Stokes::StokesTypes type, const String& name
1315
+
) const {
1316
+
ThrowIf(! imPtr, "This image does not have Stokes " + Stokes::name(type));
1317
+
LatticeExprNode node(*imPtr);
1318
+
LatticeExpr<Float> le(node);
1319
+
ImageExpr<Float> ie(le, name);
1320
+
ie.setUnits(_image->units());
1321
+
_fiddleStokesCoordinate(ie, type);
1322
+
return ie;
1323
+
}
1358
1324
1359
-
1360
-
ImageInterface<Float>* ImagePolarimetry::_makeSubImage (IPosition& blc,
1361
-
IPosition& trc,
1362
-
Int axis, Int pix) const
1363
-
{
1364
-
blc(axis) = pix;
1365
-
trc(axis) = pix;
1325
+
ImageInterface<Float>* ImagePolarimetry::_makeSubImage(
1326
+
IPosition& blc, IPosition& trc, Int axis, Int pix
1327
+
) const {
1328
+
blc[axis] = pix;
1329
+
trc[axis] = pix;
1366
1330
LCSlicer slicer(blc, trc, RegionType::Abs);
1367
1331
ImageRegion region(slicer);
1368
1332
return new SubImage<Float>(*_image, region);
1369
1333
}
1370
1334
1371
-
1372
-
LatticeExprNode ImagePolarimetry::_makePolIntNode(LogIO& os, Bool debias, Float clip, Float sigma,
1373
-
Bool doLin, Bool doCirc)
1374
-
{
1375
-
LatticeExprNode linNode, circNode, node;
1376
-
//
1377
-
Float sigma2 = 0.0;
1378
-
if (doLin) {
1379
-
if (debias) {
1380
-
if (sigma > 0.0) {
1381
-
sigma2 = sigma;
1382
-
} else {
1383
-
sigma2 = ImagePolarimetry::sigma(clip);
1384
-
}
1385
-
}
1386
-
linNode = LatticeExprNode(pow(*_stokes[ImagePolarimetry::U],2) +
1387
-
pow(*_stokes[ImagePolarimetry::Q],2));
1388
-
}
1389
-
//
1390
-
if (doCirc) {
1391
-
if (debias) {
1392
-
if (sigma > 0.0) {
1393
-
sigma2 = sigma;
1394
-
} else {
1395
-
sigma2 = ImagePolarimetry::sigma(clip);
1396
-
}
1397
-
}
1398
-
circNode = LatticeExprNode(pow(*_stokes[ImagePolarimetry::V],2));
1399
-
}
1400
-
//
1401
-
Float sigmasq = sigma2 * sigma2;
1402
-
if (doLin && doCirc) {
1403
-
if (debias) {
1404
-
node = linNode + circNode - LatticeExprNode(sigmasq);
1405
-
os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq) << LogIO::POST;
1406
-
} else {
1407
-
node = linNode + circNode;
1408
-
}
1409
-
} else if (doLin) {
1410
-
if (debias) {
1411
-
node = linNode - LatticeExprNode(sigmasq);
1412
-
os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq) << LogIO::POST;
1413
-
} else {
1414
-
node = linNode;
1415
-
}
1416
-
} else if (doCirc) {
1417
-
if (debias) {
1418
-
node = circNode - LatticeExprNode(sigmasq);
1419
-
os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq) << LogIO::POST;
1420
-
} else {
1421
-
node = circNode;
1422
-
}
1423
-
}
1424
-
//
1425
-
return LatticeExprNode(sqrt(node));
1335
+
LatticeExprNode ImagePolarimetry::_makePolIntNode(
1336
+
LogIO& os, Bool debias, Float clip, Float sigma, Bool doLin, Bool doCirc
1337
+
) {
1338
+
LatticeExprNode linNode, circNode, node;
1339
+
Float sigma2 = debias ? (sigma > 0.0 ? sigma : this->sigma(clip)) : 0.0;
1340
+
if (doLin) {
1341
+
linNode = LatticeExprNode(pow(*_stokes[U], 2) + pow(*_stokes[Q], 2));
1342
+
}
1343
+
if (doCirc) {
1344
+
circNode = LatticeExprNode(pow(*_stokes[V], 2));
1345
+
}
1346
+
Float sigmasq = sigma2 * sigma2;
1347
+
if (doLin && doCirc) {
1348
+
node = linNode + circNode;
1349
+
if (debias) {
1350
+
node = node - LatticeExprNode(sigmasq);
1351
+
os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq)
1352
+
<< LogIO::POST;
1353
+
}
1354
+
}
1355
+
else if (doLin) {
1356
+
node = linNode;
1357
+
if (debias) {
1358
+
node = node - LatticeExprNode(sigmasq);
1359
+
os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq)
1360
+
<< LogIO::POST;
1361
+
}
1362
+
}
1363
+
else if (doCirc) {
1364
+
node = circNode;
1365
+
if (debias) {
1366
+
node = node - LatticeExprNode(sigmasq);
1367
+
os << LogIO::NORMAL << "Debiasing with sigma = " << sqrt(sigmasq)
1368
+
<< LogIO::POST;
1369
+
}
1370
+
}
1371
+
return LatticeExprNode(sqrt(node));
1426
1372
}
1427
1373
1428
-
1429
1374
Bool ImagePolarimetry::_rmPrimaryFit(Float& nTurns, Float& rmFitted, Float& rmErrFitted,
1430
1375
Float& pa0Fitted, Float& pa0ErrFitted,
1431
1376
Float& rChiSqFitted, const Vector<Float>& wsq,
1432
1377
const Vector<Float>& pa, const Vector<Float>& paerr,
1433
1378
Float rmMax, const String& /*posString*/)
1434
1379
{
1435
1380
static Vector<Float> plotPA;
1436
1381
static Vector<Float> plotPAErr;
1437
1382
static Vector<Float> plotPAErrY1;
1438
1383
static Vector<Float> plotPAErrY2;