diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc
index b0c92c6db629f22973c36faf198955b51775a2b2..7ea02509383487a152f6b0dfe5bf3105f498d072 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc
@@ -73,11 +73,6 @@ namespace LOFAR
// PolarizationsToCheck.resize(NumPolarizations);
// DeterminePolarizationsToCheck(UseOnlyXpolarizations);
- TimeslotData.resize(NumPairs*NumBands);
- for (int i = 0; i < NumPairs*NumBands; i++)
- { TimeslotData[i].resize(NumPolarizations, NumChannels, WindowSize);
- }
-
int index = 0;
for (int i = 0; i < NumAntennae; i++)
{ for(int j = i; j < NumAntennae; j++)
@@ -164,66 +159,104 @@ namespace LOFAR
}
}
- //===============>>> DFTImager::ImageTimeslot <<<===============
- /* This function iterates over baseline and band and uses FlagBaselineBand() to determine
- for each one if it needs to be flagged. It treats the autocorrelations separately,
- to detect entire malfunctioning telescopes. Finally it writes the flags.
+ //===============>>> DFTImager::ImageBaseline <<<===============
+ /*
+ Cube<Complex> TimeslotData(NumBands*NumPairs, NumPolarizations, NumChannels);
+ Matrix<Double> UVWData(NumBands*NumPairs, 3);
+ vector< Cube<float> > Image(NumBands, Cube<float> (Resolution, Resolution, NumChannels));
*/
- void DFTImager::ImageTimeslot()
+ void DFTImager::ImageBaseline(Matrix<Complex>& TimeslotData,
+ Vector<Double>& UVWData,
+ Cube<float>& Image)
{
-/* Cube<bool> Image(NumPolarizations, NumChannels, NumPairs*NumBands, false);
- vector<bool> test(NumPairs*NumBands);
- for (int i = 0; i < NumBands; i++)
+ vector<double> arrayLM(Resolution+1);
+ const double pi = 3.14;
+ const double Freq = 60000000.0; // needs to be parameter
+ const double c = 299792458.0;
+ for (int i=0; i <= Resolution; i++)
+ { arrayLM[i] = -1 + i*2/Resolution;
+ }
+ for (int k = 0; k < NumChannels; k++)
{
- vector< vector<double> > CrossCorrRMS(NumPolarizations, NumPairs - NumAntennae);
- for(int j = 0; j < NumAntennae; j++)
+ Matrix<float> skymap(Resolution, Resolution);
+ for (int l = 0; l < Resolution; l++)
{
- for(int k = j+1; k < NumAntennae; k++)
+ //cout << l << ":";
+ for (int m = 0; m < Resolution; m++)
{
- int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- Matrix<Bool> image = Image.xyPlane(index);
- if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
- { //we skip the non-existent telescopes
- test = ImageBaselineBand(&image,
- &(TimeslotData[index]),
- Position);
+ //cout << m << ";";
+ double temp = (arrayLM[l]*arrayLM[l] + arrayLM[m]*arrayLM[m]);
+ if (temp < 1.0)
+ {
+ double n = std::sqrt(1 - temp);
+ Complex I = TimeslotData(k, 0) + TimeslotData(k, NumPolarizations-1);
+ Complex X = UVWData(0)*arrayLM[l] + UVWData(1)*arrayLM[m]+ UVWData(2) * (n - 1.0);
+ Complex z = Complex(n, 0)*I*exp(Complex(2,0)*Complex(pi, 0)*Complex(0,1)*Complex(Freq, 0)*X/Complex(c, 0));
+ skymap(l, m) += 2 * z.real();
}
}
}
- //this code could be separated in to different functions
- for(int j = 0; j < NumAntennae; j++) //write the data
+ Image.xyPlane(k) = Image.xyPlane(k) + skymap;
+ }
+ }
+ //===============>>> DFTImager::ImageTimeslot <<<===============
+ /*
+ Cube<Complex> TimeslotData(NumBands*NumPairs, NumPolarizations, NumChannels);
+ Matrix<Double> UVWData(NumBands*NumPairs, 3);
+ vector< Cube<float> > Image(NumBands, Cube<float> (Resolution, Resolution, NumChannels));
+ */
+ void DFTImager::ImageTimeslot(Cube<Complex>& TimeslotData,
+ Matrix<Double>& UVWData,
+ vector<Cube<float> >& Image)
+ {
+ for (int k = 0; k < NumBands; k++)
+ {
+ for(int l = 0; l < NumAntennae; l++)
{
- for(int k = j; k < NumAntennae; k++)
+ for(int m = l+1; m < NumAntennae; m++)
{
- int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- Matrix<bool> image = Image.xyPlane(index);
- (*Imagefile).WriteImageData(&image);
+ int index = k*NumPairs + BaselineIndex[pairii(l, m)];
+ if ((BaselineLengths[BaselineIndex[pairii(l, m)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
+ { //we skip the non-existent telescopes
+ if (UVWData(0, index)*UVWData(0, index) + UVWData(1, index)*UVWData(1, index) + UVWData(2, index)*UVWData(2, index) > 0.1)
+ {
+ Matrix<Complex> TD = TimeslotData.xyPlane(index);
+ Vector<Double> UD = UVWData.column(index);
+ ImageBaseline(TD, UD, Image[k]);
+ }
+ }
}
}
- }*/
+ }
}
-
//===============>>> DFTImager::UpdateTimeslotData <<<===============
/* This function reads the visibility data for one timeslot and checks if
a mosaicing mode is active*/
/* The datastructure Timeslotdata is rather complex because it flattens
half-filled antenna x antenna x bands matrix into a vector of NumPairs X bands length. */
- bool DFTImager::UpdateTimeslotData(vector<int>* OldFields,
- vector<int>* OldBands,
- int* TimeCounter,
- Table* TimeslotTable,
- double* Time)
+ bool DFTImager::UpdateTimeslotData(vector<int>& OldFields,
+ vector<int>& OldBands,
+ int* TimeCounter,
+ Table* TimeslotTable,
+ Cube<Complex>& TimeslotData,
+ Matrix<Double>& UVWData,
+ double* Time)
{
int rowcount = (*TimeslotTable).nrow();
+
ROTableVector<Int> antenna1((*TimeslotTable), "ANTENNA1");
ROTableVector<Int> antenna2((*TimeslotTable), "ANTENNA2");
ROTableVector<Int> bandnr ((*TimeslotTable), "DATA_DESC_ID");
ROTableVector<Int> fieldid ((*TimeslotTable), "FIELD_ID");
ROTableVector<Double> time ((*TimeslotTable), "TIME_CENTROID");//for testing purposes
ROArrayColumn<Complex> data ((*TimeslotTable), "DATA");
+ ROArrayColumn<Double> uvw ((*TimeslotTable), "UVW");
+
+ Matrix<Double> uvwdata(3, rowcount);
Cube<Complex> Data(NumPolarizations, NumChannels, rowcount);
data.getColumn(Data); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
+ uvw.getColumn(uvwdata);
(*Time) = time(0);//for testing purposes, might be useful in the future
bool NewField_or_Frequency = false;
@@ -234,19 +267,21 @@ namespace LOFAR
int field = fieldid(i);
int band = bandnr(i);
int index = (band % NumBands) * NumPairs + bi;
- if ((*TimeCounter) > WindowSize - 1)
+ if ((*TimeCounter) > 1)
{
- if (field != (*OldFields)[bi]) //pointing mosaicing
- { NewField_or_Frequency = true;
+ if (field != OldFields[bi]) //pointing mosaicing
+ { //NewField_or_Frequency = true;
+ cout << "OldField:" << OldFields[bi] << " NewField:" << field << endl;
}
- if (band != (*OldBands)[index]) //frequency mosaicing
- { NewField_or_Frequency = true;
+ if (band != OldBands[index]) //frequency mosaicing
+ { //NewField_or_Frequency = true;
+ cout << "OldBand:" << OldBands[bi] << " NewBand:" << band << endl;
}
}
- (*OldFields)[bi] = field;
- (*OldBands)[index] = band;
-
- TimeslotData[index].xyPlane((*TimeCounter) % WindowSize) = Data.xyPlane(i); //TimeslotData is only WindowSize size!
+ OldFields[bi] = field;
+ OldBands[index] = band;
+ TimeslotData.xyPlane(index) = Data.xyPlane(i);
+ UVWData.column(index) = uvwdata.column(i);
}
(*TimeCounter)++; //we want to reset if we detect a gap in DATA_DESC_ID or FIELD. Maybe TIME too???
return NewField_or_Frequency; //assuming everybody is changing field or frequency at the same time!
@@ -255,14 +290,17 @@ namespace LOFAR
//===============>>> DFTImager::FlagDataOrBaselines <<<===============
/* This function iterates over the data per timeslot and uses Flagtimeslot()
to actually flag datapoints (if flagDatapoints), and entire baselines (if flagRMS)*/
- void DFTImager::MakeImage(int Resolution)
+ void DFTImager::MakeImage(int resolution)
{
- TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
- int TimeCounter = 0;
- double Time = 0.0;//for testing purposes
- vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
- vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
- vector< Cube<float> > Image(NumBands);
+ Resolution = resolution;
+ TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
+ int TimeCounter = 0;
+ double Time = 0.0;//for testing purposes
+ vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
+ vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
+ Cube<Complex> TimeslotData(NumPolarizations, NumChannels, NumBands*NumPairs);
+ Matrix<Double> UVWData(3, NumBands*NumPairs);
+ vector< Cube<float> > Image(NumBands, Cube<float> (Resolution, Resolution, NumChannels, 0.0));
int step = NumTimeslots / 10 + 1; //not exact but it'll do
int row = 0;
@@ -270,10 +308,12 @@ namespace LOFAR
while (!timeslot_iter.pastEnd())
{
Table TimeslotTable = timeslot_iter.table();
- bool NewFieldorFreq = UpdateTimeslotData(&OldFields,
- &OldBands,
+ bool NewFieldorFreq = UpdateTimeslotData(OldFields,
+ OldBands,
&TimeCounter,
&TimeslotTable,
+ TimeslotData,
+ UVWData,
&Time);
if (NewFieldorFreq)
{
@@ -282,13 +322,15 @@ namespace LOFAR
}
cout << "Processing: " << MVTime(Time/(24*3600)).string(MVTime::YMD) << endl; //for testing purposes
- ImageTimeslot();
+ ImageTimeslot(TimeslotData,
+ UVWData,
+ Image);
timeslot_iter++;
if (row++ % step == 0) // to tell the user how much % we have processed,
{ cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
}
}
- Imagefile->WriteImage(&Image);
+ Imagefile->WriteImage(Image);
}
//===============>>> DFTImager <<<===============
} //namespace CS1
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h
index 6767aa4b0785cfad92a5178901824045a47edb8f..2c736d2518a5d12cb6abca16cf2a1c42d30d2058 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h
@@ -33,13 +33,9 @@ namespace LOFAR
namespace CS1
{
using casa::Cube;
- using casa::Matrix;
using casa::Complex;
- using casa::Int;
- using casa::Double;
- using casa::Bool;
- using casa::Table;
- using std::list;
+ using casa::Matrix;
+ using casa::Vector;
using std::vector;
typedef pair<int, int> pairii;
@@ -51,7 +47,7 @@ namespace LOFAR
Image_File* Imagefile);
~DFTImager();
- void MakeImage(int Resolution);
+ void MakeImage(int resolution);
protected:
int NumAntennae;
@@ -59,7 +55,6 @@ namespace LOFAR
int NumBands;
int NumChannels;
int NumPolarizations;
- int WindowSize;
int NumTimeslots;
double MinThreshold;
double MaxThreshold;
@@ -67,9 +62,7 @@ namespace LOFAR
vector<double> BaselineLengths;
vector<pairii> PairsIndex;
map<pairii, int> BaselineIndex;
- vector< Cube<Complex> > TimeslotData;
vector<casa::String> AntennaNames;
- Cube< int > Statistics;
int Resolution;
MS_File* MSfile;
@@ -77,12 +70,19 @@ namespace LOFAR
private:
void ComputeBaselineLengths();
- bool UpdateTimeslotData(vector<int>* OldFields,
- vector<int>* OldBands,
+ void ImageBaseline(Matrix<casa::Complex>& TimeslotData,
+ Vector<casa::Double>& UVWData,
+ Cube<float>& Image);
+ void ImageTimeslot(Cube<casa::Complex>& TimeslotData,
+ Matrix<casa::Double>& UVWData,
+ vector<Cube<float> >& Image);
+ bool UpdateTimeslotData(vector<int>& OldFields,
+ vector<int>& OldBands,
int* TimeCounter,
- Table* TimeslotTable,
+ casa::Table* TimeslotTable,
+ casa::Cube<casa::Complex>& TimeslotData,
+ casa::Matrix<casa::Double>& UVWData,
double* Time);
- void ImageTimeslot(void);
}; // DFTImager
}; // namespace CS1
}; // namespace LOFAR
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc
index b20362ffdabd35438d6f0db9f9a6212f0b9135ae..4cf485291caa8783d4d70e956d5575f5fe229700 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc
@@ -77,10 +77,10 @@ namespace LOFAR
//===============>>> Image_file::~WriteImage <<<===============
- void Image_File::WriteImage(std::vector< casa::Cube<float> >* image)
+ void Image_File::WriteImage(std::vector< casa::Cube<float> > image)
{
- casa::Cube<float> Bal = (*image)[0];
- init(Bal.shape(), image->size());
+ casa::Cube<float> Bal = image[0];
+ init(Bal.shape(), image.size());
Image->putSlice(Bal, casa::IPosition(3, 0, 0, 0), casa::IPosition(3, 0, 0, 0));
}
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h
index 876236ca933f935fda555ab4b69146590be033ff..976965386a05fa5a058129842d1e96b2d0d7ce74 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h
@@ -51,7 +51,7 @@ namespace LOFAR
public:
Image_File(const std::string& imagename);
~Image_File();
- void WriteImage(vector< casa::Cube<float> >* image);
+ void WriteImage(vector< casa::Cube<float> > image);
protected:
std::string ImageName;
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h
index c1a530418cc357cd5b72fb273956125745b2bd1c..02fd216ddf6c831a7a7decf0cf7bb4b44d53a017 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h
@@ -25,6 +25,7 @@
#include <string>
#include <tables/Tables.h>
#include <tables/Tables/TableIter.h>
+#include <casa/Arrays.h>
namespace LOFAR
{