diff --git a/.gitattributes b/.gitattributes
index 8a8400c3b6ae018d88eea24f4fd35dbca7df1703..de64c6841c537dfe79d8265b9933c370e7a87f0d 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -210,6 +210,7 @@ CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/combineflagresults.h -text
CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/cutareaaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/fitsimageset.h -text
CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachbaselineaction.h -text
+CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachcomplexcomponentaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachmsaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachpolarisationblock.h -text
CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/frequencyselectionaction.h -text
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/CMakeLists.txt b/CEP/DP3/AOFlagger/include/AOFlagger/CMakeLists.txt
index 689ea31709b9451753b34219cfbe1cf039edeb37..f816222fa715da6b3992b3484cc64b0d56f90686 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/CMakeLists.txt
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/CMakeLists.txt
@@ -131,6 +131,7 @@ install(FILES
rfi/strategy/combineflagresults.h
rfi/strategy/fitsimageset.h
rfi/strategy/foreachbaselineaction.h
+ rfi/strategy/foreachcomplexcomponentaction.h
rfi/strategy/foreachmsaction.h
rfi/strategy/foreachpolarisationblock.h
rfi/strategy/frequencyselectionaction.h
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h
index b1d4535f1e5c65ca2b5718358555e71a0abea3d4..5260ec3caa4bb54043e801a4e152b88dd641772b 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h
@@ -31,6 +31,7 @@
template<typename T>
struct OutputReceiver
{
+ void SetY(size_t) { }
};
template<>
struct OutputReceiver<class UVImager>
@@ -98,6 +99,7 @@ class Model {
static num_t GetWPosition(num_t delayDirectionDec, num_t delayDirectionRA, num_t frequency, num_t earthLattitudeAngle, num_t dx, num_t dy);
void loadUrsaMajor();
+ void loadUrsaMajorDistortingSource();
private:
std::vector<PointSource *> _sources;
double _noiseSigma, _sourceSigma;
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/observatorium.h b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/observatorium.h
index 6e019a21233277c981ea5b7559a726e1d02251eb..576236a4077aaed36d1a51580b646a577c999f24 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/observatorium.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/observatorium.h
@@ -14,11 +14,13 @@ class Observatorium
}
size_t AntennaCount() const { return _antennae.size(); }
const AntennaInfo &GetAntenna(size_t index) const { return _antennae[index]; }
+
void SetChannelWidthHz(double channelWidthHz)
{
_channelWidthHz = channelWidthHz;
}
double ChannelWidthHz() const { return _channelWidthHz; }
+
const struct BandInfo &BandInfo() const
{
return _bandInfo;
@@ -42,6 +44,7 @@ struct WSRTObservatorium : public Observatorium
WSRTn(i, antennas[i]);
AddAntenna(antennas[i]);
}
+ initBand();
}
WSRTObservatorium(size_t antenna1, size_t antenna2)
{
@@ -52,6 +55,7 @@ struct WSRTObservatorium : public Observatorium
WSRTn(antenna2, antennas[1]);
AddAntenna(antennas[0]);
AddAntenna(antennas[1]);
+ initBand();
}
private:
@@ -167,16 +171,16 @@ struct WSRTObservatorium : public Observatorium
}
void initBand()
{
- SetChannelWidthHz(10000.0);
GetBandInfo().windowIndex = 0;
- GetBandInfo().channelCount = 256;
- for(size_t i=0;i<256;++i)
+ GetBandInfo().channelCount = 16;
+ SetChannelWidthHz(10000.0 * 256.0 / GetBandInfo().channelCount);
+ for(size_t i=0;i<GetBandInfo().channelCount;++i)
{
ChannelInfo channel;
channel.frequencyIndex = i;
channel.channelWidthHz = ChannelWidthHz();
channel.frequencyHz = 147000000.0 + ChannelWidthHz() * i;
- GetBandInfo().channels[i] = channel;
+ GetBandInfo().channels.push_back(channel);
}
}
};
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachcomplexcomponentaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachcomplexcomponentaction.h
new file mode 100644
index 0000000000000000000000000000000000000000..7b98905c25a8d09ef62171aae2568ac975c8c2f9
--- /dev/null
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/foreachcomplexcomponentaction.h
@@ -0,0 +1,252 @@
+/***************************************************************************
+ * Copyright (C) 2008 by A.R. Offringa *
+ * offringa@astro.rug.nl *
+ * *
+ * This program 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 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program 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 this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#ifndef FOR_EACH_COMPLEX_COMPONENT_ACTION_H
+#define FOR_EACH_COMPLEX_COMPONENT_ACTION_H
+
+#include "artifactset.h"
+#include "actionblock.h"
+
+namespace rfiStrategy {
+
+ class ForEachComplexComponentAction : public ActionBlock
+ {
+ public:
+ ForEachComplexComponentAction() : ActionBlock(), _onAmplitude(false), _onPhase(false), _onReal(true), _onImaginary(true)
+ {
+ }
+ virtual std::string Description()
+ {
+ return "For each complex component";
+ }
+ virtual ActionType Type() const { return AdapterType; }
+ virtual void Perform(ArtifactSet &artifacts, class ProgressListener &listener)
+ {
+ size_t taskCount = 0;
+ if(_onAmplitude) ++taskCount;
+ if(_onPhase) ++taskCount;
+ if(_onReal) ++taskCount;
+ if(_onImaginary) ++taskCount;
+
+ size_t taskIndex = 0;
+
+ if(_onAmplitude) {
+ listener.OnStartTask(taskIndex, taskCount, "On amplitude");
+ performOnAmplitude(artifacts, listener);
+ listener.OnEndTask();
+ ++taskIndex;
+ }
+ if(_onPhase) {
+ listener.OnStartTask(taskIndex, taskCount, "On phase");
+ performOnPhaseRepresentation(artifacts, listener, TimeFrequencyData::PhasePart);
+ listener.OnEndTask();
+ ++taskIndex;
+ }
+ if(_onReal) {
+ listener.OnStartTask(taskIndex, taskCount, "On real");
+ performOnPhaseRepresentation(artifacts, listener, TimeFrequencyData::RealPart);
+ listener.OnEndTask();
+ ++taskIndex;
+ }
+ if(_onImaginary) {
+ listener.OnStartTask(taskIndex, taskCount, "On imaginary");
+ performOnPhaseRepresentation(artifacts, listener, TimeFrequencyData::ImaginaryPart);
+ listener.OnEndTask();
+ ++taskIndex;
+ }
+ }
+ void SetRestoreFromAmplitude(bool restoreFromAmplitude)
+ {
+ _restoreFromAmplitude = restoreFromAmplitude;
+ }
+ bool RestoreFromAmplitude() const
+ {
+ return _restoreFromAmplitude;
+ }
+ void SetOnAmplitude(bool onAmplitude)
+ {
+ _onAmplitude = onAmplitude;
+ }
+ bool OnAmplitude() const
+ {
+ return _onAmplitude;
+ }
+ void SetOnPhase(bool onPhase)
+ {
+ _onPhase = onPhase;
+ }
+ bool OnPhase() const
+ {
+ return _onPhase;
+ }
+ void SetOnReal(bool onReal)
+ {
+ _onReal = onReal;
+ }
+ bool OnReal() const
+ {
+ return _onReal;
+ }
+ void SetOnImaginary(bool onImaginary)
+ {
+ _onImaginary = onImaginary;
+ }
+ bool OnImaginary() const
+ {
+ return _onImaginary;
+ }
+ private:
+ void performOnAmplitude(ArtifactSet &artifacts, class ProgressListener &listener)
+ {
+ enum TimeFrequencyData::PhaseRepresentation contaminatedPhase =
+ artifacts.ContaminatedData().PhaseRepresentation();
+ enum TimeFrequencyData::PhaseRepresentation revisedPhase =
+ artifacts.RevisedData().PhaseRepresentation();
+ enum TimeFrequencyData::PhaseRepresentation originalPhase =
+ artifacts.OriginalData().PhaseRepresentation();
+
+ if(contaminatedPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newContaminatedData =
+ artifacts.ContaminatedData().CreateTFData(TimeFrequencyData::AmplitudePart);
+ artifacts.SetContaminatedData(*newContaminatedData);
+ delete newContaminatedData;
+ }
+ if(revisedPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newRevisedData =
+ artifacts.RevisedData().CreateTFData(TimeFrequencyData::AmplitudePart);
+ artifacts.SetRevisedData(*newRevisedData);
+ delete newRevisedData;
+ }
+ if(originalPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newOriginalData =
+ artifacts.OriginalData().CreateTFData(TimeFrequencyData::AmplitudePart);
+ artifacts.SetOriginalData(*newOriginalData);
+ delete newOriginalData;
+ }
+
+ ActionBlock::Perform(artifacts, listener);
+
+ if(_restoreFromAmplitude)
+ {
+ if(contaminatedPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newContaminatedData =
+ TimeFrequencyData::CreateTFDataFromComplexCombination(artifacts.ContaminatedData(), artifacts.ContaminatedData());
+ newContaminatedData->MultiplyImages(1.0L/M_SQRT2);
+ newContaminatedData->SetMask(artifacts.ContaminatedData());
+ artifacts.SetContaminatedData(*newContaminatedData);
+ delete newContaminatedData;
+ }
+ if(revisedPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newRevisedData =
+ TimeFrequencyData::CreateTFDataFromComplexCombination(artifacts.RevisedData(), artifacts.RevisedData());
+ newRevisedData->MultiplyImages(1.0L/M_SQRT2);
+ newRevisedData->SetMask(artifacts.RevisedData());
+ artifacts.SetRevisedData(*newRevisedData);
+ delete newRevisedData;
+ }
+ if(originalPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newOriginalData =
+ TimeFrequencyData::CreateTFDataFromComplexCombination(artifacts.OriginalData(), artifacts.OriginalData());
+ newOriginalData->MultiplyImages(1.0L/M_SQRT2);
+ newOriginalData->SetMask(artifacts.OriginalData());
+ artifacts.SetOriginalData(*newOriginalData);
+ delete newOriginalData;
+ }
+ }
+ }
+
+ void performOnPhaseRepresentation(ArtifactSet &artifacts, class ProgressListener &listener, enum TimeFrequencyData::PhaseRepresentation phaseRepresentation)
+ {
+ enum TimeFrequencyData::PhaseRepresentation
+ contaminatedPhase = artifacts.ContaminatedData().PhaseRepresentation(),
+ revisedPhase = artifacts.RevisedData().PhaseRepresentation(),
+ originalPhase = artifacts.OriginalData().PhaseRepresentation();
+
+ TimeFrequencyData
+ prevContaminated = artifacts.ContaminatedData(),
+ prevRevised = artifacts.RevisedData(),
+ prevOriginal = artifacts.OriginalData();
+
+ if(contaminatedPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newContaminatedData =
+ artifacts.ContaminatedData().CreateTFData(phaseRepresentation);
+ artifacts.SetContaminatedData(*newContaminatedData);
+ delete newContaminatedData;
+ }
+ if(revisedPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newRevisedData =
+ artifacts.RevisedData().CreateTFData(phaseRepresentation);
+ artifacts.SetRevisedData(*newRevisedData);
+ delete newRevisedData;
+ }
+ if(originalPhase == TimeFrequencyData::ComplexRepresentation)
+ {
+ TimeFrequencyData *newOriginalData =
+ artifacts.OriginalData().CreateTFData(phaseRepresentation);
+ artifacts.SetOriginalData(*newOriginalData);
+ delete newOriginalData;
+ }
+
+ ActionBlock::Perform(artifacts, listener);
+
+ if(phaseRepresentation != TimeFrequencyData::PhasePart)
+ {
+ if(contaminatedPhase == TimeFrequencyData::ComplexRepresentation)
+ setPart(artifacts.ContaminatedData(), prevContaminated);
+ if(revisedPhase == TimeFrequencyData::ComplexRepresentation)
+ setPart(artifacts.RevisedData(), prevRevised);
+ if(originalPhase == TimeFrequencyData::ComplexRepresentation)
+ setPart(artifacts.OriginalData(), prevOriginal);
+ }
+ }
+ void setPart(TimeFrequencyData &changedData, TimeFrequencyData &prevData)
+ {
+ TimeFrequencyData *newData, *otherPart;
+ switch(changedData.PhaseRepresentation())
+ {
+ default:
+ case TimeFrequencyData::RealPart:
+ otherPart = prevData.CreateTFData(TimeFrequencyData::ImaginaryPart);
+ newData = TimeFrequencyData::CreateTFDataFromComplexCombination(changedData, *otherPart);
+ break;
+ case TimeFrequencyData::ImaginaryPart:
+ otherPart = prevData.CreateTFData(TimeFrequencyData::RealPart);
+ newData = TimeFrequencyData::CreateTFDataFromComplexCombination(*otherPart, changedData);
+ break;
+ }
+ changedData = *newData;
+ delete newData;
+ delete otherPart;
+ }
+
+ bool _restoreFromAmplitude;
+ bool _onAmplitude, _onPhase, _onReal, _onImaginary;
+ };
+
+} // namespace
+
+#endif // FOR_EACH_COMPLEX_COMPONENT_ACTION_H
diff --git a/CEP/DP3/AOFlagger/src/gui/mswindow.cpp b/CEP/DP3/AOFlagger/src/gui/mswindow.cpp
index 1fe8f1377a2f7aea4887797327dffc740f3d8305..9a32ffd69e0d41e8adf42c0a42240fa1b5d599c2 100644
--- a/CEP/DP3/AOFlagger/src/gui/mswindow.cpp
+++ b/CEP/DP3/AOFlagger/src/gui/mswindow.cpp
@@ -1246,7 +1246,7 @@ void MSWindow::onShowImagePlane()
void MSWindow::onAddToImagePlane()
{
try {
- if(_timeFrequencyWidget.GetMetaData() != 0)
+ if(_timeFrequencyWidget.GetMetaData() != 0 && _timeFrequencyWidget.GetMetaData()->HasUVW())
_imagePlaneWindow->AddData(GetActiveData(), _timeFrequencyWidget.GetMetaData());
else if(_spatialMetaData != 0)
_imagePlaneWindow->AddData(GetActiveData(), _spatialMetaData);
@@ -1372,6 +1372,7 @@ void MSWindow::onSimulateCorrelation()
{
Model model;
model.loadUrsaMajor();
+ model.loadUrsaMajorDistortingSource();
WSRTObservatorium wsrtObservatorium;
model.SimulateObservation(*_imagePlaneWindow->GetImager(), wsrtObservatorium, -M_PIn-0.05, 0.05, 147000000.0);
@@ -1382,6 +1383,7 @@ void MSWindow::onSimulateDoubledBaselineCorrelation()
{
Model model;
model.loadUrsaMajor();
+ model.loadUrsaMajorDistortingSource();
WSRTObservatorium wsrtObservatorium(0,5);
std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, -M_PIn-0.05, 0.05, 147000000.0);
diff --git a/CEP/DP3/AOFlagger/src/imaging/model.cpp b/CEP/DP3/AOFlagger/src/imaging/model.cpp
index 743c98c753f9eab8a5e725a5331188c35fa041a6..a3549ef6df7a4ff73e55073bf1c62eaee7224afb 100644
--- a/CEP/DP3/AOFlagger/src/imaging/model.cpp
+++ b/CEP/DP3/AOFlagger/src/imaging/model.cpp
@@ -38,11 +38,11 @@ Model::~Model()
template<typename T>
void Model::SimulateObservation(struct OutputReceiver<T> &receiver, Observatorium &observatorium, num_t delayDirectionDEC, num_t delayDirectionRA, num_t frequency)
{
- size_t frequencySteps = 1;
+ size_t channelCount = observatorium.BandInfo().channelCount;
- for(size_t f=0;f<frequencySteps;++f)
+ for(size_t f=0;f<channelCount;++f)
{
- double channelFrequency = frequency + observatorium.ChannelWidthHz() * f * 256/frequencySteps;
+ double channelFrequency = frequency + observatorium.ChannelWidthHz() * f;
receiver.SetY(f);
for(size_t i=0;i<observatorium.AntennaCount();++i)
{
@@ -67,9 +67,10 @@ template void Model::SimulateObservation(struct OutputReceiver<TimeFrequencyData
std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> Model::SimulateObservation(class Observatorium &observatorium, num_t delayDirectionDEC, num_t delayDirectionRA, num_t frequency)
{
+ size_t channelCount = observatorium.BandInfo().channelCount;
OutputReceiver<TimeFrequencyData> tfOutputter;
- tfOutputter._real = Image2D::CreateZeroImagePtr(12*60*60/10, 1);
- tfOutputter._imaginary = Image2D::CreateZeroImagePtr(12*60*60/10, 1);
+ tfOutputter._real = Image2D::CreateZeroImagePtr(12*60*60/10, channelCount);
+ tfOutputter._imaginary = Image2D::CreateZeroImagePtr(12*60*60/10, channelCount);
SimulateObservation(tfOutputter, observatorium, delayDirectionDEC, delayDirectionRA, frequency);
TimeFrequencyMetaData *metaData = new TimeFrequencyMetaData();
@@ -77,7 +78,22 @@ std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> Model::SimulateObservatio
metaData->SetAntenna2(observatorium.GetAntenna(1));
metaData->SetBand(observatorium.BandInfo());
std::vector<double> times;
- for(num_t t=0.0;t<12*60*60;t+=10) times.push_back(t);
+ std::vector<UVW> uvws;
+ num_t wavelength = 1.0L / frequency;
+ double
+ dx = metaData->Antenna1().position.x - metaData->Antenna2().position.x,
+ dy = metaData->Antenna1().position.y - metaData->Antenna2().position.y,
+ dz = metaData->Antenna1().position.z - metaData->Antenna2().position.z;
+ for(num_t t=0.0;t<12*60*60;t+=10)
+ {
+ times.push_back(t);
+ num_t earthLattitudeApprox = t*(M_PIn/12.0/60.0/60.0);
+ UVW uvw;
+ GetUVPosition(uvw.u, uvw.v, earthLattitudeApprox, delayDirectionDEC, delayDirectionRA, dx, dy, dz, wavelength);
+ uvw.w = 0.0;
+ uvws.push_back(uvw);
+ }
+ metaData->SetUVW(uvws);
metaData->SetObservationTimes(times);
return std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr>(TimeFrequencyData(StokesIPolarisation, tfOutputter._real, tfOutputter._imaginary), TimeFrequencyMetaDataPtr(metaData));
@@ -247,3 +263,15 @@ void Model::loadUrsaMajor()
AddSource(cd + s*rs*2, cr + s*-131, 5.0/8.0 + fluxoffset); // Zeta
AddSource(cd + s*rs*-36, cr + s*-192, 7.0/8.0 + fluxoffset); // Alkaid
}
+
+void Model::loadUrsaMajorDistortingSource()
+{
+ double
+ s = 0.00005, //scale
+ rs = 8.0; // stretch in dec
+ double cd = -M_PIn - 0.05;
+ double cr = 0.05;
+ double fluxoffset = 0.0;
+
+ AddSource(cd + s*rs*160, cr + s*300, 4.0 + fluxoffset); // Dubhe
+}
diff --git a/CEP/DP3/AOFlagger/src/rfi/strategy/actionfactory.cpp b/CEP/DP3/AOFlagger/src/rfi/strategy/actionfactory.cpp
index c60585c5c53d0cb8dfcc126f754ee56dbb0ac7e8..b21c6216d47369e00297df62ea76a4077898fea0 100644
--- a/CEP/DP3/AOFlagger/src/rfi/strategy/actionfactory.cpp
+++ b/CEP/DP3/AOFlagger/src/rfi/strategy/actionfactory.cpp
@@ -28,6 +28,7 @@
#include <AOFlagger/rfi/strategy/cutareaaction.h>
#include <AOFlagger/rfi/strategy/frequencyselectionaction.h>
#include <AOFlagger/rfi/strategy/foreachbaselineaction.h>
+#include <AOFlagger/rfi/strategy/foreachcomplexcomponentaction.h>
#include <AOFlagger/rfi/strategy/foreachpolarisationblock.h>
#include <AOFlagger/rfi/strategy/foreachmsaction.h>
#include <AOFlagger/rfi/strategy/fringestopaction.h>
@@ -56,6 +57,7 @@ const std::vector<std::string> ActionFactory::GetActionList()
list.push_back("Combine flag results");
list.push_back("Cut area");
list.push_back("For each baseline");
+ list.push_back("For each complex component");
list.push_back("For each polarisation");
list.push_back("For each measurement set");
list.push_back("Frequency selection");
@@ -91,6 +93,8 @@ Action *ActionFactory::CreateAction(const std::string &action)
return new CutAreaAction();
else if(action == "For each baseline")
return new ForEachBaselineAction();
+ else if(action == "For each complex component")
+ return new ForEachComplexComponentAction();
else if(action == "For each measurement set")
return new ForEachMSAction();
else if(action == "For each polarisation")