Revert "Removed station response packages, which are not used in LOFAR and...

Revert "Removed station response packages, which are not used in LOFAR and externally available as https://github.com/lofar-astron/LOFARBeam". Was still needed by BBSKernel, oops.

This reverts commit 34770298.

CEP/Calibration/CMakeLists.txt

@@ -2,3 +2,7 @@
 
 lofar_add_package(BBSKernel)
 lofar_add_package(BBSControl)
+lofar_add_package(ExpIon)
+lofar_add_package(pystationresponse)
+lofar_add_package(ElementResponse)
+lofar_add_package(StationResponse)

CEP/Calibration/ElementResponse/CMakeLists.txt

+# $Id$
+
+lofar_package(ElementResponse 0.1 DEPENDS Common)
+
+# Uncomment to check for unsafe conversions (gcc), for example conversion of
+# size_t to unsigned int (truncation).
+#add_definitions(-Wconversion)
+
+add_subdirectory(include/ElementResponse)
+add_subdirectory(src)

CEP/Calibration/ElementResponse/include/ElementResponse/CMakeLists.txt

+# $Id$
+
+# Create symbolic link to include directory.
+execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink
+  ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_BINARY_DIR}/include/${PACKAGE_NAME})
+
+# Install header files.
+install(FILES
+  ElementResponse.h
+  DESTINATION include/${PACKAGE_NAME})

CEP/Calibration/ElementResponse/include/ElementResponse/ElementResponse.h

+//# ElementResponse.h: Functions to compute the (idealized) response of a LOFAR
+//# LBA or HBA dual dipole antenna.
+//#
+//# Copyright (C) 2011
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+
+#ifndef LOFAR_ELEMENTRESPONSE_H
+#define LOFAR_ELEMENTRESPONSE_H
+
+// \file
+// Functions to compute the (idealized) response of a LOFAR LBA or HBA dual
+// dipole antenna.
+
+#include <complex>
+
+namespace LOFAR
+{
+
+// \addtogroup ElementResponse
+// @{
+
+// Compute the response of an idealized LOFAR LBA dual dipole antenna to
+// radiation at frequency freq (Hz) arriving from the direction given by theta,
+// phi (rad). The antenna model is described in a spherical coordinate system
+// with coordinates theta (zenith angle) and phi (azimith). The +X dipole is at
+// azimuth zero, the +Y dipole is at azimuth PI / 2.0.
+//
+// Preconditions:
+// --------------
+// freq: Frequency in Hz in the range [10 MHz, 100 MHz].
+// theta: Zenith angle in rad in the range [0.0, PI / 2.0].
+// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
+//
+void element_response_lba(double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2]);
+
+// Compute the response of an idealized LOFAR HBA dual dipole antenna to
+// radiation at frequency freq (Hz) arriving from the direction given by theta,
+// phi (rad). The antenna model is described in a spherical coordinate system
+// with coordinates theta (zenith angle) and phi (azimith). The +X dipole is at
+// azimuth zero, the +Y dipole is at azimuth PI / 2.0.
+//
+// Preconditions:
+// --------------
+// freq: Frequency in Hz in the range [120 MHz, 240 MHz].
+// theta: Zenith angle in rad in the range [0.0, PI / 2.0].
+// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
+//
+void element_response_hba(double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2]);
+
+// Compute the response of an idealized LOFAR dual dipole antenna to radiation
+// at frequency freq (Hz) arriving from the direction given by theta, phi (rad).
+// The antenna model is described in a spherical coordinate system with
+// coordinates theta (zenith angle) and phi (azimith). The +X dipole is at
+// azimuth zero, the +Y dipole is at azimuth PI / 2.0.
+//
+// This function uses a set of user defined coefficients to evaluate the beam
+// model. The coeff_shape parameter defines the shape of the coefficient array
+// as no. of harmonics x degree in theta x degree in frequency x 2. The last
+// dimension is implicit and always equal to 2. The coeff parameter points to an
+// array of coefficients of the proper size, stored in row-major order
+// ("C"-order). The freq_center and freq_range parameters define the frequency
+// range covered by the model described by the set of coefficients.
+//
+// Preconditions:
+// --------------
+// freq: Frequency in Hz in the range [freq_center - freq_range,
+//     freq_center + freq_range].
+// theta: Zenith angle in rad in the range [0.0, PI / 2.0].
+// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
+// freq_range, freq_center: Frequency center and range in Hz, should be > 0.
+// coeff_shape: Shape of the coefficient array, all dimensions should be > 0.
+//
+void element_response(double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2], double freq_center,
+    double freq_range, const unsigned int (&coeff_shape)[3],
+    const std::complex<double> coeff[]);
+
+// @}
+
+} //# namespace LOFAR
+
+#endif

CEP/Calibration/ElementResponse/src/CMakeLists.txt

+# $Id: CMakeLists.txt 18775 2011-09-06 13:36:45Z zwieten $
+
+include(LofarPackageVersion)
+
+lofar_add_library(elementresponse
+  Package__Version.cc
+  ElementResponse.cc)

CEP/Calibration/ElementResponse/src/DefaultCoeffHBA.cc

+// Beam model coefficients converted by convert_coeff.py.
+// Conversion performed on 2011/09/14/08:23:16 UTC using: 
+//     convert_coeff.py element_beam_HBA.coeff DefaultCoeffHBA.cc default_hba
+
+#include <complex>
+
+// Center frequency, and frequency range for which the beam model coefficients
+// are valid. The beam model is parameterized in terms of a normalized
+// frequency f' in the range [-1.0, 1.0]. The appropriate conversion is:
+//
+//     f' = (f - center) / range
+//
+const double default_hba_freq_center = 180e6;
+const double default_hba_freq_range = 60e6;
+
+// Shape of the coefficient array: 2x5x5x2 (the size of the last dimension is
+// implied, and always equal to 2).
+//
+const unsigned int default_hba_coeff_shape[3] = {2, 5, 5};
+
+// The array of coefficients in row-major order ("C"-order).
+//
+const std::complex<double> default_hba_coeff[100] = {
+    std::complex<double>(0.9989322499459223, 0.0003305895124867), std::complex<double>(1.0030546028872600, 0.0002157249025076),
+    std::complex<double>(0.0003002209532403, 0.0007909077657054), std::complex<double>(0.0022051270911392, 0.0003834815341981),
+    std::complex<double>(-0.0003856663268042, 0.0008435910525861), std::complex<double>(0.0004887765294093, 0.0002777796480946),
+    std::complex<double>(-0.0000699366665322, 0.0005136144371953), std::complex<double>(0.0001520602842105, 0.0001303481681886),
+    std::complex<double>(0.0000512381993616, 0.0001550785137302), std::complex<double>(0.0000819244737818, 0.0000466470412396),
+    std::complex<double>(0.0249658150445263, -0.0122024663463393), std::complex<double>(-0.0917825091832822, -0.0062606338208358),
+    std::complex<double>(-0.0083709499453879, -0.0289759752488368), std::complex<double>(-0.0689260153643395, -0.0111348626546314),
+    std::complex<double>(0.0116296166994115, -0.0307342946951178), std::complex<double>(-0.0171249717275797, -0.0080642275561593),
+    std::complex<double>(0.0012408055399100, -0.0191295543986957), std::complex<double>(-0.0051031652662961, -0.0037143632875100),
+    std::complex<double>(-0.0022414352263751, -0.0060474723525871), std::complex<double>(-0.0024377933436567, -0.0012852163337395),
+    std::complex<double>(-0.6730977722052307, 0.0940030437973656), std::complex<double>(0.3711597596859299, 0.0557089394867947),
+    std::complex<double>(0.2119250520015808, 0.2155514942677135), std::complex<double>(0.6727380529527980, 0.0989550572104158),
+    std::complex<double>(-0.0419944347289523, 0.2355624543349744), std::complex<double>(0.1917656461134636, 0.0732470381581913),
+    std::complex<double>(0.0048918921441903, 0.1588912409502319), std::complex<double>(0.0575369727210951, 0.0344677222786687),
+    std::complex<double>(0.0241014578366618, 0.0547046570516960), std::complex<double>(0.0219986510834463, 0.0112189146988984),
+    std::complex<double>(0.0665319393516388, -0.1418009730472832), std::complex<double>(-0.7576728614553603, -0.0472040122949963),
+    std::complex<double>(-0.1017024786435272, -0.3302620837788515), std::complex<double>(-0.5600906156274197, -0.0797555201430585),
+    std::complex<double>(0.0889729243872774, -0.3406964719938829), std::complex<double>(-0.1342560801672904, -0.0515926960946038),
+    std::complex<double>(-0.0149335262655201, -0.2084962323582034), std::complex<double>(-0.0327252678958813, -0.0172371907472848),
+    std::complex<double>(-0.0362395089905272, -0.0661322227928722), std::complex<double>(-0.0141568558526096, -0.0042676979206835),
+    std::complex<double>(0.1121669548152054, 0.0504713119323919), std::complex<double>(0.1882531376700409, 0.0088411256350159),
+    std::complex<double>(0.0066968933526899, 0.1181452711088882), std::complex<double>(0.0981630367567397, 0.0129921405004959),
+    std::complex<double>(-0.0347327225501659, 0.1186585563636635), std::complex<double>(0.0102831315790362, 0.0046275244914932),
+    std::complex<double>(0.0070209144233666, 0.0689639468490938), std::complex<double>(-0.0020239346031291, -0.0025499069613344),
+    std::complex<double>(0.0132702874173192, 0.0207916487187541), std::complex<double>(0.0004387107229914, -0.0017223838914815),
+    std::complex<double>(-0.0004916757488397, 0.0000266213616248), std::complex<double>(0.0006516553273188, -0.0000433166563288),
+    std::complex<double>(-0.0004357897643121, 0.0000320567996700), std::complex<double>(0.0005818285824826, -0.0001021069650381),
+    std::complex<double>(-0.0001047488648808, -0.0000302146563592), std::complex<double>(0.0001593350153828, -0.0000879125663990),
+    std::complex<double>(-0.0000141882506567, -0.0000941521783975), std::complex<double>(-0.0000004226298134, -0.0000245060763932),
+    std::complex<double>(-0.0000177429496833, -0.0000561890408003), std::complex<double>(-0.0000018388829279, 0.0000032387726477),
+    std::complex<double>(0.0162495046881796, -0.0010736997976255), std::complex<double>(-0.0175635905033026, 0.0012997068962173),
+    std::complex<double>(0.0138897851110661, -0.0014876219938565), std::complex<double>(-0.0150211436594772, 0.0029712291209158),
+    std::complex<double>(0.0031705620225488, 0.0004838463688512), std::complex<double>(-0.0034418973689263, 0.0024603729467258),
+    std::complex<double>(0.0003028387544878, 0.0026905629457281), std::complex<double>(0.0006768121359769, 0.0005901486396051),
+    std::complex<double>(0.0004634797107989, 0.0016976603895716), std::complex<double>(0.0003344773954073, -0.0001499932789294),
+    std::complex<double>(-0.1492097398080444, 0.0123735410547393), std::complex<double>(0.1393121453502456, -0.0121117146246749),
+    std::complex<double>(-0.1217628319418324, 0.0222643129255504), std::complex<double>(0.1108579917761457, -0.0262986164183475),
+    std::complex<double>(-0.0273147374272124, 0.0098595182007132), std::complex<double>(0.0208992817013466, -0.0205929453727953),
+    std::complex<double>(-0.0002152227668601, -0.0089220757225133), std::complex<double>(-0.0074792188817697, -0.0043562231368076),
+    std::complex<double>(-0.0012019994038721, -0.0079939660050373), std::complex<double>(-0.0035807498769946, 0.0014801422733613),
+    std::complex<double>(0.1567990061437258, -0.0143275575385193), std::complex<double>(-0.1043118778001582, 0.0106756004832779),
+    std::complex<double>(0.1151024257152241, -0.0225518489392044), std::complex<double>(-0.0593437249231851, 0.0216080058910987),
+    std::complex<double>(0.0142781186223020, -0.0057037138045721), std::complex<double>(0.0151043140114779, 0.0141435752121475),
+    std::complex<double>(-0.0057143555179676, 0.0141142700941743), std::complex<double>(0.0251435557201315, -0.0005753615445942),
+    std::complex<double>(0.0004475745352473, 0.0102135659618127), std::complex<double>(0.0090474375150397, -0.0032177128650026),
+    std::complex<double>(-0.0459124372023251, 0.0044990718645418), std::complex<double>(0.0135433541303599, -0.0021789296923529),
+    std::complex<double>(-0.0306136798186735, 0.0064963361606382), std::complex<double>(-0.0046440676338940, -0.0037281688158807),
+    std::complex<double>(-0.0006372791846825, 0.0008894047150233), std::complex<double>(-0.0181611528840412, -0.0011106177431486),
+    std::complex<double>(0.0032325387394458, -0.0048123509184894), std::complex<double>(-0.0136340313457176, 0.0021185000810664),
+    std::complex<double>(0.0001287985092565, -0.0032079544559908), std::complex<double>(-0.0045503800737417, 0.0015366231416036)
+};

CEP/Calibration/ElementResponse/src/DefaultCoeffLBA.cc

+// Beam model coefficients converted by convert_coeff.py.
+// Conversion performed on 2011/09/14/08:23:09 UTC using: 
+//     convert_coeff.py element_beam_LBA.coeff DefaultCoeffLBA.cc default_lba
+
+#include <complex>
+
+// Center frequency, and frequency range for which the beam model coefficients
+// are valid. The beam model is parameterized in terms of a normalized
+// frequency f' in the range [-1.0, 1.0]. The appropriate conversion is:
+//
+//     f' = (f - center) / range
+//
+const double default_lba_freq_center = 55e6;
+const double default_lba_freq_range = 45e6;
+
+// Shape of the coefficient array: 2x5x5x2 (the size of the last dimension is
+// implied, and always equal to 2).
+//
+const unsigned int default_lba_coeff_shape[3] = {2, 5, 5};
+
+// The array of coefficients in row-major order ("C"-order).
+//
+const std::complex<double> default_lba_coeff[100] = {
+    std::complex<double>(0.9982445079290715, 0.0000650863154389), std::complex<double>(1.0006230902158257, -0.0022053287681416),
+    std::complex<double>(0.0001002692200362, 0.0006838211278268), std::complex<double>(-0.0003660049052840, -0.0008418920419220),
+    std::complex<double>(-0.0010581424498791, 0.0015237878543047), std::complex<double>(0.0007398729642721, 0.0028468649470433),
+    std::complex<double>(-0.0039458389254656, -0.0007048354913730), std::complex<double>(0.0007040177887611, 0.0007856369612188),
+    std::complex<double>(-0.0031701591723043, -0.0010521154166512), std::complex<double>(-0.0007213036752903, -0.0007227764008022),
+    std::complex<double>(0.0550606068782634, 0.0011958385659938), std::complex<double>(-0.0160912944232080, 0.0703645376267940),
+    std::complex<double>(0.0033849565901213, -0.0244636379385135), std::complex<double>(0.0234264238829944, 0.0084068836453700),
+    std::complex<double>(0.0557107413978542, -0.0634701730653090), std::complex<double>(-0.0139549526991330, -0.1175401658864208),
+    std::complex<double>(0.1336911750356096, 0.0202651327657687), std::complex<double>(-0.0113385668361727, -0.0339262369086247),
+    std::complex<double>(0.0962263571740972, 0.0440074333288440), std::complex<double>(0.0313595045238824, 0.0230763038515351),
+    std::complex<double>(-0.8624889445327827, -0.1522883072804402), std::complex<double>(-0.0386800869486029, -0.7569350701887934),
+    std::complex<double>(0.0891332399420108, 0.1876527151756476), std::complex<double>(-0.1012363483900640, -0.1975118891151966),
+    std::complex<double>(-0.6404795825927633, 0.7568775384981410), std::complex<double>(0.0767245154665722, 1.3441875993523555),
+    std::complex<double>(-0.8758406699506004, 0.3350237639226141), std::complex<double>(0.2824832769101577, 0.6821307442669313),
+    std::complex<double>(-0.3144282315609649, -0.2763869580286276), std::complex<double>(-0.1705959031354030, -0.0712085950559831),
+    std::complex<double>(0.4039567648146965, 0.0810473144253429), std::complex<double>(-0.0350803390479135, 0.5214591717801087),
+    std::complex<double>(0.2232030356124932, -0.2248154851829713), std::complex<double>(0.4704343293662089, -0.3552101485419532),
+    std::complex<double>(0.9646419509627557, -0.8095088593139815), std::complex<double>(0.1635280638865702, -1.4854352979459096),
+    std::complex<double>(1.0331569921006993, 0.0509705885336283), std::complex<double>(0.1501121326521990, -0.5193414816770609),
+    std::complex<double>(0.4715775965513117, 0.5077361528286819), std::complex<double>(0.3847391427972284, 0.1136717951238837),
+    std::complex<double>(-0.0756250564248881, 0.0056622911723172), std::complex<double>(-0.1267444401630109, -0.0349676272376008),
+    std::complex<double>(-0.1793752883639813, 0.0720222655359702), std::complex<double>(-0.2678542619793421, 0.3152115802895427),
+    std::complex<double>(-0.3718069213271066, 0.2275266747872172), std::complex<double>(-0.1372223722572021, 0.4314989948093362),
+    std::complex<double>(-0.3316657641578328, -0.1655909947939444), std::complex<double>(-0.2158100484836540, 0.0614504774034524),
+    std::complex<double>(-0.1901597954359592, -0.2294955549701665), std::complex<double>(-0.1864961465389693, -0.0486276177310768),
+    std::complex<double>(-0.0000762326746410, 0.0000118155774181), std::complex<double>(0.0000118903581604, -0.0000251324432498),
+    std::complex<double>(-0.0002204197663391, -0.0000213776348027), std::complex<double>(0.0001477083861977, 0.0000599750510518),
+    std::complex<double>(-0.0003281057522772, -0.0000770207588466), std::complex<double>(0.0003478997686964, 0.0001481982639746),
+    std::complex<double>(-0.0000625695757282, 0.0000249138990722), std::complex<double>(-0.0000960097542525, 0.0002521364065803),
+    std::complex<double>(0.0001275344578325, 0.0000652362392482), std::complex<double>(-0.0003113309221942, 0.0001956734476566),
+    std::complex<double>(0.0029807707669629, -0.0003262084082071), std::complex<double>(0.0001639620574332, -0.0000266272685197),
+    std::complex<double>(0.0076282580587895, 0.0026614359017468), std::complex<double>(-0.0044850263974801, -0.0058337192660638),
+    std::complex<double>(0.0124258438959177, 0.0067985224235178), std::complex<double>(-0.0126349778957970, -0.0100656881493938),
+    std::complex<double>(0.0059031372522229, 0.0008660479915339), std::complex<double>(0.0039660364524413, -0.0100356333791398),
+    std::complex<double>(-0.0020520685193773, -0.0028564379463666), std::complex<double>(0.0121039958869239, -0.0059701468961263),
+    std::complex<double>(-0.0229975846564195, 0.0010565261888195), std::complex<double>(0.0019573207027441, 0.0050550600926414),
+    std::complex<double>(-0.0682274156850413, -0.0758159820140411), std::complex<double>(0.0497303968865466, 0.1019681987654797),
+    std::complex<double>(-0.1757936183439326, -0.1363710820472197), std::complex<double>(0.1765450269056824, 0.1555919358121995),
+    std::complex<double>(-0.1541299429420569, -0.0281422177614844), std::complex<double>(0.0816399676454817, 0.0691599035109852),
+    std::complex<double>(-0.0235110916473515, 0.0306385386726702), std::complex<double>(-0.0474273292450285, 0.0116831908947225),
+    std::complex<double>(0.0333560984394624, -0.0009767086536162), std::complex<double>(0.0141704479374002, -0.0205386534626779),
+    std::complex<double>(0.0562541280098909, 0.0743149092143081), std::complex<double>(-0.0226634801339250, -0.1439026188572270),
+    std::complex<double>(0.1238595124159999, 0.1766108700786397), std::complex<double>(-0.1307647072780430, -0.2090615438301942),
+    std::complex<double>(0.1557916917691289, 0.0646351862895731), std::complex<double>(0.0170294191358757, -0.1027926845803498),
+    std::complex<double>(0.0543537332385954, -0.0366524906364179), std::complex<double>(0.1127180664279469, -0.0176607923511174),
+    std::complex<double>(-0.0126732549319889, 0.0002042370658763), std::complex<double>(-0.0101360135082899, 0.0114084024114141),
+    std::complex<double>(-0.0102147881225462, -0.0176848554302252), std::complex<double>(-0.0051268936720694, 0.0527621533959941),
+    std::complex<double>(-0.0110701836450407, -0.0593085026046026), std::complex<double>(0.0140598301629874, 0.0738668439833535),
+    std::complex<double>(-0.0389912915621699, -0.0301364165752433), std::complex<double>(-0.0462331759359031, 0.0405864871628086),
+    std::complex<double>(-0.0251598701859194, 0.0115712688652445), std::complex<double>(-0.0563476280247398, 0.0079787883434624)
+};

CEP/Calibration/ElementResponse/src/ElementResponse.cc

+//# ElementResponse.cc: Functions to compute the (idealized) response of a LOFAR
+//# LBA or HBA dual dipole antenna.
+//#
+//# Copyright (C) 2011
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+
+#include <lofar_config.h>
+#include <ElementResponse/ElementResponse.h>
+#include <cmath>
+
+// The coefficients are kept in an unnamed namespace which effectively makes
+// them invisible outside this translation unit.
+namespace
+{
+// PI / 2.0
+const double pi_2 = 1.570796326794896619231322;
+
+#include "DefaultCoeffLBA.cc"
+#include "DefaultCoeffHBA.cc"
+}
+
+namespace LOFAR
+{
+
+void element_response_lba(double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2])
+{
+    element_response(freq, theta, phi, response, default_lba_freq_center,
+        default_lba_freq_range, default_lba_coeff_shape, default_lba_coeff);
+}
+
+void element_response_hba(double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2])
+{
+    element_response(freq, theta, phi, response, default_hba_freq_center,
+        default_hba_freq_range, default_hba_coeff_shape, default_hba_coeff);
+}
+
+void element_response(double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2], double freq_center,
+    double freq_range, const unsigned int (&coeff_shape)[3],
+    const std::complex<double> coeff[])
+{
+    // Initialize the response to zero.
+    response[0][0] = 0.0;
+    response[0][1] = 0.0;
+    response[1][0] = 0.0;
+    response[1][1] = 0.0;
+
+    // Clip directions below the horizon.
+    if(theta >= pi_2)
+    {
+        return;
+    }
+
+    const unsigned int nHarmonics  = coeff_shape[0];
+    const unsigned int nPowerTheta = coeff_shape[1];
+    const unsigned int nPowerFreq  = coeff_shape[2];
+
+    // The model is parameterized in terms of a normalized frequency in the
+    // range [-1, 1]. The appropriate conversion is taken care of below.
+    freq = (freq - freq_center) / freq_range;
+
+    // The variables sign and kappa are used to compute the value of kappa
+    // mentioned in the description of the beam model [kappa = (-1)^k * (2 * k
+    //+ 1)] incrementally.
+    int sign = 1, kappa = 1;
+
+    std::complex<double> P[2], Pj[2];
+    for(unsigned int k = 0; k < nHarmonics; ++k)
+    {
+        // Compute the (diagonal) projection matrix P for the current harmonic.
+        // This requires the evaluation of two polynomials in theta and freq (of
+        // degree nPowerTheta in theta and nPowerFreq in freq), one for each
+        // element of P. The polynomials are evaluated using Horner's rule.
+
+        // Horner's rule requires backward iteration of the coefficients, so
+        // move the iterator to the first location past the end of the block of
+        // coefficients for the current harmonic (k).
+        coeff += nPowerTheta * nPowerFreq * 2;
+
+        // Evaluate the polynomial. Note that the iterator is always decremented
+        // before it is dereferenced, using the prefix operator--. After
+        // evaluation of the polynomial, the iterator points exactly to the
+        // beginning of the block of coefficients for the current harmonic (k),
+        // that is, all the decrements together exactly cancel the increment
+        // aplied above.
+
+        // Evaluate the highest order term. Note that the order of the
+        // assigments is important because of the iterator decrement, i.e. P[1]
+        // should be assigned first.
+        P[1] = *--coeff;
+        P[0] = *--coeff;
+
+        for(unsigned int i = 0; i < nPowerFreq - 1; ++i)
+        {
+            P[1] = P[1] * freq + *--coeff;
+            P[0] = P[0] * freq + *--coeff;
+        }
+
+        // Evaluate the remaining terms.
+        for(unsigned int j = 0; j < nPowerTheta - 1; ++j)
+        {
+            Pj[1] = *--coeff;
+            Pj[0] = *--coeff;
+
+            for(unsigned int i = 0; i < nPowerFreq - 1; ++i)
+            {
+                Pj[1] = Pj[1] * freq + *--coeff;
+                Pj[0] = Pj[0] * freq + *--coeff;
+            }
+
+            P[1] = P[1] * theta + Pj[1];
+            P[0] = P[0] * theta + Pj[0];
+        }
+
+        // Because the increment and decrements cancel, the iterator points to
+        // the same location as at the beginning of this iteration of the outer
+        // loop. The next iteration should use the coefficients for the next
+        // harmonic (k), so we move the iterator to the start of that block.
+        coeff += nPowerTheta * nPowerFreq * 2;
+
+        // Compute the Jones matrix for the current harmonic, by rotating P over
+        // kappa * az, and add it to the result.
+        const double angle = sign * kappa * phi;
+        const double caz = std::cos(angle);
+        const double saz = std::sin(angle);
+
+        response[0][0] += caz * P[0];
+        response[0][1] += -saz * P[1];
+        response[1][0] += saz * P[0];
+        response[1][1] += caz * P[1];
+
+        // Update sign and kappa.
+        sign = -sign;
+        kappa += 2;
+    }
+}
+
+} //# namespace LOFAR

CEP/Calibration/ElementResponse/src/convert_coeff.py

+#!/usr/bin/env python3
+
+# Script to convert an ASCII beam model coefficient file to a .cc file for
+# inclusion in the library. Whenever the beam model coefficients file are
+# updated, this script can be run to automatically update the corresponding .cc
+# files.
+#
+# $Id$
+
+import sys
+import time
+import re
+from functools import reduce
+
+def flat_index(shape, index):
+    """
+    Compute the flat index of the element with the provided (N-dimensional)
+    index of an N-dimensional array with the provided shape. Row-major order
+    (or "C"-order) is assumed in the computation of the flattened index. The
+    index is range checked against the provided shape.
+    """
+
+    assert len(shape) == len(index)
+
+    if len(shape) == 0:
+        return 0
+
+    assert index[0] < shape[0]
+    flat = index[0]
+
+    for i in range(1, len(shape)):
+        assert index[i] < shape[i]
+        flat *= shape[i]
+        flat += index[i]
+
+    return flat
+
+def regex(name, type, signed = True):
+    """
+    Return a regular expression to match a (possibly signed) int or float, using
+    the named group syntax. The matching group will by assigned the provided
+    name.
+    """
+
+    expr = None
+    if type == "int":
+        expr = "\d+"
+    elif type == "float":
+        expr = "(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?"
+
+    assert expr, "unknown type: %s" % type
+
+    if signed:
+        expr = "[-+]?" + expr
+
+    return "(?P<%s>%s)" % (name, expr)
+
+def main(args):
+    print("converting %s -> %s (variable name: %s)" % (args[0], args[1], args[2]))
+
+    HEADER, COEFF = list(range(2))
+    state = HEADER
+
+    shape = None
+    coeff = None
+    freqAvg = None
+    freqRange = None
+
+    line_no = 0
+    count = 0
+
+    fin = open(args[0], "r")
+    for line in fin:
+        # Remove leading and trailing whitespace.
+        line = line.strip()
+
+        # Skip empty lines and comments.
+        if len(line) == 0 or line[0] == "#":
+            line_no += 1
+            continue
+
+        # Parse header information.
+        if state == HEADER:
+            match = re.match("^d\s+%s\s+k\s+%s\s+pwrT\s+%s\s+pwrF\s+%s\s+"
+                "freqAvg\s+%s\s+freqRange\s+%s$" % (regex("d", "int", False),
+                regex("k", "int", False), regex("pwrT", "int", False),
+                regex("pwrF", "int", False), regex("freqAvg", "float", False),
+                regex("freqRange", "float", False)), line)
+            assert match, "unable to parse header: \"%s\"" % line
+
+            shape = (int(match.group("k")), int(match.group("pwrT")),
+                int(match.group("pwrF")), int(match.group("d")))
+            assert shape[3] == 2, "unsupported array shape, expected d == 2"
+
+            size = reduce(lambda x, y: x * y, shape)
+            print("coefficient array shape:", shape, "(%d total)" % size)
+
+            freqAvg = match.group("freqAvg")
+            freqRange = match.group("freqRange")
+
+            coeff = [None for x in range(size)]
+            state = COEFF
+
+        # Parse coefficients.
+        elif state == COEFF:
+            match = re.match("^%s\s+%s\s+%s\s+%s\s+%s\s+%s$"
+                % (regex("d", "int", False), regex("k", "int", False),
+                regex("pwrT", "int", False), regex("pwrF", "int", False),
+                regex("re", "float"), regex("im", "float")), line)
+            assert match, "unable to parse line #%d" % line_no
+
+            d = int(match.group("d"))
+            k = int(match.group("k"))
+            pwrT = int(match.group("pwrT"))
+            pwrF = int(match.group("pwrF"))
+
+            index = flat_index(shape, (k, pwrT, pwrF, d))
+            coeff[index] = "std::complex<double>(%s, %s)" % (match.group("re"), match.group("im"))
+
+            count += 1
+
+        line_no += 1
+
+    fin.close()
+    assert not (coeff is None) and all(coeff)
+
+    # Write the output.
+    fout = open(args[1], "w")
+
+    print("// Beam model coefficients converted by convert_coeff.py.", file=fout)
+    print("// Conversion performed on %s UTC using: " % time.strftime("%Y/%m/%d/%H:%M:%S", time.gmtime()), file=fout)
+    print("//     convert_coeff.py %s %s %s" % (args[0], args[1], args[2]), file=fout)
+    print(file=fout)
+    print("#include <complex>", file=fout)
+    print(file=fout)
+    print("// Center frequency, and frequency range for which the beam model coefficients", file=fout)
+    print("// are valid. The beam model is parameterized in terms of a normalized", file=fout)
+    print("// frequency f' in the range [-1.0, 1.0]. The appropriate conversion is:", file=fout)
+    print("//", file=fout)
+    print("//     f' = (f - center) / range", file=fout)
+    print("//", file=fout)
+    print("const double %s_freq_center = %s;" % (args[2], freqAvg), file=fout)
+    print("const double %s_freq_range = %s;" % (args[2], freqRange), file=fout)
+    print(file=fout)
+    print("// Shape of the coefficient array: %dx%dx%dx2 (the size of the last dimension is" % (shape[0], shape[1], shape[2]), file=fout)
+    print("// implied, and always equal to 2).", file=fout)
+    print("//", file=fout)
+    print("const unsigned int %s_coeff_shape[3] = {%d, %d, %d};" % (args[2], shape[0], shape[1], shape[2]), file=fout)
+    print(file=fout)
+    print("// The array of coefficients in row-major order (\"C\"-order).", file=fout)
+    print("//", file=fout)
+    print("const std::complex<double> %s_coeff[%d] = {" % (args[2], len(coeff)), file=fout)
+
+    i = 0
+    while i < len(coeff):
+        if i + 2 < len(coeff):
+            print("    %s, %s," % (coeff[i], coeff[i + 1]), file=fout)
+            i += 2
+        elif i + 2 == len(coeff):
+            print("    %s, %s" % (coeff[i], coeff[i + 1]), file=fout)
+            i += 2
+        else:
+            print("    %s" % coeff[i], file=fout)
+            i += 1
+
+    print("};", file=fout)
+
+    fout.close()
+
+if __name__ == "__main__":
+    if len(sys.argv) != 4:
+        print("convert a beam model coefficient (.coeff) file to a C++ (.cc) file.")
+        print("usage: convert_coeff.py <input-file> <output-file> <variable-name>")
+        sys.exit(1)
+
+    main(sys.argv[1:])

CEP/Calibration/ElementResponse/src/element_beam_HBA.coeff

+d 2 k 2 pwrT 5 pwrF 5 freqAvg 180e6 freqRange 60e6
+0 0 0 0 0.9989322499459223 0.0003305895124867 
+1 0 0 0 1.0030546028872600 0.0002157249025076 
+0 1 0 0 -0.0004916757488397 0.0000266213616248 
+1 1 0 0 0.0006516553273188 -0.0000433166563288 
+0 0 1 0 0.0249658150445263 -0.0122024663463393 
+1 0 1 0 -0.0917825091832822 -0.0062606338208358 
+0 1 1 0 0.0162495046881796 -0.0010736997976255 
+1 1 1 0 -0.0175635905033026 0.0012997068962173 
+0 0 2 0 -0.6730977722052307 0.0940030437973656 
+1 0 2 0 0.3711597596859299 0.0557089394867947 
+0 1 2 0 -0.1492097398080444 0.0123735410547393 
+1 1 2 0 0.1393121453502456 -0.0121117146246749 
+0 0 3 0 0.0665319393516388 -0.1418009730472832 
+1 0 3 0 -0.7576728614553603 -0.0472040122949963 
+0 1 3 0 0.1567990061437258 -0.0143275575385193 
+1 1 3 0 -0.1043118778001582 0.0106756004832779 
+0 0 4 0 0.1121669548152054 0.0504713119323919 
+1 0 4 0 0.1882531376700409 0.0088411256350159 
+0 1 4 0 -0.0459124372023251 0.0044990718645418 
+1 1 4 0 0.0135433541303599 -0.0021789296923529 
+0 0 0 1 0.0003002209532403 0.0007909077657054 
+1 0 0 1 0.0022051270911392 0.0003834815341981 
+0 1 0 1 -0.0004357897643121 0.0000320567996700 
+1 1 0 1 0.0005818285824826 -0.0001021069650381 
+0 0 1 1 -0.0083709499453879 -0.0289759752488368 
+1 0 1 1 -0.0689260153643395 -0.0111348626546314 
+0 1 1 1 0.0138897851110661 -0.0014876219938565 
+1 1 1 1 -0.0150211436594772 0.0029712291209158 
+0 0 2 1 0.2119250520015808 0.2155514942677135 
+1 0 2 1 0.6727380529527980 0.0989550572104158 
+0 1 2 1 -0.1217628319418324 0.0222643129255504 
+1 1 2 1 0.1108579917761457 -0.0262986164183475 
+0 0 3 1 -0.1017024786435272 -0.3302620837788515 
+1 0 3 1 -0.5600906156274197 -0.0797555201430585 
+0 1 3 1 0.1151024257152241 -0.0225518489392044 
+1 1 3 1 -0.0593437249231851 0.0216080058910987 
+0 0 4 1 0.0066968933526899 0.1181452711088882 
+1 0 4 1 0.0981630367567397 0.0129921405004959 
+0 1 4 1 -0.0306136798186735 0.0064963361606382 
+1 1 4 1 -0.0046440676338940 -0.0037281688158807 
+0 0 0 2 -0.0003856663268042 0.0008435910525861 
+1 0 0 2 0.0004887765294093 0.0002777796480946 
+0 1 0 2 -0.0001047488648808 -0.0000302146563592 
+1 1 0 2 0.0001593350153828 -0.0000879125663990 
+0 0 1 2 0.0116296166994115 -0.0307342946951178 
+1 0 1 2 -0.0171249717275797 -0.0080642275561593 
+0 1 1 2 0.0031705620225488 0.0004838463688512 
+1 1 1 2 -0.0034418973689263 0.0024603729467258 
+0 0 2 2 -0.0419944347289523 0.2355624543349744 
+1 0 2 2 0.1917656461134636 0.0732470381581913 
+0 1 2 2 -0.0273147374272124 0.0098595182007132 
+1 1 2 2 0.0208992817013466 -0.0205929453727953 
+0 0 3 2 0.0889729243872774 -0.3406964719938829 
+1 0 3 2 -0.1342560801672904 -0.0515926960946038 
+0 1 3 2 0.0142781186223020 -0.0057037138045721 
+1 1 3 2 0.0151043140114779 0.0141435752121475 
+0 0 4 2 -0.0347327225501659 0.1186585563636635 
+1 0 4 2 0.0102831315790362 0.0046275244914932 
+0 1 4 2 -0.0006372791846825 0.0008894047150233 
+1 1 4 2 -0.0181611528840412 -0.0011106177431486 
+0 0 0 3 -0.0000699366665322 0.0005136144371953 
+1 0 0 3 0.0001520602842105 0.0001303481681886 
+0 1 0 3 -0.0000141882506567 -0.0000941521783975 
+1 1 0 3 -0.0000004226298134 -0.0000245060763932 
+0 0 1 3 0.0012408055399100 -0.0191295543986957 
+1 0 1 3 -0.0051031652662961 -0.0037143632875100 
+0 1 1 3 0.0003028387544878 0.0026905629457281 
+1 1 1 3 0.0006768121359769 0.0005901486396051 
+0 0 2 3 0.0048918921441903 0.1588912409502319 
+1 0 2 3 0.0575369727210951 0.0344677222786687 
+0 1 2 3 -0.0002152227668601 -0.0089220757225133 
+1 1 2 3 -0.0074792188817697 -0.0043562231368076 
+0 0 3 3 -0.0149335262655201 -0.2084962323582034 
+1 0 3 3 -0.0327252678958813 -0.0172371907472848 
+0 1 3 3 -0.0057143555179676 0.0141142700941743 
+1 1 3 3 0.0251435557201315 -0.0005753615445942 
+0 0 4 3 0.0070209144233666 0.0689639468490938 
+1 0 4 3 -0.0020239346031291 -0.0025499069613344 
+0 1 4 3 0.0032325387394458 -0.0048123509184894 
+1 1 4 3 -0.0136340313457176 0.0021185000810664 
+0 0 0 4 0.0000512381993616 0.0001550785137302 
+1 0 0 4 0.0000819244737818 0.0000466470412396 
+0 1 0 4 -0.0000177429496833 -0.0000561890408003 
+1 1 0 4 -0.0000018388829279 0.0000032387726477 
+0 0 1 4 -0.0022414352263751 -0.0060474723525871 
+1 0 1 4 -0.0024377933436567 -0.0012852163337395 
+0 1 1 4 0.0004634797107989 0.0016976603895716 
+1 1 1 4 0.0003344773954073 -0.0001499932789294 
+0 0 2 4 0.0241014578366618 0.0547046570516960 
+1 0 2 4 0.0219986510834463 0.0112189146988984 
+0 1 2 4 -0.0012019994038721 -0.0079939660050373 
+1 1 2 4 -0.0035807498769946 0.0014801422733613 
+0 0 3 4 -0.0362395089905272 -0.0661322227928722 
+1 0 3 4 -0.0141568558526096 -0.0042676979206835 
+0 1 3 4 0.0004475745352473 0.0102135659618127 
+1 1 3 4 0.0090474375150397 -0.0032177128650026 
+0 0 4 4 0.0132702874173192 0.0207916487187541 
+1 0 4 4 0.0004387107229914 -0.0017223838914815 
+0 1 4 4 0.0001287985092565 -0.0032079544559908 
+1 1 4 4 -0.0045503800737417 0.0015366231416036 

CEP/Calibration/ElementResponse/src/element_beam_LBA.coeff

+d 2 k 2 pwrT 5 pwrF 5 freqAvg 55e6 freqRange 45e6
+ 0 0 0 0        0.9982445079290715     0.0000650863154389  
+ 1 0 0 0        1.0006230902158257    -0.0022053287681416  
+ 0 1 0 0       -0.0000762326746410     0.0000118155774181  
+ 1 1 0 0        0.0000118903581604    -0.0000251324432498  
+ 0 0 1 0        0.0550606068782634     0.0011958385659938  
+ 1 0 1 0       -0.0160912944232080     0.0703645376267940  
+ 0 1 1 0        0.0029807707669629    -0.0003262084082071  
+ 1 1 1 0        0.0001639620574332    -0.0000266272685197  
+ 0 0 2 0       -0.8624889445327827    -0.1522883072804402  
+ 1 0 2 0       -0.0386800869486029    -0.7569350701887934  
+ 0 1 2 0       -0.0229975846564195     0.0010565261888195  
+ 1 1 2 0        0.0019573207027441     0.0050550600926414  
+ 0 0 3 0        0.4039567648146965     0.0810473144253429  
+ 1 0 3 0       -0.0350803390479135     0.5214591717801087  
+ 0 1 3 0        0.0333560984394624    -0.0009767086536162  
+ 1 1 3 0        0.0141704479374002    -0.0205386534626779  
+ 0 0 4 0       -0.0756250564248881     0.0056622911723172  
+ 1 0 4 0       -0.1267444401630109    -0.0349676272376008  
+ 0 1 4 0       -0.0126732549319889     0.0002042370658763  
+ 1 1 4 0       -0.0101360135082899     0.0114084024114141  
+ 0 0 0 1        0.0001002692200362     0.0006838211278268  
+ 1 0 0 1       -0.0003660049052840    -0.0008418920419220  
+ 0 1 0 1       -0.0002204197663391    -0.0000213776348027  
+ 1 1 0 1        0.0001477083861977     0.0000599750510518  
+ 0 0 1 1        0.0033849565901213    -0.0244636379385135  
+ 1 0 1 1        0.0234264238829944     0.0084068836453700  
+ 0 1 1 1        0.0076282580587895     0.0026614359017468  
+ 1 1 1 1       -0.0044850263974801    -0.0058337192660638  
+ 0 0 2 1        0.0891332399420108     0.1876527151756476  
+ 1 0 2 1       -0.1012363483900640    -0.1975118891151966  
+ 0 1 2 1       -0.0682274156850413    -0.0758159820140411  
+ 1 1 2 1        0.0497303968865466     0.1019681987654797  
+ 0 0 3 1        0.2232030356124932    -0.2248154851829713  
+ 1 0 3 1        0.4704343293662089    -0.3552101485419532  
+ 0 1 3 1        0.0562541280098909     0.0743149092143081  
+ 1 1 3 1       -0.0226634801339250    -0.1439026188572270  
+ 0 0 4 1       -0.1793752883639813     0.0720222655359702  
+ 1 0 4 1       -0.2678542619793421     0.3152115802895427  
+ 0 1 4 1       -0.0102147881225462    -0.0176848554302252  
+ 1 1 4 1       -0.0051268936720694     0.0527621533959941  
+ 0 0 0 2       -0.0010581424498791     0.0015237878543047  
+ 1 0 0 2        0.0007398729642721     0.0028468649470433  
+ 0 1 0 2       -0.0003281057522772    -0.0000770207588466  
+ 1 1 0 2        0.0003478997686964     0.0001481982639746  
+ 0 0 1 2        0.0557107413978542    -0.0634701730653090  
+ 1 0 1 2       -0.0139549526991330    -0.1175401658864208  
+ 0 1 1 2        0.0124258438959177     0.0067985224235178  
+ 1 1 1 2       -0.0126349778957970    -0.0100656881493938  
+ 0 0 2 2       -0.6404795825927633     0.7568775384981410  
+ 1 0 2 2        0.0767245154665722     1.3441875993523555  
+ 0 1 2 2       -0.1757936183439326    -0.1363710820472197  
+ 1 1 2 2        0.1765450269056824     0.1555919358121995  
+ 0 0 3 2        0.9646419509627557    -0.8095088593139815  
+ 1 0 3 2        0.1635280638865702    -1.4854352979459096  
+ 0 1 3 2        0.1238595124159999     0.1766108700786397  
+ 1 1 3 2       -0.1307647072780430    -0.2090615438301942  
+ 0 0 4 2       -0.3718069213271066     0.2275266747872172  
+ 1 0 4 2       -0.1372223722572021     0.4314989948093362  
+ 0 1 4 2       -0.0110701836450407    -0.0593085026046026  
+ 1 1 4 2        0.0140598301629874     0.0738668439833535  
+ 0 0 0 3       -0.0039458389254656    -0.0007048354913730  
+ 1 0 0 3        0.0007040177887611     0.0007856369612188  
+ 0 1 0 3       -0.0000625695757282     0.0000249138990722  
+ 1 1 0 3       -0.0000960097542525     0.0002521364065803  
+ 0 0 1 3        0.1336911750356096     0.0202651327657687  
+ 1 0 1 3       -0.0113385668361727    -0.0339262369086247  
+ 0 1 1 3        0.0059031372522229     0.0008660479915339  
+ 1 1 1 3        0.0039660364524413    -0.0100356333791398  
+ 0 0 2 3       -0.8758406699506004     0.3350237639226141  
+ 1 0 2 3        0.2824832769101577     0.6821307442669313  
+ 0 1 2 3       -0.1541299429420569    -0.0281422177614844  
+ 1 1 2 3        0.0816399676454817     0.0691599035109852  
+ 0 0 3 3        1.0331569921006993     0.0509705885336283  
+ 1 0 3 3        0.1501121326521990    -0.5193414816770609  
+ 0 1 3 3        0.1557916917691289     0.0646351862895731  
+ 1 1 3 3        0.0170294191358757    -0.1027926845803498  
+ 0 0 4 3       -0.3316657641578328    -0.1655909947939444  
+ 1 0 4 3       -0.2158100484836540     0.0614504774034524  
+ 0 1 4 3       -0.0389912915621699    -0.0301364165752433  
+ 1 1 4 3       -0.0462331759359031     0.0405864871628086  
+ 0 0 0 4       -0.0031701591723043    -0.0010521154166512  
+ 1 0 0 4       -0.0007213036752903    -0.0007227764008022  
+ 0 1 0 4        0.0001275344578325     0.0000652362392482  
+ 1 1 0 4       -0.0003113309221942     0.0001956734476566  
+ 0 0 1 4        0.0962263571740972     0.0440074333288440  
+ 1 0 1 4        0.0313595045238824     0.0230763038515351  
+ 0 1 1 4       -0.0020520685193773    -0.0028564379463666  
+ 1 1 1 4        0.0121039958869239    -0.0059701468961263  
+ 0 0 2 4       -0.3144282315609649    -0.2763869580286276  
+ 1 0 2 4       -0.1705959031354030    -0.0712085950559831  
+ 0 1 2 4       -0.0235110916473515     0.0306385386726702  
+ 1 1 2 4       -0.0474273292450285     0.0116831908947225  
+ 0 0 3 4        0.4715775965513117     0.5077361528286819  
+ 1 0 3 4        0.3847391427972284     0.1136717951238837  
+ 0 1 3 4        0.0543537332385954    -0.0366524906364179  
+ 1 1 3 4        0.1127180664279469    -0.0176607923511174  
+ 0 0 4 4       -0.1901597954359592    -0.2294955549701665  
+ 1 0 4 4       -0.1864961465389693    -0.0486276177310768  
+ 0 1 4 4       -0.0251598701859194     0.0115712688652445  
+ 1 1 4 4       -0.0563476280247398     0.0079787883434624

