Skip to content
Snippets Groups Projects
Commit 3ee28f88 authored by Bas van der Tol's avatar Bas van der Tol
Browse files

Compute phases per baseline in Simulator

parent da49cbf7
No related branches found
No related tags found
No related merge requests found
Pipeline #17584 failed
...@@ -38,7 +38,7 @@ void radec2lmn(const Position& reference, const Position& position, ...@@ -38,7 +38,7 @@ void radec2lmn(const Position& reference, const Position& position,
* *
* \f[ \mathrm{phases}(p) = e^{\mathrm{stationphases}(p)} \f] * \f[ \mathrm{phases}(p) = e^{\mathrm{stationphases}(p)} \f]
* *
* @param nStation Number of stations * @param nBaseline Number of baselines
* @param nChannel Number of channels * @param nChannel Number of channels
* @param lmn LMN coordinates of source, should be length 3 * @param lmn LMN coordinates of source, should be length 3
* @param uvw Station UVW coordinates, matrix of shape (3, nSt) * @param uvw Station UVW coordinates, matrix of shape (3, nSt)
...@@ -46,7 +46,7 @@ void radec2lmn(const Position& reference, const Position& position, ...@@ -46,7 +46,7 @@ void radec2lmn(const Position& reference, const Position& position,
* @param shift Output matrix, shift per station, matrix of shape (3, nSt) * @param shift Output matrix, shift per station, matrix of shape (3, nSt)
* @param stationPhases Output vector, store per station \f$(x_1,y_1)\f$ * @param stationPhases Output vector, store per station \f$(x_1,y_1)\f$
*/ */
void phases(size_t nStation, size_t nChannel, const double* lmn, void phases(size_t nBaseline, size_t nChannel, const double* lmn,
const casacore::Matrix<double>& uvw, const casacore::Matrix<double>& uvw,
const casacore::Vector<double>& freq, const casacore::Vector<double>& freq,
Simulator::Matrix<dcomplex>& shift, Simulator::Matrix<dcomplex>& shift,
...@@ -63,7 +63,7 @@ Simulator::Simulator(const Position& reference, size_t nStation, ...@@ -63,7 +63,7 @@ Simulator::Simulator(const Position& reference, size_t nStation,
const std::vector<Baseline>& baselines, const std::vector<Baseline>& baselines,
const casacore::Vector<double>& freq, const casacore::Vector<double>& freq,
const casacore::Vector<double>& chanWidths, const casacore::Vector<double>& chanWidths,
const casacore::Matrix<double>& stationUVW, const casacore::Matrix<double>& uvw,
casacore::Cube<dcomplex>& buffer, bool correctFreqSmearing, casacore::Cube<dcomplex>& buffer, bool correctFreqSmearing,
bool stokesIOnly) bool stokesIOnly)
: itsReference(reference), : itsReference(reference),
...@@ -75,12 +75,12 @@ Simulator::Simulator(const Position& reference, size_t nStation, ...@@ -75,12 +75,12 @@ Simulator::Simulator(const Position& reference, size_t nStation,
itsBaselines(baselines), itsBaselines(baselines),
itsFreq(freq), itsFreq(freq),
itsChanWidths(chanWidths), itsChanWidths(chanWidths),
itsStationUVW(stationUVW), itsUVW(uvw),
itsBuffer(buffer), itsBuffer(buffer),
itsShiftBuffer(), itsShiftBuffer(),
itsSpectrumBuffer() { itsSpectrumBuffer() {
itsShiftBuffer.resize(itsNChannel, nStation); itsShiftBuffer.resize(itsNChannel, itsNBaseline);
itsStationPhases.resize(nStation); itsPhases.resize(itsNBaseline);
if (stokesIOnly) { if (stokesIOnly) {
itsSpectrumBuffer.resize(1, itsNChannel); itsSpectrumBuffer.resize(1, itsNChannel);
} else { } else {
...@@ -98,8 +98,8 @@ void Simulator::visit(const PointSource& component) { ...@@ -98,8 +98,8 @@ void Simulator::visit(const PointSource& component) {
radec2lmn(itsReference, component.position(), lmn); radec2lmn(itsReference, component.position(), lmn);
// Compute station phase shifts. // Compute station phase shifts.
phases(itsNStation, itsNChannel, lmn, itsStationUVW, itsFreq, itsShiftBuffer, phases(itsNBaseline, itsNChannel, lmn, itsUVW, itsFreq, itsShiftBuffer,
itsStationPhases); itsPhases);
// Compute component spectrum. // Compute component spectrum.
spectrum(component, itsNChannel, itsFreq, itsSpectrumBuffer, itsStokesIOnly); spectrum(component, itsNChannel, itsFreq, itsSpectrumBuffer, itsStokesIOnly);
...@@ -126,8 +126,7 @@ void Simulator::visit(const PointSource& component) { ...@@ -126,8 +126,7 @@ void Simulator::visit(const PointSource& component) {
if (p == q) { if (p == q) {
buffer += itsNChannel * nCorr; buffer += itsNChannel * nCorr;
} else { } else {
const dcomplex* shiftP = &(itsShiftBuffer(0, p)); const dcomplex* shift = &(itsShiftBuffer(0, bl));
const dcomplex* shiftQ = &(itsShiftBuffer(0, q));
const dcomplex* spectrum = itsSpectrumBuffer.data(); const dcomplex* spectrum = itsSpectrumBuffer.data();
double smearterm = 1.0; double smearterm = 1.0;
...@@ -135,25 +134,23 @@ void Simulator::visit(const PointSource& component) { ...@@ -135,25 +134,23 @@ void Simulator::visit(const PointSource& component) {
for (size_t ch = 0; ch < itsNChannel; ++ch) { for (size_t ch = 0; ch < itsNChannel; ++ch) {
if (itsCorrectFreqSmearing) { if (itsCorrectFreqSmearing) {
smearterm = smearterm =
computeSmearterm(itsStationPhases[q] - itsStationPhases[p], computeSmearterm(itsPhases[bl],
itsChanWidths[ch] / 2); itsChanWidths[ch] / 2);
} }
// Compute baseline phase shift. // Compute baseline phase shift.
// Compute visibilities. // Compute visibilities.
*buffer++ += (*shiftQ) * conj(*shiftP) * (*spectrum++) * smearterm; *buffer++ += (*shift) * (*spectrum++) * smearterm;
++shiftP; ++shift;
++shiftQ;
} // Channels. } // Channels.
} else { } else {
for (size_t ch = 0; ch < itsNChannel; ++ch) { for (size_t ch = 0; ch < itsNChannel; ++ch) {
if (itsCorrectFreqSmearing) if (itsCorrectFreqSmearing)
smearterm = smearterm =
computeSmearterm(itsStationPhases[q] - itsStationPhases[p], computeSmearterm(itsPhases[bl],
itsChanWidths[ch] / 2); itsChanWidths[ch] / 2);
// Compute baseline phase shift. // Compute baseline phase shift.
const dcomplex blShift = (*shiftQ) * conj(*shiftP) * smearterm; const dcomplex blShift = (*shift) * smearterm;
++shiftP; ++shift;
++shiftQ;
// Compute visibilities. // Compute visibilities.
*buffer++ += blShift * (*spectrum++); *buffer++ += blShift * (*spectrum++);
...@@ -172,8 +169,8 @@ void Simulator::visit(const GaussianSource& component) { ...@@ -172,8 +169,8 @@ void Simulator::visit(const GaussianSource& component) {
radec2lmn(itsReference, component.position(), lmn); radec2lmn(itsReference, component.position(), lmn);
// Compute station phase shifts. // Compute station phase shifts.
phases(itsNStation, itsNChannel, lmn, itsStationUVW, itsFreq, itsShiftBuffer, phases(itsNBaseline, itsNChannel, lmn, itsUVW, itsFreq, itsShiftBuffer,
itsStationPhases); itsPhases);
// Compute component spectrum. // Compute component spectrum.
spectrum(component, itsNChannel, itsFreq, itsSpectrumBuffer, itsStokesIOnly); spectrum(component, itsNChannel, itsFreq, itsSpectrumBuffer, itsStokesIOnly);
...@@ -208,11 +205,8 @@ void Simulator::visit(const GaussianSource& component) { ...@@ -208,11 +205,8 @@ void Simulator::visit(const GaussianSource& component) {
if (p == q) { if (p == q) {
buffer += itsNChannel * nCorr; buffer += itsNChannel * nCorr;
} else { } else {
double u = itsStationUVW(0, q); double u = itsUVW(0, bl);
double v = itsStationUVW(1, q); double v = itsUVW(1, bl);
u -= itsStationUVW(0, p);
v -= itsStationUVW(1, p);
// Rotate (u, v) by the position angle and scale with the major // Rotate (u, v) by the position angle and scale with the major
// and minor axis lengths (FWHM in rad). // and minor axis lengths (FWHM in rad).
...@@ -224,22 +218,20 @@ void Simulator::visit(const GaussianSource& component) { ...@@ -224,22 +218,20 @@ void Simulator::visit(const GaussianSource& component) {
const double uvPrime = (-2.0 * casacore::C::pi * casacore::C::pi) * const double uvPrime = (-2.0 * casacore::C::pi * casacore::C::pi) *
(uPrime * uPrime + vPrime * vPrime); (uPrime * uPrime + vPrime * vPrime);
const dcomplex* shiftP = &(itsShiftBuffer(0, p)); const dcomplex* shift = &(itsShiftBuffer(0, bl));
const dcomplex* shiftQ = &(itsShiftBuffer(0, q));
const dcomplex* spectrum = itsSpectrumBuffer.data(); const dcomplex* spectrum = itsSpectrumBuffer.data();
if (itsStokesIOnly) { if (itsStokesIOnly) {
for (size_t ch = 0; ch < itsNChannel; ++ch) { for (size_t ch = 0; ch < itsNChannel; ++ch) {
// Compute baseline phase shift. // Compute baseline phase shift.
dcomplex blShift = (*shiftQ) * conj(*shiftP); dcomplex blShift = (*shift);
++shiftP; ++shift;
++shiftQ;
const double ampl = exp((itsFreq[ch] * itsFreq[ch]) / const double ampl = exp((itsFreq[ch] * itsFreq[ch]) /
(casacore::C::c * casacore::C::c) * uvPrime); (casacore::C::c * casacore::C::c) * uvPrime);
if (itsCorrectFreqSmearing) if (itsCorrectFreqSmearing)
blShift *= blShift *=
computeSmearterm(itsStationPhases[q] - itsStationPhases[p], computeSmearterm(itsPhases[bl],
itsChanWidths[ch] / 2); itsChanWidths[ch] / 2);
blShift *= ampl; blShift *= ampl;
...@@ -250,16 +242,15 @@ void Simulator::visit(const GaussianSource& component) { ...@@ -250,16 +242,15 @@ void Simulator::visit(const GaussianSource& component) {
} else { } else {
for (size_t ch = 0; ch < itsNChannel; ++ch) { for (size_t ch = 0; ch < itsNChannel; ++ch) {
// Compute baseline phase shift. // Compute baseline phase shift.
dcomplex blShift = (*shiftQ) * conj(*shiftP); dcomplex blShift = (*shift);
++shiftP; ++shift;
++shiftQ;
const double ampl = exp((itsFreq[ch] * itsFreq[ch]) / const double ampl = exp((itsFreq[ch] * itsFreq[ch]) /
(casacore::C::c * casacore::C::c) * uvPrime); (casacore::C::c * casacore::C::c) * uvPrime);
if (itsCorrectFreqSmearing) if (itsCorrectFreqSmearing)
blShift *= blShift *=
computeSmearterm(itsStationPhases[q] - itsStationPhases[p], computeSmearterm(itsPhases[bl],
itsChanWidths[ch] / 2); itsChanWidths[ch] / 2);
blShift *= ampl; blShift *= ampl;
...@@ -298,22 +289,22 @@ inline float computeSmearterm(double uvw, double halfwidth) { ...@@ -298,22 +289,22 @@ inline float computeSmearterm(double uvw, double halfwidth) {
} }
// Compute station phase shifts. // Compute station phase shifts.
inline void phases(size_t nStation, size_t nChannel, const double* lmn, inline void phases(size_t nBaseline, size_t nChannel, const double* lmn,
const casacore::Matrix<double>& uvw, const casacore::Matrix<double>& uvw,
const casacore::Vector<double>& freq, const casacore::Vector<double>& freq,
Simulator::Matrix<dcomplex>& shift, Simulator::Matrix<dcomplex>& shift_out,
std::vector<double>& stationPhases) { std::vector<double>& phases_out) {
dcomplex* shiftdata = shift.data(); dcomplex* shiftdata = shift_out.data();
const double cinv = 1 / casacore::C::c; const double cinv = 1 / casacore::C::c;
for (size_t st = 0; st < nStation; ++st) { for (size_t bl = 0; bl < nBaseline; ++bl) {
const double phase = casacore::C::_2pi * cinv * const double phase = casacore::C::_2pi * cinv *
(uvw(0, st) * lmn[0] + uvw(1, st) * lmn[1] + (uvw(0, bl) * lmn[0] + uvw(1, bl) * lmn[1] +
uvw(2, st) * (lmn[2] - 1.0)); uvw(2, bl) * (lmn[2] - 1.0));
stationPhases[st] = phase; phases_out[bl] = phase;
for (size_t ch = 0; ch < nChannel; ++ch) { for (size_t ch = 0; ch < nChannel; ++ch) {
const double chPhase = stationPhases[st] * freq[ch]; const double chPhase = phases_out[bl] * freq[ch];
*shiftdata = dcomplex(cos(chPhase), sin(chPhase)); *shiftdata = dcomplex(cos(chPhase), sin(chPhase));
++shiftdata; ++shiftdata;
} // Channels. } // Channels.
......
...@@ -87,10 +87,10 @@ class Simulator : public ModelComponentVisitor { ...@@ -87,10 +87,10 @@ class Simulator : public ModelComponentVisitor {
const std::vector<Baseline> itsBaselines; const std::vector<Baseline> itsBaselines;
const casacore::Vector<double> itsFreq; const casacore::Vector<double> itsFreq;
const casacore::Vector<double> itsChanWidths; const casacore::Vector<double> itsChanWidths;
const casacore::Matrix<double> itsStationUVW; const casacore::Matrix<double> itsUVW;
casacore::Cube<dcomplex> itsBuffer; casacore::Cube<dcomplex> itsBuffer;
Matrix<dcomplex> itsShiftBuffer; Matrix<dcomplex> itsShiftBuffer;
std::vector<double> itsStationPhases; std::vector<double> itsPhases;
Matrix<dcomplex> itsSpectrumBuffer; Matrix<dcomplex> itsSpectrumBuffer;
}; };
......
...@@ -38,12 +38,16 @@ Simulator MakeSimulator(bool correct_freq_smearing, bool stokes_i_only, ...@@ -38,12 +38,16 @@ Simulator MakeSimulator(bool correct_freq_smearing, bool stokes_i_only,
} }
std::vector<double> chan_widths(kNChan, 1.0e6); std::vector<double> chan_widths(kNChan, 1.0e6);
casacore::Matrix<double> uvw(3, kNStations); size_t nBaseline = baselines.size();
for (size_t st = 0; st < kNStations; ++st) { casacore::Matrix<double> uvw(3, nBaseline);
uvw(0, st) = st * 5000;
uvw(1, st) = st * 1000; for (size_t bl = 0; bl < nBaseline; ++bl) {
uvw(2, st) = 0; size_t p = baselines[bl].first;
size_t q = baselines[bl].second;
uvw(0, bl) = (q-p) * 5000;
uvw(1, bl) = (q-p) * 1000;
uvw(2, bl) = 0;
} }
return Simulator(kReference, kNStations, baselines, chan_freqs, chan_widths, return Simulator(kReference, kNStations, baselines, chan_freqs, chan_widths,
......
...@@ -191,9 +191,9 @@ void OnePredict::initializeThreadData() { ...@@ -191,9 +191,9 @@ void OnePredict::initializeThreadData() {
const size_t nCr = itsStokesIOnly ? 1 : info().ncorr(); const size_t nCr = itsStokesIOnly ? 1 : info().ncorr();
const size_t nThreads = getInfo().nThreads(); const size_t nThreads = getInfo().nThreads();
itsStationUVW.resize(3, nSt); // itsStationUVW.resize(3, nSt);
itsUVWSplitIndex = base::nsetupSplitUVW(info().nantenna(), info().getAnt1(), // itsUVWSplitIndex = base::nsetupSplitUVW(info().nantenna(), info().getAnt1(),
info().getAnt2()); // info().getAnt2());
if (!itsPredictBuffer) { if (!itsPredictBuffer) {
itsPredictBuffer = std::make_shared<base::PredictBuffer>(); itsPredictBuffer = std::make_shared<base::PredictBuffer>();
...@@ -316,8 +316,8 @@ bool OnePredict::process(const DPBuffer& bufin) { ...@@ -316,8 +316,8 @@ bool OnePredict::process(const DPBuffer& bufin) {
itsTimerPredict.start(); itsTimerPredict.start();
base::nsplitUVW(itsUVWSplitIndex, itsBaselines, tempBuffer.getUVW(), // base::nsplitUVW(itsUVWSplitIndex, itsBaselines, tempBuffer.getUVW(),
itsStationUVW); // itsStationUVW);
double time = tempBuffer.getTime(); double time = tempBuffer.getTime();
// Set up directions for beam evaluation // Set up directions for beam evaluation
...@@ -374,7 +374,7 @@ bool OnePredict::process(const DPBuffer& bufin) { ...@@ -374,7 +374,7 @@ bool OnePredict::process(const DPBuffer& bufin) {
(itsApplyBeam ? itsPredictBuffer->GetPatchModel(thread) (itsApplyBeam ? itsPredictBuffer->GetPatchModel(thread)
: itsPredictBuffer->GetModel(thread)); : itsPredictBuffer->GetModel(thread));
simulators.emplace_back(itsPhaseRef, nSt, itsBaselines, info().chanFreqs(), simulators.emplace_back(itsPhaseRef, nSt, itsBaselines, info().chanFreqs(),
info().chanWidths(), itsStationUVW, simulatedest, info().chanWidths(), tempBuffer.getUVW(), simulatedest,
itsCorrectFreqSmearing, itsStokesIOnly); itsCorrectFreqSmearing, itsStokesIOnly);
} }
std::vector<base::Patch::ConstPtr> curPatches(pool->NThreads()); std::vector<base::Patch::ConstPtr> curPatches(pool->NThreads());
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment