diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h index 5bdf572f12b831345e277ec16316abd66abd5740..bc48dce31941190da74c1a7dae2709a93b2af31f 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h @@ -192,6 +192,8 @@ class MSWindow : public Gtk::Window { void onUnrollPhaseButtonPressed(); void showError(const std::string &description); void onSimulateCorrelation(); + + void getSetData(double &ra, double &dec, double &factor); void onSimulateSourceSetA(); void onSimulateSourceSetB(); void onSimulateSourceSetC(); @@ -215,7 +217,8 @@ class MSWindow : public Gtk::Window { Glib::RefPtr<Gtk::RadioAction> _mapLogButton, _mapBWButton, _mapColorButton, _rangeFullButton, _rangeWinsorizedButton, _rangeSpecifiedButton, - _gaussianTestSetsButton, _rayleighTestSetsButton; + _gaussianTestSetsButton, _rayleighTestSetsButton, + _ncpSetButton, _b1834SetButton; //std::vector<Gtk::Window*> _subWindows; class ImagePlaneWindow *_imagePlaneWindow; Gtk::Window diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/frequencyconvolutionframe.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/frequencyconvolutionframe.h index a043819b6919463e2ad40942d92ee4a8abdd044b..2169db2167a255b204cf1cf13b5f3922ff1ffc36 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/frequencyconvolutionframe.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/frequencyconvolutionframe.h @@ -36,21 +36,32 @@ class FrequencyConvolutionFrame : public Gtk::Frame { FrequencyConvolutionFrame(rfiStrategy::FrequencyConvolutionAction &action, EditStrategyWindow &editStrategyWindow) : Gtk::Frame("Frequency convolution"), _editStrategyWindow(editStrategyWindow), _action(action), - _sizeXLabel("Size x:"), - _sizeYLabel("Size y"), - _sizeXScale(1, 1024, 1), - _sizeYScale(1, 1024, 1), + _rectangularKernelButton("Rectangular kernel"), + _sincKernelButton("Sinc kernel"), + _convolutionSizeLabel("Convolution size:"), + _convolutionSizeScale(1, 1024, 1), + _inSamplesButton("Size in samples"), _applyButton(Gtk::Stock::APPLY) { - _box.pack_start(_sizeXLabel); + Gtk::RadioButton::Group kernelGroup; + + _rectangularKernelButton.set_group(kernelGroup); + _box.pack_start(_rectangularKernelButton); + _sincKernelButton.set_group(kernelGroup); + _box.pack_start(_sincKernelButton); + + if(_action.KernelKind() == rfiStrategy::FrequencyConvolutionAction::RectangleKernel) + _rectangularKernelButton.set_active(true); + else + _sincKernelButton.set_active(true); - _box.pack_start(_sizeXScale); - _sizeXScale.set_value(_action.SizeX()); + _box.pack_start(_convolutionSizeLabel); + + _box.pack_start(_convolutionSizeScale); + _convolutionSizeScale.set_value(_action.ConvolutionSize()); - _box.pack_start(_sizeYLabel); - - _box.pack_start(_sizeYScale); - _sizeYScale.set_value(_action.SizeY()); + _box.pack_start(_inSamplesButton); + _inSamplesButton.set_active(_action.InSamples()); _buttonBox.pack_start(_applyButton); _applyButton.signal_clicked().connect(sigc::mem_fun(*this, &FrequencyConvolutionFrame::onApplyClicked)); @@ -66,14 +77,20 @@ class FrequencyConvolutionFrame : public Gtk::Frame { Gtk::VBox _box; Gtk::HButtonBox _buttonBox; - Gtk::Label _sizeXLabel, _sizeYLabel; - Gtk::HScale _sizeXScale, _sizeYScale; + Gtk::RadioButton _rectangularKernelButton, _sincKernelButton; + Gtk::Label _convolutionSizeLabel; + Gtk::HScale _convolutionSizeScale; + Gtk::CheckButton _inSamplesButton; Gtk::Button _applyButton; void onApplyClicked() { - _action.SetSizeX((unsigned) _sizeXScale.get_value()); - _action.SetSizeY((unsigned)_sizeYScale.get_value()); + if(_rectangularKernelButton.get_active()) + _action.SetKernelKind(rfiStrategy::FrequencyConvolutionAction::RectangleKernel); + else + _action.SetKernelKind(rfiStrategy::FrequencyConvolutionAction::SincKernel); + _action.SetConvolutionSize(_convolutionSizeScale.get_value()); + _action.SetInSamples(_inSamplesButton.get_active()); _editStrategyWindow.UpdateAction(&_action); } }; diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/fringestoppingframe.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/fringestoppingframe.h index 51f9b3df52b0d1b1eaedeff54d2988ebdc7a22c8..a1a7994ba5816b58fe1c648eb3db011f08f77be0 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/fringestoppingframe.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/gui/strategyframes/fringestoppingframe.h @@ -20,6 +20,8 @@ #ifndef FRINGESTOPPINGFRAME_H #define FRINGESTOPPINGFRAME_H +#include <sstream> + #include <gtkmm/box.h> #include <gtkmm/button.h> #include <gtkmm/buttonbox.h> @@ -45,46 +47,52 @@ class FringeStoppingFrame : public Gtk::Frame { _maxWindowSizeScale(0, 2048, 16), _fitChannelsIndividuallyButton("Fit channels individually"), _onlyFringeStopButton("No fit, only fringe stop"), + _raLabel("Right ascension:"), + _decLabel("Declination:"), _applyButton(Gtk::Stock::APPLY) { _box.pack_start(_fringesToConsiderLabel); - _fringesToConsiderLabel.show(); _box.pack_start(_fringesToConsiderScale); _fringesToConsiderScale.set_value(_action.FringesToConsider()); - _fringesToConsiderScale.show(); _box.pack_start(_minWindowSizeLabel); - _minWindowSizeLabel.show(); _box.pack_start(_minWindowSizeScale); _minWindowSizeScale.set_value(_action.MinWindowSize()); - _minWindowSizeScale.show(); _box.pack_start(_maxWindowSizeLabel); - _maxWindowSizeLabel.show(); _box.pack_start(_maxWindowSizeScale); _maxWindowSizeScale.set_value(_action.MaxWindowSize()); - _maxWindowSizeScale.show(); _box.pack_start(_fitChannelsIndividuallyButton); _fitChannelsIndividuallyButton.set_active(_action.FitChannelsIndividually()); - _fitChannelsIndividuallyButton.show(); _box.pack_start(_onlyFringeStopButton); _onlyFringeStopButton.set_active(_action.OnlyFringeStop()); - _onlyFringeStopButton.show(); + + _box.pack_start(_raLabel); + + std::ostringstream raStr; + raStr << _action.NewPhaseCentreRA(); + _raEntry.set_text(raStr.str()); + _box.pack_start(_raEntry); + + _box.pack_start(_decLabel); + + std::ostringstream decStr; + decStr << _action.NewPhaseCentreDec(); + _decEntry.set_text(decStr.str()); + _box.pack_start(_decEntry); _buttonBox.pack_start(_applyButton); _applyButton.signal_clicked().connect(sigc::mem_fun(*this, &FringeStoppingFrame::onApplyClicked)); - _applyButton.show(); _box.pack_start(_buttonBox); - _buttonBox.show(); add(_box); - _box.show(); + _box.show_all(); } private: EditStrategyWindow &_editStrategyWindow; @@ -100,6 +108,9 @@ class FringeStoppingFrame : public Gtk::Frame { Gtk::HScale _maxWindowSizeScale; Gtk::CheckButton _fitChannelsIndividuallyButton; Gtk::CheckButton _onlyFringeStopButton; + Gtk::Label _raLabel; + Gtk::Label _decLabel; + Gtk::Entry _raEntry, _decEntry; Gtk::Button _applyButton; void onApplyClicked() @@ -109,6 +120,11 @@ class FringeStoppingFrame : public Gtk::Frame { _action.SetMaxWindowSize((size_t) _maxWindowSizeScale.get_value()); _action.SetFitChannelsIndividually(_fitChannelsIndividuallyButton.get_active()); _action.SetOnlyFringeStop(_onlyFringeStopButton.get_active()); + std::string + raStr = _raEntry.get_text(), + decStr = _decEntry.get_text(); + _action.SetNewPhaseCentreRA(atof(raStr.c_str())); + _action.SetNewPhaseCentreDec(atof(decStr.c_str())); _editStrategyWindow.UpdateAction(&_action); } }; diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h index 4d4dfeb683f534f7f468651c4ed26b2ad9bf2058..1fef4ffa0abe13af163870efed9396493c8160b5 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h @@ -131,11 +131,14 @@ class Model { void SimulateObservation(struct OutputReceiver<T> &receiver, class Observatorium &observatorium, num_t delayDirectionDEC, num_t delayDirectionRA, num_t frequency); static void GetUVPosition(num_t &u, num_t &v, num_t earthLattitudeAngle, num_t delayDirectionDEC, num_t delayDirectionRA, num_t dx, num_t dy, num_t dz, num_t waveLength); - static num_t GetWPosition(num_t delayDirectionDec, num_t delayDirectionRA, num_t frequency, num_t earthLattitudeAngle, num_t dx, num_t dy); + static num_t GetWPosition(num_t delayDirectionDec, num_t delayDirectionRA, num_t frequency, num_t earthLattitudeAngle, num_t dx, num_t dy) + { + return UVImager::GetWPosition(delayDirectionDec, delayDirectionRA, frequency, earthLattitudeAngle, dx, dy); + } - void loadUrsaMajor(); - void loadUrsaMajorDistortingSource(); - void loadUrsaMajorDistortingVariableSource(bool weak=false, bool slightlyMiss=false); + void loadUrsaMajor(double ra, double dec, double factor); + void loadUrsaMajorDistortingSource(double ra, double dec, double factor); + void loadUrsaMajorDistortingVariableSource(double ra, double dec, double factor, bool weak=false, bool slightlyMiss=false); private: std::vector<Source *> _sources; double _noiseSigma, _sourceSigma; diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/uvimager.h b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/uvimager.h index 6f48a97e799a8abbf52fced5cd88f3ebdbb0f9e5..f2f47e2dc8e4fe9618ae03ac411e5a03761495c7 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/uvimager.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/uvimager.h @@ -70,6 +70,29 @@ class UVImager { static num_t GetFringeStopFrequency(size_t time, const Baseline &baseline, num_t delayDirectionRA, num_t delayDirectionDec, num_t frequency, TimeFrequencyMetaDataCPtr metaData); //static double GetFringeCount(long double timeStart, long double timeEnd, const Baseline &baseline, long double delayDirectionRA, long double delayDirectionDec, long double frequency); static num_t GetFringeCount(size_t timeIndexStart, size_t timeIndexEnd, unsigned channelIndex, TimeFrequencyMetaDataCPtr metaData); + + static numl_t GetWPosition(numl_t delayDirectionDec, numl_t delayDirectionRA, numl_t frequency, numl_t earthLattitudeAngle, numl_t dx, numl_t dy) + { + numl_t wavelength = 299792458.0L / frequency; + numl_t raSinEnd = sinn(-delayDirectionRA - earthLattitudeAngle); + numl_t raCosEnd = cosn(-delayDirectionRA - earthLattitudeAngle); + numl_t decCos = cosn(delayDirectionDec); + // term "+ dz * decCos" is eliminated because of subtraction + num_t wPosition = + (dx*raCosEnd - dy*raSinEnd) * (-decCos) / wavelength; + return wPosition; + } + + static numl_t TimeToEarthLattitude(unsigned x, TimeFrequencyMetaDataCPtr metaData) + { + return TimeToEarthLattitude(metaData->ObservationTimes()[x]); + } + + static numl_t TimeToEarthLattitude(double time) + { + return time*M_PInl/(12.0*60.0*60.0); + } + void Empty(); void PerformFFT(); bool HasUV() const { return _uvReal != 0; } diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h b/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h index a614f167699177260c6d252f042fcca128ed5c0a..1d8efa09b4d83a4b305f243b48e7d1d8ef8d0c2d 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h @@ -123,13 +123,13 @@ struct Baseline { Baseline(EarthPosition _antenna1, EarthPosition _antenna2) : antenna1(_antenna1), antenna2(_antenna2) { } - num_t Distance() { + num_t Distance() const { num_t dx = antenna1.x-antenna2.x; num_t dy = antenna1.y-antenna2.y; num_t dz = antenna1.z-antenna2.z; return sqrtn(dx*dx+dy*dy+dz*dz); } - num_t Angle() { + num_t Angle() const { num_t dz = antenna1.z-antenna2.z; // baseline is either orthogonal to the earths axis, or // the length of the baseline is zero. @@ -140,6 +140,9 @@ struct Baseline { num_t length = sqrtn(dx*dx + dy*dy + 1.0); return acosn(1.0/length); } + num_t DeltaX() const { return antenna2.x-antenna1.x; } + num_t DeltaY() const { return antenna2.y-antenna1.y; } + num_t DeltaZ() const { return antenna2.z-antenna1.z; } }; struct Frequency { diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/foreachsimulatedbaselineaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/foreachsimulatedbaselineaction.h index 0059ab79a2be64bd248e94e15525bad9dc72a669..9d4878c3abad4d36eb5a9b00969d1825766f21a1 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/foreachsimulatedbaselineaction.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/foreachsimulatedbaselineaction.h @@ -43,10 +43,14 @@ namespace rfiStrategy { virtual ActionType Type() const { return ForEachSimulatedBaselineActionType; } virtual void Perform(ArtifactSet &artifacts, class ProgressListener &listener) { + double dec = 0.5*M_PI + 0.12800; + double ra = -0.03000; + double factor = 1.0; + struct Observatorium *observatorium = new WSRTObservatorium(); class Model *model = new Model(); - model->loadUrsaMajor(); - model->loadUrsaMajorDistortingSource(); + model->loadUrsaMajor(dec, ra, factor); + model->loadUrsaMajorDistortingSource(dec, ra, factor); if(observatorium != 0 && model != 0) { diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/frequencyconvolutionaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/frequencyconvolutionaction.h index e85e22dda54e9edb2573c94a5e37afef049cbd5a..515d38e45f3609642ad101880f64d130aaa9643e 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/frequencyconvolutionaction.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/frequencyconvolutionaction.h @@ -76,6 +76,9 @@ namespace rfiStrategy { enum KernelKind KernelKind() const { return _kernelKind; } void SetKernelKind(enum KernelKind kind) { _kernelKind = kind; } + + bool InSamples() const { return _inSamples; } + void SetInSamples(bool inSamples) { _inSamples = inSamples; } private: Image2DPtr sincConvolution(TimeFrequencyMetaDataCPtr metaData, Image2DCPtr source) { diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/fringestopaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/fringestopaction.h index 3e670db329a17311f58e27012403a6bed7226ca7..0da225e442f1da13a934b1236a7fc81eea3148c2 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/fringestopaction.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/fringestopaction.h @@ -21,6 +21,8 @@ #ifndef FRINGESTOPACTION_H #define FRINGESTOPACTION_H +#include <AOFlagger/msio/types.h> + #include <AOFlagger/strategy/actions/action.h> namespace rfiStrategy { @@ -28,13 +30,16 @@ namespace rfiStrategy { class FringeStopAction : public Action { public: - FringeStopAction() : _fringesToConsider(1.0L), _maxWindowSize(128), _fitChannelsIndividually(true), _onlyFringeStop(false) { } + FringeStopAction() : _fringesToConsider(1.0L), _maxWindowSize(128), _fitChannelsIndividually(true), _onlyFringeStop(false), _newPhaseCentreRA(0.0), _newPhaseCentreDec(0.5 * M_PInl) { } virtual ~FringeStopAction() { } + virtual std::string Description() { return "Fringe stop recovery"; } + virtual void Initialize() { } + virtual void Perform(class ArtifactSet &artifacts, class ProgressListener &listener); long double FringesToConsider() const { return _fringesToConsider; } @@ -52,12 +57,21 @@ namespace rfiStrategy { bool OnlyFringeStop() const { return _onlyFringeStop; } void SetOnlyFringeStop(bool onlyFringeStop) throw() { _onlyFringeStop = onlyFringeStop; } + virtual ActionType Type() const { return FringeStopActionType; } + + long double NewPhaseCentreRA() const { return _newPhaseCentreRA; } + void SetNewPhaseCentreRA(long double newPhaseCentreRA) { _newPhaseCentreRA = newPhaseCentreRA; } + + long double NewPhaseCentreDec() const { return _newPhaseCentreDec; } + void SetNewPhaseCentreDec(long double newPhaseCentreDec) { _newPhaseCentreDec = newPhaseCentreDec; } + private: long double _fringesToConsider; size_t _minWindowSize, _maxWindowSize; bool _fitChannelsIndividually; bool _onlyFringeStop; + long double _newPhaseCentreRA, _newPhaseCentreDec; }; } diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/timeconvolutionaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/timeconvolutionaction.h index 5d5c64c3e57681b308b4aa615f9406b1af79e156..c17d0045fbeb52bb5ae10bd216f18cad957121b6 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/timeconvolutionaction.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/timeconvolutionaction.h @@ -188,7 +188,7 @@ private: const num_t sincScale = ActualSincScaleInSamples(artifacts, band.channels[y].frequencyHz); if(y == image->Height()/2) { - AOLogger::Debug << "Horizontal sinc scale: " << sincScale << " (filter scale: " << ActualSincScaleAsRaDecDist(artifacts, band.channels[y].frequencyHz) * ((180.0/M_PI)*(60*60)) << " arcsec)\n"; + AOLogger::Debug << "Horizontal sinc scale: " << sincScale << " (filter scale: " << Angle::ToString(ActualSincScaleAsRaDecDist(artifacts, band.channels[y].frequencyHz)) << ")\n"; } if(sincScale > 1.0) { @@ -225,7 +225,7 @@ private: const num_t sincScale = ActualSincScaleInSamples(artifacts, band.channels[y].frequencyHz); if(y == image->Height()/2) { - AOLogger::Debug << "Horizontal sinc scale: " << sincScale << " (filter scale: " << ActualSincScaleAsRaDecDist(artifacts, band.channels[y].frequencyHz) * ((180.0/M_PI)*(60*60)) << " arcsec)\n"; + AOLogger::Debug << "Horizontal sinc scale: " << sincScale << " (filter scale: " << Angle::ToString(ActualSincScaleAsRaDecDist(artifacts, band.channels[y].frequencyHz)) << ")\n"; } if(sincScale > 1.0) { diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/algorithms/fringestoppingfitter.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/algorithms/fringestoppingfitter.h index 7a16952976d77ccedc05f127639918fa008c2313..ed71f34e4dacecbba06f8f13e323e2a4df749b1e 100644 --- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/algorithms/fringestoppingfitter.h +++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/algorithms/fringestoppingfitter.h @@ -99,6 +99,13 @@ class FringeStoppingFitter : public SurfaceFitMethod { void PerformDynamicFrequencyFit(unsigned yStart, unsigned yEnd); void PerformDynamicFrequencyFit(unsigned yStart, unsigned yEnd, unsigned windowSize); num_t GetAmplitude(unsigned yStart, unsigned yEnd); + + numl_t NewPhaseCentreRA() const { return _newPhaseCentreRA; } + void SetNewPhaseCentreRA(long double newPhaseCentreRA) { _newPhaseCentreRA = newPhaseCentreRA; } + + numl_t NewPhaseCentreDec() const { return _newPhaseCentreDec; } + void SetNewPhaseCentreDec(long double newPhaseCentreDec) { _newPhaseCentreDec = newPhaseCentreDec; } + private: num_t CalculateFitValue(const Image2D &image, size_t y); inline num_t CalculateMaskedAverage(const Image2D &image, size_t x, size_t yFrom, size_t yLength); @@ -128,6 +135,7 @@ class FringeStoppingFitter : public SurfaceFitMethod { bool _fitChannelsIndividually; bool _returnFittedValue, _returnMeanValue; bool _fringeFit; + numl_t _newPhaseCentreDec, _newPhaseCentreRA; }; #endif diff --git a/CEP/DP3/AOFlagger/src/gui/editstrategywindow.cpp b/CEP/DP3/AOFlagger/src/gui/editstrategywindow.cpp index 57011a3963dc6e66c9199cf618ac3d55bfc5280e..b224a82b6bc5d57621b5b78e54ba44a65b3051ed 100644 --- a/CEP/DP3/AOFlagger/src/gui/editstrategywindow.cpp +++ b/CEP/DP3/AOFlagger/src/gui/editstrategywindow.cpp @@ -38,6 +38,7 @@ #include <AOFlagger/gui/strategyframes/foreachmsframe.h> #include <AOFlagger/gui/strategyframes/foreachpolarisationframe.h> #include <AOFlagger/gui/strategyframes/foreachcomplexcomponentframe.h> +#include <AOFlagger/gui/strategyframes/frequencyconvolutionframe.h> #include <AOFlagger/gui/strategyframes/fringestoppingframe.h> #include <AOFlagger/gui/strategyframes/iterationframe.h> #include <AOFlagger/gui/strategyframes/plotframe.h> @@ -300,6 +301,9 @@ void EditStrategyWindow::onSelectionChanged() case ForEachPolarisationBlockType: showRight(new ForEachPolarisationFrame(*static_cast<rfiStrategy::ForEachPolarisationBlock*>(selectedAction), *this)); break; + case FrequencyConvolutionActionType: + showRight(new FrequencyConvolutionFrame(*static_cast<rfiStrategy::FrequencyConvolutionAction*>(selectedAction), *this)); + break; case PlotActionType: showRight(new StrategyPlotFrame(*static_cast<rfiStrategy::PlotAction*>(selectedAction), *this)); break; diff --git a/CEP/DP3/AOFlagger/src/gui/mswindow.cpp b/CEP/DP3/AOFlagger/src/gui/mswindow.cpp index ae36cc346bf1505ab7c35b469e12ad9504080b21..f5525b3a5c40e9bd5ba89e1f215ad0d1959f4213 100644 --- a/CEP/DP3/AOFlagger/src/gui/mswindow.cpp +++ b/CEP/DP3/AOFlagger/src/gui/mswindow.cpp @@ -596,6 +596,14 @@ void MSWindow::createToolbar() sigc::mem_fun(*this, &MSWindow::onSetAndShowImagePlane) ); _actionGroup->add( Gtk::Action::create("AddToImagePlane", "Add to _image plane"), sigc::mem_fun(*this, &MSWindow::onAddToImagePlane) ); + + Gtk::RadioButtonGroup setGroup; + _ncpSetButton = Gtk::RadioAction::create(setGroup, "NCPSet", "Use NCP set"); + _b1834SetButton = Gtk::RadioAction::create(setGroup, "B1834Set", "Use B1834 set"); + _ncpSetButton->set_active(true); + _actionGroup->add(_ncpSetButton); + _actionGroup->add(_b1834SetButton); + _actionGroup->add( Gtk::Action::create("SimulateCorrelation", "Simulate correlation"), sigc::mem_fun(*this, &MSWindow::onSimulateCorrelation) ); _actionGroup->add( Gtk::Action::create("SimulateSourceSetA", "Simulate source set A"), @@ -758,6 +766,8 @@ void MSWindow::createToolbar() " <menuitem action='SetAndShowImagePlane'/>" " <menuitem action='AddToImagePlane'/>" " <separator/>" + " <menuitem action='NCPSet'/>" + " <menuitem action='B1834Set'/>" " <menuitem action='SimulateCorrelation'/>" " <menuitem action='SimulateSourceSetA'/>" " <menuitem action='SimulateSourceSetB'/>" @@ -1490,11 +1500,28 @@ void MSWindow::showError(const std::string &description) dialog.run(); } +void MSWindow::getSetData(double &ra, double &dec, double &factor) +{ + if(_ncpSetButton->get_active()) + { + dec = 0.5*M_PI + 0.12800; + ra = -0.03000; + factor = 1.0; + } else { + dec = 1.083; + ra = 4.865; + factor = 4.0; + } +} + void MSWindow::onSimulateCorrelation() { + double ra, dec, factor; + getSetData(ra, dec, factor); + Model model; - model.loadUrsaMajor(); - model.loadUrsaMajorDistortingSource(); + model.loadUrsaMajor(ra, dec, factor); + model.loadUrsaMajorDistortingSource(ra, dec, factor); WSRTObservatorium wsrtObservatorium; model.SimulateObservation(*_imagePlaneWindow->GetImager(), wsrtObservatorium, -0.5*M_PIn-0.05, 0.05, 147000000.0); @@ -1503,12 +1530,15 @@ void MSWindow::onSimulateCorrelation() void MSWindow::onSimulateSourceSetA() { + double ra, dec, factor; + getSetData(ra, dec, factor); + Model model; - model.loadUrsaMajor(); - model.loadUrsaMajorDistortingSource(); + model.loadUrsaMajor(ra, dec, factor); + model.loadUrsaMajorDistortingSource(ra, dec, factor); WSRTObservatorium wsrtObservatorium; - std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, 0.5*M_PI + 0.12800, -0.03000, 147000000.0, 0, 5); + std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, dec, ra, 147000000.0, 0, 5); TimeFrequencyData data = pair.first; TimeFrequencyMetaDataCPtr metaData = pair.second; @@ -1518,12 +1548,15 @@ void MSWindow::onSimulateSourceSetA() void MSWindow::onSimulateSourceSetB() { + double ra, dec, factor; + getSetData(ra, dec, factor); + Model model; - model.loadUrsaMajor(); - model.loadUrsaMajorDistortingVariableSource(false, false); + model.loadUrsaMajor(ra, dec, factor); + model.loadUrsaMajorDistortingVariableSource(ra, dec, factor, false, false); WSRTObservatorium wsrtObservatorium; - std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, 0.5*M_PI + 0.12800, -0.03000, 147000000.0, 0, 5); + std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, dec, ra, 147000000.0, 0, 5); TimeFrequencyData data = pair.first; TimeFrequencyMetaDataCPtr metaData = pair.second; @@ -1533,12 +1566,15 @@ void MSWindow::onSimulateSourceSetB() void MSWindow::onSimulateSourceSetC() { + double ra, dec, factor; + getSetData(ra, dec, factor); + Model model; - model.loadUrsaMajor(); - model.loadUrsaMajorDistortingVariableSource(true, false); + model.loadUrsaMajor(ra, dec, factor); + model.loadUrsaMajorDistortingVariableSource(ra, dec, factor, true, false); WSRTObservatorium wsrtObservatorium; - std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, 0.5*M_PI + 0.12800, -0.03000, 147000000.0, 0, 5); + std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, dec, ra, 147000000.0, 0, 5); TimeFrequencyData data = pair.first; TimeFrequencyMetaDataCPtr metaData = pair.second; @@ -1548,12 +1584,15 @@ void MSWindow::onSimulateSourceSetC() void MSWindow::onSimulateSourceSetD() { + double ra, dec, factor; + getSetData(ra, dec, factor); + Model model; - model.loadUrsaMajor(); - model.loadUrsaMajorDistortingVariableSource(false, true); + model.loadUrsaMajor(ra, dec, factor); + model.loadUrsaMajorDistortingVariableSource(ra, dec, factor, false, true); WSRTObservatorium wsrtObservatorium; - std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, 0.5*M_PI + 0.12800, -0.03000, 147000000.0, 0, 5); + std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> pair = model.SimulateObservation(wsrtObservatorium, dec, ra, 147000000.0, 0, 5); TimeFrequencyData data = pair.first; TimeFrequencyMetaDataCPtr metaData = pair.second; @@ -1563,12 +1602,15 @@ void MSWindow::onSimulateSourceSetD() void MSWindow::onSimulateFourProductCorrelation() { + double ra, dec, factor; + getSetData(ra, dec, factor); + Model model; - model.loadUrsaMajor(); + model.loadUrsaMajor(ra, dec, factor); WSRTObservatorium wsrtObservatorium; FourProductCorrelatorTester fpcTester(model, *_imagePlaneWindow->GetImager(), wsrtObservatorium); - fpcTester.SimulateObservation(0.5*M_PI + 0.12800, -0.03000, 147000000.0); + fpcTester.SimulateObservation(dec, ra, 147000000.0); _imagePlaneWindow->Update(); } diff --git a/CEP/DP3/AOFlagger/src/imaging/model.cpp b/CEP/DP3/AOFlagger/src/imaging/model.cpp index 670bd607d395c9fc1704cb28e2b22e36d4134907..3a1c62dd832e40ad81106c5d6577cbad0fa9dc81 100644 --- a/CEP/DP3/AOFlagger/src/imaging/model.cpp +++ b/CEP/DP3/AOFlagger/src/imaging/model.cpp @@ -208,52 +208,34 @@ void Model::GetUVPosition(num_t &u, num_t &v, num_t earthLattitudeAngle, num_t d v = -sinn(baselineAngle)*baselineLength; } -num_t Model::GetWPosition(num_t delayDirectionDec, num_t delayDirectionRA, num_t frequency, num_t earthLattitudeAngle, num_t dx, num_t dy) -{ - num_t wavelength = 299792458.0L / frequency; - num_t raSinEnd = sinn(-delayDirectionRA - earthLattitudeAngle); - num_t raCosEnd = cosn(-delayDirectionRA - earthLattitudeAngle); - num_t decSin = sinn(delayDirectionDec); - // term "+ dz * decCos" is eliminated because of subtraction - //num_t wPosition = - // (dx*raCosEnd - dy*raSinEnd) * (-decCos) / wavelength; - num_t wPosition = - (dx*raCosEnd - dy*raSinEnd) * (-decSin) / wavelength; - return wPosition; -} - -void Model::loadUrsaMajor() +void Model::loadUrsaMajor(double ra, double dec, double factor) { double - s = 0.00005, //scale - rs = 8.0; // stretch in dec - double cd = 0.5*M_PI + 0.12800; - double cr = -0.03000; + s = 0.00005 * factor, //scale + rs = 6.0 + 2.0 * factor; // stretch in dec double fluxoffset = 0.0; - AddSource(cd + s*rs*40, cr + s*72, 8.0/8.0 + fluxoffset); // Dubhe - AddSource(cd + s*rs*-16, cr + s*81, 4.0/8.0 + fluxoffset); // Beta - AddSource(cd + s*rs*-45, cr + s*2, 3.0/8.0 + fluxoffset); // Gamma - AddSource(cd + s*rs*-6, cr + s*-27, 2.0/8.0 + fluxoffset); // Delta - AddSource(cd + s*rs*-4, cr + s*-85, 6.0/8.0 + fluxoffset); // Alioth - 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 + AddSource(dec + s*rs*40, ra + s*72, 8.0/8.0 + fluxoffset); // Dubhe + AddSource(dec + s*rs*-16, ra + s*81, 4.0/8.0 + fluxoffset); // Beta + AddSource(dec + s*rs*-45, ra + s*2, 3.0/8.0 + fluxoffset); // Gamma + AddSource(dec + s*rs*-6, ra + s*-27, 2.0/8.0 + fluxoffset); // Delta + AddSource(dec + s*rs*-4, ra + s*-85, 6.0/8.0 + fluxoffset); // Alioth + AddSource(dec + s*rs*2, ra + s*-131, 5.0/8.0 + fluxoffset); // Zeta + AddSource(dec + s*rs*-36, ra + s*-192, 7.0/8.0 + fluxoffset); // Alkaid //AddSource(cd, cr - M_PI, 4.0); } -void Model::loadUrsaMajorDistortingSource() +void Model::loadUrsaMajorDistortingSource(double ra, double dec, double factor) { - double fluxoffset = 0.0; - - AddSource(0.5*M_PI, 0, 4.0 + fluxoffset); // NCP + AddSource(dec - 0.12800 * factor, ra + 0.015 + 0.015 * factor, 4.0); } -void Model::loadUrsaMajorDistortingVariableSource(bool weak, bool slightlyMiss) +void Model::loadUrsaMajorDistortingVariableSource(double ra, double dec, double factor, bool weak, bool slightlyMiss) { double flux = 4.0; - double dec = 0.5*M_PI; - double ra = 0.0; + dec = dec - 0.12800 * factor; + ra = ra + 0.015 + 0.015 * factor; if(slightlyMiss) { dec += 0.005; @@ -263,6 +245,6 @@ void Model::loadUrsaMajorDistortingVariableSource(bool weak, bool slightlyMiss) { flux /= 100.0; } - AddVariableSource(dec, ra, flux); // NCP + AddVariableSource(dec, ra, flux); } diff --git a/CEP/DP3/AOFlagger/src/imaging/uvimager.cpp b/CEP/DP3/AOFlagger/src/imaging/uvimager.cpp index b9706f076a3cf9ef1a779543b34be2bcd83162c2..2231e0ce76681076b032a4d8c96b1afe0382572e 100644 --- a/CEP/DP3/AOFlagger/src/imaging/uvimager.cpp +++ b/CEP/DP3/AOFlagger/src/imaging/uvimager.cpp @@ -508,26 +508,11 @@ num_t UVImager::GetFringeStopFrequency(size_t timeIndex, const Baseline &baselin const num_t earthSpeed = 2.0L * M_PIn / (24.0L * 60.0L * 60.0L); num_t earthLattitudeAngle = Date::JDToHourOfDay(Date::AipsMJDToJD(metaData->ObservationTimes()[timeIndex]))*M_PIn/12.0L; - //long double u, v; - //GetUVPosition(u, v, baseline, time, delayDirectionRA, delayDirectionDec, frequency); - //return - // earthSpeed * (u * sinn(delayDirectionRA-earthLattitudeAngle) + v * cosn(delayDirectionRA-earthLattitudeAngle)) * cosn(delayDirectionDec); num_t raSin = sinn(-delayDirectionRA - earthLattitudeAngle); num_t raCos = cosn(-delayDirectionRA - earthLattitudeAngle); num_t dx = baseline.antenna2.x - baseline.antenna1.x; num_t dy = baseline.antenna2.y - baseline.antenna1.y; - //num_t dz = baseline.antenna2.z - baseline.antenna1.z; num_t wavelength = 299792458.0L / frequency; - /*std::cout << "Angle=" << - 180.0L / M_PIn * acosn(((dx * raCos - dy * raSin) * cosn(delayDirectionDec) + dz*sinn(delayDirectionDec)) / sqrtn(dx*dx + dy*dy + dz*dz)) - << std::endl; - std::cout << "delay=" << - ((dx * raCos - dy * raSin) * cosn(delayDirectionDec) + dz*sinn(delayDirectionDec)) << "m" - << std::endl; - std::cout << "ddelay/dt=" << - (earthSpeed * (dx*raSin + dy*raCos) * cosn(delayDirectionDec)) - << "m/s" << std::endl;*/ - //return (earthSpeed * (dx*raSin + dy*raCos) * cosn(delayDirectionDec)) / wavelength; return -earthSpeed * metaData->UVW()[timeIndex].u * cosn(delayDirectionDec); } @@ -537,29 +522,6 @@ num_t UVImager::GetFringeCount(size_t timeIndexStart, size_t timeIndexEnd, unsig // with the fringe stop frequency returned above otherwise; probably because of a // mismatch in the signs of u,v,w somewhere... return -(metaData->UVW()[timeIndexEnd].w - metaData->UVW()[timeIndexStart].w) * metaData->Band().channels[channelIndex].frequencyHz / 299792458.0L; - /*double - timeStart = metaData->ObservationTimes()[timeIndexStart], - timeEnd = metaData->ObservationTimes()[timeIndexEnd]; - num_t earthLattitudeAngleStart = - Date::JDToHourOfDay(Date::AipsMJDToJD(timeStart))*M_PIn/12.0L; - num_t earthLattitudeAngleEnd = - Date::JDToHourOfDay(Date::AipsMJDToJD(timeEnd))*M_PIn/12.0L; - num_t wavelength = 299792458.0L / metaData->Band().channels[channelIndex].frequencyHz; - num_t dx = metaData->Antenna2().position.x - metaData->Antenna1().position.x; - num_t dy = metaData->Antenna2().position.y - metaData->Antenna1().position.y; - num_t delayDirectionRA = metaData->Field().delayDirectionRA; - num_t raSinStart = sinn(-delayDirectionRA - earthLattitudeAngleStart); - num_t raCosStart = cosn(-delayDirectionRA - earthLattitudeAngleStart); - num_t raSinEnd = sinn(-delayDirectionRA - earthLattitudeAngleEnd); - num_t raCosEnd = cosn(-delayDirectionRA - earthLattitudeAngleEnd); - num_t decCos = cosn(metaData->Field().delayDirectionDec); - // term "+ dz * decCos" is eliminated because of subtraction - num_t fringeCount = - ( (dx*raCosStart - dy*raSinStart) - - - (dx*raCosEnd - dy*raSinEnd) ) * (-decCos) / wavelength; - //std::cout << "Fringecount = " << fringeCount << " within " <<timeStart << "-" << timeEnd << "=" << (timeEnd-timeStart) << std::endl; - return fringeCount;*/ } void UVImager::InverseImage(class MeasurementSet &prototype, unsigned /*band*/, const Image2D &/*uvReal*/, const Image2D &/*uvImaginary*/, unsigned antenna1Index, unsigned antenna2Index) diff --git a/CEP/DP3/AOFlagger/src/strategy/actions/fringestopaction.cpp b/CEP/DP3/AOFlagger/src/strategy/actions/fringestopaction.cpp index 957caba1fc56f5d72d460415b151720171cc6f4b..07df84180ac4c12a233404b4b7b61dad5dfdddb6 100644 --- a/CEP/DP3/AOFlagger/src/strategy/actions/fringestopaction.cpp +++ b/CEP/DP3/AOFlagger/src/strategy/actions/fringestopaction.cpp @@ -38,13 +38,15 @@ namespace rfiStrategy { throw BadUsageException("Baseline antenna info not set"); if(!artifacts.MetaData()->HasObservationTimes()) throw BadUsageException("Baseline observation times not set"); - + FringeStoppingFitter fitter; fitter.SetFringesToConsider(_fringesToConsider); fitter.SetMinWindowSize(_minWindowSize); fitter.SetMaxWindowSize(_maxWindowSize); fitter.SetFitChannelsIndividually(_fitChannelsIndividually); fitter.SetMetaData(artifacts.MetaData()); + fitter.SetNewPhaseCentreRA(_newPhaseCentreRA); + fitter.SetNewPhaseCentreDec(_newPhaseCentreDec); fitter.Initialize(artifacts.ContaminatedData()); if(_onlyFringeStop) fitter.PerformFringeStop(); diff --git a/CEP/DP3/AOFlagger/src/strategy/algorithms/fringestoppingfitter.cpp b/CEP/DP3/AOFlagger/src/strategy/algorithms/fringestoppingfitter.cpp index 09df49ae792b8b92e36ab8294907f650008616db..fd19352802c0b2bcfa9519b3ef5141b0355f1ba9 100644 --- a/CEP/DP3/AOFlagger/src/strategy/algorithms/fringestoppingfitter.cpp +++ b/CEP/DP3/AOFlagger/src/strategy/algorithms/fringestoppingfitter.cpp @@ -28,7 +28,7 @@ #include <AOFlagger/strategy/algorithms/sinusfitter.h> FringeStoppingFitter::FringeStoppingFitter() : - _originalData(0), _fringesToConsider(1.0), _minWindowSize(32), _maxWindowSize(128), _returnFittedValue(false), _returnMeanValue(false), _fringeFit(true) + _originalData(0), _fringesToConsider(1.0), _minWindowSize(32), _maxWindowSize(128), _returnFittedValue(false), _returnMeanValue(false), _fringeFit(true), _newPhaseCentreDec(M_PInl*0.5), _newPhaseCentreRA(0.0) { } @@ -270,8 +270,13 @@ num_t FringeStoppingFitter::GetFringeFrequency(size_t x, size_t y) void FringeStoppingFitter::GetRFIValue(num_t &r, num_t &i, int x, int y, num_t rfiPhase, num_t rfiStrength) { - num_t rotations = - UVImager::GetFringeCount(0, x, y, _metaData); + const numl_t earthRotation = UVImager::TimeToEarthLattitude(x, _metaData); + const Baseline baseline = _metaData->Baseline(); + const numl_t newWPos = + UVImager::GetWPosition(_newPhaseCentreDec, _newPhaseCentreRA, _metaData->Band().channels[y].frequencyHz, earthRotation, baseline.DeltaX(), baseline.DeltaY()); + const num_t rotations = + UVImager::GetFringeCount(0, x, y, _metaData) - newWPos; + r = cosn(rotations * 2.0 * M_PIn + rfiPhase) * rfiStrength; i = -sinn(rotations * 2.0 * M_PIn + rfiPhase) * rfiStrength; } @@ -298,6 +303,10 @@ void FringeStoppingFitter::MinimizeRFIFitError(num_t &phase, num_t &litude, S num_t sumR = 0.0, sumI = 0.0; size_t n = 0; + + const num_t dx = _metaData->Antenna2().position.x - _metaData->Antenna1().position.x; + const num_t dy = _metaData->Antenna2().position.y - _metaData->Antenna1().position.y; + for(unsigned t=xStart;t<xEnd;++t) { const num_t vR = real->Value(t); @@ -305,13 +314,16 @@ void FringeStoppingFitter::MinimizeRFIFitError(num_t &phase, num_t &litude, S if(std::isfinite(vR) && std::isfinite(vI)) { + const num_t tRotation = UVImager::TimeToEarthLattitude(_metaData->ObservationTimes()[t]); const num_t tauge = UVImager::GetFringeCount(0, t, y, _metaData); + const num_t taugeNew = UVImager::GetWPosition(_newPhaseCentreDec, _newPhaseCentreRA, _metaData->Band().channels[y].frequencyHz, tRotation, dx, dy); + const num_t phaseShift = tauge - taugeNew; - sumR += vR * cosn(-2.0 * M_PIn * tauge); - sumR += vI * sinn(-2.0 * M_PIn * tauge); + sumR += vR * cosn(-2.0 * M_PIn * phaseShift); + sumR += vI * sinn(-2.0 * M_PIn * phaseShift); - sumI += vR * sinn(-2.0 * M_PIn * tauge); - sumI -= vI * cosn(-2.0 * M_PIn * tauge); + sumI += vR * sinn(-2.0 * M_PIn * phaseShift); + sumI -= vI * cosn(-2.0 * M_PIn * phaseShift); ++n; } } @@ -376,7 +388,7 @@ void FringeStoppingFitter::PerformDynamicFrequencyFitOnOneRow(SampleRowCPtr real else windowEnd = real->Size(); num_t windowPhase, windowStrength; - MinimizeRFIFitError(windowPhase, windowStrength, real, imaginary, windowStart, windowEnd, y); + MinimizeRFIFitError(windowPhase, windowStrength, real, imaginary, windowStart, windowEnd, y); num_t rfiR, rfiI; GetRFIValue(rfiR, rfiI, x, y, windowPhase, windowStrength);