CEP/Calibration/ExpIon/CMakeLists.txt

+# $Id: CMakeLists.txt 14609 2009-12-07 09:10:48Z loose $
+
+# Do not split the following line, otherwise makeversion will fail!
+lofar_package(ExpIon 1.0 DEPENDS pyparameterset pyparmdb)
+
+include(LofarFindPackage)
+lofar_find_package(Boost REQUIRED COMPONENTS python thread)
+lofar_find_package(Casacore REQUIRED COMPONENTS python scimath)
+add_subdirectory(src)

CEP/Calibration/ExpIon/src/CMakeLists.txt

+# $Id: CMakeLists.txt 14508 2010-03-24 13:18:18Z vdtol $
+
+include(PythonInstall)
+
+lofar_add_bin_scripts(
+  calibrate-ion
+  readms-part.sh
+  readms-part.py
+  parmdbwriter.py)
+
+lofar_add_library(_baselinefitting MODULE baselinefitting.cc)
+
+set_target_properties(_baselinefitting PROPERTIES
+  PREFIX ""
+  LIBRARY_OUTPUT_DIRECTORY ${PYTHON_BUILD_DIR}/lofar/expion)
+
+# This is a quick-and-dirty fix to install the Python binding module in the
+# right place. It will now be installed twice, because lofar_add_library()
+# will install it in $prefix/$libdir
+install(TARGETS _baselinefitting
+  DESTINATION ${PYTHON_INSTALL_DIR}/lofar/expion)
+
+# Python modules.
+python_install(
+  __init__.py
+  io.py
+  format.py
+  ionosphere.py
+  acalc.py
+  client.py
+  sphere.py
+  error.py
+  mpfit.py
+  readms.py
+  parmdbmain.py
+  baselinefitting.py
+  repairGlobaldb.py
+  MMionosphere.py
+  fitClockTEC.py
+  PosTools.py
+  read_sagecal.py
+  DESTINATION lofar/expion)
+

CEP/Calibration/ExpIon/src/PosTools.py

