Skip to content

Commit

Permalink
Merge pull request #1144 from kettenis/maxant
Browse files Browse the repository at this point in the history 
Add UVFITS support for large antenna numbers
  • Loading branch information
@tammojan
tammojan committed Dec 3, 2021
2 parents abff44a + f72e8f0 commit a30f2ed
Show file tree
Hide file tree
Showing 3 changed files with 133 additions and 51 deletions.
77 changes: 64 additions & 13 deletions msfits/MSFits/MSFitsInput.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1138,6 +1138,10 @@ void MSFitsInput::fillMSMainTableColWise(Int& nField, Int& nSpW) {
}
// get index for baseline
Int iBsln = getIndex(pType, "BASELINE");
// get indices for subarray, antenna1 and antenna2
Int iSubarr = getIndex(pType, "SUBARRAY");
Int iAnt1 = getIndex(pType, "ANTENNA1");
Int iAnt2 = getIndex(pType, "ANTENNA2");
// get indices for time
Int iTime0 = getIndex(pType, "DATE", 0);
Int iTime1 = getIndex(pType, "DATE", 1);
Expand Down Expand Up @@ -1231,8 +1235,18 @@ void MSFitsInput::fillMSMainTableColWise(Int& nField, Int& nSpW) {
// make 0-based
fieldId = (Int) _priGroup.parm(iSource) - 1;
}
Float baseline = _priGroup.parm(iBsln);
Int arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
Int arrayId = 0;
std::pair<Int, Int> ants;
if (iBsln >= 0) {
Float baseline = _priGroup.parm(iBsln);
ants = _extractAntennas(baseline);
arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
} else {
Int antenna1 = _priGroup.parm(iAnt1);
Int antenna2 = _priGroup.parm(iAnt2);
ants = _extractAntennas(antenna1, antenna2);
arrayId = _priGroup.parm(iSubarr);
}
_nArray = max(_nArray, arrayId + 1);
for (Int k = 0; k < nif; ++k) {
Int index = group * nif + k;
Expand All @@ -1241,7 +1255,6 @@ void MSFitsInput::fillMSMainTableColWise(Int& nField, Int& nSpW) {
uvw(1, index) = _priGroup.parm(iV) * C::c;
uvw(2, index) = _priGroup.parm(iW) * C::c;
// Convert from units of seconds to meters
std::pair<Int, Int> ants = _extractAntennas(baseline);
ant1[index] = ants.first;
ant2[index] = ants.second;
}
Expand Down Expand Up @@ -1461,6 +1474,10 @@ void MSFitsInput::fillMSMainTable(Int& nField, Int& nSpW) {
}
// get index for baseline
Int iBsln = getIndex(pType, "BASELINE");
// get indices for subarray, antenna1 and antenna2
Int iSubarr = getIndex(pType, "SUBARRAY");
Int iAnt1 = getIndex(pType, "ANTENNA1");
Int iAnt2 = getIndex(pType, "ANTENNA2");
// get indices for time
Int iTime0 = getIndex(pType, "DATE", 0);
Int iTime1 = getIndex(pType, "DATE", 1);
Expand Down Expand Up @@ -1556,10 +1573,19 @@ void MSFitsInput::fillMSMainTable(Int& nField, Int& nSpW) {
uvw *= C::c;

// Extract array/baseline/antenna info
Float baseline = _priGroup.parm(iBsln);
Int arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
Int arrayId = 0;
std::pair<Int, Int> ants;
if (iBsln >= 0) {
Float baseline = _priGroup.parm(iBsln);
ants = _extractAntennas(baseline);
arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
} else {
Int antenna1 = _priGroup.parm(iAnt1);
Int antenna2 = _priGroup.parm(iAnt2);
ants = _extractAntennas(antenna1, antenna2);
arrayId = _priGroup.parm(iSubarr);
}
_nArray = max(_nArray, arrayId + 1);
std::pair<Int, Int> ants = _extractAntennas(baseline);
Int ant1 = ants.first;
Int ant2 = ants.second;
// Ensure arrayId-specific params are of correct length:
Expand Down Expand Up @@ -3196,6 +3222,9 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
}

Int iBsln = getIndex(TType, "BASELINE");
Int iSubarr = getIndex(TType, "SUBARRAY");
Int iAnt1 = getIndex(TType, "ANTENNA1");
Int iAnt2 = getIndex(TType, "ANTENNA2");
Int iTime0 = getIndex(TType, "DATE");
Int iSource = getIndex(TType, "SOURCE");
Int iFreq = getIndex(TType, "FREQSEL");
Expand Down Expand Up @@ -3239,9 +3268,18 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
ScalarColumn<Float> colUU(tb, TType(iU));
ScalarColumn<Float> colVV(tb, TType(iV));
ScalarColumn<Float> colWW(tb, TType(iW));
ScalarColumn<Float> colBL(tb, TType(iBsln));
ScalarColumn<Float> colBL;
ScalarColumn<Float> colSubarr, colAnt1, colAnt2;
ArrayColumn<Float> colVis(tb, TType(iVis));

if (iBsln >= 0) {
colBL = ScalarColumn<Float>(tb, TType(iBsln));
} else {
colSubarr = ScalarColumn<Float>(tb, TType(iSubarr));
colAnt1 = ScalarColumn<Float>(tb, TType(iAnt1));
colAnt2 = ScalarColumn<Float>(tb, TType(iAnt2));
}

Int visL = 1;
for (uInt i = 0; i < _nPixel.nelements(); i++)
visL *= _nPixel(i);
Expand Down Expand Up @@ -3275,10 +3313,19 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
uvw *= C::c;

// Extract array/baseline/antenna info
Float baseline = colBL.asfloat(0);
Int arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
Int arrayId;
std::pair<Int, Int> ants;
if (iBsln >= 0) {
Float baseline = colBL.asfloat(0);
ants = _extractAntennas(baseline);
arrayId = Int(100.0 * (baseline - Int(baseline) + 0.001));
} else {
Int antenna1 = _priGroup.parm(iAnt1);
Int antenna2 = _priGroup.parm(iAnt2);
ants = _extractAntennas(antenna1, antenna2);
arrayId = _priGroup.parm(iSubarr) - 1;
}
_nArray = max(_nArray, arrayId + 1);
std::pair<Int, Int> ants = _extractAntennas(baseline);
Int ant1 = ants.first;
Int ant2 = ants.second;

Expand Down Expand Up @@ -3459,9 +3506,7 @@ void MSFitsInput::fillMSMainTable(BinaryTable& bt) {
}
}

