From a347be86b12c186b78bf5d378e14c20073d2849d Mon Sep 17 00:00:00 2001 From: Pieter Donker <donker@astron.nl> Date: Mon, 18 Jan 2016 10:50:06 +0000 Subject: [PATCH] Task #8966: backup --- .../include/APL/CAL_Protocol/AntennaGains.h | 169 ++++---- .../include/APL/CAL_Protocol/CMakeLists.txt | 11 +- .../include/APL/CAL_Protocol/SpectralWindow.h | 126 +++--- .../include/APL/CAL_Protocol/SubArray.h | 218 +++++----- MAC/APL/PAC/CAL_Protocol/src/AntennaGains.cc | 135 +++--- .../PAC/CAL_Protocol/src/CAL_Protocol.prot | 275 ++++++------ MAC/APL/PAC/CAL_Protocol/src/CMakeLists.txt | 5 +- .../PAC/CAL_Protocol/src/SpectralWindow.cc | 193 ++++----- MAC/APL/PAC/CAL_Protocol/src/SubArray.cc | 312 +++++++------- MAC/APL/PAC/CMakeLists.txt | 2 +- MAC/APL/PAC/Cal_Server/src/ACMProxy.h | 4 +- MAC/APL/PAC/Cal_Server/src/CMakeLists.txt | 10 +- MAC/APL/PAC/Cal_Server/src/CalServer.cc | 245 +++++------ MAC/APL/PAC/Cal_Server/src/CalServer.h | 19 +- .../PAC/Cal_Server/src/CalibrationAlgorithm.h | 8 +- .../PAC/Cal_Server/src/CalibrationThread.cc | 20 +- .../PAC/Cal_Server/src/CalibrationThread.h | 13 +- .../src/RemoteStationCalibration.cc | 395 +++++++++--------- .../Cal_Server/src/SubArraySubscription.cc | 22 +- MAC/APL/PAC/Cal_Server/src/SubArrays.cc | 12 +- MAC/APL/PAC/Cal_Server/src/SubArrays.h | 9 +- .../CalibrationControl/CalibrationControl.cc | 60 ++- 22 files changed, 1087 insertions(+), 1176 deletions(-) diff --git a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/AntennaGains.h b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/AntennaGains.h index 57c362ee419..1d15a7c4e89 100644 --- a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/AntennaGains.h +++ b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/AntennaGains.h @@ -30,93 +30,94 @@ #include <pthread.h> namespace LOFAR { - namespace CAL { + namespace CAL { + +/** + * This class holds the results of a remote station calibration: the + * calibrated antenna gains. Along with the gains a quality measure + * is computed by the calibration algorithm. This quality measure indicates + * the confidence in the computed antenna gain. + */ +class AntennaGains +{ +public: + + /*@{*/ + /** + * Constructors + */ + AntennaGains(); + AntennaGains(uint nRCUs, uint nsubbands); + /*@}*/ + + /** + * Destructor + */ + virtual ~AntennaGains(); + // Make a copy of this class + AntennaGains* clone() const; + + /** + * Get reference to the array with calibrated antenna gains. + * @return a reference to the calibrated gains. A three dimensional array of + * complex doubles with dimensions: nantennas x npol x nsubbands + */ + const blitz::Array<std::complex<double>, 2>& getGains() const { return m_gains; } /** - * This class holds the results of a remote station calibration: the - * calibrated antenna gains. Along with the gains a quality measure - * is computed by the calibration algorithm. This quality measure indicates - * the confidence in the computed antenna gain. + * Get reference to the array with quality measure. + * @return a reference to the quality measure array. A 3-dimensional array + * of doubles with nantennas x npol x nsubbands elements. */ - class AntennaGains - { - public: - - /*@{*/ - /** - * Constructors - */ - AntennaGains(); - AntennaGains(int nantennas, int npol, int nsubbands); - /*@}*/ - - /** - * Destructor - */ - virtual ~AntennaGains(); - - /** - * Get reference to the array with calibrated antenna gains. - * @return a reference to the calibrated gains. A three dimensional array of - * complex doubles with dimensions: nantennas x npol x nsubbands - */ - const blitz::Array<std::complex<double>, 3>& getGains() const { return m_gains; } - - /** - * Get reference to the array with quality measure. - * @return a reference to the quality measure array. A 3-dimensional array - * of doubles with nantennas x npol x nsubbands elements. - */ - const blitz::Array<double, 3>& getQuality() const { return m_quality; } - - /** - * has the calibration algorithm producing this result completed? - */ - inline bool isDone() const { bool done; lock(); done = m_done; unlock(); return done; } - - /** - * set the complete status. - */ - inline void setDone(bool value = true) { lock(); m_done = value; unlock(); } - - /** - * assignment operator - * @param rhs Right-hand-side of assignment. It is not const - * because the rhs must be locked during the assignment. - */ - AntennaGains& operator=(const AntennaGains& rhs); - - /** - * lock/unlock - */ - inline int lock() const { return pthread_mutex_lock((pthread_mutex_t*)m_mutex); } - inline int unlock() const { return pthread_mutex_unlock((pthread_mutex_t*)m_mutex); } - - public: - /*@{*/ - /** - * marshalling methods - */ - size_t getSize() const; - size_t pack (char* buffer) const; - size_t unpack (const char* buffer); - /*@}*/ - - private: - /** - * Prevent copy constructor - */ - AntennaGains(const AntennaGains & copy); - - private: - blitz::Array<std::complex<double>, 3> m_gains; - blitz::Array<double, 3> m_quality; - bool m_done; // has the calibration finished - - const pthread_mutex_t* m_mutex; // control access to the m_done flag - }; - - }; // namespace CAL + const blitz::Array<double, 2>& getQuality() const { return m_quality; } + + /** + * has the calibration algorithm producing this result completed? + */ + inline bool isDone() const { bool done; lock(); done = m_done; unlock(); return done; } + + /** + * set the complete status. + */ + inline void setDone(bool value = true) { lock(); m_done = value; unlock(); } + + /** + * assignment operator + * @param rhs Right-hand-side of assignment. It is not const + * because the rhs must be locked during the assignment. + */ + AntennaGains& operator=(const AntennaGains& rhs); + + /** + * lock/unlock + */ + inline int lock() const { return pthread_mutex_lock((pthread_mutex_t*)m_mutex); } + inline int unlock() const { return pthread_mutex_unlock((pthread_mutex_t*)m_mutex); } + +public: + /*@{*/ + /** + * marshalling methods + */ + size_t getSize() const; + size_t pack (char* buffer) const; + size_t unpack (const char* buffer); + /*@}*/ + +private: + /** + * Prevent copy constructor + */ + AntennaGains(const AntennaGains & copy); + + blitz::Array<std::complex<double>, 2> m_gains; + blitz::Array<double, 2> m_quality; + bool m_done; // has the calibration finished + + const pthread_mutex_t* m_mutex; // control access to the m_done flag +}; + + }; // namespace CAL }; // namespace LOFAR #endif /* CALIBRATIONRESULT_H_ */ diff --git a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/CMakeLists.txt b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/CMakeLists.txt index 274dd0a17f5..97f9999210e 100644 --- a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/CMakeLists.txt +++ b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/CMakeLists.txt @@ -2,23 +2,18 @@ # List of header files that will be installed. set(inst_HEADERS - AntennaArray.h - AntennaArrayData.h AntennaGains.h - SharedResource.h SpectralWindow.h - SubArray.h - CalibrationInterface.h - ACC.h) + SubArray.h) # Add definition of include path suffix definition, because the header files -# reside in a non-standard location. +# reside in a non-standard location. set(${PACKAGE_NAME}_INCLUDE_PATH_SUFFIX APL/CAL_Protocol CACHE INTERNAL "Include path suffix for package \"${PACKAGE_NAME}\"") # Create symbolic link to include directory, create directory first, if needed. file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/include/APL) -execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink +execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_BINARY_DIR}/include/APL/CAL_Protocol) diff --git a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SpectralWindow.h b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SpectralWindow.h index 6b0fa59972c..81cf671a0ec 100644 --- a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SpectralWindow.h +++ b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SpectralWindow.h @@ -32,81 +32,65 @@ namespace LOFAR { namespace CAL { - /** - * Class which represents a window on the electromagnetic spectrum. - * - * A window is defined by three parameters: - * - sampling frequency - * - nyquist_zone - * - nsubbands - * - * The band starting at frequency ((nyquist_zone - 1) * (sampling_frequency / 2.0)) and - * ending at (nyquist_zone * (sampling_frequency / 2.0)) is filtered into - * nsubband equal frequency bins. - * - * The method getSubbandFreq(subband) returns the center frequency of a subband. - */ - class SpectralWindow - { - public: - // Constructors - SpectralWindow(); - SpectralWindow(std::string name, double sampling_freq, - int nyquist_zone, int numsubbands, uint32 rcucontrol); - virtual ~SpectralWindow(); - - // Return the name of the spectral window. - std::string getName() const { return m_name; } - - // Return the sampling frequency for this window - double getSamplingFrequency() const { return m_sampling_freq; } - - // Return the nyquist zone for this window. - int getNyquistZone() const { return m_nyquist_zone; } - - // Return the number of subbands for the spectral window. - int getNumSubbands() const { return m_numsubbands; } - - // Return the width of the subbands. - double getSubbandWidth() const { return m_sampling_freq / (2.0 * MAX_SUBBANDS); } - - // Return frequency of a specific subband. - double getSubbandFreq(int subband) const; - - // Returns try if spectralWindow is ment for the HBA antennas. - bool isForHBA() const; - - // Returns the rcumode of this SPW - int rcumode() const; - - // Returns the rcumode of this SPW - uint32 raw_rcumode() const; - - // Output function - ostream& print (ostream& os) const; - - public: - /*@{*/ - // marshalling methods - size_t getSize() const; - size_t pack (char* buffer) const; - size_t unpack(const char *buffer); - /*@}*/ - - private: - std::string m_name; // name of the spectral window - double m_sampling_freq; // sampling frequency - uint16 m_nyquist_zone; // defines the window - uint16 m_numsubbands; // number of subbands - uint32 m_rcucontrol; // RCU control setting - }; - +// Class which represents a window on the electromagnetic spectrum. // -// operator<< +// A window is defined by two parameters: +// - sampling frequency +// - nyquist_zone // +// The band starting at frequency ((nyquist_zone - 1) * (sampling_frequency / 2.0)) and +// ending at (nyquist_zone * (sampling_frequency / 2.0)) is split into MAX_SUBBANDS frequency bins +class SpectralWindow +{ +public: + // Constructors + SpectralWindow(); + explicit SpectralWindow(uint rcumode); + SpectralWindow(const string& name, double sampling_freq, int nyquist_zone, bool LBAfilterOn); + ~SpectralWindow(); + + // Return the name of the spectral window. + string name() const { return itsName; } + + // Return the sampling frequency for this window + double samplingFrequency() const { return itsSamplingFreq; } + + // Return the nyquist zone for this window. + int nyquistZone() const { return itsNyquistZone; } + + // Return the LBA filter setting + int LBAfilterOn() const { return itsLBAfilterOn; } + + // Return the rcumode of SPW (only defined for HBA SPW's). + uint rcumodeHBA() const; + + // Return the width of the subbands. + double subbandWidth() const { return itsSamplingFreq / (2.0 * MAX_SUBBANDS); } + + // Return centre frequency of a specific subband. + double subbandFreq(int subband) const; + + /*@{*/ + // marshalling methods + size_t getSize() const; + size_t pack (char* buffer) const; + size_t unpack(const char *buffer); + /*@}*/ + + // call for operator<< + ostream& print (ostream& os) const; + +private: + string itsName; // name of the spectral window + double itsSamplingFreq; // sampling frequency + uint16 itsNyquistZone; // defines the window + bool itsLBAfilterOn; // 10-30Mhz filter switch on/off +}; + +// operator<< inline ostream& operator<< (ostream& os, const SpectralWindow& spw) { - return (spw.print(os)); + return (spw.print(os)); } }; // namespace CAL diff --git a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SubArray.h b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SubArray.h index b50abb153ae..eb2cbec7cc2 100644 --- a/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SubArray.h +++ b/MAC/APL/PAC/CAL_Protocol/include/APL/CAL_Protocol/SubArray.h @@ -21,117 +21,139 @@ //# //# $Id$ -#ifndef SUBARRAY_H_ -#define SUBARRAY_H_ +#ifndef CALPROTOCOL_SUBARRAY_H_ +#define CALPROTOCOL_SUBARRAY_H_ + +#include <blitz/array.h> #include <Common/lofar_map.h> -#include <Common/lofar_list.h> +//#include <Common/lofar_list.h> #include <Common/lofar_string.h> #include <Common/lofar_bitset.h> +#include <ApplCommon/StationDatatypes.h> -#include <APL/RTCCommon/Subject.h> -#include <APL/RSP_Protocol/MEPHeader.h> #include <APL/CAL_Protocol/SpectralWindow.h> -#include <APL/CAL_Protocol/AntennaArray.h> -/*#include "ACC.h"*/ -#include "SharedResource.h" #include <APL/CAL_Protocol/AntennaGains.h> +//#include <APL/RTCCommon/Subject.h> +//#include <APL/RSP_Protocol/MEPHeader.h> +//#include <APL/CAL_Protocol/AntennaArray.h> +//#include "SharedResource.h" +/*#include "ACC.h"*/ namespace LOFAR { - using EPA_Protocol::MEPHeader; + //using EPA_Protocol::MEPHeader; namespace CAL { // forward declarations -class ACC; -class CalibrationInterface; +//class ACC; +//class CalibrationInterface; -class SubArray : public AntennaArray, public RTC::Subject +//class SubArray : public AntennaArray, public RTC::Subject +class SubArray { public: - // Default constructor. - SubArray(); - - // Construct a subarray. - // @param name The name of the subarray. - // @param geoloc The geographical location of the subarray. - // @param pos The antenna positions of the parent array elements (nantennas x npolarizations x 3-coordinates). - // @param select Select for each polarization dipole of each antenna whether it is included (true) in the subarray. - // @param sampling_frequency The sampling frequency this runs at. - // @param nyquist_zone The nyquist zone in which we wish to measure. - // @param nsubbands The number of subbands of the spectral window. - // @param rcucontrol The RCU control setting (LB, HBL, HBH, etc). - SubArray(string name, - const blitz::Array<double, 1>& geoloc, - const blitz::Array<double, 3>& pos, - const blitz::Array<bool, 2>& select, - double sampling_frequency, - int nyquist_zone, - int nsubbands, - uint32 rcucontrol); - SubArray(string name); // used to return unknown subarray - virtual ~SubArray(); - - // Start (background) calibration of the subarray - // using the specified algorithm and ACC as input. - // @param cal The calibration algorithm to use. - // @param acc The Array Correlation Cube on which to calibrate. - void calibrate(CalibrationInterface* cal, ACC& acc); - - // Get calibration result (if available). - // @param cal Calibration result - bool getGains(AntennaGains*& cal, int buffer = FRONT); - - // get bitset containing the rcu's of the subArray. - typedef bitset<MEPHeader::MAX_N_RCUS> RCUmask_t; - RCUmask_t getRCUMask() const; - - // Abort background calibration. - void abortCalibration(); - - // Check whether calibration has completed. - bool isDone(); - - // Used to clear the 'done' flag after updating all subscriptions. - void clearDone(); - - // Get a reference to the spectral window for this subarray. - const SpectralWindow& getSPW() const; - - // Assignement operator. - SubArray& operator=(const SubArray& rhs); - - // Enumeration of buffer positions. - enum { - FRONT = 0, - BACK = 1 - }; - - //@{ - // marshalling methods - size_t getSize() const; - size_t pack (char* buffer) const; - size_t unpack(const char *buffer); - //@} + // default constructor + SubArray(); // needed for packing/unpacking + + // Construct a subarray. + // @param name The name of the subarray. + // @param geoloc The geographical location of the subarray. + // @param pos The antenna positions of the parent array elements (nantennas x npolarizations x 3-coordinates). + // @param select Select for each polarization dipole of each antenna whether it is included (true) in the subarray. + // @param sampling_frequency The sampling frequency this runs at. + // @param nyquist_zone The nyquist zone in which we wish to measure. + // @param nsubbands The number of subbands of the spectral window. + // @param rcucontrol The RCU control setting (LB, HBL, HBH, etc). + SubArray(const string& name, + const string& antennaSet, + RCUmask_t RCUmask, + bool LBAfilterOn, + double sampling_frequency, + int nyquist_zone); + SubArray(const string& name); // used to return unknown subarray + ~SubArray(); + + // Start (background) calibration of the subarray + // using the specified algorithm and ACC as input. + // @param cal The calibration algorithm to use. + // @param acc The Array Correlation Cube on which to calibrate. + void updateGains(const AntennaGains& newGains); + //void calibrate(CalibrationInterface* cal, ACC& acc); + + // Get calibration result (if available). + // @param cal Calibration result + AntennaGains& getGains() + { return (*itsGains); } + + // Write gains to a file + void writeGains(); + + // Abort background calibration. + //void abortCalibration(); + + // get name of Subarray + const string& name() const + { return (itsName); } + + // get antennaSetname of Subarray + const string& antennaSet() const + { return (itsAntennaSet); } + + // get bitset containing the rcu's of the subArray. + RCUmask_t RCUMask(uint rcumode) const; + RCUmask_t RCUMask() const + { return (itsRCUmask); } + + // Check whether calibration has completed. + //bool isDone(); + + // Used to clear the 'done' flag after updating all subscriptions. + //void clearDone(); + + // Get a reference to the spectral window for this subarray. + const SpectralWindow& SPW() const + { return (itsSPW); } + + // Get a reference to the RCUmodes + const blitz::Array<uint,1>& RCUmodes() const + { return (itsRCUmodes); } + + // Test if RCUmode is used by this array + bool usesRCUmode(int rcumode) const; + + // Assignment operator. + SubArray& operator=(const SubArray& rhs); + // Enumeration of buffer positions. + + //@{ + // marshalling methods + size_t getSize() const; + size_t pack (char* buffer) const; + size_t unpack(const char *buffer); + //@} + + // call for operator<< + ostream& print (ostream& os) const; private: - // prevent copy - SubArray(const SubArray& other); // no implementation - - int m_antenna_count; // number of seleted antennas - blitz::Array<bool, 2> m_antenna_selection; // antenna selection dimensions: - // (nantennas x npol) - RCUmask_t itsRCUmask; - - SpectralWindow m_spw; // the spectral window for this subarray - AntennaGains* m_result[BACK + 1]; // two calibration result records + // prevent copy + SubArray(const SubArray& other); // no implementation + + string itsName; // unique name choosen by user. + string itsAntennaSet; // name of the used antennaSet + SpectralWindow itsSPW; // the spectral window for this subarray + RCUmask_t itsRCUmask; // used RCU's + + blitz::Array<uint,1> itsRCUmodes; // redundant info (=AntSet+SPW for each rcu in RCUmask) + bitset<NR_RCU_MODES+1> itsRCUuseFlags; // which modes are used. +// bool itsLBAfilterOn; + AntennaGains* itsGains; // calibration result record }; -// -// getRCUMask() -// -inline SubArray::RCUmask_t SubArray::getRCUMask() const +// operator<< +inline ostream& operator<< (ostream& os, const SubArray& anSA) { - return (itsRCUmask); + return (anSA.print(os)); } // ------------------- SubArraymap ------------------- @@ -140,16 +162,16 @@ inline SubArray::RCUmask_t SubArray::getRCUMask() const class SubArrayMap : public map<string, SubArray*> { public: - //@{ - // marshalling methods - size_t getSize() const; - size_t pack (char* buffer) const; - size_t unpack(const char *buffer); - //@} + //@{ + // marshalling methods + size_t getSize() const; + size_t pack (char* buffer) const; + size_t unpack(const char *buffer); + //@} }; }; // namespace CAL }; // namespace LOFAR -#endif +#endif diff --git a/MAC/APL/PAC/CAL_Protocol/src/AntennaGains.cc b/MAC/APL/PAC/CAL_Protocol/src/AntennaGains.cc index 2b6459b740d..2828ba94220 100644 --- a/MAC/APL/PAC/CAL_Protocol/src/AntennaGains.cc +++ b/MAC/APL/PAC/CAL_Protocol/src/AntennaGains.cc @@ -30,101 +30,112 @@ using namespace std; using namespace blitz; -using namespace LOFAR; -using namespace CAL; +namespace LOFAR { + namespace CAL { AntennaGains::AntennaGains() : m_mutex(new pthread_mutex_t) { - ASSERT(m_mutex); - pthread_mutex_init((pthread_mutex_t*)m_mutex, 0); - lock(); + ASSERT(m_mutex); + pthread_mutex_init((pthread_mutex_t*)m_mutex, 0); + lock(); - m_done = false; + m_done = false; - m_gains.resize(1,1,1); - m_gains = 0; - - m_quality.resize(1,1,1); - m_quality = 0; + m_gains.resize(1,1,1); + m_gains = 0; - unlock(); + m_quality.resize(1,1,1); + m_quality = 0; + + unlock(); } -AntennaGains::AntennaGains(int nantennas, int npol, int nsubbands) : m_mutex(new pthread_mutex_t) +AntennaGains::AntennaGains(uint nrRCUs, uint nsubbands) : m_mutex(new pthread_mutex_t) { - ASSERT(m_mutex); - pthread_mutex_init((pthread_mutex_t*)m_mutex, 0); - lock(); m_done = false; unlock(); - - if (nantennas < 0 || npol < 0 || nsubbands < 0) { - nantennas = 0; - npol = 0; - nsubbands = 0; - } - - m_gains.resize(nantennas, npol, nsubbands); - m_gains = 1; + ASSERT(m_mutex); + pthread_mutex_init((pthread_mutex_t*)m_mutex, 0); + lock(); m_done = false; unlock(); + + m_gains.resize(nrRCUs, nsubbands); + m_gains = 1; - m_quality.resize(nantennas, npol, nsubbands); - m_quality = 1; + m_quality.resize(nrRCUs, nsubbands); + m_quality = 1; } AntennaGains::~AntennaGains() { - if (m_mutex) - delete m_mutex; + if (m_mutex) + delete m_mutex; +} + +// +// clone() +// +AntennaGains* AntennaGains::clone() const +{ + AntennaGains* theClone = new AntennaGains(m_gains.extent(firstDim), m_gains.extent(secondDim)); + ASSERTSTR(theClone, "Could not clone the AntennaGains class"); + + theClone->m_gains = m_gains; + theClone->m_quality = m_quality; + theClone->m_done = m_done; + + return (theClone); } size_t AntennaGains::getSize() const { - return - MSH_size(m_gains) - + MSH_size(m_quality) - + sizeof(bool); + return + MSH_size(m_gains) + + MSH_size(m_quality) + + sizeof(bool); } size_t AntennaGains::pack(char* buffer) const { - size_t offset = 0; + size_t offset = 0; - lock(); - offset = MSH_pack(buffer, offset, m_gains); - offset = MSH_pack(buffer, offset, m_quality); - memcpy((char*)buffer + offset, &m_done, sizeof(bool)); - offset += sizeof(bool); - unlock(); + lock(); + offset = MSH_pack(buffer, offset, m_gains); + offset = MSH_pack(buffer, offset, m_quality); + memcpy((char*)buffer + offset, &m_done, sizeof(bool)); + offset += sizeof(bool); + unlock(); - return offset; + return offset; } size_t AntennaGains::unpack(const char* buffer) { - size_t offset = 0; + size_t offset = 0; - lock(); - offset = MSH_unpack(buffer, offset, m_gains); - offset = MSH_unpack(buffer, offset, m_quality); - memcpy(&m_done, (char*)buffer + offset, sizeof(bool)); - offset += sizeof(bool); - unlock(); + lock(); + offset = MSH_unpack(buffer, offset, m_gains); + offset = MSH_unpack(buffer, offset, m_quality); + memcpy(&m_done, (char*)buffer + offset, sizeof(bool)); + offset += sizeof(bool); + unlock(); - return offset; + return offset; } AntennaGains& AntennaGains::operator=(const AntennaGains& rhs) { - if (this != &rhs) { - lock(); rhs.lock(); - - m_gains.resize(rhs.m_gains.shape()); - m_gains = rhs.m_gains.copy(); - - m_quality.resize(rhs.m_quality.shape()); - m_quality = rhs.m_quality.copy(); - - m_done = rhs.m_done; - rhs.unlock(); unlock(); - } - return *this; + if (this != &rhs) { + lock(); rhs.lock(); + + m_gains.resize(rhs.m_gains.shape()); + m_gains = rhs.m_gains.copy(); + + m_quality.resize(rhs.m_quality.shape()); + m_quality = rhs.m_quality.copy(); + + m_done = rhs.m_done; + rhs.unlock(); unlock(); + } + return *this; } + } // namespace CAL +} // namespace LOFAR diff --git a/MAC/APL/PAC/CAL_Protocol/src/CAL_Protocol.prot b/MAC/APL/PAC/CAL_Protocol/src/CAL_Protocol.prot index 0c06263a0b8..f334d1b1fb5 100644 --- a/MAC/APL/PAC/CAL_Protocol/src/CAL_Protocol.prot +++ b/MAC/APL/PAC/CAL_Protocol/src/CAL_Protocol.prot @@ -34,12 +34,11 @@ enum { ERR_ALREADY_REGISTERED, // subarray already registered ERR_ONLY_ONE_SUBARRAY, // only 1 subarray per client allowed ERR_NO_SUBARRAY_NAME, // name of subarray is empty - ERR_NO_ANTENNAS, // antenna selection is empty - ERR_CONFLICT // array conflicts with other array + ERR_NO_ANTENNAS // antenna selection is empty }; // -// The following events are defined in the CAL protocol +// The following events are defined in the ICAL protocol // // START (subarrayname, antennaSet, receiverset, RCUsettings); // STARTACK (subarrayname, status); @@ -65,173 +64,161 @@ PRELUDE_END; // A "param" has a "name" and a "type". // event = { - signal = START; // start calibration of a subarray - dir = IN; - param = { // name of the subarray (used as nodeid in key-value logger) - name = "name"; - type = "string"; - }; - param = { // name of the array of which it is a subarray - name = "parent"; - type = "string"; - }; - param = { // subset of the receivers (2 receivers per antenna) - name = "subset"; - type = "bitset<LOFAR::MAX_RCUS>"; - }; - param = { - // RCU control byte controls the setting of all - // RCU's for this subarray - name = "rcumode"; - type = "RSP_Protocol::RCUSettings"; - userdefined; - }; - // TODO: add this parameter for faster calibration in the future -// param = { // set of subbands -// name = "subbandset"; -// type = "bitset<LOFAR::MAX_SUBBANDS>"; -// }; - // nyquist_zone is deducted from RCUSettings - // sampling_frequency and nyquist_zone determine the spectral window + signal = START; // start calibration of a subarray + dir = IN; + param = { // name of the subarray (used as nodeid in key-value logger) + name = "name"; + type = "string"; + }; + param = { // name of the array of which it is a subarray + name = "antennaSet"; + type = "string"; + }; + param = { // subset of the receivers (2 receivers per antenna) + name = "rcuMask"; + type = "LOFAR::bitset<LOFAR::MAX_RCUS>"; + }; + param = { + // RCU control byte controls the setting of all + // RCU's for this subarray + name = "rcumode"; + type = "int32"; +// type = "RSP_Protocol::RCUSettings"; +// userdefined; + }; + // TODO: add this parameter for faster calibration in the future + // param = { // set of subbands + // name = "subbandset"; + // type = "LOFAR::bitset<LOFAR::MAX_SUBBANDS>"; + // }; + // nyquist_zone is deducted from RCUSettings + // sampling_frequency and nyquist_zone determine the spectral window }; event = { - signal = STARTACK; - dir = OUT; - param = { - name = "name"; - type = "string"; - }; - param = { - name = "status"; - type = "int16"; - }; + signal = STARTACK; + dir = OUT; + param = { + name = "name"; + type = "string"; + }; + param = { + name = "status"; + type = "int16"; + }; }; event = { - signal = STOP; // stop calibration of a subarray - dir = IN; - param = { - name = "name"; - type = "string"; - }; + signal = STOP; // stop calibration of a subarray + dir = IN; + param = { + name = "name"; + type = "string"; + }; }; event = { - signal = STOPACK; - dir = OUT; - param = { - name = "name"; - type = "string"; - }; - param = { - name = "status"; - type = "int16"; - }; + signal = STOPACK; + dir = OUT; + param = { + name = "name"; + type = "string"; + }; + param = { + name = "status"; + type = "int16"; + }; }; event = { - signal = SUBSCRIBE; // subscribe to antenna gain updates - dir = IN; - param = { // name of the subarray - name = "name"; - type = "string"; - }; - param = { // set of subbands - name = "subbandset"; - type = "bitset<LOFAR::MAX_SUBBANDS>"; - }; + signal = SUBSCRIBE; // subscribe to antenna gain updates + dir = IN; + param = { // name of the subarray + name = "name"; + type = "string"; + }; + param = { // set of subbands + name = "subbandset"; + type = "LOFAR::bitset<LOFAR::MAX_SUBBANDS>"; + }; }; event = { - signal = SUBSCRIBEACK; - dir = OUT; - // TODO: return name too - // param = { // name of the subarray - // name = "name"; - // type = "string"; - // }; - param = { - name = "status"; - type = "int16"; - }; - param = { - name = "handle"; - type = "memptr_t"; - }; - param = { // return the subarray class - name = "subarray"; - type = "CAL::SubArray"; - userdefined; - }; + signal = SUBSCRIBEACK; + dir = OUT; + param = { // name of the subarray + name = "name"; + type = "string"; + }; + param = { + name = "status"; + type = "int16"; + }; + param = { // return the subarray class + name = "subarray"; + type = "CAL::SubArray"; + userdefined; + }; }; event = { - // TODO: rename to UPDATEGAINS - signal = UPDATE; // antenna gains update - dir = OUT; - param = { // time of the gains - name = "timestamp"; - type = "RTC::Timestamp"; - userdefined; - }; - - // status == FAILURE indicates that the subarray was stoppe d - // This is not yet implemented. - param = { // status - name = "status"; - type = "int16"; - }; - param = { // handle - name = "handle"; - type = "memptr_t"; - }; - // TODO: AntennaGains is not good name, since the gains are - // computed for the complete signal path for a receiver. - param = { // gains - name = "gains"; - type = "CAL::AntennaGains"; - userdefined; - }; + signal = UPDATE; // antenna gains update + dir = OUT; + param = { // name of the subarray + name = "name"; + type = "string"; + }; + // status == FAILURE indicates that the subarray was stopped + // This is not yet implemented: stopped subarrays are not updated anymore + param = { // status + name = "status"; + type = "int16"; + }; + param = { // time of the gains + name = "timestamp"; + type = "RTC::Timestamp"; + userdefined; + }; + param = { // gains + name = "gains"; + type = "CAL::AntennaGains"; + userdefined; + }; }; event = { - signal = UNSUBSCRIBE; // unsubscribe from updates - dir = IN; - param = { - name = "name"; - type = "string"; - }; - param = { - name = "handle"; - type = "memptr_t"; - }; + signal = UNSUBSCRIBE; // unsubscribe from updates + dir = IN; + param = { + name = "name"; + type = "string"; + }; }; event = { - signal = UNSUBSCRIBEACK; - dir = OUT; - param = { - name = "name"; - type = "string"; - }; - param = { - name = "status"; - type = "int16"; - }; - param = { - name = "handle"; - type = "memptr_t"; - }; + signal = UNSUBSCRIBEACK; + dir = OUT; + param = { + name = "name"; + type = "string"; + }; + param = { + name = "status"; + type = "int16"; + }; + param = { + name = "handle"; + type = "memptr_t"; + }; }; event = { - signal = GETSUBARRAY; // Ask info about allocations - dir = IN; - param = { - name = "subarrayname"; // may be empty - type = "string"; - }; + signal = GETSUBARRAY; // Ask info about allocations + dir = IN; + param = { + name = "subarrayname"; // may be empty + type = "string"; + }; }; event = { @@ -244,7 +231,7 @@ event = { param = { name = "subarraymap"; type = "CAL::SubArrayMap"; - userdefined; + userdefined; }; }; diff --git a/MAC/APL/PAC/CAL_Protocol/src/CMakeLists.txt b/MAC/APL/PAC/CAL_Protocol/src/CMakeLists.txt index 1b432a64a2d..c02b51e6fb9 100644 --- a/MAC/APL/PAC/CAL_Protocol/src/CMakeLists.txt +++ b/MAC/APL/PAC/CAL_Protocol/src/CMakeLists.txt @@ -12,12 +12,9 @@ include_directories(${CMAKE_BINARY_DIR}/include/MAC) lofar_add_library(cal_protocol Package__Version.cc CAL_Protocol.cc - AntennaArray.cc - AntennaArrayData.cc AntennaGains.cc SpectralWindow.cc - SubArray.cc - ACC.cc) + SubArray.cc) add_dependencies(cal_protocol CAL_Protocol-CAL_Protocol) lofar_add_bin_program(versioncal_protocol versioncal_protocol.cc) diff --git a/MAC/APL/PAC/CAL_Protocol/src/SpectralWindow.cc b/MAC/APL/PAC/CAL_Protocol/src/SpectralWindow.cc index b019f5ae298..e6660f21538 100644 --- a/MAC/APL/PAC/CAL_Protocol/src/SpectralWindow.cc +++ b/MAC/APL/PAC/CAL_Protocol/src/SpectralWindow.cc @@ -27,161 +27,120 @@ #include <APL/CAL_Protocol/SpectralWindow.h> -#include <blitz/array.h> -#include <sstream> - #include <MACIO/Marshalling.tcc> #include <APL/RTCCommon/MarshallBlitz.h> -#include <math.h> - -using namespace std; -using namespace blitz; -using namespace LOFAR; -using namespace CAL; +namespace LOFAR { + namespace CAL { SpectralWindow::SpectralWindow() : - m_name("undefined"), m_sampling_freq(0), m_nyquist_zone(0), m_numsubbands(0), m_rcucontrol(0) + itsName("undefined"), itsSamplingFreq(0), itsNyquistZone(0), itsLBAfilterOn(0) { - LOG_TRACE_OBJ("SpectralWindow()"); + LOG_TRACE_OBJ("SpectralWindow()"); } -SpectralWindow::SpectralWindow(std::string name, double sampling_freq, - int nyquist_zone, int numsubbands, uint32 rcucontrol) : - m_name(name), - m_sampling_freq(sampling_freq), - m_nyquist_zone(nyquist_zone), - m_numsubbands(numsubbands), - m_rcucontrol(rcucontrol) -{ - LOG_TRACE_OBJ(formatString("SpectralWindow(%s,%f,%d,%d,%08X)", - name.c_str(), sampling_freq, nyquist_zone, numsubbands, rcucontrol)); -} - -SpectralWindow::~SpectralWindow() +SpectralWindow::SpectralWindow(uint rcumode) { - LOG_TRACE_OBJ("~SpectralWindow()"); + switch (rcumode) { + case 1: *this = SpectralWindow("rcumode1", 200.0e6, 1, false); // 1 + case 2: *this = SpectralWindow("rcumode2", 200.0e6, 1, true); // 2 + case 3: *this = SpectralWindow("rcumode3", 200.0e6, 1, false); // 1 + case 4: *this = SpectralWindow("rcumode4", 200.0e6, 1, true); // 2 + case 5: *this = SpectralWindow("rcumode5", 200.0e6, 2, false); // 3 + case 6: *this = SpectralWindow("rcumode6", 160.0e6, 3, false); // 4 + case 7: *this = SpectralWindow("rcumode7", 200.0e6, 3, false); // 5 + default: + ASSERTSTR(rcumode >= 1 && rcumode <= (uint)NR_RCU_MODES, "rcumode must have value: 1.." << NR_RCU_MODES); + } } -double SpectralWindow::getSubbandFreq(int subband) const +SpectralWindow::SpectralWindow(const string& name, double sampling_freq, + int nyquist_zone, bool LBAfilterOn) : + itsName (name), + itsSamplingFreq(sampling_freq), + itsNyquistZone (nyquist_zone), + itsLBAfilterOn (LBAfilterOn) { -// LOG_TRACE_OBJ_STR("SpectralWindow::getSubbandFreq(" << subband << " of " << m_numsubbands << ")"); - - ASSERT(m_numsubbands); - ASSERT(subband >= 0 && subband <= m_numsubbands); - - bool is160 = ::fabs(160e6 - m_sampling_freq) < 1e-4; - float freqOffset = (is160 ? 80.0e6 : 100.0e6) * (m_nyquist_zone - 1); - - return (freqOffset + ((subband % m_numsubbands) * getSubbandWidth())); + LOG_TRACE_OBJ(formatString("SpectralWindow(%s,%f,%d,%s)", + name.c_str(), sampling_freq, nyquist_zone, (LBAfilterOn ? "ON" : " OFF"))); } - -// -// isForHBA(): bool -// -bool SpectralWindow::isForHBA() const +SpectralWindow::~SpectralWindow() { - LOG_DEBUG (formatString("isForHBA(%06X)", m_rcucontrol)); - - switch (m_rcucontrol) { - case 0x017900: // LB_10_90 - case 0x057900: // LB_30_80 - case 0x037A00: // LBH_10_90 - case 0x077A00: // LBH_30_80 - return (false); - break; - - case 0x07A400: // HB_110_190 - case 0x079400: // HB_170_230 - case 0x078400: // HB_210_250 - return (true); - break; - - default: - LOG_WARN(formatString("Unknown RCUcontrol setting (0x%X), assuming LBA array", - m_rcucontrol)); - break; - } - - return (false); // assume LBA + LOG_TRACE_OBJ("~SpectralWindow()"); } // -// rcumode(): rcumode +// rcumodeHBA() // -int SpectralWindow::rcumode() const +uint SpectralWindow::rcumodeHBA() const { - LOG_DEBUG (formatString("rcumode(%06X)", m_rcucontrol)); - - switch (m_rcucontrol) { - case 0x017900: return (1); - case 0x057900: return (2); - case 0x037A00: return (3); - case 0x077A00: return (4); - case 0x07A400: return (5); - case 0x079400: return (6); - case 0x078400: return (7); - default: return (0); - } + switch (itsNyquistZone) { + case 1: return (0); // LBA + case 2: return (5); + case 3: return (itsSamplingFreq == 200.0e6 ? 7 : 6); + default: ASSERTSTR(false, "rcuMode cannot be determined, illegal nyquistzone."); + } } -uint32 SpectralWindow::raw_rcumode() const +double SpectralWindow::subbandFreq(int subband) const { - LOG_DEBUG (formatString("rcumode(%06X)", m_rcucontrol)); - return (m_rcucontrol); + ASSERT(subband >= 0 && subband < MAX_SUBBANDS); + + bool is160 = ::fabs(160e6 - itsSamplingFreq) < 1e-4; + float freqOffset = (is160 ? 80.0e6 : 100.0e6) * (itsNyquistZone - 1); + + return (freqOffset + (subband * subbandWidth())); } -// -// print -// +// print function for operator<< ostream& SpectralWindow::print(ostream& os) const { - os << formatString("SpectralWindow(name=%s,f=%f,nyq=%d,#sub=%d,rcuctrl=%08X)", - m_name.c_str(), m_sampling_freq, m_nyquist_zone, m_numsubbands, m_rcucontrol); - return (os); + os << "SpectralWindow " << itsName << ":sampleFreq=" << itsSamplingFreq << ", nyquistzone=" << itsNyquistZone; + os << ", LBAfilter=" << (itsLBAfilterOn ? "ON" : "OFF"); + return (os); } +// +// ---------- pack and unpack functions ---------- +// size_t SpectralWindow::getSize() const { - return MSH_size(m_name) + - sizeof(m_sampling_freq) + - sizeof(m_nyquist_zone) + - sizeof(m_numsubbands) + - sizeof(m_rcucontrol); + return MSH_size(itsName) + + sizeof(itsSamplingFreq) + + sizeof(itsNyquistZone) + + sizeof(itsLBAfilterOn); } size_t SpectralWindow::pack(char* buffer) const { - size_t offset = 0; - - offset = MSH_pack(buffer, offset, m_name); - memcpy(buffer + offset, &m_sampling_freq, sizeof(m_sampling_freq)); - offset += sizeof(m_sampling_freq); - memcpy(buffer + offset, &m_nyquist_zone, sizeof(m_nyquist_zone)); - offset += sizeof(m_nyquist_zone); - memcpy(buffer + offset, &m_numsubbands, sizeof(m_numsubbands)); - offset += sizeof(m_numsubbands); - memcpy(buffer + offset, &m_rcucontrol, sizeof(m_rcucontrol)); - offset += sizeof(m_rcucontrol); - - return offset; + size_t offset = 0; + + offset = MSH_pack(buffer, offset, itsName); + memcpy(buffer + offset, &itsSamplingFreq, sizeof(itsSamplingFreq)); + offset += sizeof(itsSamplingFreq); + memcpy(buffer + offset, &itsNyquistZone, sizeof(itsNyquistZone)); + offset += sizeof(itsNyquistZone); + memcpy(buffer + offset, &itsLBAfilterOn, sizeof(itsLBAfilterOn)); + offset += sizeof(itsLBAfilterOn); + + return offset; } size_t SpectralWindow::unpack(const char* buffer) { - size_t offset = 0; - - offset = MSH_unpack(buffer, offset, m_name); - memcpy(&m_sampling_freq, buffer + offset, sizeof(m_sampling_freq)); - offset += sizeof(m_sampling_freq); - memcpy(&m_nyquist_zone, buffer + offset, sizeof(m_nyquist_zone)); - offset += sizeof(m_nyquist_zone); - memcpy(&m_numsubbands, buffer + offset, sizeof(m_numsubbands)); - offset += sizeof(m_numsubbands); - memcpy(&m_rcucontrol, buffer + offset, sizeof(m_rcucontrol)); - offset += sizeof(m_rcucontrol); - - return offset; + size_t offset = 0; + + offset = MSH_unpack(buffer, offset, itsName); + memcpy(&itsSamplingFreq, buffer + offset, sizeof(itsSamplingFreq)); + offset += sizeof(itsSamplingFreq); + memcpy(&itsNyquistZone, buffer + offset, sizeof(itsNyquistZone)); + offset += sizeof(itsNyquistZone); + memcpy(&itsLBAfilterOn, buffer + offset, sizeof(itsLBAfilterOn)); + offset += sizeof(itsLBAfilterOn); + + return offset; } + } // namespace ICAL +} // namespace LOFAR diff --git a/MAC/APL/PAC/CAL_Protocol/src/SubArray.cc b/MAC/APL/PAC/CAL_Protocol/src/SubArray.cc index a2585d2b7d9..0f0e8d74076 100644 --- a/MAC/APL/PAC/CAL_Protocol/src/SubArray.cc +++ b/MAC/APL/PAC/CAL_Protocol/src/SubArray.cc @@ -21,237 +21,231 @@ //# //# $Id$ -#include <APL/CAL_Protocol/SubArray.h> -#include <APL/CAL_Protocol/CalibrationInterface.h> +#include <lofar_config.h> #include <Common/LofarLogger.h> #include <Common/StringUtil.h> #include <Common/hexdump.h> +#include <ApplCommon/AntennaSets.h> + #include <MACIO/Marshalling.tcc> #include <APL/RTCCommon/MarshallBlitz.h> +#include <APL/CAL_Protocol/SubArray.h> +//#include <APL/CAL_Protocol/CalibrationInterface.h> -using namespace std; -using namespace blitz; -using namespace LOFAR; -using namespace CAL; +namespace LOFAR { + namespace CAL { + using namespace blitz; // forward declaration class CalibrationInterface; - + // // SubArray() // -SubArray::SubArray() : - AntennaArray(), - m_antenna_selection(), - m_spw("undefined", 0, 0, 0, 0) +SubArray::SubArray() : + itsSPW("undefined", 0, 0, false), + itsGains(new AntennaGains) { - LOG_TRACE_OBJ("SubArray()"); - - m_result[FRONT] = 0; - m_result[BACK] = 0; + LOG_TRACE_OBJ("SubArray()"); } // -// SubArray(name, geoloc, pos, select, freq, nyquist, nrSubbands, rcuControl) +// SubArray(name, antennaSet, select, freq, nyquist, nrSubbands, rcuControl) // -SubArray::SubArray (string name, - const Array<double,1>& geoloc, - const Array<double,3>& pos, - const Array<bool, 2>& select, - double sampling_frequency, - int nyquist_zone, - int nsubbands, - uint32 rcucontrol) : - AntennaArray(name, geoloc, pos), - m_antenna_selection(select), - m_spw(name + "_spw", sampling_frequency, nyquist_zone, nsubbands, rcucontrol) +SubArray::SubArray (const string& name, + const string& antennaSet, + RCUmask_t RCUmask, + bool LBAfilterOn, + double sampling_frequency, + int nyquist_zone) : + itsName(name), + itsAntennaSet (antennaSet), + itsSPW (name + "_spw", sampling_frequency, nyquist_zone, LBAfilterOn), + itsRCUmask (RCUmask & globalAntennaSets()->RCUallocation(antennaSet)) { - LOG_TRACE_OBJ(formatString("SubArray(%s,%f,%d,%d,%08X)", - name.c_str(), sampling_frequency, nyquist_zone, nsubbands, rcucontrol)); - - // assert sizes - ASSERT(m_antenna_selection.extent(firstDim) == m_pos.extent(firstDim) - && m_antenna_selection.extent(secondDim) == m_pos.extent(secondDim) - && m_pos.extent(thirdDim) == 3); - - LOG_DEBUG_STR("name=" << name); - LOG_DEBUG_STR("select=" << select); - LOG_DEBUG_STR("m_pos=" << m_pos); - - // construct a new pos(ition) array containing only the coordinates of the - // the subarray. Remember in m_rcuindex of each of these elements to which - // of the global rcus the newpos element refers. - - // make array at least big enough - m_rcuindex.resize(m_pos.extent(firstDim), m_pos.extent(secondDim)); - m_rcuindex = -1; - - // will contain positions of antenna's that contribute to subarray - Array<double, 3> newpos(m_pos.shape()); - newpos = 0.0; - - // clear RCUbitset - itsRCUmask.reset(); - - // loop over inputs and update our admin. - int sel = 0; - for (int ant = 0; ant < m_pos.extent(firstDim); ant++) { - if (sum(m_antenna_selection(ant, Range::all())) > 0) { - for (int pol = 0; pol < m_pos.extent(secondDim); pol++) { - if (m_antenna_selection(ant, pol)) { - newpos(sel, pol, Range::all()) = m_pos(ant, pol, Range::all()); - m_rcuindex(sel, pol) = ant * m_pos.extent(secondDim) + pol; - itsRCUmask.set(ant * N_POL + pol); - } - } // for each pol - sel++; - } // antenna in subarray? - } // for each antenna - m_antenna_count = sel; - m_pos = newpos; // overwrite original positions - - ASSERT(m_antenna_count > 0); - - // resize the arrays - m_pos.resizeAndPreserve(m_antenna_count, m_pos.extent(secondDim), m_pos.extent(thirdDim)); - m_rcuindex.resizeAndPreserve(m_antenna_count, m_rcuindex.extent(secondDim)); - - LOG_DEBUG_STR("m_pos=" << m_pos); - LOG_DEBUG_STR("m_rcuindex=" << m_rcuindex); - LOG_DEBUG_STR("itsRCUmask=" << itsRCUmask); - - // create calibration result objects - m_result[FRONT] = new AntennaGains(m_pos.extent(firstDim), m_pos.extent(secondDim), m_spw.getNumSubbands()); - m_result[BACK] = new AntennaGains(m_pos.extent(firstDim), m_pos.extent(secondDim), m_spw.getNumSubbands()); - ASSERT(m_result[FRONT] && m_result[BACK]); + LOG_DEBUG(formatString("SubArray(%s,%f,%d,%s)", + name.c_str(), sampling_frequency, nyquist_zone, (LBAfilterOn ? "ON" : "OFF"))); + + // create calibration result objects [ant x pol x subbands] + itsGains = new AntennaGains(itsRCUmask.count(), MAX_SUBBANDS); // TODO: does this work with non contiguous RCUmasks???? + ASSERT(itsGains); + + // fill rcumode array + string RCUinputs(globalAntennaSets()->RCUinputs(itsAntennaSet)); + itsRCUmodes.resize(MAX_RCUS); + itsRCUmodes = 0; + for (int rcu = 0; rcu < MAX_RCUS; rcu++) { + if (itsRCUmask.test(rcu)) { + switch (RCUinputs[rcu]) { + case 'l': itsRCUmodes(rcu) = (LBAfilterOn ? 2 : 1); break; + case 'h': itsRCUmodes(rcu) = (LBAfilterOn ? 4 : 3); break; + case 'H': itsRCUmodes(rcu) = itsSPW.rcumodeHBA(); break; + case '.': itsRCUmodes(rcu) = 0; break; + default: ASSERTSTR(false, "RCUinput #" << rcu << " contains illegal specification"); + } // switch + itsRCUuseFlags.set(itsRCUmodes(rcu)); + } // if + } // for } SubArray::~SubArray() { - LOG_DEBUG_STR("SubArray destructor"); - if (m_result[FRONT]) { - delete m_result[FRONT]; - } - if (m_result[BACK]) { - delete m_result[BACK]; - } + LOG_DEBUG_STR("SubArray destructor: " << itsName); + if (itsGains) { + delete itsGains; + } } -void SubArray::calibrate(CalibrationInterface* cal, ACC& acc) +SubArray& SubArray::operator=(const SubArray& that) { - ASSERT(m_result[FRONT]); - - if (cal) { - acc.setSelection(m_antenna_selection); - cal->calibrate(*this, acc, *m_result[FRONT]); - } - m_result[FRONT]->setDone(); + LOG_DEBUG_STR("SubArray operator= : " << that.itsName); + + if (this != &that) { + itsName = that.itsName; + itsAntennaSet = that.itsAntennaSet; + itsSPW = that.itsSPW; + itsRCUmask = that.itsRCUmask; + itsRCUmodes.resize(that.itsRCUmodes.shape()); + itsRCUmodes = that.itsRCUmodes.copy(); + itsGains = that.itsGains->clone(); + } + + return (*this); } -bool SubArray::getGains(AntennaGains*& cal, int buffer) +bool SubArray::usesRCUmode(int rcumode) const { - cal = 0; - ASSERT(buffer >= FRONT && buffer <= BACK && m_result[buffer]); - - cal = m_result[buffer]; + ASSERTSTR(rcumode >= 0 && rcumode <= NR_RCU_MODES, "RCUmode must be in the range 0..7"); - return m_result[buffer]->isDone(); + return (itsRCUuseFlags.test(rcumode)); } -void SubArray::abortCalibration() -{} - -const SpectralWindow& SubArray::getSPW() const +// +// RCUMask(rcumode) +// +RCUmask_t SubArray::RCUMask(uint rcumode) const { - return m_spw; + RCUmask_t result; + if (!itsRCUuseFlags.test(rcumode)) { // is this mode used anywhere? + return (result); + } + + // mode is used somewhere, is this the only used mode? + if (itsRCUuseFlags.count() == 1) { + return (itsRCUmask); + } + + // several modes are used, compose the right mask. + for (int rcu = 0; rcu < MAX_RCUS; rcu++) { + if (itsRCUmodes(rcu) == rcumode) { + result.set(rcu); + } + } // for + return (result); } -SubArray& SubArray::operator=(const SubArray& rhs) +void SubArray::updateGains(const AntennaGains& newGains) { - if (this != &rhs) { - // base-class assignment - AntennaArray::operator=(rhs); - - // assign spectral window - m_spw = rhs.getSPW(); - - // clear m_result pointers - this->m_result[FRONT] = 0; - this->m_result[BACK] = 0; - - itsRCUmask = rhs.itsRCUmask; - } - - return *this; + itsGains = newGains.clone(); } -bool SubArray::isDone() +// +// writeGains() +// +void SubArray::writeGains() { - ASSERT(m_result[FRONT]); - return m_result[FRONT]->isDone(); + time_t now = time(0); + struct tm* t = gmtime(&now); + char filename[PATH_MAX]; + + snprintf(filename, PATH_MAX, "%s_%04d%02d%02d_%02d%02d%02d_gain_%dx%dx%d.dat", + itsName.c_str(), + t->tm_year + 1900, t->tm_mon + 1, t->tm_mday, + t->tm_hour, t->tm_min, t->tm_sec, + itsGains->getGains().extent(firstDim), + itsGains->getGains().extent(secondDim), + itsGains->getGains().extent(thirdDim)); + LOG_DEBUG_STR("writeGains(" << name() << ") to " << filename); + + FILE* gainFile = fopen(filename, "w"); + + if (!gainFile) { + LOG_ERROR_STR("failed to open file: " << filename); + return; + } + + if (fwrite(itsGains->getGains().data(), sizeof(complex<double>), itsGains->getGains().size(), gainFile) != + (size_t)itsGains->getGains().size()) { + LOG_ERROR_STR("failed to write to file: " << filename); + } + + (void)fclose(gainFile); } -void SubArray::clearDone() +// +// print function for operator<< +// +ostream& SubArray::print (ostream& os) const { - ASSERT(m_result[FRONT]); - m_result[FRONT]->setDone(false); + os << "SubArray " << itsName << ":AntennaSet=" << itsAntennaSet << ", SPW={" << itsSPW; + os << "}, RCUmask=" << itsRCUmask << ", RCUmodeFlags=" << itsRCUuseFlags; + return (os); } +// +// ---------- pack and unpack routines ---------- +// size_t SubArray::getSize() const { return - MSH_size(m_name) - + MSH_size(m_geoloc) - + MSH_size(m_pos) - + MSH_size(m_rcuindex) + MSH_size(itsName) + + MSH_size(itsAntennaSet) + MSH_size(itsRCUmask) - + m_spw.getSize(); + + itsSPW.getSize(); } size_t SubArray::pack(char* buffer) const { - size_t offset = 0; + size_t offset = 0; - offset = MSH_pack(buffer, offset, m_name); - offset = MSH_pack(buffer, offset, m_geoloc); - offset = MSH_pack(buffer, offset, m_pos); - offset = MSH_pack(buffer, offset, m_rcuindex); - offset = MSH_pack(buffer, offset, itsRCUmask); - offset += m_spw.pack(buffer + offset); + offset = MSH_pack(buffer, offset, itsName); + offset = MSH_pack(buffer, offset, itsAntennaSet); + offset = MSH_pack(buffer, offset, itsRCUmask); + offset += itsSPW.pack(buffer + offset); - return offset; + return offset; } size_t SubArray::unpack(const char* buffer) { - size_t offset = 0; + size_t offset = 0; - offset = MSH_unpack(buffer, offset, m_name); - offset = MSH_unpack(buffer, offset, m_geoloc); - offset = MSH_unpack(buffer, offset, m_pos); - offset = MSH_unpack(buffer, offset, m_rcuindex); - offset = MSH_unpack(buffer, offset, itsRCUmask); - offset += m_spw.unpack(buffer + offset); + offset = MSH_unpack(buffer, offset, itsName); + offset = MSH_unpack(buffer, offset, itsAntennaSet); + offset = MSH_unpack(buffer, offset, itsRCUmask); + offset += itsSPW.unpack(buffer + offset); - return offset; + return offset; } // -------------------- SubArrayMap -------------------- size_t SubArrayMap::getSize() const { - return (MSH_size(*this)); + return MSH_size(*this); } size_t SubArrayMap::pack(char* buffer) const { - size_t offset = 0; - return MSH_pack(buffer, offset, (*this)); + size_t offset = 0; + return MSH_pack(buffer, offset, (*this)); } size_t SubArrayMap::unpack(const char* buffer) { - size_t offset = 0; - return MSH_unpack(buffer, offset, (*this)); + size_t offset = 0; + return MSH_unpack(buffer, offset, (*this)); } + } // namespace CAL +} // namespace LOFAR diff --git a/MAC/APL/PAC/CMakeLists.txt b/MAC/APL/PAC/CMakeLists.txt index 55b1004e3a6..16d6a0b9fa4 100644 --- a/MAC/APL/PAC/CMakeLists.txt +++ b/MAC/APL/PAC/CMakeLists.txt @@ -5,7 +5,7 @@ lofar_add_package(IBS_Protocol) lofar_add_package(Cal_Server) lofar_add_package(ITRFBeamServer) lofar_add_package(SHMInfo_Server) -lofar_add_package(ICAL_Protocol) +#lofar_add_package(ICAL_Protocol) # THIS DOES NOT BUILD ON KIS001 YET!!! #lofar_add_package(LBA_Calibration) #lofar_add_package(ITRFCalServer) diff --git a/MAC/APL/PAC/Cal_Server/src/ACMProxy.h b/MAC/APL/PAC/Cal_Server/src/ACMProxy.h index d6d9f560a72..693a86fc572 100644 --- a/MAC/APL/PAC/Cal_Server/src/ACMProxy.h +++ b/MAC/APL/PAC/Cal_Server/src/ACMProxy.h @@ -24,7 +24,7 @@ #ifndef ACMPROXY_H_ #define ACMPROXY_H_ -#include <APL/CAL_Protocol/ACC.h> +#include "ACC.h" #include <APL/RSP_Protocol/XCStatistics.h> #include <APL/RTCCommon/Timestamp.h> @@ -90,5 +90,5 @@ namespace LOFAR { }; // namespace CAL }; // namespace LOFAR - + #endif /* ACMPROXY_H_ */ diff --git a/MAC/APL/PAC/Cal_Server/src/CMakeLists.txt b/MAC/APL/PAC/Cal_Server/src/CMakeLists.txt index 131b47d927f..febfa6153c4 100644 --- a/MAC/APL/PAC/Cal_Server/src/CMakeLists.txt +++ b/MAC/APL/PAC/Cal_Server/src/CMakeLists.txt @@ -15,16 +15,16 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR}) lofar_add_library(caldevel Package__Version.cc - CalibrationAlgorithm.cc + AntennaArray.cc + AntennaArrayData.cc + ACC.cc CalibrationThread.cc DipoleModel.cc DipoleModelData.cc - RemoteStationCalibration.cc Source.cc SourceData.cc - SubArrays.cc - SubArraySubscription.cc) - + SubArrayMgr.cc) + lofar_add_bin_program(versioncal_server versioncal_server.cc) lofar_add_bin_program(CalServer CalServer.cc ACMProxy.cc) lofar_add_bin_program(calinfo calinfo.cc) diff --git a/MAC/APL/PAC/Cal_Server/src/CalServer.cc b/MAC/APL/PAC/Cal_Server/src/CalServer.cc index b9c9fef137c..62f1b7e1af5 100644 --- a/MAC/APL/PAC/Cal_Server/src/CalServer.cc +++ b/MAC/APL/PAC/Cal_Server/src/CalServer.cc @@ -29,9 +29,14 @@ #include <Common/lofar_bitset.h> #include <Common/Exception.h> #include <Common/Version.h> + +#include <ApplCommon/LofarDirs.h> + #include <ApplCommon/StationConfig.h> #include <ApplCommon/StationDatatypes.h> +#include <ApplCommon/AntennaSets.h> +#include <APL/APLCommon/AntennaField.h> #include <APL/RTCCommon/daemonize.h> #include <APL/CAL_Protocol/CAL_Protocol.ph> #include <APL/RSP_Protocol/RSP_Protocol.ph> @@ -39,10 +44,9 @@ #include <MACIO/MACServiceInfo.h> #include "CalServer.h" -#include "SubArrays.h" +#include "SubArrayMgr.h" #include "SubArraySubscription.h" -#include "RemoteStationCalibration.h" -#include "CalibrationAlgorithm.h" +//#include "CalibrationAlgorithm.h" #include <Cal_Server/Package__Version.h> #ifdef USE_CAL_THREAD @@ -65,12 +69,13 @@ using namespace std; using namespace blitz; -using namespace LOFAR; -using namespace CAL; -using namespace RTC; -using namespace RSP_Protocol; -using namespace CAL_Protocol; -using namespace GCF::TM; + +namespace LOFAR { + using namespace RTC; + using namespace RSP_Protocol; + using namespace CAL_Protocol; + using namespace GCF::TM; + namespace CAL { #define NPOL 2 @@ -81,16 +86,15 @@ Exception::TerminateHandler t(Exception::terminate); // static pointer used for signal handler. static CalServer* thisCalServer = 0; - + // // CalServer constructor // CalServer::CalServer(const string& name, ACCs& accs) : GCFTask((State)&CalServer::initial, name), m_accs(accs), - m_cal(0), m_converter(0), - m_sampling_frequency(0.0), + itsClockSetting(0.0), m_n_rspboards(0), m_n_rcus(0), itsHasSecondRing (0), @@ -137,7 +141,6 @@ CalServer::CalServer(const string& name, ACCs& accs) // CalServer::~CalServer() { - if (m_cal) delete m_cal; if (m_converter) delete m_converter; #ifdef USE_CAL_THREAD if (m_calthread) delete m_calthread; @@ -178,7 +181,7 @@ GCFEvent::TResult CalServer::finishing_state(GCFEvent& event, GCFPortInterface case F_ENTRY: { // turn off 'old' timers itsCheckTimer->cancelAllTimers(); - SubArray::RCUmask_t rcus2switchOff; + RCUmask_t rcus2switchOff; rcus2switchOff.reset(); for (uint i = 0; i < m_n_rcus; i++) { rcus2switchOff.set(i); @@ -190,11 +193,11 @@ GCFEvent::TResult CalServer::finishing_state(GCFEvent& event, GCFPortInterface case F_TIMER: { GCFScheduler::instance()->stop(); } break; - + default: LOG_DEBUG("finishing_state, default"); break; - } + } return (GCFEvent::HANDLED); } @@ -223,10 +226,10 @@ void CalServer::remove_client(GCFPortInterface* port) while (iter != end) { if (iter->second == port) { // stop subarray if it is still there. - SubArray* subarray = m_subarrays.getByName(iter->first); + SubArray* subarray = itsSubArrays.getByName(iter->first); if (subarray) { _disableRCUs(subarray); - m_subarrays.schedule_remove(subarray); + itsSubArrays.scheduleRemove(subarray); } // add to dead list. @@ -282,11 +285,9 @@ GCFEvent::TResult CalServer::initial(GCFEvent& e, GCFPortInterface& port) itsPowerOffDelay = globalParameterSet()->getTime("CalServer.PowerOffDelay", 180); LOG_INFO_STR("PowerOff delay = " << itsPowerOffDelay << " seconds."); - // Setup calibration algorithm - m_cal = new RemoteStationCalibration(m_sources, m_dipolemodels, *m_converter); #ifdef USE_CAL_THREAD // Setup calibration thread - m_calthread = new CalibrationThread(&m_subarrays, m_cal, m_globallock, itsDataDir); + m_calthread = new CalibrationThread(&itsSubArrays, m_globallock, itsDataDir); pthread_mutex_unlock(&m_globallock); // unlock global lock #endif @@ -346,10 +347,10 @@ GCFEvent::TResult CalServer::initial(GCFEvent& e, GCFPortInterface& port) exit(EXIT_FAILURE); } - // get clock value and convert to MHz - m_sampling_frequency = getclockack.clock * (uint32)1.0e6; + // get clock value in MHz + itsClockSetting = getclockack.clock; - LOG_INFO_STR("Initial sampling frequency: " << m_sampling_frequency); + LOG_INFO_STR("Initial sampling frequency: " << itsClockSetting << " MHz"); // subscribe to clock change updates RSPSubclockEvent subclock; @@ -454,8 +455,8 @@ GCFEvent::TResult CalServer::enabled(GCFEvent& e, GCFPortInterface& port) case RSP_UPDCLOCK: { RSPUpdclockEvent updclock(e); // use new sampling frequency - m_sampling_frequency = updclock.clock * (uint32)1.0e6; - LOG_INFO_STR("New sampling frequency: " << m_sampling_frequency); + itsClockSetting = updclock.clock * (uint32)1.0e6; + LOG_INFO_STR("New sampling frequency: " << itsClockSetting); } break; @@ -489,7 +490,7 @@ GCFEvent::TResult CalServer::enabled(GCFEvent& e, GCFPortInterface& port) } } break; - + case F_CONNECTED: { LOG_INFO(formatString("CONNECTED: port '%s'", port.getName().c_str())); } @@ -497,7 +498,7 @@ GCFEvent::TResult CalServer::enabled(GCFEvent& e, GCFPortInterface& port) case F_TIMER: { // first check Poweroff timer - SubArray::RCUmask_t rcus2switchOff; + RCUmask_t rcus2switchOff; rcus2switchOff.reset(); time_t now(time(0L)); for (uint i = 0; i < m_n_rcus; i++) { // loop over all rcus @@ -532,11 +533,11 @@ GCFEvent::TResult CalServer::enabled(GCFEvent& e, GCFPortInterface& port) (void)m_calthread->join(); #endif - m_subarrays.mutex_lock(); + itsSubArrays.mutex_lock(); undertaker(); // destroy dead clients, done here to prevent possible use of closed port - m_subarrays.undertaker(); // remove subarrays scheduled for deletion - m_subarrays.creator(); // bring new subarrays to life - m_subarrays.mutex_unlock(); + itsSubArrays.removeDeadArrays(); // remove subarrays scheduled for deletion + itsSubArrays.activateArrays(); // bring new subarrays to life + itsSubArrays.mutex_unlock(); if (GET_CONFIG("CalServer.WriteACCToFile", i)) { write_acc(); @@ -548,17 +549,13 @@ GCFEvent::TResult CalServer::enabled(GCFEvent& e, GCFPortInterface& port) m_calthread->run(); #else bool writeGains(GET_CONFIG("CalServer.WriteGainsToFile", i) == 1); - if (GET_CONFIG("CalServer.DisableCalibration", i)) { - m_subarrays.calibrate(0, m_accs.getFront(), writeGains, itsDataDir); - } else { - m_subarrays.calibrate(m_cal, m_accs.getFront(), writeGains, itsDataDir); - } - m_subarrays.updateAll(); + itsSubArrays.calibrate(m_accs.getFront(), writeGains, itsDataDir); + itsSubArrays.updateAll(); #endif } #ifdef USE_CAL_THREAD - m_subarrays.updateAll(); + itsSubArrays.updateAll(); #endif } break; @@ -624,9 +621,9 @@ GCFEvent::TResult CalServer::handle_cal_start(GCFEvent& e, GCFPortInterface &por ack.name = start.name; // find parent AntennaArray - const AntennaArray* parent = m_arrays.getByName(start.parent); + const AntennaArray* parent = m_arrays.getByName(start.antennaSet); - if (m_subarrays.getByName(start.name)) { + if (itsSubArrays.getByName(start.name)) { LOG_ERROR_STR("A subarray with name='" << start.name << "' has already been registered."); ack.status = ERR_ALREADY_REGISTERED; @@ -636,10 +633,10 @@ GCFEvent::TResult CalServer::handle_cal_start(GCFEvent& e, GCFPortInterface &por } else if (!parent) { // parent not found, set error status - LOG_ERROR_STR("Parent array '" << start.parent << "' not found."); + LOG_ERROR_STR("Parent array '" << start.antennaSet << "' not found."); ack.status = ERR_NO_PARENT; - } else if (start.subset.count() == 0) { + } else if (start.rcuMask.count() == 0) { // empty selection LOG_ERROR("Empty antenna selection not allowed."); ack.status = ERR_NO_ANTENNAS; @@ -662,7 +659,7 @@ GCFEvent::TResult CalServer::handle_cal_start(GCFEvent& e, GCFPortInterface &por int maxI(positions.extent(firstDim)*positions.extent(secondDim)); for (int i = 0; i < maxI; i++) { // subset is one-dimensional (receiver), select is two-dimensional (antenna, polarization) - if (start.subset[i]) { + if (start.rcuMask[i]) { select(i/2,i%2) = true; } } @@ -689,54 +686,50 @@ GCFEvent::TResult CalServer::handle_cal_start(GCFEvent& e, GCFPortInterface &por { LOG_INFO("ACC shape and parent array positions shape don't match."); LOG_ERROR_STR("ACC.shape=" << m_accs.getFront().getACC().shape()); - LOG_ERROR_STR("'" << start.parent << "'.shape=" << positions.shape()); + LOG_ERROR_STR("'" << start.antennaSet << "'.shape=" << positions.shape()); LOG_ERROR_STR("Expecting AntenneArray with " << m_accs.getFront().getACC().extent(fourthDim) << " antennas."); ack.status = ERR_RANGE; } - else + else #endif - if ((start.subset & invalidmask).any()) { + if ((start.rcuMask & invalidmask).any()) { LOG_INFO("CAL_START: Invalid receiver subset."); ack.status = ERR_RANGE; } else { - // create subarray to calibrate - SubArray* subarray = new SubArray(start.name, - parent->getGeoLoc(), - positions.copy(), - select, - m_sampling_frequency, - start.rcumode()(0).getNyquistZone(), - GET_CONFIG("CalServer.N_SUBBANDS", i), - start.rcumode()(0).getRaw()); - - if (m_subarrays.conflicting(subarray)) { - // error is already printed - ack.status = ERR_NO_PARENT; - } - else { - m_subarrays.schedule_add(subarray); - - // calibration will start within one second - // set the spectral inversion right - RSPSetbypassEvent specInvCmd; - bool SIon(start.rcumode()(0).getNyquistZone() == 2);// on or off? - specInvCmd.timestamp = Timestamp(0,0); - specInvCmd.rcumask = start.subset; - specInvCmd.settings().resize(1); - specInvCmd.settings()(0).setXSI(SIon); - specInvCmd.settings()(0).setYSI(SIon); - LOG_DEBUG_STR("NyquistZone = " << start.rcumode()(0).getNyquistZone() - << " setting spectral inversion " << ((SIon) ? "ON" : "OFF")); - itsRSPDriver->send(specInvCmd); - - RCUSettings rcu_settings; - rcu_settings().resize(1); - rcu_settings()(0) = start.rcumode()(0); - - _enableRCUs(subarray, rcu_settings, SCHEDULING_DELAY + 4); - } // conflict? + RCUmask_t validRCUs(start.rcuMask & globalAntennaSets()->RCUallocation(start.antennaSet)); + + // create subarray to calibrate + RCUSettings::Control RCUcontrol; // for conversion from rcumode to LBAfilter setting + RCUcontrol.setMode((RCUSettings::Control::RCUMode)start.rcumode); + SubArray* subarray = new SubArray(start.name, + start.antennaSet, + validRCUs, + RCUcontrol.LBAfilter(), + itsClockSetting * 1.0e6, + RCUcontrol.getNyquistZone()); + + itsSubArrays.scheduleAdd(subarray); + + // calibration will start within one second + // set the spectral inversion right + RSPSetbypassEvent specInvCmd; + bool SIon(RCUcontrol.getNyquistZone() == 2);// on or off? + specInvCmd.timestamp = Timestamp(0,0); + specInvCmd.rcumask = start.rcuMask; + specInvCmd.settings().resize(1); + specInvCmd.settings()(0).setXSI(SIon); + specInvCmd.settings()(0).setYSI(SIon); + LOG_DEBUG_STR("NyquistZone = " << RCUcontrol.getNyquistZone() + << " setting spectral inversion " << ((SIon) ? "ON" : "OFF")); + itsRSPDriver->send(specInvCmd); + + RCUSettings rcu_settings; + rcu_settings().resize(1); + rcu_settings()(0).setMode(start.rcumode); + + _enableRCUs(subarray, rcu_settings, SCHEDULING_DELAY + 4); } } port.send(ack); // send ack @@ -757,9 +750,9 @@ GCFEvent::TResult CalServer::handle_cal_stop(GCFEvent& e, GCFPortInterface &port ack.status = CAL_Protocol::CAL_SUCCESS; // return success: don't bother client with our admin port.send(ack); - _disableRCUs(m_subarrays.getByName(stop.name)); // switch off unused RCUs + _disableRCUs(itsSubArrays.getByName(stop.name)); // switch off unused RCUs - m_subarrays.schedule_remove(stop.name); // stop calibration + itsSubArrays.scheduleRemove(stop.name); // stop calibration // remove subarray-port entry from the map. map<string, GCFPortInterface*>::iterator iter = m_clients.begin(); @@ -778,78 +771,55 @@ GCFEvent::TResult CalServer::handle_cal_stop(GCFEvent& e, GCFPortInterface &port // // handle_cal_subscribe(event, port) // -GCFEvent::TResult CalServer::handle_cal_subscribe(GCFEvent& e, GCFPortInterface &port) +GCFEvent::TResult CalServer::handle_cal_subscribe(GCFEvent& event, GCFPortInterface &port) { - GCFEvent::TResult status = GCFEvent::HANDLED; - - CALSubscribeEvent subscribe(e); + CALSubscribeEvent subscribe(event); CALSubscribeackEvent ack; - ack.status = CAL_Protocol::CAL_SUCCESS; + ack.name = subscribe.name; // get subarray by name - SubArray* subarray = m_subarrays.getByName(subscribe.name); - + SubArray* subarray = itsSubArrays.getByName(subscribe.name); if (subarray) { - // create subscription - SubArraySubscription* subscription = new SubArraySubscription(subarray, - subscribe.subbandset, - port); - - ack.handle = (CAL_Protocol::memptr_t)subscription; - subarray->attach(subscription); // attach subscription to the subarray + ack.status = CAL_SUCCESS; ack.subarray = *subarray; // return subarray positions - // don't register subarray in cal_subscribe - // it has already been registerd in cal_start - LOG_DEBUG_STR("Subscription succeeded: " << subscribe.name); + // create subscription + itsSubArrays.addSubscription(subscribe.name, &port); + + LOG_INFO_STR("Subscription succeeded: " << subscribe.name); } else { ack.status = ERR_NO_SUBARRAY; - ack.handle = 0; - //memset(&ack.subarray, 0, sizeof(SubArray)); - // doesn't work with gcc-3.4 ack.subarray = SubArray(); - (void)new((void*)&ack.subarray) SubArray(); - - LOG_WARN_STR("Subarray not found: " << subscribe.name); + LOG_INFO_STR("Subarray not found: " << subscribe.name); } port.send(ack); - return status; + return (GCFEvent::HANDLED); } // // handle_cal_unsubscribe(event,port) // -GCFEvent::TResult CalServer::handle_cal_unsubscribe(GCFEvent& e, GCFPortInterface &port) +GCFEvent::TResult CalServer::handle_cal_unsubscribe(GCFEvent& event, GCFPortInterface &port) { - GCFEvent::TResult status = GCFEvent::HANDLED; - - CALUnsubscribeEvent unsubscribe(e); - + CALUnsubscribeEvent unsubscribe(event); // create ack CALUnsubscribeackEvent ack; - ack.name = unsubscribe.name; - ack.handle = unsubscribe.handle; - ack.status = CAL_Protocol::CAL_SUCCESS; - - // find associated subarray - SubArray* subarray = m_subarrays.getByName(unsubscribe.name); - if (subarray) { - // detach subscription, this destroys the subscription - subarray->detach((SubArraySubscription*)unsubscribe.handle); + ack.name = unsubscribe.name; - // handle is no longer valid + if (itsSubArrays.removeSubscription(unsubscribe.name)) { + ack.status = CAL_SUCCESS; LOG_INFO_STR("Subscription deleted: " << unsubscribe.name); } else { ack.status = ERR_NO_SUBARRAY; - LOG_WARN_STR("Subscription failed. Subbarray not found: " << unsubscribe.name); + LOG_INFO_STR("Unsubscribe failed. Subbarray not found: " << unsubscribe.name); } port.send(ack); - return status; + return (GCFEvent::HANDLED); } // @@ -862,7 +832,7 @@ GCFEvent::TResult CalServer::handle_cal_getsubarray(GCFEvent& e, GCFPortInterfac CALGetsubarrayEvent request(e); CALGetsubarrayackEvent ack; ack.status = CAL_Protocol::CAL_SUCCESS; - ack.subarraymap = m_subarrays.getSubArrays(request.subarrayname); // let SubArrays do the job + ack.subarraymap = itsSubArrays.getSubArrays(request.subarrayname); // let SubArrays do the job // correct status is name was given but nothing was found. if (!request.subarrayname.empty() && ack.subarraymap.size() == 0) { @@ -877,11 +847,11 @@ GCFEvent::TResult CalServer::handle_cal_getsubarray(GCFEvent& e, GCFPortInterfac // // (subarray*) // -void CalServer::_enableRCUs(SubArray* subarray, RCUSettings rcu_settings, int delay) +void CalServer::_enableRCUs(SubArray* subarray, int delay) { // increment the usecount of the receivers - SubArray::RCUmask_t rcuMask = subarray->getRCUMask(); - SubArray::RCUmask_t rcus2switchOn; + RCUmask_t rcuMask = subarray->RCUMask(); + RCUmask_t rcus2switchOn; rcus2switchOn.reset(); Timestamp timeStamp; for (uint r = 0; r < m_n_rcus; r++) { @@ -903,7 +873,7 @@ void CalServer::_enableRCUs(SubArray* subarray, RCUSettings rcu_settings, int de enableCmd.timestamp = timeStamp; enableCmd.rcumask = rcus2switchOn; enableCmd.settings().resize(1); - enableCmd.settings()(0) = rcu_settings()(0); + //enableCmd.settings()(0) = rcu_settings()(0); enableCmd.settings()(0).setEnable(true); sleep (1); LOG_INFO("Enabling some rcu's because they are used for the first time"); @@ -911,7 +881,7 @@ void CalServer::_enableRCUs(SubArray* subarray, RCUSettings rcu_settings, int de } // when the lbl inputs are selected swap the X and the Y. - int rcumode(subarray->getSPW().rcumode()); + int rcumode(subarray->SPW().rcumode()); if (rcumode == 1 || rcumode == 2) { // LBLinput used? LOG_INFO("LBL inputs are used, swapping X and Y RCU's"); RSPSetswapxyEvent swapCmd; @@ -929,13 +899,13 @@ void CalServer::_disableRCUs(SubArray* subarray) { // decrement the usecount of the receivers and build mask of receiver that may be disabled. bool allSwitchedOff(true); - SubArray::RCUmask_t rcus2switchOff; + RCUmask_t rcus2switchOff; rcus2switchOff.reset(); Timestamp powerOffTime; powerOffTime.setNow(itsPowerOffDelay); if (subarray) { // when no subarray is defined skip this loop: switch all rcus off. - SubArray::RCUmask_t rcuMask = subarray->getRCUMask(); + RCUmask_t rcuMask = subarray->RCUMask(); for (uint r = 0; r < m_n_rcus; r++) { if (rcuMask.test(r)) { if (--itsRCUcounts[r] == 0) { @@ -979,7 +949,7 @@ void CalServer::_disableRCUs(SubArray* subarray) // // _powerdownRCUs(rcus2switchOff) // -void CalServer::_powerdownRCUs(SubArray::RCUmask_t rcus2switchOff) +void CalServer::_powerdownRCUs(RCUmask_t rcus2switchOff) { // anything to disable? Tell the RSPDriver. if (rcus2switchOff.any()) { @@ -1051,7 +1021,7 @@ GCFEvent::TResult CalServer::handle_cal_getsubarray(GCFEvent& e, GCFPortInterfac ack.status = CAL_Protocol::CAL_SUCCESS; // find associated subarray - SubArray* subarray = m_subarrays.getByName(getsubarray.name); + SubArray* subarray = itsSubArrays.getByName(getsubarray.name); if (subarray) { // return antenna positions @@ -1170,3 +1140,6 @@ int main(int argc, char** argv) return 0; } + + }; // CAL +}; / LOFAR diff --git a/MAC/APL/PAC/Cal_Server/src/CalServer.h b/MAC/APL/PAC/Cal_Server/src/CalServer.h index c2439fa9236..fbbf1c6f10c 100644 --- a/MAC/APL/PAC/Cal_Server/src/CalServer.h +++ b/MAC/APL/PAC/Cal_Server/src/CalServer.h @@ -27,12 +27,14 @@ #include <Common/lofar_list.h> #include <Common/lofar_map.h> #include <Common/lofar_string.h> -#include <APL/CAL_Protocol/ACC.h> +#include <ApplCommon/StationDatatypes.h> +#include "ACC.h" #include <APL/CAL_Protocol/SubArray.h> #include <APL/RSP_Protocol/RCUSettings.h> #include "Source.h" +#include "AntennaArray.h" #include "DipoleModel.h" -#include "SubArrays.h" +#include "SubArrayMgr.h" #include <GCF/TM/GCF_Control.h> #include <AMCBase/ConverterClient.h> @@ -43,7 +45,6 @@ namespace LOFAR { using GCF::TM::GCFTCPPort; using GCF::TM::GCFPortInterface; using GCF::TM::GCFTimerPort; - namespace CAL { // forward declarations @@ -71,7 +72,7 @@ public: void remove_client(GCFPortInterface* port); // increment RCU usagecounters and enable newly used RCUs - void _enableRCUs(SubArray* subarray, RSP_Protocol::RCUSettings rcu_settings, int delay); + void _enableRCUs(SubArray* subarray, int delay); // decrement RCU usagecounters and disable unused RCUs void _disableRCUs(SubArray* subarray); @@ -85,7 +86,7 @@ public: static void sigintHandler(int signum); void finish(); GCFEvent::TResult finishing_state(GCFEvent& event, GCFPortInterface& port); - + /*@}*/ /*@{*/ @@ -96,14 +97,14 @@ public: GCFEvent::TResult handle_cal_unsubscribe(GCFEvent& e, GCFPortInterface &port); GCFEvent::TResult handle_cal_getsubarray(GCFEvent& e, GCFPortInterface &port); /*@}*/ - + // Write ACC to file if configured to do so. void write_acc(); // Helper functions bool _dataOnRing (uint ringNr) const; void _updateDataStream(uint delay); - void _powerdownRCUs (SubArray::RCUmask_t rcus2switchOff); + void _powerdownRCUs (RCUmask_t rcus2switchOff); private: // ----- DATA MEMBERS ----- @@ -113,7 +114,7 @@ private: Sources m_sources; // source catalog (read from file) DipoleModels m_dipolemodels; // dipole model (read from file) - SubArrays m_subarrays; // the subarrays (created by clients) + SubArrayMgr itsSubArrays; // the subarrays (created by clients) ACCs& m_accs; // front and back ACC buffers (received from ACMServer) CalibrationAlgorithm* m_cal; // pointer to the calibration algorithm to use @@ -121,7 +122,7 @@ private: AMC::ConverterClient* m_converter; // interface for coordinate conversion (Astronomical Measures Conversion) // Current sampling frequency of the system. - double m_sampling_frequency; + int itsClockSetting; // remember number of RSP boards and number of rcus uint m_n_rspboards; diff --git a/MAC/APL/PAC/Cal_Server/src/CalibrationAlgorithm.h b/MAC/APL/PAC/Cal_Server/src/CalibrationAlgorithm.h index ad4ebfd137f..c807238c606 100644 --- a/MAC/APL/PAC/Cal_Server/src/CalibrationAlgorithm.h +++ b/MAC/APL/PAC/Cal_Server/src/CalibrationAlgorithm.h @@ -24,7 +24,7 @@ #ifndef CALIBRATIONALGORITHM_H_ #define CALIBRATIONALGORITHM_H_ -#include <APL/CAL_Protocol/CalibrationInterface.h> +//#include "CalibrationInterface.h" #include "Source.h" #include "DipoleModel.h" #include <AMCBase/Converter.h> @@ -32,16 +32,16 @@ namespace LOFAR { namespace CAL { - class CalibrationAlgorithm : public CalibrationInterface + class CalibrationAlgorithm //: public CalibrationInterface { public: CalibrationAlgorithm(const Sources& sources, DipoleModels& dipolemodels, AMC::Converter& converter); virtual ~CalibrationAlgorithm(); - + virtual const Sources& getSources() { return m_sources; } virtual DipoleModels& getDipoleModels() { return m_dipolemodels; } virtual AMC::Converter& getConverter() { return m_converter; } - + private: const Sources& m_sources; DipoleModels& m_dipolemodels; diff --git a/MAC/APL/PAC/Cal_Server/src/CalibrationThread.cc b/MAC/APL/PAC/Cal_Server/src/CalibrationThread.cc index bd6e5bc8630..6b356647621 100644 --- a/MAC/APL/PAC/Cal_Server/src/CalibrationThread.cc +++ b/MAC/APL/PAC/Cal_Server/src/CalibrationThread.cc @@ -28,21 +28,19 @@ #include <APL/RTCCommon/PSAccess.h> #include <pthread.h> -#include "SubArrays.h" +#include "SubArrayMgr.h" #include "CalibrationThread.h" using namespace LOFAR; using namespace CAL; -CalibrationThread::CalibrationThread(SubArrays* subarrays, - CalibrationInterface* cal, +CalibrationThread::CalibrationThread(SubArrayMgr* subarrays, pthread_mutex_t& globallock, const string& dataDir) : - m_subarrays (subarrays), - m_cal (cal), - m_acc (0), + m_subarrays (subarrays), + m_acc (0), itsDataDir (dataDir), - m_thread (0), + m_thread (0), m_globallock(globallock) { } @@ -84,11 +82,7 @@ void* CalibrationThread::thread_main(void* thisthread) if (thread->m_acc) { bool writeGains(GET_CONFIG("CalServer.WriteGainsToFile", i) == 1); - if (GET_CONFIG("CalServer.DisableCalibration", i)) { - thread->m_subarrays->calibrate(0, *thread->m_acc, writeGains, thread->itsDataDir); - } else { - thread->m_subarrays->calibrate(thread->m_cal, *thread->m_acc, writeGains, thread->itsDataDir); - } + thread->m_subarrays->calibrate(*thread->m_acc, writeGains); } thread->m_acc->readUnlock(); @@ -111,6 +105,6 @@ int CalibrationThread::join() LOG_WARN_STR("pthread_join returned " << returncode); } return returncode; -} +} #endif diff --git a/MAC/APL/PAC/Cal_Server/src/CalibrationThread.h b/MAC/APL/PAC/Cal_Server/src/CalibrationThread.h index a8a163f70dd..717fbe25f4f 100644 --- a/MAC/APL/PAC/Cal_Server/src/CalibrationThread.h +++ b/MAC/APL/PAC/Cal_Server/src/CalibrationThread.h @@ -26,19 +26,17 @@ #ifdef USE_CAL_THREAD #include <pthread.h> -#include <APL/CAL_Protocol/ACC.h> +#include "ACC.h" namespace LOFAR { namespace CAL { - class SubArrays; - class CalibrationInterface; + class SubArrayMgr; class CalibrationThread { public: - CalibrationThread(SubArrays* subarrays, - CalibrationInterface* cal, + CalibrationThread(SubArrayMgr* subarrays, pthread_mutex_t& globallock, const string& dataDir); virtual ~CalibrationThread(); @@ -48,10 +46,9 @@ namespace LOFAR { void run(); static void* thread_main(void* thisthread); int join(); - + private: - SubArrays* m_subarrays; - CalibrationInterface* m_cal; + SubArrayMgr* m_subarrays; ACC* m_acc; string itsDataDir; diff --git a/MAC/APL/PAC/Cal_Server/src/RemoteStationCalibration.cc b/MAC/APL/PAC/Cal_Server/src/RemoteStationCalibration.cc index bb9fc698002..91bd02ca0f1 100644 --- a/MAC/APL/PAC/Cal_Server/src/RemoteStationCalibration.cc +++ b/MAC/APL/PAC/Cal_Server/src/RemoteStationCalibration.cc @@ -60,7 +60,6 @@ RemoteStationCalibration::RemoteStationCalibration(Sources& sources, DipoleModel itsBaseLineRestriction = 4.0; // baseline restriction in wavelengths itsMaxBaseLine = 20.0; // max base line in meters itsSpeedOfLight = 2.99792e8; - } //--------------------------------------------------------------------------------------------------------------------- @@ -83,17 +82,17 @@ void RemoteStationCalibration::calibrate(const SubArray& subarray, ACC& acc, Ant // double start_freq = 10e6; double stop_freq = 90e6; - + const SpectralWindow& spw = subarray.getSPW(); // get spectral window DipoleModels& dipolemodels = getDipoleModels(); // get dipole models const Sources& sources = getSources(); // get sky model const DipoleModel* dipolemodel = dipolemodels.getByName("LBAntenna"); const Array<double, 1> geoLoc( subarray.getGeoLoc() ); const Array<double, 3> antennaPos( subarray.getAntennaPos() ); - + time_t obstime; struct tm* timeinfo; - + time(&obstime); timeinfo = gmtime(&obstime); logfile << "filling timeinfo" << endl; @@ -102,11 +101,11 @@ void RemoteStationCalibration::calibrate(const SubArray& subarray, ACC& acc, Ant timeinfo->tm_mday = 18; timeinfo->tm_hour = 8; timeinfo->tm_min = 32; - timeinfo->tm_sec = 29; + timeinfo->tm_sec = 29; logfile << "setting obstime " << asctime(timeinfo) << endl; obstime = mktime( timeinfo ); - logfile << "setting obstime done " << obstime << endl; - + logfile << "setting obstime done " << obstime << endl; + if (!dipolemodel) { LOG_FATAL("Failed to load dipolemodel 'LBAntenna'"); exit(EXIT_FAILURE); @@ -118,51 +117,51 @@ void RemoteStationCalibration::calibrate(const SubArray& subarray, ACC& acc, Ant LOG_INFO_STR("calibrate: subarray name=" << subarray.getName()); LOG_INFO_STR("calibrate: num_antennas=" << subarray.getNumAntennas()); - + int npolarizations = 2; c_start1 = clock(); for (int pol = 0; pol < npolarizations; pol++) { logfile << "calibrate: working on polarization " << pol << " of " << npolarizations << endl; - + blitz::Array<bool, 1> isClean(detectRFI(acc, pol, spw.getSubbandWidth()) ); logfile << "detectRFI done" << " time=" << getClock() << endl; - + blitz::Array<bool, 1> isInBand(setPassBand(start_freq, stop_freq, spw) ); logfile << "setPassBand done" << " time=" << getClock() << endl; - - + + blitz::Array<bool, 1> isSelected( where(isClean && isInBand, true, false)); logfile << "isSelected done" << " time=" << getClock() << endl; - //bmToFile(isSelected,"isSelected",fileFormat); - - blitz::Array<bool, 1> isGood(detectBadElements(acc, pol, isSelected)); + //bmToFile(isSelected,"isSelected",fileFormat); + + blitz::Array<bool, 1> isGood(detectBadElements(acc, pol, isSelected)); logfile << "detectBadElements done" << " time=" << getClock() << endl; - + for (int sb = 0; sb < spw.getNumSubbands(); sb++) { logfile << "sb = " << sb << " time=" << getClock() << endl; // process only if sb = clean and in band if (isSelected(sb)) { //logfile << "calibrate: working on subband " << sb + 1 << " of " << spw.getNumSubbands() << endl; - + double freq = sb * spw.getSubbandWidth() + spw.getSamplingFrequency() * (spw.getNyquistZone() - 1); logfile << "freq = " << freq << endl; - + Timestamp acmtime; Array<complex<double>, 2> acm = acc.getACM(sb, pol, pol, acmtime); - + //obstime += sb; logfile << "obstime=" << obstime+sb << endl; acmtime = Timestamp(obstime+sb, 0); - + // construct local sky model (LSM) vector<Source> LSM = make_local_sky_model(sources, geoLoc, acmtime); logfile << "LSM done" << " time=" << getClock() << endl; - + // geef alleen geldige elementen door // antennePos, acm int nelem = acm.extent(firstDim); - int gelem = sum(isGood); + int gelem = sum(isGood); Array<double, 3> goodAntPos(gelem,antennaPos.extent(secondDim),antennaPos.extent(thirdDim)); Array<complex<double>, 2> goodAcm(gelem,gelem); int id1 = 0; @@ -173,115 +172,115 @@ void RemoteStationCalibration::calibrate(const SubArray& subarray, ACC& acc, Ant for (int idx2 = 0; idx2 < nelem; idx2++) { if (isGood(idx2)) { goodAcm(id1,id2) = acm(idx1,idx2); - id2++; + id2++; } } id1++; - } - } - + } + } + Array<bool, 2> restriction( setRestriction(goodAntPos, freq) ); logfile << "restriction done" << " time=" << getClock() << endl; - + Array<complex<double>, 2> asrc(computeA(LSM, goodAntPos, freq)); logfile << "asrc done" << " time=" << getClock() << endl; - + Array<double, 1> flux( computeFlux(goodAcm, asrc, restriction)); - logfile << "flux done" << " time=" << getClock() << endl; - + logfile << "flux done" << " time=" << getClock() << endl; + int maxIter = 5; double diffStop = 0.001; Array<complex<double>, 1> gain( computeGain(goodAcm, asrc, flux, restriction, maxIter, diffStop)); logfile << "gain done" << " time=" << getClock() << endl; - + Array<std::complex<double>, 1> cal( conj(1.0 / gain)); //bmToFile(cal,"cal",fileFormat); - + // checkCal() ?? //Array<complex<double>, 1> cal0( checkCal(cal, isGood, spw.getNumSubbands(), goodAcm.rows())); - + // convert computed gain to full array, and set bad elemts to zero int ngains = gains.getGains().extent(firstDim); Array<complex<double>, 1> gain_full(ngains); gain_full = 0.0; id1 = 0; for (int idx1 = 0;idx1 < ngains; idx1++) { - if (isGood(idx1)) { + if (isGood(idx1)) { gain_full(idx1) = gain(id1); id1++; } } - + ASSERT(gains.getGains().extent(firstDim) == gain_full.extent(firstDim)); gains.getGains()(Range::all(), pol, sb) = gain_full; - + logfile << "gains[" << sb << "]=" << gain_full << endl; - + } else { //logfile << "calibrate: subband " << sb+1 << " contains RFI and is not calibrated" << endl; } } } //interpolate_bad_subbands(); -#endif +#endif gains.setDone(true); // when finished } //--------------------------------------------------------------------------------------------------------------------- #if 0 // detects RFI on each channel // returns true if clean -// -blitz::Array<bool,1> RemoteStationCalibration::detectRFI(ACC& acc, int pol, double N) +// +blitz::Array<bool,1> RemoteStationCalibration::detectRFI(ACC& acc, int pol, double N) { // TODO put in conf file double rfiThresholdFactor = 10.0; - + int nelem = acc.getNAntennas(); int nch = acc.getNSubbands(); Timestamp acmtime; - + blitz::Array<double,1> teststat(nch); - + teststat = 0; logfile << "nelem=" << nelem << " N=" << N << endl; double threshold = sqrt(2.0 * nelem / N); logfile << "threshold=" << threshold << endl; - + for (int ch = 0; ch < nch; ch++) { blitz::Array<std::complex<double>,2> Rhat(acc.getACM(ch, pol, pol, acmtime)); blitz::Array<std::complex<double>,1> diagRhat( bmDiag(Rhat) ); blitz::Array<std::complex<double>,2> Rwhite(Rhat(tensor::i, tensor::j) / sqrt( diagRhat(tensor::i) * diagRhat(tensor::j))); - + teststat(ch) = pow2( bmNorm(Rwhite, fro)); - } + } //bmToFile(teststat,"teststat",fileFormat); - + blitz::Array<double,1> teststatL(nch); teststatL = 0; teststatL(Range(0,nch-2)) = teststat(Range(1,nch-1)); - + blitz::Array<double,1> teststatR(nch); teststatR = 0; - teststatR(Range(1,nch-1)) = teststat(Range(0,nch-2)); - + teststatR(Range(1,nch-1)) = teststat(Range(0,nch-2)); + blitz::Array<double,2> absLR(nch,2); absLR(Range::all(),0) = abs(teststat(Range::all()) - teststatL(Range::all())); - absLR(Range::all(),1) = abs(teststat(Range::all()) - teststatR(Range::all())); + absLR(Range::all(),1) = abs(teststat(Range::all()) - teststatR(Range::all())); blitz::Array<double,1> maxLR( max(absLR(tensor::i,tensor::j),tensor::j) ); //bmToFile(maxLR,"maxLR",fileFormat); - + blitz::Array<bool,1> clean( where(maxLR < (rfiThresholdFactor * threshold), true, false)); //bmToFile(clean,"cleanBands",fileFormat); - return (clean); + return (clean); } //--------------------------------------------------------------------------------------------------------------------- -blitz::Array<bool,1> RemoteStationCalibration::setPassBand(double start_freq, double stop_freq , const SpectralWindow& spw) +blitz::Array<bool,1> RemoteStationCalibration::setPassBand(double start_freq, double stop_freq , const SpectralWindow& spw) { double freq; blitz::Array<bool,1> inpassband(spw.getNumSubbands()); - + inpassband = true; for (int sb = 0; sb < spw.getNumSubbands(); sb++) { freq = sb * spw.getSubbandWidth() + spw.getSamplingFrequency() * (spw.getNyquistZone() - 1); @@ -290,7 +289,7 @@ blitz::Array<bool,1> RemoteStationCalibration::setPassBand(double start_freq, do } } //bmToFile(inpassband,"inPassBand",fileFormat); - return(inpassband); + return(inpassband); } //--------------------------------------------------------------------------------------------------------------------- @@ -302,43 +301,43 @@ blitz::Array<bool,1> RemoteStationCalibration::detectBadElements(ACC& acc, int p logfile << "nsel=" << nsel << endl; blitz::Array<double,2> ac(nelem,nsel); Timestamp acmtime; - + int sel = 0; for (int ch = 0; ch < nch; ch++) { if (selected(ch) == true) { ac(Range::all(),sel) = abs( bmDiag(acc.getACM(ch, pol, pol, acmtime).copy())); - sel++; + sel++; } } - + blitz::Array<double,2> med2(nelem,nsel); blitz::Array<double,1> med1(bmMedian(ac)); - med2 = med1(tensor::j); + med2 = med1(tensor::j); blitz::Array<double,2> diff( abs(ac - med2) ); - + blitz::Array<double,2> st2(nelem,nsel); blitz::Array<double,1> st1(bmStd(ac)); - st2 = st1(tensor::j); + st2 = st1(tensor::j); blitz::Array<double,2> threshold( 5 * st2); - + blitz::Array<int,2> inband( where(diff < threshold, 1, 0) ); blitz::Array<int,1> sum_ch( sum(inband(tensor::i,tensor::j),tensor::j) ); blitz::Array<bool,1> good( where(sum_ch == nsel, true, false)); - - //bmToFile(good,"goodElements",fileFormat); + + //bmToFile(good,"goodElements",fileFormat); return (good); } //--------------------------------------------------------------------------------------------------------------------- vector<Source> RemoteStationCalibration::make_local_sky_model( - const Sources& sources, - const Array<double, 1>& geoloc, + const Sources& sources, + const Array<double, 1>& geoloc, Timestamp& acmtime ) { double obstime = acmtime.sec() + acmtime.usec() / 1e6; - logfile << "obstime=" << obstime << endl; - + logfile << "obstime=" << obstime << endl; + std::vector<Source> skymodel = sources.getSources(); #if 0 AMC::Converter& converter = getConverter(); // can be used for coordinate conversions @@ -352,23 +351,23 @@ vector<Source> RemoteStationCalibration::make_local_sky_model( // number of seconds in a day: 86400 double JulianDay = obstime / 86400 + 2440587.5; logfile << "JulianDay=" << JulianDay << endl; - + LOG_INFO_STR("calibrate: time of observation in Julian Days: " << JulianDay - 2453500); // geographical location of station needed - //double geolon = geoloc(0); - //double geolat = geoloc(1); - + //double geolon = geoloc(0); + //double geolat = geoloc(1); + // used fixed values to compare with MatLab code double geolon = 6 + 52.0/60 + 3.17/3600; double geolat = 52 + 54.0/60 + 42.6/3600; - + logfile << "geolon=" << geolon << logfile << " , geolat=" << geolat << endl; // conversion from JD to Greenwich Star Time (GST) in seconds double TU = (JulianDay - 2451545) / 36525, - GST = (JulianDay + 0.5)*86400 + 24110.54841 + 8640184.812866 * TU + + GST = (JulianDay + 0.5)*86400 + 24110.54841 + 8640184.812866 * TU + 0.093104 * TU * TU + -6.2e-6 * TU * TU * TU; //logfile << "GST=" << GST << endl; @@ -396,9 +395,9 @@ vector<Source> RemoteStationCalibration::make_local_sky_model( if (el > 0) local_skymodel.push_back(Source((skymodel.begin() + idx)->getName(), -::cos(el) * ::sin(az), ::cos(el) * ::cos(az), const_cast<blitz::Array<double, 2>&>((skymodel.begin() + idx)->getFluxes()))); - + /* - logfile << "local_skymodel-" << idx << "=" << (skymodel.begin() + idx)->getName() << ", " + logfile << "local_skymodel-" << idx << "=" << (skymodel.begin() + idx)->getName() << ", " << (-::cos(el) * ::sin(az)) << ", " << (::cos(el) * ::cos(az)) << ", " << const_cast<blitz::Array<double, 2>&>((skymodel.begin() + idx)->getFluxes()) << endl; @@ -415,19 +414,19 @@ Array<complex<double>, 2> RemoteStationCalibration::computeA( { int nelem = AntennaPos.extent(firstDim); int nsrc = LSM.size(); - - + + Array<complex<double>, 2> asrc(nelem, nsrc); asrc = 0.0; - + Array<double, 1> xpos(AntennaPos(Range::all(), 0, 0)); Array<double, 1> ypos(AntennaPos(Range::all(), 0, 1)); //logfile << "freq=" << freq << endl; //logfile << "xpos=" << xpos << endl; //logfile << "ypos=" << ypos << endl; - + double k = 2 * M_PI * freq / itsSpeedOfLight; - + for (int idx = 0; idx < nsrc; idx++) { double l = (LSM.begin() + idx)->getRA(); double m = (LSM.begin() + idx)->getDEC(); @@ -442,7 +441,7 @@ Array<complex<double>, 2> RemoteStationCalibration::computeA( //--------------------------------------------------------------------------------------------------------------------- blitz::Array<bool, 2> RemoteStationCalibration::setRestriction( - const Array<double, 3>& AntennaPos, + const Array<double, 3>& AntennaPos, double freq ) { // three dimensions, Nantennas x Npol x 3 (x,y,z) @@ -453,12 +452,12 @@ blitz::Array<bool, 2> RemoteStationCalibration::setRestriction( Array<double, 2> u( xpos(tensor::j) - xpos(tensor::i) ); Array<double, 2> v( ypos(tensor::j) - ypos(tensor::i) ); - Array<double, 2> uvdist( sqrt(u * u + v * v) ); - - double minBaseLine = min( (itsBaseLineRestriction * (itsSpeedOfLight / freq)), itsMaxBaseLine); - + Array<double, 2> uvdist( sqrt(u * u + v * v) ); + + double minBaseLine = min( (itsBaseLineRestriction * (itsSpeedOfLight / freq)), itsMaxBaseLine); + blitz::Array<bool, 2> restriction( where(uvdist < minBaseLine, true, false)); - + // diagonal is always 1 for (int idx1 = 0; idx1 < uvdist.rows(); idx1++) { restriction(idx1, idx1) = 1; @@ -477,74 +476,74 @@ void RemoteStationCalibration::setSelection( { int dim1, dim2; blitz::Array<double, 2> isDiag(nelem*nelem, nelem); - blitz::Array<double, 2> eye( bmEye<double>(nelem)); + blitz::Array<double, 2> eye( bmEye<double>(nelem)); isDiag = bmKhatrirao( eye, eye); //bmToFile(isDiag,"isDiag",fileFormat); dim1 = isDiag.extent(firstDim); dim2 = isDiag.extent(secondDim); - + blitz::Array<int, 2> sel1(nelem, nelem); blitz::Array<int, 1> sel2(nelem*nelem); blitz::Array<int, 2> sel3(nelem*nelem, nelem*nelem); int selectedPoints; int idx; - + sel1 = where( bmTriu(restriction, 1) == true, 1, 0); sel2 = bmMatrixToVector(sel1); sel3 = bmDiag(sel2); selectedPoints = sum(sel2); //logfile << "Selected points Re=" << selectedPoints << endl; - + dim2 += selectedPoints; blitz::Array<double, 2> seltriu(nelem*nelem, selectedPoints); idx = 0; for (int idx1 = 0; idx1 < selectedPoints; idx1++) { while (sel2(idx) == 0) { idx++; } seltriu(Range::all(), idx1) = cast<double>(sel3(Range::all(),idx)); - idx++; + idx++; } //bmToFile(seltriu,"seltriu",fileFormat); - + sel1 = where(bmTril(restriction, -1) == true, 1, 0); sel2 = bmMatrixToVector(sel1); sel3 = bmDiag(sel2); - + selectedPoints = sum(sel2); dim2 += selectedPoints; //logfile << "Selected points Im=" << selectedPoints << endl; - + blitz::Array<double, 2> seltril(nelem*nelem, selectedPoints); idx = 0; for (int idx1 = 0; idx1 < selectedPoints; idx1++) { while (sel2(idx) == 0) { idx++; } seltril(Range::all(), idx1) = cast<double>(sel3(Range::all(),idx)); - idx++; + idx++; } //bmToFile(seltril,"seltril",fileFormat); - + blitz::Array<double, 2> isRe(seltriu + seltril); blitz::Array<double, 2> isIm(seltriu - seltril); //bmToFile(isRe,"isRe",fileFormat); //bmToFile(isIm,"isIm",fileFormat); - + //logfile << "set Isel and pinvIsel" << endl; int t1,t2; - Isel.resize(dim1, dim2); + Isel.resize(dim1, dim2); t1 = 0; t2 = isDiag.cols() - 1; Isel(Range::all(), Range(t1, t2)) = isDiag(Range::all(), Range::all()).copy(); - t1 = isDiag.cols(); t2 = isDiag.cols() + isRe.cols() - 1; + t1 = isDiag.cols(); t2 = isDiag.cols() + isRe.cols() - 1; Isel(Range::all(), Range(t1, t2)) = isRe(Range::all(), Range::all()).copy(); - t1 = isDiag.cols() + isRe.cols(); t2 = isDiag.cols() + isRe.cols() + isIm.cols() - 1; - Isel(Range::all(), Range(t1, t2)) = isIm(Range::all(), Range::all()).copy(); - + t1 = isDiag.cols() + isRe.cols(); t2 = isDiag.cols() + isRe.cols() + isIm.cols() - 1; + Isel(Range::all(), Range(t1, t2)) = isIm(Range::all(), Range::all()).copy(); + blitz::Array<double, 2> pinv_isel(dim1, dim2); isRe *= 0.5; isIm *= 0.5; - t1 = 0; t2 = isDiag.cols() - 1; + t1 = 0; t2 = isDiag.cols() - 1; pinv_isel(Range::all(), Range(t1, t2)) = isDiag(Range::all(), Range::all()).copy(); - t1 = isDiag.cols(); t2 = isDiag.cols() + isRe.cols() - 1; + t1 = isDiag.cols(); t2 = isDiag.cols() + isRe.cols() - 1; pinv_isel(Range::all(), Range(t1, t2)) = isRe(Range::all(), Range::all()).copy(); - t1 = isDiag.cols() + isRe.cols(); t2 = isDiag.cols() + isRe.cols() + isIm.cols() - 1; + t1 = isDiag.cols() + isRe.cols(); t2 = isDiag.cols() + isRe.cols() + isIm.cols() - 1; pinv_isel(Range::all(), Range(t1, t2)) = isIm(Range::all(), Range::all()).copy(); pinvIsel.resize(dim2,dim1); pinvIsel = bmTrans(pinv_isel); @@ -556,50 +555,50 @@ void RemoteStationCalibration::setSelection_2( blitz::Array<int,2>& posIsel, blitz::Array<double,2>& valIsel, blitz::Array<int,2>& posInvIsel, - blitz::Array<double,2>& valInvIsel, + blitz::Array<double,2>& valInvIsel, blitz::Array<bool, 2>& restriction, int nelem ) { int dim1, dim2; blitz::Array<int, 2> isDiag(nelem*nelem, nelem); - blitz::Array<int, 2> eye( bmEye<int>(nelem)); + blitz::Array<int, 2> eye( bmEye<int>(nelem)); isDiag = bmKhatrirao( eye, eye); dim1 = isDiag.extent(firstDim); dim2 = isDiag.extent(secondDim); - + blitz::Array<int, 2> sel1(where(bmTriu(restriction, 1) == true, 1, 0)); blitz::Array<int, 1> sel2(bmMatrixToVector(sel1)); blitz::Array<int, 2> sel3(bmDiag(sel2)); int selectedPoints; int idx; - + selectedPoints = sum(sel2); dim2 += selectedPoints; logfile << "Selected points Re=" << selectedPoints << endl; - - + + blitz::Array<int, 2> seltriu(nelem*nelem, selectedPoints); idx = 0; for (int idx1 = 0; idx1 < selectedPoints; idx1++) { while (sel2(idx) == 0) { idx++; } seltriu(Range::all(), idx1) = sel3(Range::all(),idx); } - + sel1 = where(bmTril(restriction, -1) == true, 1, 0); sel2 = bmMatrixToVector(sel1); sel3 = bmDiag(sel2); - + selectedPoints = sum(sel2); dim2 += selectedPoints; logfile << "Selected points Im=" << selectedPoints << endl; - + blitz::Array<int, 2> seltril(nelem*nelem, selectedPoints); idx = 0; for (int idx1 = 0; idx1 < selectedPoints; idx1++) { while (sel2(idx) == 0) { idx++; } seltril(Range::all(), idx1) = sel3(Range::all(),idx); } - + blitz::Array<int, 2> isRe(seltriu + seltril); blitz::Array<int, 2> isIm(seltriu - seltril); @@ -675,7 +674,7 @@ void RemoteStationCalibration::setSelection_2( } } logfile << "InvIsel done" << endl; - + } //--------------------------------------------------------------------------------------------------------------------- @@ -697,23 +696,23 @@ blitz::Array<std::complex<double>,1> RemoteStationCalibration::multIsel( blitz::Array<double,1> RemoteStationCalibration::computeFlux( Array<std::complex<double>, 2>& acm, - Array<std::complex<double>, 2>& asrc, + Array<std::complex<double>, 2>& asrc, Array<bool, 2>& restriction ) { Array<std::complex<double>, 2> tmp0( cast<std::complex<double> >( bmSqr( bmAbs( bmMult( bmTransH(asrc), asrc))))); Array<std::complex<double>, 2> tmp1( bmInv(tmp0)); - + Array<std::complex<double>, 2> tmp2( bmTransH( bmKhatrirao( bmConj(asrc), asrc))); Array<std::complex<double>, 2> tmp3(acm.rows()*acm.cols(),1); tmp3(Range::all(),0) = where( bmMatrixToVector(restriction) == false, bmMatrixToVector(acm), 0.0); Array<double,2> tmp4( bmReal( bmMult(tmp1, tmp2, tmp3))); - + Array<double,1> flux( tmp4(Range::all(),0)); - + flux = flux / flux(0); //flux = where(restriction, 0, flux); flux = where(flux < 0.0, 0.0 , flux); - + //bmToFile(flux,"flux",fileFormat); return (flux); } @@ -723,38 +722,38 @@ blitz::Array<double,1> RemoteStationCalibration::computeFlux( Array<std::complex<double>, 2> RemoteStationCalibration::computeAlpha( Array<std::complex<double>, 2>& sigma, Array<std::complex<double>, 2>& acm, - Array<std::complex<double>, 2>& asrc, + Array<std::complex<double>, 2>& asrc, Array<bool, 2>& restriction ) { int nelem = acm.extent(firstDim); Array<std::complex<double>, 2> alpha(nelem, nelem); alpha = 0.0; - + Array<std::complex<double>, 2> Rhat(acm.copy()); - + // .* (1 - mask) - Rhat = where(restriction == false, Rhat, 0.0); - - // R0 --> asrc * flux * asrc' + Rhat = where(restriction == false, Rhat, 0.0); + + // R0 --> asrc * flux * asrc' Array<std::complex<double>, 2> R0( bmMult(asrc, sigma, bmTransH(asrc))); // .* (1 - mask) R0 = where(restriction == false, R0, 0.0); - + Array<std::complex<double>, 1> w(nelem); - + for (int idx1 = 0; idx1 < nelem; idx1++) { for (int idx2 = 0; idx2 < nelem; idx2++) { w = 0.0; for (int idx3 = 0; idx3 < nelem; idx3++) { - + if ((R0(idx1,idx3) != 0.0) && (R0(idx2,idx3) != 0.0)) { if ((idx1 != idx3) || (idx2 != idx3)) { w(idx3) = abs( pow2( R0(idx1,idx3) * Rhat(idx2,idx3))); - alpha(idx1,idx2) += w(idx3) * conj((Rhat(idx1,idx3) * R0(idx2,idx3)) / (Rhat(idx2,idx3) * R0(idx1,idx3))); + alpha(idx1,idx2) += w(idx3) * conj((Rhat(idx1,idx3) * R0(idx2,idx3)) / (Rhat(idx2,idx3) * R0(idx1,idx3))); } } } - + if (sum(w) != 0.0) { alpha(idx1, idx2) /= sum(w); } @@ -767,9 +766,9 @@ Array<std::complex<double>, 2> RemoteStationCalibration::computeAlpha( Array<std::complex<double>, 1> RemoteStationCalibration::computeGain( - Array<std::complex<double>, 2>& acm, + Array<std::complex<double>, 2>& acm, Array<std::complex<double>, 2>& asrc, - Array<double, 1>& flux, + Array<double, 1>& flux, Array<bool, 2>& restriction, int maxiter, double diffstop ) @@ -777,56 +776,56 @@ Array<std::complex<double>, 1> RemoteStationCalibration::computeGain( //bmToFile(acm,"acm",fileFormat); //bmToFile(asrc,"asrc",fileFormat); //bmToFile(flux,"flux",fileFormat); - + logfile << "start of computeGain" << " time=" << getClock() << endl; - + int nelem = asrc.extent(firstDim); int nsrc = asrc.extent(secondDim); int nrestriction = sum(restriction); logfile << "asrc=" << asrc.shape() << endl; - - + + // for experimental Isel blitz::Array<int,2> posIsel; blitz::Array<double,2> valIsel; blitz::Array<int,2> posInvIsel; blitz::Array<double,2> valInvIsel; - setSelection_2(posIsel, valIsel, posInvIsel ,valInvIsel , restriction, nelem); + setSelection_2(posIsel, valIsel, posInvIsel ,valInvIsel , restriction, nelem); bmToFile(posIsel,"posIsel",fileFormat); bmToFile(valIsel,"valIsel",fileFormat); bmToFile(posInvIsel,"posInvIsel",fileFormat); bmToFile(valInvIsel,"valInvIsel",fileFormat); - - - blitz::Array<double, 2> Isel; + + + blitz::Array<double, 2> Isel; blitz::Array<double, 2> pinvIsel; setSelection(Isel, pinvIsel, restriction, nelem); bmToFile(Isel,"Isel",fileFormat); - bmToFile(pinvIsel,"pinvIsel",fileFormat); - + bmToFile(pinvIsel,"pinvIsel",fileFormat); + logfile << "selection done" << endl; - + Array<std::complex<double>, 2> gHat(nelem, maxiter+1); Array<std::complex<double>, 2> sigmaHat(nsrc, maxiter+1); Array<std::complex<double>, 2> sigmanHat(nrestriction, maxiter+1); Array<std::complex<double>, 2> sigman(nelem, nelem); gHat = 0.0; sigmaHat(Range::all(), 0) = cast<std::complex<double> >(flux); - sigmanHat(Range(0, nelem-1),0) = 1.0; + sigmanHat(Range(0, nelem-1),0) = 1.0; sigman = 0.0; - + // implementation using WALS int iter; for (iter = 1; iter < (maxiter+1); iter++) { logfile << "iteration=" << iter << " time=" << getClock() << endl; - + // estimate alpha = g * (1 ./ g) Array<std::complex<double>, 2> sigma( bmDiag( sigmaHat(Range::all(),iter-1))); Array<complex<double>, 2> alpha( computeAlpha( sigma, acm, asrc, restriction)); logfile << "alpha done" << " time=" << getClock() << endl; - - - // Note: the result reduced by this EVD algorithm + + + // Note: the result reduced by this EVD algorithm // differ from MatLab's eig-function for the lower eigenvalues Array<std::complex<double>,2> eigenvec; // V Array<std::complex<double>,2> eigenval; // D @@ -835,17 +834,17 @@ Array<std::complex<double>, 1> RemoteStationCalibration::computeGain( //bmToFile(eigenvec,"eigenvec",fileFormat); //bmToFile(eigenval,"eigenval",fileFormat); - + TinyVector<int, 1> maxidx( maxIndex( abs( bmDiag(eigenval)))); logfile << "maxidx=" << maxidx << endl; - + gHat(Range::all(), iter) = bmConj( eigenvec(Range::all(), maxidx(0))); - //bmToFile(gHat,"gHat1",fileFormat); + //bmToFile(gHat,"gHat1",fileFormat); Array<std::complex<double>,2> GA( bmMult( bmDiag( gHat(Range::all(),iter)), asrc)); Array<std::complex<double>,2> Rest(bmMult( GA, bmDiag( sigmaHat(Range::all(), iter-1)), bmTransH(GA))); //bmToFile(GA,"GA",fileFormat); //bmToFile(Rest,"Rest",fileFormat); - + // convert Matrix Rest and acm into a Vector and leave out diagonal Array<std::complex<double>,2> RestV((Rest.rows()-1)*Rest.cols(),1); Array<std::complex<double>,2> RhatV((acm.rows()-1)*acm.cols(),1); @@ -856,89 +855,89 @@ Array<std::complex<double>, 1> RemoteStationCalibration::computeGain( RestV(idx,0) = Rest(idx1, idx2); RhatV(idx,0) = acm(idx1, idx2); idx++; - } + } } } logfile << "MatrixToVector done" << " time=" << getClock() << endl; - - + + Array<double,2> normg( abs( sqrt( bmMult( bmPinv(RestV), RhatV)))); //bmToFile(normg,"normg",fileFormat); - + logfile << "normg done" << " time=" << getClock() << endl; gHat(Range::all(),iter) = normg(0,0) * gHat(Range::all(), iter) / (gHat(0, iter) / abs( gHat(0, iter))); //bmToFile(gHat,"gHat2",fileFormat); - + Array<std::complex<double>,2> R0( pow2(normg(0,0)) * Rest); Array<std::complex<double>,2> sigmanHatM(sigmanHat.rows(),1); - sigmanHatM(Range::all(),0) = sigmanHat(Range::all(), iter-1); - + sigmanHatM(Range::all(),0) = sigmanHat(Range::all(), iter-1); + Array<std::complex<double>,2> tmp1( bmMult(Isel, sigmanHatM)); sigman = bmVectorToMatrix(tmp1(Range::all(),0)); // speed experiment //Array<std::complex<double>,1> tmp01( sigmanHat(Range::all(), iter-1)); - //Array<std::complex<double>,1> tmp1( multIsel(posIsel, valIsel, tmp01)); + //Array<std::complex<double>,1> tmp1( multIsel(posIsel, valIsel, tmp01)); //sigman = bmVectorToMatrix(tmp1); - - + + //bmToFile(sigman,"sigman",fileFormat); logfile << "sigman tmp1 done" << " time=" << getClock() << endl; - + Array<std::complex<double>,2> tmp2(acm.rows()*acm.cols(),1); - tmp2(Range::all(),0) = (bmMatrixToVector(acm) - bmMatrixToVector(R0)); + tmp2(Range::all(),0) = (bmMatrixToVector(acm) - bmMatrixToVector(R0)); logfile << "tmp2 done" << " time=" << getClock() << endl; - - Array<std::complex<double>,2> tmp3( bmMult( pinvIsel, tmp2)); + + Array<std::complex<double>,2> tmp3( bmMult( pinvIsel, tmp2)); sigmanHat(Range::all(), iter) = tmp3(Range::all(),0); // speed experiment //Array<std::complex<double>,1> tmp03( bmMatrixToVector(acm) - bmMatrixToVector(R0)); //Array<std::complex<double>,1> tmp3( multIsel(posInvIsel, valInvIsel, tmp03)); //sigmanHat(Range::all(), iter) = tmp3(Range::all()); - + //bmToFile(sigmanHat,"sigmanHat",fileFormat); logfile << "tmp3 done" << " time=" << getClock() << endl; - + Array<std::complex<double>,2> tmp4(sigmanHat.rows(),1); - tmp4(Range::all(),0) = sigmanHat(Range::all(), iter); + tmp4(Range::all(),0) = sigmanHat(Range::all(), iter); logfile << "tmp4 done" << " time=" << getClock() << endl; - + Array<std::complex<double>,2> tmp5( bmMult( Isel, tmp4)); sigman = bmVectorToMatrix(tmp5(Range::all(),0)); // speed experiment //Array<std::complex<double>,1> tmp05( sigmanHat(Range::all(), iter)); //Array<std::complex<double>,1> tmp5( multIsel(posIsel, valIsel, tmp05)); //sigman = bmVectorToMatrix(tmp5(Range::all())); - - + + //bmToFile(sigman,"sigman",fileFormat); - + logfile << "sigman tmp5 done" << " time=" << getClock() << endl; Array<std::complex<double>,2> R(R0 + sigman); //bmToFile(R,"R",fileFormat); - Array<std::complex<double>,2> r1( cast<std::complex<double> >( bmAbs( bmSqr( bmConj( bmMult(bmTransH(GA), bmInv(R), GA)))))); + Array<std::complex<double>,2> r1( cast<std::complex<double> >( bmAbs( bmSqr( bmConj( bmMult(bmTransH(GA), bmInv(R), GA)))))); Array<std::complex<double>,2> r1a( bmInv( r1)); - logfile << "r1a done" << " time=" << getClock() << endl; - Array<std::complex<double>,2> r2(acm - sigman); + logfile << "r1a done" << " time=" << getClock() << endl; + Array<std::complex<double>,2> r2(acm - sigman); Array<std::complex<double>,2> r2a( bmMult(bmTransH(GA), bmInv(R), r2)); - logfile << "r2a done" << " time=" << getClock() << endl; + logfile << "r2a done" << " time=" << getClock() << endl; Array<std::complex<double>,1> r2b( bmDiag( bmMult(r2a, bmInv(R), GA))); Array<std::complex<double>,2> r2c(r2b.rows(),1); - r2c(Range::all(),0) = r2b(Range::all()); - logfile << "r2b done" << " time=" << getClock() << endl; - Array<double,2> sh( bmReal( bmMult(r1a, r2c))); + r2c(Range::all(),0) = r2b(Range::all()); + logfile << "r2b done" << " time=" << getClock() << endl; + Array<double,2> sh( bmReal( bmMult(r1a, r2c))); //bmToFile(sh,"sigmaHat1",fileFormat); - + Array<double,1> sh1( sh(Range::all(),0) / sh(0,0)); sh1 = where(sh1 > 0.0, sh1, 0.0); - + sigmaHat(Range::all(), iter) = cast<std::complex<double> >(sh1(Range::all())); - + //bmToFile(sigmaHat,"sigmaHat2",fileFormat); logfile << "sigmaHat done" << " time=" << getClock() << endl; Array<std::complex<double>,2> thetaPrev(gHat.rows()+sigmaHat.rows()+sigmanHat.rows(),1); - Array<std::complex<double>,2> theta(gHat.rows()+sigmaHat.rows()+sigmanHat.rows(),1); + Array<std::complex<double>,2> theta(gHat.rows()+sigmaHat.rows()+sigmanHat.rows(),1); int start, stop; start = 0; stop = gHat.rows() - 1; @@ -949,12 +948,12 @@ Array<std::complex<double>, 1> RemoteStationCalibration::computeGain( thetaPrev(Range(start, stop),0) = sigmaHat(Range::all(), iter-1); theta(Range(start, stop),0) = sigmaHat(Range::all(), iter); start = stop + 1; - stop = start + sigmanHat.rows() - 1; + stop = start + sigmanHat.rows() - 1; thetaPrev(Range(start, stop),0) = sigmanHat(Range::all(), iter-1); theta(Range(start, stop),0) = sigmanHat(Range::all(), iter); //bmToFile(thetaPrev,"thetaPrev",fileFormat); //bmToFile(theta,"theta",fileFormat); - + Array<double,2> difftest( abs( bmMult( bmPinv(thetaPrev), theta) - 1.0)); //bmToFile(difftest,"difftest",fileFormat); //logfile << "difftest done" << " time=" << getClock() << endl; @@ -962,30 +961,30 @@ Array<std::complex<double>, 1> RemoteStationCalibration::computeGain( if (difftest(0,0) < diffstop) { logfile << "(diff < diffstop) --> break" << endl; break; } } // end wals if (iter == (maxiter+1)) { iter--; } - + //bmToFile(gHat,"gHat",fileFormat); //bmToFile(sigmaHat,"sigmaHat",fileFormat); - + // only for printing sigman //Array<std::complex<double>,2> tmp4(sigmanHat.rows(),1); - //tmp4(Range::all(),0) = sigmanHat(Range::all(), iter); + //tmp4(Range::all(),0) = sigmanHat(Range::all(), iter); //Array<std::complex<double>,2> tmp5( bmMult( Isel, tmp4)); - //sigman = bmVectorToMatrix(tmp5(Range::all(),0)); + //sigman = bmVectorToMatrix(tmp5(Range::all(),0)); //bmToFile(sigman,"sigman",fileFormat); - + return (gHat(Range::all(), iter)); } //--------------------------------------------------------------------------------------------------------------------- /* Array<std::complex<double>, 1> RemoteStationCalibration::checkCal( - Array<std::complex<double>, 1>& cal, + Array<std::complex<double>, 1>& cal, blitz::Array<bool, 1>& good, int nch, int nelem ) { Array<double, 1> cal0( cal / bmRepmat( bmMedian( bmAbs(cal)))); cal0 = where(cal0 == NAN, 0.0, cal0); - + } */ //--------------------------------------------------------------------------------------------------------------------- diff --git a/MAC/APL/PAC/Cal_Server/src/SubArraySubscription.cc b/MAC/APL/PAC/Cal_Server/src/SubArraySubscription.cc index a55f3cd306e..1dabbc95b0b 100644 --- a/MAC/APL/PAC/Cal_Server/src/SubArraySubscription.cc +++ b/MAC/APL/PAC/Cal_Server/src/SubArraySubscription.cc @@ -36,19 +36,17 @@ void SubArraySubscription::update(Subject* subject) { ASSERT(subject == static_cast<Subject*>(m_subarray)); - AntennaGains* calibratedGains = 0; + //AntennaGains calibratedGains; // get gains from the FRONT buffer - if (m_subarray->getGains(calibratedGains, SubArray::FRONT)) { - - CALUpdateEvent update; - update.timestamp.setNow(0); - update.status = CAL_SUCCESS; - update.handle = (memptr_t)this; - - update.gains = *calibratedGains; - - if (m_port.isConnected()) m_port.send(update); - } + //calibratedGains = m_subarray->getGains(); + CALUpdateEvent update; + update.name = m_subarray->name(); + update.timestamp.setNow(0); + update.status = CAL_SUCCESS; + //update.handle = (memptr_t)this; + update.gains = m_subarray->getGains(); + + if (m_port.isConnected()) m_port.send(update); } diff --git a/MAC/APL/PAC/Cal_Server/src/SubArrays.cc b/MAC/APL/PAC/Cal_Server/src/SubArrays.cc index ab63cf7014e..62d71a40336 100644 --- a/MAC/APL/PAC/Cal_Server/src/SubArrays.cc +++ b/MAC/APL/PAC/Cal_Server/src/SubArrays.cc @@ -23,7 +23,7 @@ #include <lofar_config.h> #include <Common/LofarLogger.h> -#include <APL/CAL_Protocol/ACC.h> +#include "ACC.h" #include "SubArrays.h" //using namespace std; @@ -34,7 +34,7 @@ namespace LOFAR { // forward declaration class CalibrationInterface; - + // // SubArrays() // @@ -83,7 +83,7 @@ bool SubArrays::conflicting(SubArray* newArray) const if (activeRCUmode != newRCUmode) { // diffent rcumode? check rcus SubArray::RCUmask_t activeRCUs(it->second->getRCUMask()); if ((newRCUs & activeRCUs).any()) { - LOG_INFO_STR ("New subarray " << newArray->getName() << " would conflict with " << + LOG_INFO_STR ("New subarray " << newArray->getName() << " would conflict with " << it->second->getName() << ". Refused."); return (true); // we have a conflict } @@ -97,7 +97,7 @@ bool SubArrays::conflicting(SubArray* newArray) const SubArray::RCUmask_t activeRCUs(it->second->getRCUMask()); // Note: inDeadList is a relatively expensive call, call it as last-check. if ((newRCUs & activeRCUs).any() && !inDeadList(newArray->getName())) { - LOG_INFO_STR ("New subarray " << newArray->getName() << " would conflict with " << + LOG_INFO_STR ("New subarray " << newArray->getName() << " would conflict with " << it->second->getName() << ". Refused."); return (true); // we have a conflict } @@ -353,7 +353,7 @@ SubArrayMap SubArrays::getSubArrays(const string& optionalName) { else { // get subarray by name SubArrayMap::const_iterator subarray; - if (((subarray = m_arrays.find(optionalName)) != m_arrays.end()) || + if (((subarray = m_arrays.find(optionalName)) != m_arrays.end()) || ((subarray = m_new_arrays.find(optionalName)) != m_new_arrays.end())) { answer[optionalName] = subarray->second; } @@ -406,7 +406,7 @@ void SubArrays::writeGains(SubArray* anSubArr, const string& dataDir) return; } - if (fwrite(gains->getGains().data(), sizeof(complex<double>), gains->getGains().size(), gainFile) != + if (fwrite(gains->getGains().data(), sizeof(complex<double>), gains->getGains().size(), gainFile) != (size_t)gains->getGains().size()) { LOG_ERROR_STR("failed to write to file: " << filename); } diff --git a/MAC/APL/PAC/Cal_Server/src/SubArrays.h b/MAC/APL/PAC/Cal_Server/src/SubArrays.h index 8d39b88f693..0015cecc6e2 100644 --- a/MAC/APL/PAC/Cal_Server/src/SubArrays.h +++ b/MAC/APL/PAC/Cal_Server/src/SubArrays.h @@ -25,9 +25,10 @@ #define SUBARRAYS_H_ #include <APL/CAL_Protocol/SubArray.h> +#include "SharedResource.h" namespace LOFAR { - using EPA_Protocol::MEPHeader; + //using EPA_Protocol::MEPHeader; namespace CAL { // forward declarations @@ -98,8 +99,8 @@ private: bool remove(const string& name); bool remove(SubArray*& subarray); - SubArrayMap m_new_arrays; // subarrays that should be added - // for the next run of the + SubArrayMap m_new_arrays; // subarrays that should be added + // for the next run of the // calibration algorithm SubArrayMap m_arrays; list<SubArray*> m_dead_arrays; // subarrays that have been stopped @@ -108,5 +109,5 @@ private: }; // namespace CAL }; // namespace LOFAR -#endif +#endif diff --git a/MAC/APL/StationCU/src/CalibrationControl/CalibrationControl.cc b/MAC/APL/StationCU/src/CalibrationControl/CalibrationControl.cc index a482c6044dd..f370429824c 100644 --- a/MAC/APL/StationCU/src/CalibrationControl/CalibrationControl.cc +++ b/MAC/APL/StationCU/src/CalibrationControl/CalibrationControl.cc @@ -53,7 +53,7 @@ namespace LOFAR { using namespace Controller_Protocol; using namespace CAL_Protocol; namespace StationCU { - + // static pointer to this object for signal handler static CalibrationControl* thisCalibrationControl = 0; static uint16 gResult = CT_RESULT_NO_ERROR; @@ -145,13 +145,13 @@ void CalibrationControl::setState(CTState::CTstateNr newState) CTState cts; itsPropertySet->setValue(PN_FSM_CURRENT_ACTION, GCFPVString(cts.name(newState))); } -} +} // // convertFilterSelection(string) : uint8 // -int32 CalibrationControl::convertFilterSelection(const string& filterselection, const string& antennaSet) +int32 CalibrationControl::convertFilterSelection(const string& filterselection, const string& antennaSet) { // international stations don't have the LBL's connected, force them to use the LBH inputs. string tmpAntennaSet(antennaSet); // modifyable copy @@ -179,7 +179,7 @@ int32 CalibrationControl::convertFilterSelection(const string& filterselection, if (filterselection == "HBA_170_230") { return(6); } // 160 Mhz if (filterselection == "HBA_210_250") { return(7); } // 200 Mhz - LOG_WARN_STR ("filterselection value '" << filterselection << + LOG_WARN_STR ("filterselection value '" << filterselection << "' not recognized, using LBA_10_70"); return (1); } @@ -189,13 +189,13 @@ int32 CalibrationControl::convertFilterSelection(const string& filterselection, // // Connect to PVSS and report state back to StartDaemon // -GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, +GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, GCFPortInterface& port) { LOG_DEBUG_STR ("initial:" << eventName(event) << "@" << port.getName()); GCFEvent::TResult status = GCFEvent::HANDLED; - + switch (event.signal) { case F_INIT: break; @@ -227,7 +227,7 @@ GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, itsTimerPort->setTimer(0.0); } break; - + case F_TIMER: if (!itsPropertySetInitialized) { itsPropertySetInitialized = true; @@ -249,7 +249,7 @@ GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, itsPropertySet->setValue(PN_CC_FILTER, GCFPVString ("")); itsPropertySet->setValue(PN_CC_NYQUISTZONE, GCFPVInteger(0)); itsPropertySet->setValue(PN_CC_RCUS, GCFPVString ("")); - + // Start ParentControl task LOG_DEBUG ("Enabling ParentControl task and wait for my name"); itsParentPort = itsParentControl->registerTask(this); @@ -271,7 +271,7 @@ GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, LOG_DEBUG_STR("Received CONNECT(" << msg.cntlrName << ")"); setState(CTState::CONNECTED); sendControlResult(port, CONTROL_CONNECTED, msg.cntlrName, CT_RESULT_NO_ERROR); - + // let ParentControl watch over the start and stop times for extra safety. ptime startTime = time_from_string(globalParameterSet()-> getString("Observation.startTime")); @@ -287,7 +287,7 @@ GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, default: status = _defaultEventHandler(event, port); break; - } + } return (status); } @@ -297,7 +297,7 @@ GCFEvent::TResult CalibrationControl::initial_state(GCFEvent& event, // // wait for CLAIM event // -GCFEvent::TResult CalibrationControl::started_state(GCFEvent& event, +GCFEvent::TResult CalibrationControl::started_state(GCFEvent& event, GCFPortInterface& port) { LOG_DEBUG_STR ("started:" << eventName(event) << "@" << port.getName()); @@ -316,7 +316,7 @@ GCFEvent::TResult CalibrationControl::started_state(GCFEvent& event, break; case F_CONNECTED: // CONNECT must be from Calserver. - ASSERTSTR (&port == itsCalServer, + ASSERTSTR (&port == itsCalServer, "F_CONNECTED event from port " << port.getName()); itsTimerPort->cancelAllTimers(); itsPropertySet->setValue(PN_CC_CONNECTED,GCFPVBool(true)); @@ -328,7 +328,7 @@ GCFEvent::TResult CalibrationControl::started_state(GCFEvent& event, case F_DISCONNECTED: port.close(); - ASSERTSTR (&port == itsCalServer, + ASSERTSTR (&port == itsCalServer, "F_DISCONNECTED event from port " << port.getName()); itsPropertySet->setValue(PN_CC_CONNECTED,GCFPVBool(false)); LOG_WARN("Connection with CalServer failed, retry in 2 seconds"); @@ -361,14 +361,14 @@ GCFEvent::TResult CalibrationControl::started_state(GCFEvent& event, return (status); } - + // // claimed_state(event, port) // // wait for PREPARE event // -GCFEvent::TResult CalibrationControl::claimed_state(GCFEvent& event, +GCFEvent::TResult CalibrationControl::claimed_state(GCFEvent& event, GCFPortInterface& port) { LOG_DEBUG_STR ("claimed:" << eventName(event) << "@" << port.getName()); @@ -394,7 +394,7 @@ GCFEvent::TResult CalibrationControl::claimed_state(GCFEvent& event, case F_DISCONNECTED: port.close(); - ASSERTSTR (&port == itsCalServer, + ASSERTSTR (&port == itsCalServer, "F_DISCONNECTED event from port " << port.getName()); LOG_WARN("Connection with CalServer lost, going to reconnect state."); TRAN(CalibrationControl::started_state); @@ -458,12 +458,12 @@ GCFEvent::TResult CalibrationControl::claimed_state(GCFEvent& event, return (status); } - + // // active_state(event, port) // -// Normal operation state. +// Normal operation state. // GCFEvent::TResult CalibrationControl::active_state(GCFEvent& event, GCFPortInterface& port) { @@ -484,7 +484,7 @@ GCFEvent::TResult CalibrationControl::active_state(GCFEvent& event, GCFPortInter case F_DISCONNECTED: port.close(); - ASSERTSTR (&port == itsCalServer, + ASSERTSTR (&port == itsCalServer, "F_DISCONNECTED event from port " << port.getName()); LOG_DEBUG("Connection with CalServer lost, going to reconnect"); TRAN(CalibrationControl::started_state); @@ -535,7 +535,7 @@ GCFEvent::TResult CalibrationControl::active_state(GCFEvent& event, GCFPortInter LOG_INFO_STR ("Calibration of beam " << ack.name << " successfully stopped"); } else { - LOG_WARN_STR ("Calibration of beam " << ack.name << " stopped with ERROR:" + LOG_WARN_STR ("Calibration of beam " << ack.name << " stopped with ERROR:" << errorName(ack.status)); } @@ -562,7 +562,7 @@ GCFEvent::TResult CalibrationControl::active_state(GCFEvent& event, GCFPortInter else { LOG_ERROR("Stop of some calibrations failed, staying in SUSPENDED mode"); setObjectState("Cannot stop the calibration", itsPropertySet->getFullScope(), RTDB_OBJ_STATE_BROKEN); - sendControlResult(*itsParentPort, CONTROL_RELEASED, getName(), + sendControlResult(*itsParentPort, CONTROL_RELEASED, getName(), CT_RESULT_CALSTOP_FAILED); setState(CTState::SUSPENDED); } @@ -585,7 +585,7 @@ GCFEvent::TResult CalibrationControl::active_state(GCFEvent& event, GCFPortInter // // Quiting: send QUITED, wait 1 second and stop // -GCFEvent::TResult CalibrationControl::quiting_state(GCFEvent& event, +GCFEvent::TResult CalibrationControl::quiting_state(GCFEvent& event, GCFPortInterface& port) { LOG_DEBUG_STR ("quiting:" << eventName(event) << "@" << port.getName()); @@ -631,7 +631,7 @@ GCFEvent::TResult CalibrationControl::quiting_state(GCFEvent& event, // // startCalibration() // -bool CalibrationControl::startCalibration() +bool CalibrationControl::startCalibration() { itsNrBeams = itsObsPar->beams.size(); LOG_DEBUG_STR("Calibrating " << itsNrBeams << " beams."); @@ -650,17 +650,15 @@ bool CalibrationControl::startCalibration() StationConfig config; // TODO: As long as the AntennaArray.conf uses different names as SAS we have to use this dirty hack. // calStartEvent.parent = itsObsPar->getAntennaArrayName(config.hasSplitters); - calStartEvent.parent = AS->antennaField(itsObsPar->beams[i].antennaSet); - calStartEvent.rcumode().resize(1); - calStartEvent.rcumode()(0).setMode((RSP_Protocol::RCUSettings::Control::RCUMode) - convertFilterSelection(itsObsPar->filter, itsObsPar->beams[i].antennaSet)); - calStartEvent.subset = itsObsPar->getRCUbitset(0, 0, itsObsPar->beams[i].antennaSet) & + calStartEvent.antennaSet = AS->antennaField(itsObsPar->beams[i].antennaSet); + calStartEvent.rcumode = convertFilterSelection(itsObsPar->filter, itsObsPar->beams[i].antennaSet); + calStartEvent.rcuMask = itsObsPar->getRCUbitset(0, 0, itsObsPar->beams[i].antennaSet) & AS->RCUallocation(itsObsPar->beams[i].antennaSet); // Note: when HBA_DUAL is selected we should set up a calibration on both HBA_0 and HBA_1 field. - LOG_DEBUG(formatString("Sending CALSTART(%s,%s,%08X)", - calStartEvent.name.c_str(), calStartEvent.parent.c_str(), - calStartEvent.rcumode()(0).getRaw())); + LOG_DEBUG(formatString("Sending CALSTART(%s,%s,%d)", + calStartEvent.name.c_str(), calStartEvent.antennaSet.c_str(), + calStartEvent.rcumode)); itsCalServer->send(calStartEvent); beamNameArr.push_back(new GCFPVString(itsObsPar->beams[i].name)); // update array for PVSS } // for all beams -- GitLab