+from pyrap.measures import measures
+from math import *
+import pyrap.quanta as qa;
+import os
+from pyrap import tables as tab
+import numpy as np
+
+R_earth=6364.62e3;
+earth_ellipsoid_a = 6378137.0;
+earth_ellipsoid_a2 = earth_ellipsoid_a*earth_ellipsoid_a;
+earth_ellipsoid_b = 6356752.3142;
+earth_ellipsoid_b2 = earth_ellipsoid_b*earth_ellipsoid_b;
+earth_ellipsoid_e2 = (earth_ellipsoid_a2 - earth_ellipsoid_b2) / earth_ellipsoid_a2;
+posCS002=[3826577.1095  ,461022.900196, 5064892.758]
+
+def getMSinfo(MS=None):
+    print("getting info for",MS)
+    if MS is None:
+        print("No measurement set given")
+        return
+    if os.path.isdir(MS):
+        myMS=tab.table(MS)
+    else:
+        print("Do not understand the format of MS",MS,"bailing out")
+        return;
+    print("opened table",MS)
+    timerange=[np.amin(myMS.getcol('TIME_CENTROID')),np.amax(myMS.getcol('TIME_CENTROID'))]
+    timestep=myMS.getcell('INTERVAL',0)
+    
+    pointing= tab.table(myMS.getkeyword('FIELD')).getcell('PHASE_DIR',0);    
+    stations = tab.table(myMS.getkeyword('ANTENNA')).getcol('NAME')
+    station_pos = tab.table(myMS.getkeyword('ANTENNA')).getcol('POSITION')
+
+    return (timerange,timestep,pointing.flatten(),stations,station_pos)
+
+def getPPsimple(height=[450.e3,],mPosition=[0.,0.,0.],direction=[0.,0.,0.]):
+    '''get piercepoints for antenna position mPosition in m, direction ITRF in m on unit sphere and for array of heights, assuming a spherical Earth'''
+    height=np.array(height)
+    stX=mPosition[0]
+    stY=mPosition[1]
+    stZ=mPosition[2]
+    x=np.divide(stX,(R_earth+height))
+    y=np.divide(stY,(R_earth+height))
+    z=np.divide(stZ,(R_earth+height))
+    c = x*x + y*y + z*z - 1.0;
+   
+    dx=np.divide(direction[0],(R_earth+height))
+    dy=np.divide(direction[1],(R_earth+height))
+    dz=np.divide(direction[2],(R_earth+height))
+
+    a = dx*dx + dy*dy + dz*dz;
+    b = x*dx + y*dy  + z*dz;
+
+    alpha = (-b + np.sqrt(b*b - a*c))/a;
+
+    pp=np.zeros(height.shape+(3,))
+    pp[:,0]=stX+alpha*direction[0]
+    pp[:,1]=stY+alpha*direction[1]
+    pp[:,2]=stZ+alpha*direction[2]
+
+    am=np.divide(1.,pp[:,0]*dx+pp[:,1]*dy+pp[:,2]*dz)
+    return pp,am
+
+def getPPsimpleAngle(height=[450.e3,],mPosition=[0.,0.,0.],direction=[0.,0.,0.]):
+    '''get (lon,lat,h values) of piercepoints for antenna position mPosition in m, direction ITRF in m on unit sphere and for array of heights, assuming a spherical Earth'''
+    height=np.array(height)
+    stX=mPosition[0]
+    stY=mPosition[1]
+    stZ=mPosition[2]
+    x=np.divide(stX,(R_earth+height))
+    y=np.divide(stY,(R_earth+height))
+    z=np.divide(stZ,(R_earth+height))
+    
+    c = x*x + y*y + z*z - 1.0;
+   
+    dx=np.divide(direction[0],(R_earth+height))
+    dy=np.divide(direction[1],(R_earth+height))
+    dz=np.divide(direction[2],(R_earth+height))
+
+    a = dx*dx + dy*dy + dz*dz;
+    b = x*dx + y*dy  + z*dz;
+
+    alpha = (-b + np.sqrt(b*b - a*c))/a;
+
+
+    pp=np.zeros(height.shape+(3,))
+    pp[:,0]=stX+alpha*direction[0]
+    pp[:,1]=stY+alpha*direction[1]
+    pp[:,2]=stZ+alpha*direction[2]
+
+    am=np.divide(1.,pp[:,0]*dx+pp[:,1]*dy+pp[:,2]*dz)
+
+    ppl=np.zeros(height.shape+(3,))
+    ppl[:,0]=np.atan2(pp[:,1],pp[:0])
+    ppl[:,1]=np.atan2(pp[:,2],np.sqrt(pp[:0]*pp[:,0]+pp[:1]*pp[:,1]))
+    ppl[:,2]=heigth
+
+    return ppl,am
+
+def getPP(h=450e3,mPosition=[0.,0.,0.],direction=[0.,0.,0.]):
+    stationX = mPosition[0];
+    stationY = mPosition[1];
+    stationZ = mPosition[2];
+
+    ion_ellipsoid_a = earth_ellipsoid_a + h;
+    ion_ellipsoid_a2_inv = 1.0 / (ion_ellipsoid_a * ion_ellipsoid_a);
+    ion_ellipsoid_b = earth_ellipsoid_b + h;
+    ion_ellipsoid_b2_inv = 1.0 / (ion_ellipsoid_b * ion_ellipsoid_b);
+    
+    x = stationX/ion_ellipsoid_a;
+    y = stationY/ion_ellipsoid_a;
+    z = stationZ/ion_ellipsoid_b;
+    c = x*x + y*y + z*z - 1.0;
+
+    dx = direction [0]/ ion_ellipsoid_a;
+    dy = direction [1] / ion_ellipsoid_a;
+    dz = direction [2] / ion_ellipsoid_b;
+
+    a = dx*dx + dy*dy + dz*dz;
+    b = x*dx + y*dy  + z*dz;
+    alpha = (-b + sqrt(b*b - a*c))/a;
+    pp_x = stationX + alpha*direction[0];
+    pp_y = stationY + alpha*direction[1]
+    pp_z = stationZ + alpha*direction[2];
+
+    normal_x = pp_x * ion_ellipsoid_a2_inv;
+    normal_y = pp_y * ion_ellipsoid_a2_inv;
+    normal_z = pp_z * ion_ellipsoid_b2_inv;
+    norm_normal2 = normal_x*normal_x + normal_y*normal_y + normal_z*normal_z;
+    norm_normal = sqrt(norm_normal2);
+    sin_lat2 = normal_z*normal_z / norm_normal2;
+
+ 
+    g = 1.0 - earth_ellipsoid_e2*sin_lat2;
+    sqrt_g = sqrt(g);
+
+    M = earth_ellipsoid_b2 / ( earth_ellipsoid_a * g * sqrt_g );
+    N = earth_ellipsoid_a / sqrt_g;
+
+    local_ion_ellipsoid_e2 = (M-N) / ((M+h)*sin_lat2 - N - h);
+    local_ion_ellipsoid_a = (N+h) * sqrt(1.0 - local_ion_ellipsoid_e2*sin_lat2);
+    local_ion_ellipsoid_b = local_ion_ellipsoid_a*sqrt(1.0 - local_ion_ellipsoid_e2);
+
+    z_offset = ((1.0-earth_ellipsoid_e2)*N + h - (1.0-local_ion_ellipsoid_e2)*(N+h)) * sqrt(sin_lat2);
+
+    x1 = stationX/local_ion_ellipsoid_a;
+    y1 = stationY/local_ion_ellipsoid_a;
+    z1 = (stationZ-z_offset)/local_ion_ellipsoid_b;
+    c1 = x1*x1 + y1*y1 + z1*z1 - 1.0;
+
+    dx = direction[0] / local_ion_ellipsoid_a;
+    dy = direction[1] / local_ion_ellipsoid_a;
+    dz = direction[2] / local_ion_ellipsoid_b;
+    a = dx*dx + dy*dy + dz*dz;
+    b = x1*dx + y1*dy  + z1*dz;
+    alpha = (-b + sqrt(b*b - a*c1))/a;
+
+    pp_x = stationX + alpha*direction[0];
+    pp_y = stationY + alpha*direction[1]
+    pp_z = stationZ + alpha*direction[2];
+
+    normal_x = pp_x / (local_ion_ellipsoid_a * local_ion_ellipsoid_a);
+    normal_y = pp_y / (local_ion_ellipsoid_a * local_ion_ellipsoid_a);
+    normal_z = (pp_z-z_offset) / (local_ion_ellipsoid_b * local_ion_ellipsoid_b);
+
+    norm_normal2 = normal_x*normal_x + normal_y*normal_y + normal_z*normal_z;
+    norm_normal = sqrt(norm_normal2);
+    
+    pp_airmass = norm_normal / (direction[0]*normal_x + direction[1]*normal_y + direction[2]*normal_z);
+
+    return (np.array([[pp_x,pp_y,pp_z]]),pp_airmass)
+
+class PPdummy:
+    pass
+
+def getPiercePoints(time,source_positions,station_positions,height=450.e3):
+    me=measures()
+    piercepoints=PPdummy()
+    piercepoints.positions_xyz=np.zeros((source_positions.shape[0],station_positions.shape[0],3),dtype=np.float64)
+    piercepoints.positions=np.zeros((source_positions.shape[0],station_positions.shape[0],2),dtype=np.float64)
+    piercepoints.zenith_angles=np.zeros((source_positions.shape[0],station_positions.shape[0]),dtype=np.float64)
+    for isrc in range(source_positions.shape[0]) :
+        radec=source_positions[isrc]
+        for istat in range(station_positions.shape[0]):
+            stat_pos=station_positions[istat]
+            azel=radec2azel(radec[0],radec[1],time=str(time)+'s',pos=stat_pos)
+            az=azel['m0']['value'];
+            el=azel['m1']['value'];
+            lonlat=getLonLatStation(az,el,pos=stat_pos);
+
+            lon=lonlat['m0']['value'];
+            lat=lonlat['m1']['value'];
+
+            # convert to itrf coordinates on sphere with radius 1
+            diritrf=[cos(lat)*cos(lon),cos(lat)*sin(lon),sin(lat)]
+            (pp,am)=getPP(h=height,mPosition=stat_pos,direction=diritrf)
+            piercepoints.positions_xyz[isrc,istat]=np.array(pp[0])
+            piercepoints.zenith_angles[isrc,istat]=np.arccos(1./am)
+            pp1position=me.position("ITRF",str(pp[0,0])+'m',str(pp[0,1])+'m',str(pp[0,2])+'m')
+            # print "pp",  degrees(pp1position['m0']['value']),degrees(pp1position['m1']['value'])
+            piercepoints.positions[isrc,istat,0] = pp1position['m0']['value'];
+            piercepoints.positions[isrc,istat,1] = pp1position['m1']['value'];
+    return piercepoints
+
+def getLonLat(pos):
+    #converts ITRF pos in xyz to lon lat 
+    me=measures()
+    a=me.measure(me.position('ITRF',str(pos[0])+'m',str(pos[1])+'m',str(pos[2])+'m'),"ITRF");
+    return (a['m0']['value'],a['m1']['value'])
+
+def getLonLatStation(az=0,el=0,pos=posCS002):
+    #gets converts local station direction to ITRF lon/lat
+    if not isinstance(az,str):
+        az=str(az)+'rad';
+    if not isinstance(el,str):
+        el=str(el)+'rad';
+    me=measures()
+    me.do_frame(me.position('ITRF',str(pos[0])+'m',str(pos[1])+'m',str(pos[2])+'m'))
+    #me.do_frame(me.epoch('utc', 'today'))
+    direction=me.direction("AZEL",az,el)
+    return me.measure(direction,"ITRF");
+
+
+def radec2azel(ra,dec,time, pos):
+    me=measures();
+    if type(ra)!=str:
+        ra=str(ra)+'rad';
+    if type(dec)!=str:
+        dec=str(dec)+'rad';
+    phasedir=me.direction('J2000',ra,dec)
+    t=me.epoch("UTC",qa.quantity(time));
+    me.do_frame(t);
+
+    p = me.position('ITRF',str(pos[0])+'m',str(pos[1])+'m',str(pos[2])+'m')
+    me.do_frame(p);
+
+    azel = me.measure(phasedir,'azel');
+    return azel;
+
+def azel2radec(az,el,time, pos):
+    me=measures();
+    if type(az)!=str:
+        az=str(az)+'rad';
+    if type(el)!=str:
+        el=str(el)+'rad';
+    phasedir=me.direction('AZEL',az,el)
+    t=me.epoch("UTC",qa.quantity(time));
+    me.do_frame(t);
+
+    p = me.position('ITRF',str(pos[0])+'m',str(pos[1])+'m',str(pos[2])+'m')
+    me.do_frame(p);
+
+    radec = me.measure(phasedir,'RADEC');
+    return radec;
+
+
+
+def getStatPos(stations,AFPath='/opt/lofar/etc/StaticMetaData/',Field='LBA'):
+    StatPos=[];
+    for st in stations:
+        antFile=open(AFPath+st+"-AntennaField.conf")
+        line=antFile.readline()
+        while len(line)>0 and line[:len(Field)]!=Field:
+        
+            line=antFile.readline()
+        if line[:len(Field)]==Field:
+            StatPos.append([float(stp) for stp in antFile.readline().split()[2:5]])
+        else:
+            print("Field",Field,"for",st,"not found,putting zeros")
+            StatPos.append([0,0,0])
+        antFile.close()
+    return StatPos

CEP/Calibration/ExpIon/src/__init__.py

+# -*- coding: iso-8859-1 -*-
+# __init__.py: Top level .py file for python solution analysis tools.
+#
+# Copyright (C) 2010
+# ASTRON (Netherlands Institute for Radio Astronomy)
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite.
+# The LOFAR software suite is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+#
+# $Id: __init__.py 12729 2010-03-24 13:39:59Z vdtol $
+
+#from ionosphere import *
+
+__all__ = ['ionosphere', 'parmdbmain']
\ No newline at end of file

CEP/Calibration/ExpIon/src/acalc.py

