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
|
#ifndef MODEL_H
#define MODEL_H
#include <vector>
#include <cmath>
#include <utility>
#include "../structures/image2d.h"
#include "../structures/types.h"
#include "uvimager.h"
template <typename T>
struct OutputReceiver {
void SetY(size_t) {}
};
template <>
struct OutputReceiver<UVImager> {
UVImager* _imager;
void SetUVValue(size_t, double u, double v, double r, double i, double w) {
_imager->SetUVValue(u, v, r, i, w);
_imager->SetUVValue(-u, -v, r, -i, w);
}
void SetY(size_t) {}
};
template <>
struct OutputReceiver<TimeFrequencyData> {
Image2DPtr _real, _imaginary;
size_t _y;
void SetUVValue(size_t x, double, double, double r, double i, double) {
_real->SetValue(x, _y, r);
_imaginary->SetValue(x, _y, i);
}
void SetY(size_t y) { _y = y; }
};
class Model {
struct Source {
virtual ~Source() {}
virtual numl_t Dec(num_t t) const = 0;
virtual numl_t Ra(num_t t) const = 0;
virtual numl_t FluxIntensity(num_t t) const = 0;
virtual numl_t SqrtFluxIntensity(num_t t) const {
return std::sqrt(FluxIntensity(t));
}
};
struct StablePointSource final : public Source {
long double dec, ra, fluxIntensity, sqrtFluxIntensity;
numl_t Dec(num_t) const override { return dec; }
numl_t Ra(num_t) const override { return ra; }
numl_t FluxIntensity(num_t) const override { return fluxIntensity; }
numl_t SqrtFluxIntensity(num_t) const override { return sqrtFluxIntensity; }
};
struct VariablePointSource final : public Source {
long double dec, ra, fluxIntensity;
double peakTime, oneOverSigmaSq;
numl_t Dec(num_t) const override { return dec; }
numl_t Ra(num_t) const override { return ra; }
numl_t FluxIntensity(num_t t) const override {
numl_t mu = std::fmod(std::fabs(t - peakTime), 1.0);
if (mu > 0.5) mu = 1.0 - mu;
return fluxIntensity * (1.0 + std::exp(mu * mu * oneOverSigmaSq)) *
(1.0 + std::fmod(t * 1007.0, 13.0) / 26.0);
}
};
public:
Model();
void AddSource(long double dec, long double ra, long double fluxIntensity) {
StablePointSource* source = new StablePointSource();
source->dec = dec;
source->ra = ra;
source->fluxIntensity = fluxIntensity;
source->sqrtFluxIntensity = sqrt(fluxIntensity);
_sources.push_back(source);
}
void AddVariableSource(long double dec, long double ra,
long double fluxIntensity) {
VariablePointSource* source = new VariablePointSource();
source->dec = dec;
source->ra = ra;
source->fluxIntensity = fluxIntensity;
source->peakTime = 0.2;
source->oneOverSigmaSq = 1.0 / (0.3 * 0.3);
_sources.push_back(source);
}
void SimulateAntenna(double time, num_t delayDirectionDEC,
num_t delayDirectionRA, num_t dx, num_t dy,
num_t frequency, num_t earthLattitude, num_t& r,
num_t& i);
void SimulateUncoherentAntenna(double time, num_t delayDirectionDEC,
num_t delayDirectionRA, num_t dx, num_t dy,
num_t frequency, num_t earthLattitude,
num_t& r, num_t& i, size_t index);
template <typename T>
void SimulateCorrelation(struct OutputReceiver<T>& receiver,
num_t delayDirectionDEC, num_t delayDirectionRA,
num_t dx, num_t dy, num_t dz, num_t frequency,
num_t channelWidth, size_t nTimes,
double integrationTime);
void SimulateObservation(class UVImager& imager,
class Observatorium& observatorium,
num_t delayDirectionDEC, num_t delayDirectionRA) {
srand(1);
OutputReceiver<UVImager> imagerOutputter;
imagerOutputter._imager = &imager;
SimulateObservation(imagerOutputter, observatorium, delayDirectionDEC,
delayDirectionRA);
}
std::pair<TimeFrequencyData, TimeFrequencyMetaDataPtr> SimulateObservation(
size_t nTimes, class Observatorium& observatorium,
num_t delayDirectionDEC, num_t delayDirectionRA, size_t a1, size_t a2);
template <typename T>
void SimulateObservation(struct OutputReceiver<T>& receiver,
class Observatorium& observatorium,
num_t delayDirectionDEC, num_t delayDirectionRA);
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) {
return UVImager::GetWPosition(delayDirectionDec, delayDirectionRA,
frequency, earthLattitudeAngle, dx, dy);
}
void loadUrsaMajor(double ra, double dec, double factor);
void loadUrsaMajorDistortingSource(double ra, double dec, double factor,
bool slightlyMiss = false);
void loadUrsaMajorDistortingVariableSource(double ra, double dec,
double factor, bool weak = false,
bool slightlyMiss = false);
double NoiseSigma() const { return _noiseSigma; }
void SetNoiseSigma(double noiseSigma) { _noiseSigma = noiseSigma; }
private:
std::vector<Source*> _sources;
double _noiseSigma, _sourceSigma;
double _integrationTime;
};
#endif
|