std::pair<Int, Int> MSFitsInput::_extractAntennas(Float baseline) {
Int ant1 = Int(baseline) / 256;
Int ant2 = Int(baseline) - ant1 * 256;
std::pair<Int, Int> MSFitsInput::_extractAntennas(Int ant1, Int ant2) {
_nAntRow = max(_nAntRow, ant1);
_nAntRow = max(_nAntRow, ant2);
_uniqueAnts.insert(ant1);
Expand All @@ -3472,6 +3517,12 @@ std::pair<Int, Int> MSFitsInput::_extractAntennas(Float baseline) {
return make_pair(ant1, ant2);
}

std::pair<Int, Int> MSFitsInput::_extractAntennas(Float baseline) {
Int ant1 = Int(baseline) / 256;
Int ant2 = Int(baseline) - ant1 * 256;
return _extractAntennas(ant1, ant2);
}

void MSFitsInput::fillObservationTable(ConstFitsKeywordList& kwl) {
const FitsKeyword* kw;
const Regex trailing(" *$"); // trailing blanks
Expand Down
3 changes: 2 additions & 1 deletion msfits/MSFits/MSFitsInput.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -421,8 +421,9 @@ class MSFitsInput
void readRandomGroupUVFits(Int obsType);
void readPrimaryTableUVFits(Int obsType);

std::pair<Int, Int> _extractAntennas(Int antenna1, Int antenna2);
std::pair<Int, Int> _extractAntennas(Float baseline);

void _fillSysPowerTable(BinaryTable& bt);

void _doFillSysPowerSingleIF(
Expand Down
104 changes: 67 additions & 37 deletions msfits/MSFits/MSFitsOutput.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,6 +366,23 @@ uInt MSFitsOutput::get_tbf_end(const uInt rownr, const uInt nrow,
return tbfend;
}

// Define a FITS random group parameter with default scaling and offset
static void defineRandomParam(Record& ek, Int n, const String& name) {
String ptype = "ptype" + String::toString(n);
String pscal = "pscal" + String::toString(n);
String pzero = "pzero" + String::toString(n);
ek.define(ptype, name);
ek.define(pscal, 1.0);
ek.define(pzero, 0.0);
}

static void defineRandomParam(Record& ek, Int n, const String& name,
const String& comment) {
String ptype = "ptype" + String::toString(n);
defineRandomParam(ek, n, name);
ek.setComment(ptype, comment);
}

std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double& refFreq,
Double& chanbw, const String &outFITSFile,
const Block<Int>& spwidMap, Int nrspw,
Expand Down Expand Up @@ -477,6 +494,10 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
}
}

Vector<Int> antnumbers;
_handleAntNumbers(_ms, antnumbers);
Int maxant = max(antnumbers);

// Also find out what the Stokes are and make sure that they are the same
// throughout the MS. In principle we could handle the same stokes in
// different order by transposing, but this may well never happen.
Expand Down Expand Up @@ -678,8 +699,16 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
// DATE
desc.addField("date1", TpFloat);
desc.addField("date2", TpFloat);
// BASELINE
desc.addField("baseline", TpFloat);
// BASELINE if maximum antenna number < 256
// SUBARRAY, ANTENNA1 and ANTENNA2 otherwise
// (see AIPS memo 117, section 3.1.2)
if (maxant < 256) {
desc.addField("baseline", TpFloat);
} else {
desc.addField("subarray", TpFloat);
desc.addField("antenna1", TpFloat);
desc.addField("antenna2", TpFloat);
}
// FREQSEL
ScalarColumn<Int> inddid(_ms, MS::columnName(MS::DATA_DESC_ID));
desc.addField("freqsel", TpFloat);
Expand Down Expand Up @@ -748,36 +777,23 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
ek.define("crota7", 0.0);

// PTYPE PSCALE PZERO
ek.define("ptype1", "UU");
ek.define("pscal1", 1.0);
ek.define("pzero1", 0.0);
ek.define("ptype2", "VV");
ek.define("pscal2", 1.0);
ek.define("pzero2", 0.0);
ek.define("ptype3", "WW");
ek.define("pscal3", 1.0);
ek.define("pzero3", 0.0);
ek.define("ptype4", "DATE");
ek.define("pscal4", 1.0);
ek.define("pzero4", 0.0);
ek.setComment("ptype4", "Day number");
ek.define("ptype5", "DATE");
ek.define("pscal5", 1.0);
ek.define("pzero5", 0.0);
ek.setComment("ptype5", "Day fraction");
ek.define("ptype6", "BASELINE");
ek.define("pscal6", 1.0);
ek.define("pzero6", 0.0);
ek.define("ptype7", "FREQSEL");
ek.define("pscal7", 1.0);
ek.define("pzero7", 0.0);
Int idx = 1;
defineRandomParam(ek, idx++, "UU");
defineRandomParam(ek, idx++, "VV");
defineRandomParam(ek, idx++, "WW");
defineRandomParam(ek, idx++, "DATE", "Day number");
defineRandomParam(ek, idx++, "DATE", "Day fraction");
if (maxant < 256) {
defineRandomParam(ek, idx++, "BASELINE");
} else {
defineRandomParam(ek, idx++, "SUBARRAY");
defineRandomParam(ek, idx++, "ANTENNA1");
defineRandomParam(ek, idx++, "ANTENNA2");
}
defineRandomParam(ek, idx++, "FREQSEL");
if (asMultiSource) {
ek.define("ptype8", "SOURCE");
ek.define("pscal8", 1.0);
ek.define("pzero8", 0.0);
ek.define("ptype9", "INTTIM");
ek.define("pscal9", 1.0);
ek.define("pzero9", 0.0);
defineRandomParam(ek, idx++, "SOURCE");
defineRandomParam(ek, idx++, "INTTIM");
}

// EXTEND - already written by FITSGroupWriter
Expand Down Expand Up @@ -1046,7 +1062,17 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
RecordFieldPtr<Float> oww(writer.row(), "w");
RecordFieldPtr<Float> odate1(writer.row(), "date1");
RecordFieldPtr<Float> odate2(writer.row(), "date2");
RecordFieldPtr<Float> obaseline(writer.row(), "baseline");
RecordFieldPtr<Float> obaseline;
RecordFieldPtr<Float> osubarray;
RecordFieldPtr<Float> oantenna1;
RecordFieldPtr<Float> oantenna2;
if (maxant < 256) {
obaseline = RecordFieldPtr<Float> (writer.row(), "baseline");
} else {
osubarray = RecordFieldPtr<Float> (writer.row(), "subarray");
oantenna1 = RecordFieldPtr<Float> (writer.row(), "antenna1");
oantenna2 = RecordFieldPtr<Float> (writer.row(), "antenna2");
}
RecordFieldPtr<Float> ofreqsel(writer.row(), "freqsel");
RecordFieldPtr<Float> osource;
RecordFieldPtr<Float> ointtim;
Expand All @@ -1065,9 +1091,6 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
}
}

Vector<Int> antnumbers;
_handleAntNumbers(_ms, antnumbers);

// Loop through all rows.
ProgressMeter meter(0.0, nOutRow * 1.0, "UVFITS Writer", "Rows copied", "",
"", True, nOutRow / 100);
Expand Down Expand Up @@ -1343,8 +1366,15 @@ std::shared_ptr<FitsOutput> MSFitsOutput::_writeMain(Int& refPixelFreq, Double&
*odate2 = dayFraction;

// BASELINE
*obaseline = antnumbers(inant1(tbfrownr)) * 256 + antnumbers(inant2(
tbfrownr)) + inarray(tbfrownr) * 0.01;
if (maxant < 256) {
*obaseline = antnumbers(inant1(tbfrownr)) * 256 +
antnumbers(inant2(tbfrownr)) +
inarray(tbfrownr) * 0.01;
} else {
*osubarray = inarray(tbfrownr) + 1;
*oantenna1 = antnumbers(inant1(tbfrownr));
*oantenna2 = antnumbers(inant2(tbfrownr));
}

// FREQSEL (in the future it might be FREQ_GRP+1)
// *ofreqsel = inddid(i) + 1;
Expand Down

0 comments on commit a30f2ed

Please sign in to comment.