+###############################################################################
+
+# import Python modules
+from math import pi, atan2, degrees, radians
+
+# import 3rd party modules
+from numpy import *
+
+# import user modules
+from .error import *
+try:
+  import _acalc
+except:
+  __acalc = False
+else:
+  __acalc = True
+
+###############################################################################
+# NOTE: there are several cases where one needs to preserve the sign of zero,
+# e.g. +0 and -0. 
+###############################################################################
+
+def factorial( n ):
+  nn = int( n )
+  if ( nn < 0 ):
+    fac = None
+  else:
+    if __acalc:
+      fac = _acalc.factorial( nn )
+    else:
+      fac = int( 1 )
+      while ( nn > 0 ):
+        fac = fac * int( nn )
+        nn = nn - 1
+  return fac
+
+###############################################################################
+
+def binomial( n, k ):
+  nn = int( n )
+  kk = int( k )
+  if ( kk < 0 ) or ( kk > nn ):
+    binom = None
+  else:
+    if __acalc:
+      binom = _acalc.binomial( nn, kk )
+    else:
+      binom = factorial( nn ) / ( factorial( kk ) * factorial( nn - kk ) )
+  return binom
+
+###############################################################################
+
+def complex_to_r_phi( c ):
+  if __acalc:
+    [ r, phi ] = _acalc.complex_to_r_phi( [ c.real, c.imag ] )
+  else:
+    r = abs( c )
+    phi = degrees( log( c ).imag )
+  return [ r, phi ]
+
+###############################################################################
+
+def r_phi_to_complex( rp ):
+  if __acalc:
+    [ cr, ci ] = _acalc.r_phi_to_complex( rp )
+    c = complex( cr, ci )
+  else:
+    [ r, phi ] = rp
+    c = r * exp( complex( 0., 1. ) * radians( phi ) )
+  return c
+
+###############################################################################
+
+def is_array( a ):
+  return isinstance( a, type( array( [ 1 ] ) ) )
+
+###############################################################################
+
+def azeros( x ):
+  if ( len( shape( x ) ) == 0 ):
+    zero = 0.
+  else:
+    zero = zeros( shape = x.shape, dtype = x.dtype )
+  return zero
+
+###############################################################################
+
+def aones( x ):
+  if ( len( x.shape ) == 0 ):
+    one = 1.
+  else:
+    one = ones( shape = x.shape, dtype = x.dtype )
+  return one
+
+###############################################################################
+
+def aatan2( y, x ):
+  if ( shape( x ) != shape( y ) ):
+    raise error( 'x and y have different shapes' )
+  if ( len( shape( x ) ) == 0 ):
+    z = atan2( y, x )
+  else:
+    xx = x.ravel()
+    yy = y.ravel()
+    zz = array( [ atan2( yy[ i ], xx[ i ] ) for i in range( len( xx ) ) ], dtype = x.dtype )
+    z = zz.reshape( x.shape )
+  return z
+
+###############################################################################
+def asign( x ):
+# this function also separates between -0 and +0
+  if ( not is_array( x ) ):
+    s = ( - 2. * float( aatan2( x, x ) < 0. ) + 1. )
+  else:
+    s = ( - 2. * array( aatan2( x, x ) < 0., dtype = x.dtype ) + 1. )
+  return s
+
+###############################################################################
+
+def amodulo( x, y ):
+  if ( not is_array( x ) ):
+    if ( not is_array( y ) ):
+      m = x - y * floor( x / ( y + float( y == 0. ) ) )
+    else:
+      xx = x * aones( y )
+      m = xx - y * floor( x / ( y + array( y == 0., dtype = y.dtype ) ) )
+  else:
+    if ( not is_array( y ) ):
+      yy = y * aones( x )
+      m = x - yy * floor( x / ( yy + array( yy == 0., dtype = yy.dtype ) ) )
+    else:
+      m = x - y * floor( x / ( y + array( y == 0., dtype = y.dtype ) ) )
+  return m
+
+###############################################################################
+
+def aradians( x ):
+  r = x * ( pi / 180. )
+  return r
+
+###############################################################################
+
+def adegrees( x ):
+  r = x * ( 180. / pi )
+  return r
+
+###############################################################################
+
+def awhere( a ):
+   return transpose( array( where( a ) ) )
+
+###############################################################################
+
+def aput( data, sel, sub ):
+# data, sel and sub must be arrays
+
+  # check input dimensions
+  if ( len( sel ) == 0 ):
+    return data.copy()
+  if ( len( sel.ravel() ) == 0 ):
+    return data.copy()
+  if ( sel.shape[ 1 ] > len( data.shape ) ):
+    raise error( 'number of dimensions of index array is higher than that of data array' )
+  asub = array( sub )
+  if ( len( asub.shape ) == len( data.shape ) - sel.shape[ 1 ] ):
+    if ( asub.shape != data.shape[ sel.shape[ 1 ] : ] ):
+      raise error( 'shape of subarray does not match selected data' )
+    asub = resize( asub, [ sel.shape[ 0 ] ] + list( data.shape[ sel.shape[ 1 ] : ] ) )
+  elif ( len( asub.shape ) == len( data.shape ) - sel.shape[ 1 ] + 1 ):
+    if ( list( asub.shape ) != [ sel.shape[ 0 ] ] + list( data.shape[ sel.shape[ 1 ] : ] ) ):
+      raise error( 'shape of subarray does not match selected data' )
+
+  # collapse and write data
+  coffset = [ int( product( data.shape[ i : sel.shape[ 1 ] ] ) ) for i in range( 1, 1 + sel.shape[ 1 ] ) ]
+  coffset[ - 1 ] = 1
+  csel = add.reduce( sel * array( coffset ), 1 )
+  subsize = int( product( data.shape[ sel.shape[ 1 ] : ] ) )
+  suboffset = resize( arange( subsize ), [ sel.shape[ 0 ], subsize ] )
+  csel = ( transpose( resize( csel * subsize, [ subsize, sel.shape[ 0 ] ] ) ) + suboffset ).ravel()
+  cdata = data.copy().ravel()
+  csub = asub.ravel()
+  put( cdata, csel, csub )
+
+  return cdata.reshape( data.shape )
+
+###############################################################################
+
+def aget( data, sel ):
+# data and sel must be arrays
+
+  # check input dimensions
+  if ( len( sel ) == 0 ):
+    return array( [], dtype = data.dtype )
+  if ( len( sel.ravel() ) == 0 ):
+    return array( [], dtype = data.dtype )
+  if ( sel.shape[ 1 ] > len( data.shape ) ):
+    raise error( 'number of dimensions of index array is higher than that of data array' )
+
+  # collapse data along sel axes
+  cdata_len = int( product( data.shape[ 0 : sel.shape[ 1 ] ] ) )
+  cdata_shape = [ cdata_len ] + list( data.shape[ sel.shape[ 1 ] : ] )
+  cdata = data.reshape( cdata_shape )
+  coffset = [ int( product( data.shape[ i : sel.shape[ 1 ] ] ) ) for i in range( 1, 1 + sel.shape[ 1 ] ) ]
+  coffset[ - 1 ] = 1
+  csel = add.reduce( sel * array( coffset ), 1 )
+
+  return take( cdata, csel, axis = 0 ) 
+
+###############################################################################
+
+def amean_phase( data ):
+  # determine phase average (with trick to get around possible phase wrap problems)
+  phases = array( [ data ], dtype = float64 ).ravel()
+  offsets = [ 0., 90., 180., 270. ]
+  mean_phases = [ ( amodulo( ( phases + offsets[ j ] ) + 180., 360. ) - 180. ).mean()
+      for j in range( len( offsets ) ) ]
+  var_phases = [ ( ( ( amodulo( ( phases + offsets[ j ] ) + 180., 360. ) - 180. ) - 
+      mean_phases[ j ] )**2 ).mean() for j in range( len( offsets ) ) ]
+  j = var_phases.index( min( var_phases ) )
+  mean_phase = mean_phases[ j ] - offsets[ j ]
+  return float( amodulo( mean_phase + 180, 360. ) - 180. )
+
+###############################################################################
+
+def amedian_phase( data ):
+  # determine phase average (with trick to get around possible phase wrap problems)
+  mean_phase = amean_phase( data )
+  phases = array( [ data ], dtype = float64 ).ravel()
+  median_phase = median( amodulo( ( phases - mean_phase ) + 180., 360. ) - 180. )
+  median_phase = amodulo( ( median_phase + mean_phase ) + 180., 360. ) - 180.
+  return float( median_phase )
+
+###############################################################################
+

CEP/Calibration/ExpIon/src/baselinefitting.cc

+//# fitting.cc: Clock and TEC fitting using cascore 
+//# 
+//#
+//# Copyright (C) 2012
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id: $
+
+#include <lofar_config.h>
+#include <Common/OpenMP.h>
+
+#include <casacore/casa/Containers/ValueHolder.h>
+#include <casacore/casa/Containers/Record.h>
+#include <casacore/casa/Utilities/CountedPtr.h>
+#include <casacore/scimath/Fitting/LSQFit.h>
+
+#if defined(HAVE_CASACORE)
+#include <casacore/python/Converters/PycExcp.h>
+#include <casacore/python/Converters/PycBasicData.h>
+#include <casacore/python/Converters/PycValueHolder.h>
+#include <casacore/python/Converters/PycRecord.h>
+#define pyrap python
+#else
+#include <pyrap/Converters/PycExcp.h>
+#include <pyrap/Converters/PycBasicData.h>
+#include <pyrap/Converters/PycValueHolder.h>
+#include <pyrap/Converters/PycRecord.h>
+#endif
+
+#include <boost/python.hpp>
+#include <boost/python/args.hpp>
+
+using namespace casacore;
+using namespace boost::python;
+
+namespace LOFAR
+{
+namespace ExpIon
+{
+    
+ValueHolder fit(const ValueHolder &phases_vh, const ValueHolder &A_vh, const ValueHolder &init_vh, 
+                const ValueHolder &flags_vh, const ValueHolder &constant_parm_vh ) 
+{
+    // Arrays are passed as ValueHolder
+    // They are Array is extracted by the asArray<Type> methods
+    // They pointer to the actual data is obtained by the Array.data() method.
+    
+    // Get the phases 
+    Array<Float> phases = phases_vh.asArrayFloat();
+    Float *phases_data = phases.data();
+
+    // Get the flags
+    Array<Bool> flags(flags_vh.asArrayBool());
+    Bool noflags = (flags.ndim() == 0);
+    Bool *flags_data;
+    if (!noflags)
+    {
+        flags_data = flags.data();
+    }
+    
+    IPosition s = phases.shape();
+    int N_station = s[0];
+    int N_freq = s[1];
+    
+    // Get the matrix with basis functions
+    Array<Float> A = A_vh.asArrayFloat();
+    Float *A_data = A.data();
+    
+    IPosition s1 = A.shape();
+    
+    int N_coeff = s1[0];
+    int N_parm = N_station * N_coeff;
+    Float sol[N_parm];
+    
+    Array<Float> init(init_vh.asArrayFloat());
+    if (init.ndim() == 0)
+    {
+        for(int i = 0; i < N_parm; i++) sol[i] = 0.0;
+    }
+    else    
+    {
+        for(int i = 0; i < N_parm; i++) sol[i] = init.data()[i];
+    }
+    
+    // Get the flags indicating which parameters are constant_parm
+    // i.e. are not a free parameter in the minimization problem
+    Array<Bool> constant_parm(constant_parm_vh.asArrayBool());
+    Bool no_constant_parm = (constant_parm.ndim() == 0);
+    Bool *constant_parm_data;
+    if (!no_constant_parm)
+    {
+        constant_parm_data = constant_parm.data();
+    }
+    
+    Float cEq[N_parm];
+    for(int i=0; i<N_parm; ++i) cEq[i] = 0.0;
+    
+    int N_thread = OpenMP::maxThreads();
+    std::vector<LSQFit> lnl(N_thread);
+    
+    uInt nr = 0;
+    
+    for (int iter = 0; iter<1000; iter++)
+    {
+        for(int i = 0; i<N_thread; i++) {
+            lnl[i] = LSQFit(N_parm);
+        }
+        #pragma omp parallel
+        {
+            int threadNum = OpenMP::threadNum();
+            Float derivatives_re[2*N_coeff];
+            Float derivatives_im[2*N_coeff];
+            uInt idx[2*N_coeff];
+            #pragma omp for
+            for(int k = 0; k<N_freq; k++)
+            {
+                Float *A_data_k = A_data + k*N_coeff;
+                Float *phases_data_k = phases_data + k*N_station;
+                Bool *flags_data_k = flags_data + k*N_station;
+                for(int i = 1; i<N_station; i++) 
+                {
+                    Float phases_data_k_i = phases_data_k[i];
+                    Bool flags_data_k_i = flags_data_k[i];
+                    for(int j = 0; j<i; j++)
+                    {
+                        if (noflags || !(flags_data_k_i || flags_data_k[j]))
+                        {
+                            Float phase_ij_obs = phases_data_k_i - phases_data_k[j];
+                            Float phase_ij_model = 0.0;
+                            
+                            for (int l = 0; l<N_coeff; l++) 
+                            {
+                                Float coeff = A_data_k[l];
+                                phase_ij_model += (sol[i + l*N_station] - sol[j + l*N_station]) * coeff ;
+                            }
+                            Float sin_dphase, cos_dphase;
+#if defined(_LIBCPP_VERSION)
+#define sincosf __sincosf
+#endif
+                            sincosf(phase_ij_obs - phase_ij_model, &sin_dphase, &cos_dphase);
+                            Float residual_re = cos_dphase - 1.0;
+                            Float residual_im = sin_dphase;
+                            Float derivative_re = -sin_dphase;
+                            Float derivative_im = cos_dphase;
+ 
+                            
+                            for (int l = 0; l<N_coeff; l++)
+                            {
+                                Float coeff = A_data_k[l];
+                                Float a = derivative_re * coeff;
+                                Float b = derivative_im * coeff;
+                                derivatives_re[l] = a;
+                                derivatives_re[N_coeff + l] = -a;
+                                derivatives_im[l] = b;
+                                derivatives_im[N_coeff + l] = -b;
+                                idx[l] = i+l*N_station;
+                                idx[N_coeff + l] = j+l*N_station;
+                            }
+                            lnl[threadNum].makeNorm(uInt(2*N_coeff), (uInt*) idx, (Float*) derivatives_re, Float(1.0), residual_re);
+                            lnl[threadNum].makeNorm(uInt(2*N_coeff), (uInt*) idx, (Float*) derivatives_im, Float(1.0), residual_im);
+                        }
+                    }
+                }
+            }
+        }
+        
+        for(int i = 1; i<N_thread; i++)
+        {
+            lnl[0].merge(lnl[i]);
+        }
+        
+        if ((!no_constant_parm) )
+        {
+            for (int i = 0; i<N_parm; i++) 
+            {
+                if (constant_parm_data[i])
+                {
+                    cEq[i] = 1.0;
+                    lnl[0].addConstraint( (Float*) cEq, 0.0);
+                    cEq[i] = 0.0;
+                }
+            }
+        }
+        
+        if (!lnl[0].solveLoop(nr, sol, True)) 
+        {
+            cout << "Error in loop: " << nr << endl;
+            break;
+        }
+        if (lnl[0].isReady())
+        {
+            break;
+        }
+    }
+    
+    Array<Float> solutions(IPosition(2, N_station, N_coeff), sol);
+    ValueHolder result(solutions);
+    return result;
+};
+
+
+} // namespace ExpIon
+} // namespace LOFAR
+
+BOOST_PYTHON_MODULE(_baselinefitting)
+{
+    casacore::pyrap::register_convert_excp();
+    casacore::pyrap::register_convert_basicdata();
+    casacore::pyrap::register_convert_casa_valueholder();
+    casacore::pyrap::register_convert_casa_record();
+    
+    def("fit", LOFAR::ExpIon::fit);
+}

