diff --git a/base/Simulator.cc b/base/Simulator.cc index ad77bf620f3077b14191f082f7628845232368ef..649bd96f8f872b5c54259212092197c2713dc358 100644 --- a/base/Simulator.cc +++ b/base/Simulator.cc @@ -38,7 +38,7 @@ void radec2lmn(const Position& reference, const Position& position, * * \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 lmn LMN coordinates of source, should be length 3 * @param uvw Station UVW coordinates, matrix of shape (3, nSt) @@ -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 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::Vector<double>& freq, Simulator::Matrix<dcomplex>& shift, @@ -63,7 +63,7 @@ Simulator::Simulator(const Position& reference, size_t nStation, const std::vector<Baseline>& baselines, const casacore::Vector<double>& freq, const casacore::Vector<double>& chanWidths, - const casacore::Matrix<double>& stationUVW, + const casacore::Matrix<double>& uvw, casacore::Cube<dcomplex>& buffer, bool correctFreqSmearing, bool stokesIOnly) : itsReference(reference), @@ -75,12 +75,12 @@ Simulator::Simulator(const Position& reference, size_t nStation, itsBaselines(baselines), itsFreq(freq), itsChanWidths(chanWidths), - itsStationUVW(stationUVW), + itsUVW(uvw), itsBuffer(buffer), itsShiftBuffer(), itsSpectrumBuffer() { - itsShiftBuffer.resize(itsNChannel, nStation); - itsStationPhases.resize(nStation); + itsShiftBuffer.resize(itsNChannel, itsNBaseline); + itsPhases.resize(itsNBaseline); if (stokesIOnly) { itsSpectrumBuffer.resize(1, itsNChannel); } else { @@ -98,8 +98,8 @@ void Simulator::visit(const PointSource& component) { radec2lmn(itsReference, component.position(), lmn); // Compute station phase shifts. - phases(itsNStation, itsNChannel, lmn, itsStationUVW, itsFreq, itsShiftBuffer, - itsStationPhases); + phases(itsNBaseline, itsNChannel, lmn, itsUVW, itsFreq, itsShiftBuffer, + itsPhases); // Compute component spectrum. spectrum(component, itsNChannel, itsFreq, itsSpectrumBuffer, itsStokesIOnly); @@ -126,8 +126,7 @@ void Simulator::visit(const PointSource& component) { if (p == q) { buffer += itsNChannel * nCorr; } else { - const dcomplex* shiftP = &(itsShiftBuffer(0, p)); - const dcomplex* shiftQ = &(itsShiftBuffer(0, q)); + const dcomplex* shift = &(itsShiftBuffer(0, bl)); const dcomplex* spectrum = itsSpectrumBuffer.data(); double smearterm = 1.0; @@ -135,25 +134,23 @@ void Simulator::visit(const PointSource& component) { for (size_t ch = 0; ch < itsNChannel; ++ch) { if (itsCorrectFreqSmearing) { smearterm = - computeSmearterm(itsStationPhases[q] - itsStationPhases[p], + computeSmearterm(itsPhases[bl], itsChanWidths[ch] / 2); } // Compute baseline phase shift. // Compute visibilities. - *buffer++ += (*shiftQ) * conj(*shiftP) * (*spectrum++) * smearterm; - ++shiftP; - ++shiftQ; + *buffer++ += (*shift) * (*spectrum++) * smearterm; + ++shift; } // Channels. } else { for (size_t ch = 0; ch < itsNChannel; ++ch) { if (itsCorrectFreqSmearing) smearterm = - computeSmearterm(itsStationPhases[q] - itsStationPhases[p], + computeSmearterm(itsPhases[bl], itsChanWidths[ch] / 2); // Compute baseline phase shift. - const dcomplex blShift = (*shiftQ) * conj(*shiftP) * smearterm; - ++shiftP; - ++shiftQ; + const dcomplex blShift = (*shift) * smearterm; + ++shift; // Compute visibilities. *buffer++ += blShift * (*spectrum++); @@ -172,8 +169,8 @@ void Simulator::visit(const GaussianSource& component) { radec2lmn(itsReference, component.position(), lmn); // Compute station phase shifts. - phases(itsNStation, itsNChannel, lmn, itsStationUVW, itsFreq, itsShiftBuffer, - itsStationPhases); + phases(itsNBaseline, itsNChannel, lmn, itsUVW, itsFreq, itsShiftBuffer, + itsPhases); // Compute component spectrum. spectrum(component, itsNChannel, itsFreq, itsSpectrumBuffer, itsStokesIOnly); @@ -208,11 +205,8 @@ void Simulator::visit(const GaussianSource& component) { if (p == q) { buffer += itsNChannel * nCorr; } else { - double u = itsStationUVW(0, q); - double v = itsStationUVW(1, q); - - u -= itsStationUVW(0, p); - v -= itsStationUVW(1, p); + double u = itsUVW(0, bl); + double v = itsUVW(1, bl); // Rotate (u, v) by the position angle and scale with the major // and minor axis lengths (FWHM in rad). @@ -224,22 +218,20 @@ void Simulator::visit(const GaussianSource& component) { const double uvPrime = (-2.0 * casacore::C::pi * casacore::C::pi) * (uPrime * uPrime + vPrime * vPrime); - const dcomplex* shiftP = &(itsShiftBuffer(0, p)); - const dcomplex* shiftQ = &(itsShiftBuffer(0, q)); + const dcomplex* shift = &(itsShiftBuffer(0, bl)); const dcomplex* spectrum = itsSpectrumBuffer.data(); if (itsStokesIOnly) { for (size_t ch = 0; ch < itsNChannel; ++ch) { // Compute baseline phase shift. - dcomplex blShift = (*shiftQ) * conj(*shiftP); - ++shiftP; - ++shiftQ; + dcomplex blShift = (*shift); + ++shift; const double ampl = exp((itsFreq[ch] * itsFreq[ch]) / (casacore::C::c * casacore::C::c) * uvPrime); if (itsCorrectFreqSmearing) blShift *= - computeSmearterm(itsStationPhases[q] - itsStationPhases[p], + computeSmearterm(itsPhases[bl], itsChanWidths[ch] / 2); blShift *= ampl; @@ -250,16 +242,15 @@ void Simulator::visit(const GaussianSource& component) { } else { for (size_t ch = 0; ch < itsNChannel; ++ch) { // Compute baseline phase shift. - dcomplex blShift = (*shiftQ) * conj(*shiftP); - ++shiftP; - ++shiftQ; + dcomplex blShift = (*shift); + ++shift; const double ampl = exp((itsFreq[ch] * itsFreq[ch]) / (casacore::C::c * casacore::C::c) * uvPrime); if (itsCorrectFreqSmearing) blShift *= - computeSmearterm(itsStationPhases[q] - itsStationPhases[p], + computeSmearterm(itsPhases[bl], itsChanWidths[ch] / 2); blShift *= ampl; @@ -298,22 +289,22 @@ inline float computeSmearterm(double uvw, double halfwidth) { } // 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::Vector<double>& freq, - Simulator::Matrix<dcomplex>& shift, - std::vector<double>& stationPhases) { - dcomplex* shiftdata = shift.data(); + Simulator::Matrix<dcomplex>& shift_out, + std::vector<double>& phases_out) { + dcomplex* shiftdata = shift_out.data(); 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 * - (uvw(0, st) * lmn[0] + uvw(1, st) * lmn[1] + - uvw(2, st) * (lmn[2] - 1.0)); + (uvw(0, bl) * lmn[0] + uvw(1, bl) * lmn[1] + + uvw(2, bl) * (lmn[2] - 1.0)); - stationPhases[st] = phase; + phases_out[bl] = phase; 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; } // Channels. diff --git a/base/Simulator.h b/base/Simulator.h index bab285b593a2901b450fab1284845633a9cc1d35..0d23d20031af66630f8539f0dcdeceb5b7750962 100644 --- a/base/Simulator.h +++ b/base/Simulator.h @@ -87,10 +87,10 @@ class Simulator : public ModelComponentVisitor { const std::vector<Baseline> itsBaselines; const casacore::Vector<double> itsFreq; const casacore::Vector<double> itsChanWidths; - const casacore::Matrix<double> itsStationUVW; + const casacore::Matrix<double> itsUVW; casacore::Cube<dcomplex> itsBuffer; Matrix<dcomplex> itsShiftBuffer; - std::vector<double> itsStationPhases; + std::vector<double> itsPhases; Matrix<dcomplex> itsSpectrumBuffer; }; diff --git a/base/test/unit/tSimulator.cc b/base/test/unit/tSimulator.cc index 87d0f97bd93fc0b584ac6793f99ee7e3abec119a..68d806f5844dc7fd6660f5ee85f3a4f8bcdf103a 100644 --- a/base/test/unit/tSimulator.cc +++ b/base/test/unit/tSimulator.cc @@ -38,12 +38,16 @@ Simulator MakeSimulator(bool correct_freq_smearing, bool stokes_i_only, } 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) { - uvw(0, st) = st * 5000; - uvw(1, st) = st * 1000; - uvw(2, st) = 0; + casacore::Matrix<double> uvw(3, nBaseline); + + for (size_t bl = 0; bl < nBaseline; ++bl) { + 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, diff --git a/steps/OnePredict.cc b/steps/OnePredict.cc index 562b3768c775aa923754290d65b82236945f037b..af5fca30463825ab42381231242749328ace906c 100644 --- a/steps/OnePredict.cc +++ b/steps/OnePredict.cc @@ -191,9 +191,9 @@ void OnePredict::initializeThreadData() { const size_t nCr = itsStokesIOnly ? 1 : info().ncorr(); const size_t nThreads = getInfo().nThreads(); - itsStationUVW.resize(3, nSt); - itsUVWSplitIndex = base::nsetupSplitUVW(info().nantenna(), info().getAnt1(), - info().getAnt2()); + // itsStationUVW.resize(3, nSt); + // itsUVWSplitIndex = base::nsetupSplitUVW(info().nantenna(), info().getAnt1(), + // info().getAnt2()); if (!itsPredictBuffer) { itsPredictBuffer = std::make_shared<base::PredictBuffer>(); @@ -316,8 +316,8 @@ bool OnePredict::process(const DPBuffer& bufin) { itsTimerPredict.start(); - base::nsplitUVW(itsUVWSplitIndex, itsBaselines, tempBuffer.getUVW(), - itsStationUVW); + // base::nsplitUVW(itsUVWSplitIndex, itsBaselines, tempBuffer.getUVW(), + // itsStationUVW); double time = tempBuffer.getTime(); // Set up directions for beam evaluation @@ -374,7 +374,7 @@ bool OnePredict::process(const DPBuffer& bufin) { (itsApplyBeam ? itsPredictBuffer->GetPatchModel(thread) : itsPredictBuffer->GetModel(thread)); simulators.emplace_back(itsPhaseRef, nSt, itsBaselines, info().chanFreqs(), - info().chanWidths(), itsStationUVW, simulatedest, + info().chanWidths(), tempBuffer.getUVW(), simulatedest, itsCorrectFreqSmearing, itsStokesIOnly); } std::vector<base::Patch::ConstPtr> curPatches(pool->NThreads());