diff --git a/.gitattributes b/.gitattributes index 68a9983f6a4e5395cab163e250055976209bedab..48131e4da9edc112c70ba21e788cad017de8dfdf 100644 --- a/.gitattributes +++ b/.gitattributes @@ -876,6 +876,8 @@ RTCP/CNProc/bootstrap -text RTCP/CNProc/src/BeamFormer.cc -text RTCP/CNProc/src/BeamFormer.h -text RTCP/CNProc/src/CN_Processing.machinefile -text +RTCP/CNProc/src/PencilBeams.cc -text +RTCP/CNProc/src/PencilBeams.h -text RTCP/FCNP/bootstrap -text RTCP/IONProc/bootstrap -text RTCP/Interface/bootstrap -text diff --git a/RTCP/CNProc/src/CN_Processing.cc b/RTCP/CNProc/src/CN_Processing.cc index 9a03283ec066a64903f06023e34a1d49c96aa1c6..974de0830c627a86253e33c14ff95725026f343d 100644 --- a/RTCP/CNProc/src/CN_Processing.cc +++ b/RTCP/CNProc/src/CN_Processing.cc @@ -97,6 +97,7 @@ template <typename SAMPLE_TYPE> CN_Processing<SAMPLE_TYPE>::CN_Processing(Stream #endif itsPPF(0), itsBeamFormer(0), + itsPencilBeamFormer(0), itsCorrelator(0) { @@ -307,6 +308,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::preprocess(CN_C if (itsIsTransposeOutput) { unsigned nrSubbandsPerPset = configuration.nrSubbandsPerPset(); unsigned logicalNode = usedCoresPerPset * (outputPsetIndex - outputPsets.begin()) + myCore; + PencilRings pencilCoordinates( configuration.nrPencilRings(), configuration.pencilRingSize() ); // TODO: logicalNode assumes output psets are consecutively numbered itsCenterFrequencies = configuration.refFreqs(); @@ -325,6 +327,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::preprocess(CN_C itsPPF = new PPF<SAMPLE_TYPE>(itsNrStations, nrChannels, nrSamplesPerIntegration, configuration.sampleRate() / nrChannels, configuration.delayCompensation()); + itsPencilBeamFormer = new PencilBeams(pencilCoordinates, itsNrStations, nrChannels, nrSamplesPerIntegration, itsCenterFrequencies[itsCurrentSubband], configuration.sampleRate() / nrChannels, configuration.refPhaseCentre(), configuration.phaseCentres()); itsCorrelator = new Correlator(nrBeamFormedStations, itsBeamFormer->getStationMapping(), nrChannels, nrSamplesPerIntegration, configuration.correctBandPass()); @@ -396,6 +399,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::process() #endif } // itsIsTransposeInput + #if defined HAVE_MPI if(!itsDoAsyncCommunication) { if (itsIsTransposeInput || itsIsTransposeOutput) { @@ -453,6 +457,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::process() computeTimer.start(); itsBeamFormer->formBeams(itsFilteredData); + //itsPencilBeamFormer->formPencilBeams(itsFilteredData); itsCorrelator->computeFlagsAndCentroids(itsFilteredData, itsCorrelatedData); itsCorrelator->correlate(itsFilteredData, itsCorrelatedData); computeTimer.stop(); @@ -523,6 +528,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::postprocess() delete itsPPF; delete itsFilteredData; delete itsBeamFormer; + delete itsPencilBeamFormer; delete itsCorrelator; delete itsCorrelatedData; } diff --git a/RTCP/CNProc/src/CN_Processing.h b/RTCP/CNProc/src/CN_Processing.h index f0acb4abbe5023ab5b2951e4b2b6efb093a9302b..e01c4c015e136e5b18b27549c374633f06afabd5 100644 --- a/RTCP/CNProc/src/CN_Processing.h +++ b/RTCP/CNProc/src/CN_Processing.h @@ -45,6 +45,7 @@ #include <PPF.h> #include <Correlator.h> #include <BeamFormer.h> +#include <PencilBeams.h> #include <LocationInfo.h> @@ -120,6 +121,7 @@ template <typename SAMPLE_TYPE> class CN_Processing : public CN_Processing_Base #endif PPF<SAMPLE_TYPE> *itsPPF; BeamFormer *itsBeamFormer; + PencilBeams *itsPencilBeamFormer; Correlator *itsCorrelator; #if defined HAVE_BGL diff --git a/RTCP/CNProc/src/FilteredData.h b/RTCP/CNProc/src/FilteredData.h index 7823228cf62e61a89d91d0e7a1edbf760de50016..514ab3b2b83e9d32b7a0f9453a9c0675a5af802a 100644 --- a/RTCP/CNProc/src/FilteredData.h +++ b/RTCP/CNProc/src/FilteredData.h @@ -6,6 +6,7 @@ #include <Interface/Config.h> #include <Interface/MultiDimArray.h> #include <Interface/SparseSet.h> +#include <Interface/SubbandMetaData.h> namespace LOFAR { @@ -24,19 +25,22 @@ class FilteredData // is not used. MultiDimArray<fcomplex, 4> samples; //[itsNrChannels][itsNrStations][itsNrSamplesPerIntegration | 2][NR_POLARIZATIONS] CACHE_ALIGNED SparseSet<unsigned> *flags; //[itsNrStations] + Vector<SubbandMetaData> metaData; //[itsNrStations] }; inline size_t FilteredData::requiredSize(unsigned nrStations, unsigned nrChannels, unsigned nrSamplesPerIntegration) { - return align(sizeof(fcomplex) * nrChannels * nrStations * (nrSamplesPerIntegration | 2) * NR_POLARIZATIONS, 32); + return align(sizeof(fcomplex) * nrChannels * nrStations * (nrSamplesPerIntegration | 2) * NR_POLARIZATIONS, 32) + + align(sizeof(SubbandMetaData) * nrStations, 32); } inline FilteredData::FilteredData(unsigned nrStations, unsigned nrChannels, unsigned nrSamplesPerIntegration, Allocator &allocator) : samples(boost::extents[nrChannels][nrStations][nrSamplesPerIntegration | 2][NR_POLARIZATIONS], 32, allocator), - flags(new SparseSet<unsigned>[nrStations]) + flags(new SparseSet<unsigned>[nrStations]), + metaData(nrStations, 32, allocator) { } diff --git a/RTCP/CNProc/src/Makefile.am b/RTCP/CNProc/src/Makefile.am index 90a289e5d66eea84f8e1efcaafac3081f7e8b85c..7f4d7324f9101feda9713f8bfc978d9c7a049a1b 100644 --- a/RTCP/CNProc/src/Makefile.am +++ b/RTCP/CNProc/src/Makefile.am @@ -16,7 +16,8 @@ BeamFormer.h \ Correlator.h \ CN_Processing.h \ FCNP_ClientStream.h \ -AsyncCommunication.h +AsyncCommunication.h \ +PencilBeams.h NOINSTHDRS = @@ -50,7 +51,8 @@ Correlator.cc \ CN_Processing_main.cc \ CN_Processing.cc \ FCNP_ClientStream.cc \ -AsyncCommunication.cc +AsyncCommunication.cc \ +PencilBeams.cc configfilesdir=$(bindir) configfiles_DATA = \ diff --git a/RTCP/CNProc/src/PPF.cc b/RTCP/CNProc/src/PPF.cc index 44252bbbdac278bee1a9c7a6c94d6ae19e7e36ac..c26f9e1ab09992af41d966426b0faff70b6c00d6 100644 --- a/RTCP/CNProc/src/PPF.cc +++ b/RTCP/CNProc/src/PPF.cc @@ -337,6 +337,13 @@ template <typename SAMPLE_TYPE> void PPF<SAMPLE_TYPE>::filter(unsigned stat, dou if (itsDelayCompensation) { computePhaseShifts(phaseShifts, transposedData->metaData[stat].delayAtBegin, transposedData->metaData[stat].delayAfterEnd, baseFrequency); + + // forward (pencil) beam forming information + for( unsigned i = 0; i < 3; i++ ) { + filteredData->metaData[stat].beamDirectionAtBegin[i] = transposedData->metaData[stat].beamDirectionAtBegin[i]; + filteredData->metaData[stat].beamDirectionAfterEnd[i] = transposedData->metaData[stat].beamDirectionAfterEnd[i]; + } + } const SparseSet<unsigned>::Ranges &ranges = filteredData->flags[stat].getRanges(); diff --git a/RTCP/CNProc/src/PencilBeams.cc b/RTCP/CNProc/src/PencilBeams.cc new file mode 100644 index 0000000000000000000000000000000000000000..9072de04e9f46cf3d67f19d0ed5150954836d19b --- /dev/null +++ b/RTCP/CNProc/src/PencilBeams.cc @@ -0,0 +1,312 @@ +//# Always #include <lofar_config.h> first! +#include <lofar_config.h> + +#include <PencilBeams.h> +#include <Interface/MultiDimArray.h> +#include <iostream> +#include <cmath> +#include <vector> + +#ifndef M_SQRT3 + #define M_SQRT3 1.73205080756887719000 +#endif + +namespace LOFAR { +namespace RTCP { + +PencilRings::PencilRings(unsigned nrRings, double ringWidth): + itsNrRings(nrRings), + itsRingWidth(ringWidth) +{ + computeBeamCoordinates(); +} + +unsigned PencilRings::nrPencils() const +{ + // the centered hexagonal number + return 3 * itsNrRings * (itsNrRings + 1) + 1; +} + +double PencilRings::pencilEdge() const +{ + return itsRingWidth / M_SQRT3; +} + +double PencilRings::pencilWidth() const +{ + // _ // + // / \ // + // \_/ // + //|...| // + return 2.0 * pencilEdge(); +} + +double PencilRings::pencilHeight() const +{ + // _ _ // + // / \ : // + // \_/ _ // + // // + return itsRingWidth; +} + +double PencilRings::pencilWidthDelta() const +{ + // _ // + // / \_ // + // \_/ \ // + // \_/ // + // |.| // + return 1.5 * pencilEdge(); +} + +double PencilRings::pencilHeightDelta() const +{ + // _ // + // / \_ - // + // \_/ \ - // + // \_/ // + return 0.5 * itsRingWidth; +} + +void PencilRings::computeBeamCoordinates() +{ + std::vector<PencilCoord3D> &coordinates = getCoordinates(); + double dl[6], dm[6]; + + // stride for each side, starting from the top, clock-wise + + // _ // + // / \_ // + // \_/ \ // + // \_/ // + dl[0] = pencilWidthDelta(); + dm[0] = -pencilHeightDelta(); + + // _ // + // / \ // + // \_/ // + // / \ // + // \_/ // + dl[1] = 0; + dm[1] = -pencilHeight(); + + // _ // + // _/ \ // + // / \_/ // + // \_/ // + dl[2] = -pencilWidthDelta(); + dm[2] = -pencilHeightDelta(); + + // _ // + // / \_ // + // \_/ \ // + // \_/ // + dl[3] = -pencilWidthDelta(); + dm[3] = pencilHeightDelta(); + + // _ // + // / \ // + // \_/ // + // / \ // + // \_/ // + dl[4] = 0; + dm[4] = pencilHeight(); + + // _ // + // _/ \ // + // / \_/ // + // \_/ // + dl[5] = pencilWidthDelta(); + dm[5] = pencilHeightDelta(); + + // ring 0: the center pencil beam + coordinates.reserve(nrPencils()); + coordinates.push_back( PencilCoord3D( 0, 0 ) ); + + // ring 1-n: create the pencil beams from the inner ring outwards + for( unsigned ring = 1; ring <= itsNrRings; ring++ ) { + // start from the top + double l = 0; + double m = pencilHeight() * ring; + + for( unsigned side = 0; side < 6; side++ ) { + for( unsigned pencil = 0; pencil < ring; pencil++ ) { + coordinates.push_back( PencilCoord3D( l, m ) ); + l += dl[side]; m += dm[side]; + } + } + } +} + +PencilBeams::PencilBeams(PencilCoordinates coordinates, unsigned nrStations, unsigned nrChannels, unsigned nrSamplesPerIntegration, double centerFrequency, double channelBandwidth, std::vector<double> refPhaseCentre, Matrix<double> phaseCentres ) +: + itsCoordinates(coordinates.getCoordinates()), + itsNrStations(nrStations), + itsNrChannels(nrChannels), + itsNrSamplesPerIntegration(nrSamplesPerIntegration), + itsCenterFrequency(centerFrequency), + itsChannelBandwidth(channelBandwidth), + itsRefPhaseCentre(refPhaseCentre), + itsPencilBeamData( boost::extents[1][1][1][1], 32 ) + //itsPencilBeamData( boost::extents[nrChannels][coordinates.getCoordinates().size()][itsNrSamplesPerIntegration | 2][NR_POLARIZATIONS], 32 ) +{ + // derived constants + itsBaseFrequency = centerFrequency - (itsNrChannels/2) * channelBandwidth; + + // copy all phase centres and their derived constants + itsPhaseCentres.reserve( nrStations ); + itsBaselines.reserve( nrStations ); + itsBaselinesSeconds.reserve( nrStations ); + itsDelays.resize( nrStations, itsCoordinates.size() ); + itsDelayOffsets.resize( nrStations, itsCoordinates.size() ); + for( unsigned stat = 0; stat < nrStations; stat++ ) { + double x = phaseCentres[stat][0]; + double y = phaseCentres[stat][1]; + double z = phaseCentres[stat][2]; + + PencilCoord3D phaseCentre( x, y, z ); + PencilCoord3D baseLine = phaseCentre - refPhaseCentre; + PencilCoord3D baseLineSeconds = baseLine * (1.0/speedOfLight); + + itsPhaseCentres.push_back( phaseCentre ); + itsBaselines.push_back( baseLine ); + itsBaselinesSeconds.push_back( baseLineSeconds ); + + for( unsigned beam = 0; beam < itsCoordinates.size(); beam++ ) { + itsDelayOffsets[stat][beam] = baseLine * itsCoordinates[beam] * (1.0/speedOfLight); + } + } +} + +void PencilBeams::calculateDelays( unsigned stat, const PencilCoord3D &beamDir ) +{ + double compensatedDelay = itsDelayOffsets[stat][0] - itsBaselinesSeconds[stat] * beamDir; + + // centre beam does not need compensation + itsDelays[stat][0] = 0.0; + + for( unsigned i = 1; i < itsCoordinates.size(); i++ ) { + // delay[i] = baseline * (coordinate - beamdir) / c + // = (baseline * coordinate / c) - (baseline * beamdir) / c + // = delayoffset - baselinesec * beamdir + // + // further reduced by the delay we already compensate for when doing regular beam forming (the centre of the beam) + itsDelays[stat][i] = itsDelayOffsets[stat][i] - itsBaselinesSeconds[stat] * beamDir - compensatedDelay; + } +} + +fcomplex PencilBeams::phaseShift( double frequency, double delay ) +{ + double phaseShift = delay * frequency; + double phi = -2 * M_PI * phaseShift; + + return makefcomplex( std::sin(phi), std::cos(phi) ); +} + +void PencilBeams::formPartialCenterBeam( MultiDimArray<fcomplex,4> &samples, std::vector<unsigned> stations, bool add, float factor ) +{ + // the center beam has a delay of 0 for all stations, since we already + // corrected for it + + for (unsigned ch = 0; ch < itsNrChannels; ch ++) { + for (unsigned time = 0; time < itsNrSamplesPerIntegration; time ++) { + for (unsigned pol = 0; pol < NR_POLARIZATIONS; pol ++) { + //fcomplex &dest = itsPencilBeamData[ch][0][time][pol]; + fcomplex &dest = itsPencilBeamData[0][0][0][0]; + fcomplex sample = makefcomplex( 0, 0 ); + + for( unsigned stat = 0; stat < stations.size(); stat++ ) { + sample += samples[ch][stations[stat]][time][pol]; + } + + dest = (add ? dest : 0) + sample * factor; + } + } + } +} + +void PencilBeams::formPartialBeams( MultiDimArray<fcomplex,4> &samples, std::vector<unsigned> stations, std::vector<unsigned> beams, bool add, float factor ) +{ + double frequency = itsBaseFrequency; + + for (unsigned ch = 0; ch < itsNrChannels; ch ++) { + for (unsigned time = 0; time < itsNrSamplesPerIntegration; time ++) { + for (unsigned pol = 0; pol < NR_POLARIZATIONS; pol ++) { + for( unsigned beam = 0; beam < beams.size(); beam++ ) { + //fcomplex &dest = itsPencilBeamData[ch][beams[beam]][time][pol]; + fcomplex &dest = itsPencilBeamData[0][0][0][0]; + fcomplex sample = makefcomplex( 0, 0 ); + + for( unsigned stat = 0; stat < stations.size(); stat++ ) { + fcomplex shift = phaseShift( frequency, itsDelays[stations[stat]][beams[beam]] ); + sample += samples[ch][stations[stat]][time][pol] * shift; + } + + dest = (add ? dest : 0) + sample * factor; + } + } + } + + frequency += itsChannelBandwidth; + } +} + +std::ostream& operator<<(std::ostream &str, std::vector<unsigned> &v) +{ + for(unsigned i = 0; i < v.size(); i++ ) { + if( i > 0 ) str << ", "; + str << v[i]; + } + + return str; +} + +void PencilBeams::formPencilBeams( FilteredData *filteredData ) +{ + float factor = 1.0/itsNrStations; + + // TODO: fetch a list of stations to beam form. for now, we assume + // we use all stations + + // calculate the delays for each station for this integration period + for( unsigned stat = 0; stat < itsNrStations; stat++ ) { + // todo: interpolate per time step? + PencilCoord3D beamDirBegin = filteredData->metaData[stat].beamDirectionAtBegin; + PencilCoord3D beamDirEnd = filteredData->metaData[stat].beamDirectionAfterEnd; + PencilCoord3D beamDirAverage = (beamDirBegin + beamDirEnd) * 0.5; + + calculateDelays( stat, beamDirAverage ); + } + + // combine 6 stations at a time + for( unsigned stat = 0; stat < itsNrStations; stat += 6 ) { + std::vector<unsigned> stations; + + stations.reserve( 6 ); + for( unsigned i = 0; i < 6 && stat+i < itsNrStations; i++ ) { + stations.push_back( stat+i ); + } + + // form the central pencil beam (beam #0) + std::clog << "Forming central beam for stations " << stations << std::endl; + formPartialCenterBeam( filteredData->samples, stations, stat > 0, factor ); + + // form the other beams, 6 at a time + for( unsigned beam = 1; beam < itsCoordinates.size(); beam += 6 ) { + std::vector<unsigned> beams; + + for( unsigned i = 0; i < 6 && beam+i < itsCoordinates.size(); i++ ) { + beams.push_back( beam+i ); + } + + std::clog << "Forming partial beams " << beams << " for stations " << stations << std::endl; + formPartialBeams( filteredData->samples, stations, beams, stat > 0, factor ); + } + } +} + +} // namespace RTCP +} // namespace LOFAR + + diff --git a/RTCP/CNProc/src/PencilBeams.h b/RTCP/CNProc/src/PencilBeams.h new file mode 100644 index 0000000000000000000000000000000000000000..744597ff296a4cfc9104e549a5ce958df6acdf51 --- /dev/null +++ b/RTCP/CNProc/src/PencilBeams.h @@ -0,0 +1,201 @@ +#ifndef LOFAR_CNPROC_PENCIL_BEAMS_H +#define LOFAR_CNPROC_PENCIL_BEAMS_H + +#include <vector> +#include <cmath> + +#include <FilteredData.h> + +namespace LOFAR { +namespace RTCP { + +const double speedOfLight = 299792458; + +class PencilCoord3D { + public: + PencilCoord3D(double x, double y) { + itsXYZ[0] = x; + itsXYZ[1] = y; + itsXYZ[2] = sqrt( 1.0 - x*x - y*y ); + } + + PencilCoord3D(double x, double y, double z) { + itsXYZ[0] = x; + itsXYZ[1] = y; + itsXYZ[2] = z; + } + + PencilCoord3D(double xyz[3]) { + itsXYZ[0] = xyz[0]; + itsXYZ[1] = xyz[1]; + itsXYZ[2] = xyz[2]; + } + + PencilCoord3D(std::vector<double> xyz) { + itsXYZ[0] = xyz[0]; + itsXYZ[1] = xyz[1]; + itsXYZ[2] = xyz[2]; + } + + inline PencilCoord3D& operator-=( const PencilCoord3D &rhs ) + { + itsXYZ[0] -= rhs.itsXYZ[0]; + itsXYZ[1] -= rhs.itsXYZ[1]; + itsXYZ[2] -= rhs.itsXYZ[2]; + + return *this; + } + + inline PencilCoord3D& operator+=( const PencilCoord3D &rhs ) + { + itsXYZ[0] += rhs.itsXYZ[0]; + itsXYZ[1] += rhs.itsXYZ[1]; + itsXYZ[2] += rhs.itsXYZ[2]; + + return *this; + } + + inline PencilCoord3D& operator*=( const double a ) + { + itsXYZ[0] *= a; + itsXYZ[1] *= a; + itsXYZ[2] *= a; + + return *this; + } + + friend double operator*( const PencilCoord3D &lhs, const PencilCoord3D &rhs ); + friend ostream& operator<<(ostream& os, const PencilCoord3D &c); + + private: + double itsXYZ[3]; +}; + + +// PencilCoordinates are coordinates of the pencil beams that need to +// be formed. Each coordinate is a normalised vector, relative to the +// center beam. +// +// The center beam has to be included as the first coordinate of (0,0,1). +class PencilCoordinates +{ + public: + PencilCoordinates() {} + PencilCoordinates( std::vector<PencilCoord3D> coordinates ): itsCoordinates(coordinates) {} + + std::vector<PencilCoord3D>& getCoordinates() + { return itsCoordinates; } + + private: + std::vector<PencilCoord3D> itsCoordinates; +}; + +// PencilRings are rings of regular hexagons around a center beam: +// _ +// _/ \_ // +// / \_/ \ // +// \_/ \_/ // +// / \_/ \ // +// \_/ \_/ // +// \_/ // +// +// the width of each ring is defined as the distance between pencil centers +// +// assumes l horizontal (incr left), m vertical (incr up) +class PencilRings: public PencilCoordinates +{ + public: + PencilRings(unsigned nrRings, double ringWidth); + + unsigned nrPencils() const; + + // length of a pencil edge + double pencilEdge() const; + + // width of a pencil beam + double pencilWidth() const; + + // height of a pencil beam + double pencilHeight() const; + + // horizontal difference for neighbouring beams, in case of partial overlap + double pencilWidthDelta() const; + + // vertical difference for neighbouring beams, in case of partial overlap + double pencilHeightDelta() const; + + private: + void computeBeamCoordinates(); + + unsigned itsNrRings; + double itsRingWidth; +}; + +class PencilBeams +{ + public: + PencilBeams(PencilCoordinates coordinates, unsigned nrStations, unsigned nrChannels, unsigned nrSamplesPerIntegration, double centerFrequency, double channelBandwidth, std::vector<double> refPhaseCentre, Matrix<double> phaseCentres ); + + void formPencilBeams( FilteredData *filteredData ); + + private: + fcomplex phaseShift( double frequency, double delay ); + void formPartialCenterBeam( MultiDimArray<fcomplex,4> &samples, std::vector<unsigned> stations, bool add, float factor ); + void formPartialBeams( MultiDimArray<fcomplex,4> &samples, std::vector<unsigned> stations, std::vector<unsigned> beams, bool add, float factor ); + void calculateDelays( unsigned stat, const PencilCoord3D &beamDir ); + + MultiDimArray<fcomplex,4> itsPencilBeamData; // [nrChannels][itsCoordinates.size()][nrSamplesPerIntegration | 2][NR_POLARIZATIONS] // has a size of [1][1][1][1] for now, until a better storage method is devised + std::vector<PencilCoord3D> itsCoordinates; + unsigned itsNrStations; + unsigned itsNrChannels; + unsigned itsNrSamplesPerIntegration; + double itsCenterFrequency; + double itsChannelBandwidth; + double itsBaseFrequency; + Matrix<double> itsDelays; // [itsNrStations][itsCoordinates.size()] + Matrix<double> itsDelayOffsets; // [itsNrStations][itsCoordinates.size()] + PencilCoord3D itsRefPhaseCentre; + std::vector<PencilCoord3D> itsPhaseCentres; + std::vector<PencilCoord3D> itsBaselines; + std::vector<PencilCoord3D> itsBaselinesSeconds; +}; + +inline double operator*( const PencilCoord3D &lhs, const PencilCoord3D &rhs ) +{ + double sum = 0; + + sum += lhs.itsXYZ[0] * rhs.itsXYZ[0]; + sum += lhs.itsXYZ[1] * rhs.itsXYZ[1]; + sum += lhs.itsXYZ[2] * rhs.itsXYZ[2]; + return sum; +} + +inline PencilCoord3D& operator-( const PencilCoord3D &lhs, const PencilCoord3D &rhs ) +{ + return PencilCoord3D( lhs ) -= rhs; +} + +inline PencilCoord3D& operator+( const PencilCoord3D &lhs, const PencilCoord3D &rhs ) +{ + return PencilCoord3D( lhs ) += rhs; +} + +inline PencilCoord3D& operator*( const double a, const PencilCoord3D &rhs ) +{ + return PencilCoord3D( rhs ) *= a; +} + +inline PencilCoord3D& operator*( const PencilCoord3D &lhs, const double a ) +{ + return PencilCoord3D( lhs ) *= a; +} + +inline ostream& operator<<(ostream& os, const PencilCoord3D &c) +{ + return os << "(" << c.itsXYZ[0] << "," << c.itsXYZ[1] << "," << c.itsXYZ[2] << ")"; +} + +} // namespace RTCP +} // namespace LOFAR + +#endif diff --git a/RTCP/IONProc/src/ION_main.cc b/RTCP/IONProc/src/ION_main.cc index 24d13cfe6605b1776bf343c556e3ddd7a0bda228..4021b08bfa5137ba8feb9a1a5ca061c2972dd3cc 100644 --- a/RTCP/IONProc/src/ION_main.cc +++ b/RTCP/IONProc/src/ION_main.cc @@ -141,6 +141,9 @@ static void configureCNs(const Parset &parset) { CN_Command command(CN_Command::PREPROCESS); CN_Configuration configuration; + std::vector<Parset::StationRSPpair> inputs = parset.getStationNamesAndRSPboardNumbers(myPsetNumber); + Matrix<double> &phaseCentres = configuration.phaseCentres(); + configuration.nrStations() = parset.nrStations(); configuration.nrBitsPerSample() = parset.nrBitsPerSample(); @@ -156,6 +159,18 @@ static void configureCNs(const Parset &parset) configuration.outputPsets() = parset.getUint32Vector("OLAP.CNProc.outputPsets"); configuration.tabList() = parset.getUint32Vector("OLAP.CNProc.tabList"); configuration.refFreqs() = parset.subbandToFrequencyMapping(); + configuration.nrPencilRings() = parset.getUint32("OLAP.DelayComp.nrPencilRings"); + configuration.pencilRingSize() = parset.getDouble("OLAP.DelayComp.pencilRingSize"); + configuration.refPhaseCentre() = parset.getRefPhaseCentres(); + + phaseCentres.resize( inputs.size(), 3 ); + for( unsigned stat = 0; stat < inputs.size(); stat++ ) { + std::vector<double> phaseCentre = parset.getPhaseCentresOf( inputs[stat].station ); + + for( unsigned dim = 0; dim < 3; dim++ ) { + phaseCentres[stat][dim] = phaseCentre[dim]; + } + } std::clog << "configuring " << nrCoresPerPset << " cores ..."; std::clog.flush(); diff --git a/RTCP/IONProc/src/InputSection.cc b/RTCP/IONProc/src/InputSection.cc index 54f96cb341ed83184a40f9c1f1e04d5518dac2ec..25386af060ae0efe462a15d14ba1d401d4acc7d7 100644 --- a/RTCP/IONProc/src/InputSection.cc +++ b/RTCP/IONProc/src/InputSection.cc @@ -173,23 +173,29 @@ template<typename SAMPLE_TYPE> void InputSection<SAMPLE_TYPE>::preprocess(const if (itsDelayCompensation) { itsDelaysAtBegin.resize(itsNrBeams); itsDelaysAfterEnd.resize(itsNrBeams); - itsNrCalcDelaysForEachTimeNrDirections.resize(itsNrCalcDelays * itsNrBeams); + itsBeamDirectionsAtBegin.resize(itsNrBeams); + itsBeamDirectionsAfterEnd.resize(itsNrBeams); + itsAllBeamDirections.resize(itsNrCalcDelays * itsNrBeams); itsDelayComp = new WH_DelayCompensation(ps, inputs[0].station); // TODO: support more than one station std::vector<double> startTimes(itsNrCalcDelays); - for (uint i = 0; i < itsNrCalcDelays; i ++) + for (uint i = 0; i < itsNrCalcDelays; i ++) { startTimes[i] = static_cast<int64>(itsSyncedStamp + i * itsNSamplesPerSec) * itsSampleDuration; + } - itsNrCalcDelaysForEachTimeNrDirections = itsDelayComp->calcDelaysForEachTimeNrDirections(startTimes); + itsAllBeamDirections = itsDelayComp->calculateAllBeamDirections(startTimes); itsCounter = 0; - for (unsigned beam = 0; beam < itsNrBeams; beam ++) - itsDelaysAfterEnd[beam] = itsNrCalcDelaysForEachTimeNrDirections[beam]; - + for (unsigned beam = 0; beam < itsNrBeams; beam ++) { + itsBeamDirectionsAfterEnd[beam] = itsAllBeamDirections[beam]; + itsDelaysAfterEnd[beam] = itsDelayComp->getDelay( itsAllBeamDirections[beam] ); + } + itsDelaysAtBegin = itsDelaysAfterEnd; + itsBeamDirectionsAtBegin = itsBeamDirectionsAfterEnd; } itsIsRealTime = ps->realTime(); @@ -236,6 +242,7 @@ template<typename SAMPLE_TYPE> void InputSection<SAMPLE_TYPE>::computeDelays() if (itsDelayCompensation) { itsCounter ++; itsDelaysAtBegin = itsDelaysAfterEnd; // from previous integration + itsBeamDirectionsAtBegin = itsBeamDirectionsAfterEnd; // from previous integration if (itsCounter > itsNrCalcDelays - 1) { std::vector<double> stopTimes(itsNrCalcDelays); @@ -244,18 +251,22 @@ template<typename SAMPLE_TYPE> void InputSection<SAMPLE_TYPE>::computeDelays() stopTimes[i] = static_cast<int64>(itsSyncedStamp + (i + 1) * itsNSamplesPerSec) * itsSampleDuration; itsDelayTimer.start(); - itsNrCalcDelaysForEachTimeNrDirections = itsDelayComp->calcDelaysForEachTimeNrDirections(stopTimes); + itsAllBeamDirections = itsDelayComp->calculateAllBeamDirections(stopTimes); itsDelayTimer.stop(); itsCounter = 0; for (unsigned beam = 0; beam < itsNrBeams; beam ++) { - itsDelaysAfterEnd[beam] = itsNrCalcDelaysForEachTimeNrDirections[beam]; + itsBeamDirectionsAfterEnd[beam] = itsAllBeamDirections[beam]; + itsDelaysAfterEnd[beam] = itsDelayComp->getDelay( itsAllBeamDirections[beam] ); } } else { unsigned firstBeam = itsCounter * itsNrBeams; for (unsigned beam = 0; beam < itsNrBeams; beam ++) { - itsDelaysAfterEnd[beam] = itsNrCalcDelaysForEachTimeNrDirections[firstBeam++]; + itsBeamDirectionsAfterEnd[beam] = itsAllBeamDirections[firstBeam]; + itsDelaysAfterEnd[beam] = itsDelayComp->getDelay( itsAllBeamDirections[firstBeam] ); + + firstBeam++; } } @@ -340,15 +351,26 @@ template<typename SAMPLE_TYPE> void InputSection<SAMPLE_TYPE>::toComputeNode() // create and send all metadata in one "large" message std::vector<SubbandMetaData, AlignedStdAllocator<SubbandMetaData, 16> > metaData(itsNrPsets); - for (unsigned pset = 0; pset < itsNrPsets; pset ++) { - unsigned subband = itsNSubbandsPerPset * pset + subbandBase; - unsigned rspBoard = itsSubbandToRSPboardMapping[subband]; - unsigned beam = itsSubbandToBeamMapping[subband]; + if( itsDelayCompensation ) { + for (unsigned pset = 0; pset < itsNrPsets; pset ++) { + unsigned subband = itsNSubbandsPerPset * pset + subbandBase; + unsigned rspBoard = itsSubbandToRSPboardMapping[subband]; + unsigned beam = itsSubbandToBeamMapping[subband]; + const vector<double> &beamDirBegin = itsBeamDirectionsAtBegin[beam].coord().get(); + const vector<double> &beamDirEnd = itsBeamDirectionsAfterEnd[beam].coord().get(); + + metaData[pset].delayAtBegin = itsFineDelaysAtBegin[beam]; + metaData[pset].delayAfterEnd = itsFineDelaysAfterEnd[beam]; - metaData[pset].delayAtBegin = itsFineDelaysAtBegin[beam]; - metaData[pset].delayAfterEnd = itsFineDelaysAfterEnd[beam]; - metaData[pset].alignmentShift = itsBBuffers[rspBoard]->alignmentShift(beam); - metaData[pset].setFlags(itsFlags[rspBoard][beam]); + for( unsigned i = 0; i < 3; i++ ) { + metaData[pset].beamDirectionAtBegin[i] = beamDirBegin[i]; + metaData[pset].beamDirectionAfterEnd[i] = beamDirEnd[i]; + + } + + metaData[pset].alignmentShift = itsBBuffers[rspBoard]->alignmentShift(beam); + metaData[pset].setFlags(itsFlags[rspBoard][beam]); + } } stream->write(&metaData[0], metaData.size() * sizeof(SubbandMetaData)); diff --git a/RTCP/IONProc/src/InputSection.h b/RTCP/IONProc/src/InputSection.h index efb7a778ea1d262b3701448a5a152256d59b5f51..291724dba0f08088bedf56afa83b99c91f469ac9 100644 --- a/RTCP/IONProc/src/InputSection.h +++ b/RTCP/IONProc/src/InputSection.h @@ -36,6 +36,7 @@ #include <WH_DelayCompensation.h> #include <InputThread.h> #include <LogThread.h> +#include <AMCBase/Direction.h> #include <boost/multi_array.hpp> #include <pthread.h> @@ -86,7 +87,9 @@ template <typename SAMPLE_TYPE> class InputSection { std::vector<double> itsDelaysAtBegin; std::vector<double> itsDelaysAfterEnd; - std::vector<double> itsNrCalcDelaysForEachTimeNrDirections; + std::vector<AMC::Direction> itsBeamDirectionsAtBegin; + std::vector<AMC::Direction> itsBeamDirectionsAfterEnd; + std::vector<AMC::Direction> itsAllBeamDirections; unsigned itsNrCalcDelays; unsigned itsCounter; unsigned itsNrPsets; diff --git a/RTCP/IONProc/src/WH_DelayCompensation.cc b/RTCP/IONProc/src/WH_DelayCompensation.cc index e8cddccc9d61b52d74ffec35a85c7d60268fa86b..1262ca1a591054f0e1beb85d3a56cbdf5633614a 100644 --- a/RTCP/IONProc/src/WH_DelayCompensation.cc +++ b/RTCP/IONProc/src/WH_DelayCompensation.cc @@ -57,8 +57,8 @@ namespace LOFAR itsNrCalcDelays = ps->getUint32("OLAP.DelayComp.nrCalcDelays"); - // Pre-allocate and initialize storage for the delay vectors. - itsDelaysAfterEnd.resize(itsNrCalcDelays*itsNrBeams); + // Pre-allocate and initialize storage for the beam direction vectors. + itsBeamDirectionsAfterEnd.resize(itsNrCalcDelays*itsNrBeams); itsObservationEpoch.resize(itsNrCalcDelays); getBeamDirections(ps); @@ -80,7 +80,7 @@ namespace LOFAR itsConverter = 0; } - vector<double> WH_DelayCompensation::calcDelaysForEachTimeNrDirections(vector<double> &startIntegrationTime) + vector<Direction> WH_DelayCompensation::calculateAllBeamDirections(vector<double> &startIntegrationTime) { LOG_TRACE_FLOW(AUTO_FUNCTION_NAME); for (uint i = 0; i < startIntegrationTime.size(); i++) { @@ -88,15 +88,25 @@ namespace LOFAR } // Calculate the delays for the epoch after the end of the current time - // interval. Put the results in itsDelaysAtBegin and itsDelaysAfterEnd. - calculateDelays(); + // interval. Put the results in itsBeamDirectionsAfterEnd. + calculateDirections(); - return itsDelaysAfterEnd; + return itsBeamDirectionsAfterEnd; + } + + const Position& WH_DelayCompensation::getPositionDiffs() const + { + return itsPhasePositionDiffs; + } + + double WH_DelayCompensation::getDelay( Direction &dir ) const + { + return dir * itsPhasePositionDiffs * (1.0 / speedOfLight); } //##---------------- Private methods ----------------##// - void WH_DelayCompensation::getBeamDirections(const Parset *ps) + void WH_DelayCompensation::getBeamDirections(const Parset *ps) { LOG_TRACE_FLOW(AUTO_FUNCTION_NAME); @@ -156,8 +166,7 @@ namespace LOFAR return new ConverterImpl(); } - - void WH_DelayCompensation::calculateDelays() + void WH_DelayCompensation::calculateDirections() { LOG_TRACE_FLOW(AUTO_FUNCTION_NAME); @@ -175,9 +184,8 @@ namespace LOFAR ASSERTSTR(result.direction.size() == itsObservationEpoch.size() * itsNrBeams, result.direction.size() << " == " << itsObservationEpoch.size() * itsNrBeams); - LOG_TRACE_CALC("Beamlet directions:"); for (uint i = 0; i < result.direction.size(); ++i) { - LOG_TRACE_CALC_STR(" [" << i << "] = " << result.direction[i]); + itsBeamDirectionsAfterEnd[i] = result.direction[i]; } // From the source coordinates in ITRF, calculate the geometrical @@ -188,8 +196,7 @@ namespace LOFAR // etc. LOG_TRACE_CALC("Beamlet geometrical delays:"); for (uint i = 0; i < itsObservationEpoch.size() * itsNrBeams; ++i) { - itsDelaysAfterEnd[i] = (result.direction[i] * itsPhasePositionDiffs) * (1.0 / speedOfLight); - LOG_TRACE_CALC_STR(" [" << i << "]: " << itsDelaysAfterEnd[i]); + LOG_TRACE_CALC_STR(" [" << i << "]: " << getDelay( itsBeamDirectionsAfterEnd[i]) ); } } diff --git a/RTCP/IONProc/src/WH_DelayCompensation.h b/RTCP/IONProc/src/WH_DelayCompensation.h index 6a16f5c7904d2161fd69d9e12e157e386e22ce68..41a8652d572d2d3190c58961d257e85950d8d277 100644 --- a/RTCP/IONProc/src/WH_DelayCompensation.h +++ b/RTCP/IONProc/src/WH_DelayCompensation.h @@ -87,9 +87,13 @@ namespace LOFAR ~WH_DelayCompensation(); - vector<double> calcDelaysForEachTimeNrDirections(vector<double> &startIntegrationTime); + vector<AMC::Direction> calculateAllBeamDirections(vector<double> &startIntegrationTime); + + const AMC::Position& getPositionDiffs() const; + + double getDelay( AMC::Direction &direction ) const; - private: + private: // Copying is not allowed WH_DelayCompensation (const WH_DelayCompensation& that); WH_DelayCompensation& operator= (const WH_DelayCompensation& that); @@ -121,11 +125,9 @@ namespace LOFAR // itsResultAfterEnd. void doConvert(); - // Calculate the delays for all stations for the epoch after the end of - // the current time interval and store the results in itsDelaysAfterEnd. - // \post itsDelaysAtBegin contain the data previously contained by - // itsDelaysAtEnd. - void calculateDelays(); + // Calculate the beam directions for all stations for the epoch after the end of + // the current time interval and store the results in itsBeamDirectionsAfterEnd. + void calculateDirections(); // Number of beams. const uint itsNrBeams; @@ -152,9 +154,9 @@ namespace LOFAR // Station to reference station position difference vectors. AMC::Position itsPhasePositionDiffs; - // Delays for all stations for the epoch after the end of the current - // time interval. - vector<double> itsDelaysAfterEnd; + // Beam directions for all stations, for the epoch after the end of + // the current time interval. + vector<AMC::Direction> itsBeamDirectionsAfterEnd; unsigned itsNrCalcDelays; diff --git a/RTCP/Interface/include/Interface/CN_Configuration.h b/RTCP/Interface/include/Interface/CN_Configuration.h index 63b2d258ffbf98c33067bea8fb68021b3d3ba0ce..71cd2cb2ad6e07fd647e31b2f5b5489c6bc59244 100644 --- a/RTCP/Interface/include/Interface/CN_Configuration.h +++ b/RTCP/Interface/include/Interface/CN_Configuration.h @@ -24,6 +24,7 @@ #define LOFAR_INTERFACE_CN_CONFIGURATION_H #include <Stream/Stream.h> +#include <Interface/MultiDimArray.h> #include <vector> @@ -46,16 +47,23 @@ class CN_Configuration double &sampleRate(); std::vector<unsigned> &inputPsets(), &outputPsets(), &tabList(); std::vector<double> &refFreqs(); + unsigned &nrPencilRings(); + double &pencilRingSize(); + std::vector<double> &refPhaseCentre(); + Matrix<double> &phaseCentres(); void read(Stream *); void write(Stream *); static const unsigned MAX_PSETS = 64; static const unsigned MAX_SUBBANDS = 248; + static const unsigned MAX_STATIONS = 100; private: std::vector<unsigned> itsInputPsets, itsOutputPsets, itsTabList; std::vector<double> itsRefFreqs; + std::vector<double> itsRefPhaseCentre; + Matrix<double> itsPhaseCentres; struct MarshalledData { @@ -73,6 +81,10 @@ class CN_Configuration unsigned itsRefFreqsSize; unsigned itsInputPsets[MAX_PSETS], itsOutputPsets[MAX_PSETS], itsTabList[MAX_PSETS]; double itsRefFreqs[MAX_SUBBANDS]; + unsigned itsNrPencilRings; + double itsPencilRingSize; + double itsRefPhaseCentre[3]; + double itsPhaseCentres[MAX_STATIONS * 3]; } itsMarshalledData; }; @@ -147,6 +159,26 @@ inline std::vector<double> & CN_Configuration::refFreqs() return itsRefFreqs; } +inline unsigned &CN_Configuration::nrPencilRings() +{ + return itsMarshalledData.itsNrPencilRings; +} + +inline double &CN_Configuration::pencilRingSize() +{ + return itsMarshalledData.itsPencilRingSize; +} + +inline std::vector<double> &CN_Configuration::refPhaseCentre() +{ + return itsRefPhaseCentre; +} + +inline Matrix<double> &CN_Configuration::phaseCentres() +{ + return itsPhaseCentres; +} + } // namespace RTCP } // namespace LOFAR diff --git a/RTCP/Interface/include/Interface/SubbandMetaData.h b/RTCP/Interface/include/Interface/SubbandMetaData.h index fa92687012bdbc267adb8b54462bf014e6dd78c5..ab7ff4290e52e30141df5e4fd6b057474bc44198 100644 --- a/RTCP/Interface/include/Interface/SubbandMetaData.h +++ b/RTCP/Interface/include/Interface/SubbandMetaData.h @@ -39,6 +39,8 @@ struct SubbandMetaData void setFlags(const SparseSet<unsigned> &); float delayAtBegin, delayAfterEnd; + double beamDirectionAtBegin[3], beamDirectionAfterEnd[3]; + unsigned alignmentShift; private: diff --git a/RTCP/Interface/src/CN_Configuration.cc b/RTCP/Interface/src/CN_Configuration.cc index 2838da0fb195975bd146871653fa1a81e6d39768..ddb8e5a9d0fd4331268232c46f46fde3de700f45 100644 --- a/RTCP/Interface/src/CN_Configuration.cc +++ b/RTCP/Interface/src/CN_Configuration.cc @@ -43,6 +43,14 @@ void CN_Configuration::read(Stream *str) itsRefFreqs.resize(itsMarshalledData.itsRefFreqsSize); memcpy(&itsRefFreqs[0], itsMarshalledData.itsRefFreqs, itsMarshalledData.itsRefFreqsSize * sizeof(double)); + + itsRefPhaseCentre.resize(3); + memcpy(&itsRefPhaseCentre[0], itsMarshalledData.itsRefPhaseCentre, 3 * sizeof(double)); + + itsPhaseCentres.resize(nrStations(),3); + for( unsigned stat = 0; stat < nrStations(); stat++ ) { + memcpy(&itsPhaseCentres[stat][0], &itsMarshalledData.itsPhaseCentres[stat*3], 3 * sizeof(double)); + } } @@ -63,7 +71,12 @@ void CN_Configuration::write(Stream *str) itsMarshalledData.itsRefFreqsSize = itsRefFreqs.size(); assert(itsMarshalledData.itsRefFreqsSize <= MAX_SUBBANDS); memcpy(itsMarshalledData.itsRefFreqs, &itsRefFreqs[0], itsMarshalledData.itsRefFreqsSize * sizeof(double)); - + + memcpy(itsMarshalledData.itsRefPhaseCentre, &itsRefPhaseCentre[0], 3 * sizeof(double)); + for( unsigned stat = 0; stat < nrStations(); stat++ ) { + memcpy(&itsMarshalledData.itsPhaseCentres[stat*3], &itsPhaseCentres[stat][0], 3 * sizeof(double)); + } + str->write(&itsMarshalledData, sizeof itsMarshalledData); } diff --git a/RTCP/Run/src/RTCP.parset b/RTCP/Run/src/RTCP.parset index a373a594ec572778cfce734af0c344502f1b65d2..32e1fd424c2f8f31ddfe1bf44cfc21c37aa3c185 100644 --- a/RTCP/Run/src/RTCP.parset +++ b/RTCP/Run/src/RTCP.parset @@ -20,6 +20,8 @@ Observation.MSNameMask = /data/L${YEAR}_${MSNUMBER}/SB${SUBBAND}.MS # Variables for the DelayCompensation OLAP.DelayComp.positionType = ITRF # should be ITRF OLAP.DelayComp.nrCalcDelays = 16 +OLAP.DelayComp.nrPencilRings = 0 +OLAP.DelayComp.pencilRingSize = 0.01 OLAP.IPHeaderSize = 32 OLAP.EPAHeaderSize = 16