CEP/Calibration/ExpIon/src/baselinefitting.py

+"""Fit basefunctions to phases using all baselines
+
+Application: Clock-TEC separation
+
+
+The phase can be described by a linear model
+
+.. math::
+   phases = A p
+   
+where the columns of matrix A contain the basefunctions and p are the parameters
+
+The same equation holds for the phases of multiple stations
+The dimension of phases is then N_freq x N_station and 
+the dimension of p is N_freq x N_param
+
+The cost function that will be minimized is 
+
+.. math::
+   \\sum_{i=1}^{N_{station}-1}\\sum_{j=0}^{i-1} \\sum_{k=0}^{N_{freq}} \\| \\exp(\imath(\\Delta phase_{ijk} - \\Delta modelphase_{ijk}))  \\|^2
+where
+.. math::
+   \\Delta phase_{ijk} =
+   \\Delta modelphase_{ijk} = 
+
+
+
+"""
+  
+ 
+import _baselinefitting
+ 
+def fit(phases, A, p_0 = None, flags = None, constant_parms = None):
+    """see module description for detailed info"""
+    return _baselinefitting.fit(phases, A, p_0, flags, constant_parms)
\ No newline at end of file

CEP/Calibration/ExpIon/src/calibrate-ion

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Copyright (C) 2007
+# ASTRON (Netherlands Institute for Radio Astronomy)
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite.
+# The LOFAR software suite is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+#
+# $Id$
+
+import sys
+import os
+import re
+import lofar.expion.ionosphere as ionosphere
+import lofar.parameterset
+
+from pylab import *
+
+arguments = sys.argv.__iter__()
+scriptname = os.path.basename(arguments.next())
+
+def print_error(msg):
+   print "%s: %s" % (scriptname, msg)
+
+sky_name = 'sky'
+instrument_name = 'instrument'
+clusterdescfile = "~/CEP.clusterdesc"
+globaldb = 'globaldb'
+
+def display_help_and_exit(dummy = None):
+   print '''usage:
+   %s [options] <parsetfile> <gdsfiles> ''' % scriptname
+   print
+   print '''arguments:
+   <gdsfiles>     gds file(s)'''
+   print
+   print '''options:
+   -h, --help        display this help text
+   --cluster-desc <cluster description file>
+                     Define cluster description file
+                     (default: ~/CEP.clusterdesc)
+   --sky-name <skyname>
+                     Define basename for the sky model parmdb
+                     (default: sky)
+   --instrument-name <instrumentname>
+                     Define basename of the instrument parmdb
+                     (default: instrument)
+   '''
+   exit()
+
+def set_clusterdesc(arguments):
+   global clusterdescfile
+   try:
+      clusterdescfile = arguments.next()
+      if clusterdescfile[0] == '-':
+         raise ValueError
+   except (KeyError, ValueError):
+      print_error( "--cluster-desc should be followed by the name of the cluster description file")
+      exit()
+
+def set_instrument_name(arguments):
+   global instrument_name
+   try:
+      instrument_name = arguments.next()
+      if instrument_name[0] == '-':
+         raise ValueError
+   except (KeyError, ValueError):
+      print_error( "--instrument-name should be followed by the basename of the instrument parmdb")
+      exit()
+
+def set_sky_name(arguments):
+   global sky_name
+   try:
+      sky_name = arguments.next()
+      if sky_name[0] == '-':
+         raise ValueError
+   except (KeyError, ValueError):
+      print_error( "--sky-name should be followed by the basename of the sky parmdb")
+      exit()
+
+def set_globaldb(arguments):
+   global globaldb
+   try:
+      globaldb = arguments.next()
+      if globaldb[0] == '-':
+         raise ValueError
+   except (KeyError, ValueError):
+      print_error( "--globaldb should be followed by the global")
+      exit()
+
+options = { '-h':                display_help_and_exit,
+            '--help':            display_help_and_exit,
+            '--cluster-desc':    set_clusterdesc,
+            '--sky':             set_sky_name,
+            '--instrument-name': set_instrument_name,
+            '--globaldb':              set_globaldb}
+
+while True:
+   try: 
+      argument = arguments.next()
+   except StopIteration:
+      print_error( "No parset file and no gds file(s) specified" )
+      display_help_and_exit()
+   if argument[0] == '-':
+      try:
+         options[argument](arguments)
+      except KeyError:
+         print_error( "Unknown option: " + argument )
+         display_help_and_exit()
+   else:
+      break
+
+parsetfile = argument
+parset = lofar.parameterset.parameterset( parsetfile )
+
+print clusterdescfile
+clusterdescfile = os.path.expanduser( clusterdescfile )
+print clusterdescfile
+
+clusterdesc = lofar.parameterset.parameterset( clusterdescfile )
+   
+gdsfiles = []
+gdsfiles.extend(arguments)
+
+if len(gdsfiles) == 0 :
+   print_error( "No gds file(s) or globaldb specified" )
+   display_help_and_exit()
+
+
+print "parset-file: " + repr(parsetfile)
+print "gds-files: " + repr(gdsfiles)
+print "instrument-name: " + repr(instrument_name)
+print "sky-name: " + repr(sky_name)
+
+stations = parset.getStringVector("ExpIon.Stations", [])
+sources = parset.getStringVector("ExpIon.Sources", [])
+
+d = {'XX' : 0, 'YY' : 1}
+l = parset.getStringVector("ExpIon.Polarizations", ['XX', 'YY'])
+polarizations = [d[key] for key in l]
+
+DirectionalGainEnable = parset.getBool( "ExpIon.DirectionalGain.Enable", False ) 
+GainEnable = parset.getBool( "ExpIon.Gain.Enable", False ) 
+PhasorsEnable = parset.getBool( "ExpIon.Phasors.Enable", False )
+RotationEnable = parset.getBool( "ExpIon.Rotation.Enable", False )
+print "RotationEnable:", RotationEnable
+
+print repr(stations)
+print repr(sources)
+print repr(DirectionalGainEnable)
+
+ion_model = ionosphere.IonosphericModel(gdsfiles, clusterdescfile,
+                                        GainEnable = GainEnable,
+                                        DirectionalGainEnable = DirectionalGainEnable,
+                                        RotationEnable = RotationEnable,
+                                        PhasorsEnable = PhasorsEnable,
+                                        stations = stations,
+                                        sources = sources,
+                                        polarizations = polarizations,
+                                        instrument_name = instrument_name,
+                                        globaldb = globaldb)
+
+def operation_clocktec ( step ):
+   ClockEnable = parset.getBool('.'.join(["ExpIon.Steps", step, "Clock.Enable"]), True )
+   TECEnable = parset.getBool('.'.join(["ExpIon.Steps", step, "TEC.Enable"]), True )
+   print '.'.join(["ExpIon.Steps", step, "Stations"])
+   stations = parset.getStringVector('.'.join(["ExpIon.Steps", step, "Stations"]), [])
+   print "stations: ", stations
+   if ClockEnable or TECEnable :
+      ion_model.ClockTEC( ClockEnable = ClockEnable, TECEnable = TECEnable, stations = stations )
+   
+
+def operation_fitmodel( step ) :
+   height = parset.getFloat('.'.join(["ExpIon.Steps", step, "height"]), 200.0e3 )
+   ion_model.calculate_piercepoints(height = height)
+
+   order = parset.getInt('.'.join(["ExpIon.Steps", step, "order"]), 2 )
+   ion_model.calculate_basevectors( order = order )
+
+   ion_model.fit_model()
+
+def operation_makemovie( step ) :
+   npixels = parset.getInt('.'.join(["ExpIon.Steps", step, "npixels"]), 100 )
+   extent = parset.getFloat('.'.join(["ExpIon.Steps", step, "extent"]), 0 )
+   clim = parset.getFloatVector('.'.join(["ExpIon.Steps", step, "clim"]), [] )
+   if len(clim) == 2:
+      vmin = clim[0]
+      vmax = clim[1]
+   else:
+      vmin = 0
+      vmax = 0
+   ion_model.make_movie(extent = extent, npixels = npixels, vmin = vmin, vmax = vmax)
+
+def operation_store ( step ):
+   ClockTEC_parmdbname = parset.getString('.'.join(["ExpIon.Steps", step, "ClockTEC.parmdb"]), "ClockTEC.parmdb")
+   phases_parmdbname = parset.getString('.'.join(["ExpIon.Steps", step, "Phases.parmdb"]), "ionosphere")
+   ion_model.write_to_parmdb(  )
+   #ion_model.write_phases_to_parmdb( ClockTEC_parmdbname, phases_parmdbname )
+   
+def operation_plot ( step ):
+   print ion_model.TEC.shape
+   figure(1)
+   plot(ion_model.Clock[0,:,:])
+   figure(2)
+   plot(ion_model.TEC[0,:,:])
+   show()
+
+Operations = { "CLOCKTEC": operation_clocktec ,
+               "FITMODEL": operation_fitmodel,
+               "MAKEMOVIE": operation_makemovie,
+               "STORE":  operation_store,
+               "PLOT":  operation_plot }
+
+steps = parset.getStringVector("ExpIon.Steps", [] )
+
+for step in steps:
+   operation = parset.getString( '.'.join( [ "ExpIon.Steps", step, "Operation" ] ) )
+   Operations[ operation ] ( step )
+