File: uvimager.h

package info (click to toggle)
aoflagger 2.13.0-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 4,232 kB
  • sloc: cpp: 61,805; python: 60; sh: 23; makefile: 8
file content (171 lines) | stat: -rw-r--r-- 6,660 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
#ifndef UVIMAGER_H
#define UVIMAGER_H

#include "../structures/timefrequencymetadata.h"
#include "../structures/measurementset.h"
#include "../structures/date.h"

#include "../structures/timefrequencydata.h"

struct SingleFrequencySingleBaselineData {
	casacore::Complex data;
	bool flag;
	bool available;
	double time;
	unsigned field;
};

class UVImager {
	public:
		enum ImageKind { Homogeneous, Flagging };
		UVImager(unsigned long xRes, unsigned long yRes, ImageKind imageKind=Homogeneous);
		~UVImager();
		void Image(class MeasurementSet &measurementSet, unsigned band);
		void Image(class MeasurementSet &measurementSet, unsigned band, const class IntegerDomain &frequencies);
		void Image(const class TimeFrequencyData &data, TimeFrequencyMetaDataCPtr metaData, unsigned frequencyIndex);
		void Image(const class TimeFrequencyData &data, TimeFrequencyMetaDataCPtr metaData)
		{
			for(unsigned y=0;y<data.ImageHeight();++y)
				Image(data, metaData, y);
		}
		void Image(const class TimeFrequencyData &data, class SpatialMatrixMetaData *metaData);
		void InverseImage(class MeasurementSet &prototype, unsigned band, const class Image2D &uvReal, const class Image2D &uvImaginary, unsigned antenna1, unsigned antenna2);
		const class Image2D &WeightImage() const { return _uvWeights; }
		const class Image2D &RealUVImage() const { return _uvReal; }
		const class Image2D &ImaginaryUVImage() const { return _uvImaginary; }
		void SetInvertFlagging(bool newValue) { _invertFlagging = newValue; }
		void SetDirectFT(bool directFT) { _directFT = directFT; }

		/**
		 * This function calculates the uv position, but it's not optimized for speed, so it's not to be used in an imager.
		 * @param [out] u the u position (in the uv-plane domain)
		 * @param [out] v the v position (in the uv-plane domain)
		 * @param [in] timeIndex the time index to calculate the u,v position for
		 * @param [in] frequencyIndex the frequency index to calculate the u,v position for
		 * @param [in] metaData information about the baseline
		 */
		static void GetUVPosition(num_t &u, num_t &v, size_t timeIndex, size_t frequencyIndex, TimeFrequencyMetaDataCPtr metaData);

		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.Empty(); }
		bool HasFFT() const { return !_uvFTReal.Empty(); }
		const class Image2D &FTReal() const { return _uvFTReal; }
		const class Image2D &FTImaginary() const { return _uvFTImaginary; }
		class Image2D &FTReal() { return _uvFTReal; }
		class Image2D &FTImaginary() { return _uvFTImaginary; }
		void SetUVScaling(num_t newScale)
		{
			_uvScaling = newScale;
		}
		num_t UVScaling() const {
			return _uvScaling;
		}
		void ApplyWeightsToUV();
		void SetUVValue(num_t u, num_t v, num_t r, num_t i, num_t weight);

		template<typename T>
		static T FrequencyToWavelength(const T frequency)
		{
			return SpeedOfLight() / frequency; 
		}
		static long double SpeedOfLight()
		{
			return 299792458.0L;
		}
		numl_t ImageDistanceToDecRaDistance(numl_t imageDistance) const
		{
			return imageDistance * _uvScaling;
		}
		static numl_t AverageUVDistance(TimeFrequencyMetaDataCPtr metaData, const double frequencyHz)
		{
			const std::vector<UVW> &uvw = metaData->UVW();
			numl_t avgDist = 0.0;
			for(std::vector<UVW>::const_iterator i=uvw.begin();i!=uvw.end();++i)
			{
				numl_t dist = i->u*i->u + i->v*i->v;
				avgDist += sqrtnl(dist);
			}
			return avgDist * frequencyHz / (SpeedOfLight() * (numl_t) uvw.size());
		}
		static numl_t UVTrackLength(TimeFrequencyMetaDataCPtr metaData, const double frequencyHz)
		{
			const std::vector<UVW> &uvw = metaData->UVW();
			numl_t length = 0.0;
			std::vector<UVW>::const_iterator i=uvw.begin();
			if(i == uvw.end()) return 0.0;
			while((i+1)!=uvw.end())
			{
				std::vector<UVW>::const_iterator n=i;
				++n;
				const numl_t
					du = n->u - i->u,
					dv = n->v - i->v;
				length += sqrtnl(du*du + dv*dv);
				i=n;
			}
			return length * frequencyHz / SpeedOfLight();
		}
		numl_t ImageDistanceToFringeSpeedInSamples(numl_t imageDistance, double frequencyHz, TimeFrequencyMetaDataCPtr metaData) const
		{
			return ImageDistanceToDecRaDistance(imageDistance) * AverageUVDistance(metaData, frequencyHz) / (0.5 * (numl_t) metaData->UVW().size());
		}
	private:
		void Clear();
		struct AntennaCache {
			num_t wavelength;
			num_t dx, dy, dz;
		};
		void Image(const class IntegerDomain &frequencies);
		void Image(const IntegerDomain &frequencies, const IntegerDomain &antenna1Domain, const IntegerDomain &antenna2Domain);
		void Image(unsigned frequencyIndex, class AntennaInfo &antenna1, class AntennaInfo &antenna2, SingleFrequencySingleBaselineData *data);

		// This is the fast variant.
		void GetUVPosition(num_t &u, num_t &v, const SingleFrequencySingleBaselineData &data, const AntennaCache &cache);
		void SetUVFTValue(num_t u, num_t v, num_t r, num_t i, num_t weight);


		unsigned long _xRes, _yRes;
		unsigned long _xResFT, _yResFT;
		num_t _uvScaling;
		class Image2D _uvReal, _uvImaginary, _uvWeights;
		class Image2D _uvFTReal, _uvFTImaginary;
		class Image2D _timeFreq;
		MeasurementSet *_measurementSet;
		unsigned _antennaCount, _fieldCount;
		AntennaInfo *_antennas;
		BandInfo _band;
		FieldInfo *_fields;
		size_t _scanCount;
		ImageKind _imageKind;
		bool _invertFlagging, _directFT;
		bool _ignoreBoundWarnings;
};

#endif