ARTS 2.5.11 (git: 725533f0)
m_raw.cc
Go to the documentation of this file.
1
10#include "artstime.h"
11#include "raw.h"
12
13
14void yColdAtmHot(Vector& y,
15 const Vector& cold,
16 const Vector& atm,
17 const Vector& hot,
18 const Numeric& cold_temp,
19 const Numeric& hot_temp,
20 const Index& calib,
21 const Verbosity&)
22{
23 ARTS_USER_ERROR_IF (cold.nelem() not_eq atm.nelem() or atm.nelem() not_eq hot.nelem(),
24 "Length of vectors must be correct");
25
26 y.resize(atm.nelem());
27 if (calib) {
28 y = Raw::Calibration::calibration(cold, atm, hot, cold_temp, hot_temp);
29 } else {
30 y = Raw::Calibration::systemtemp(cold, hot, cold_temp, hot_temp);
31 }
32}
33
34
35void ybatchColdAtmHotAtmCycle(ArrayOfVector& ybatch,
36 ArrayOfTime& sensor_time,
37 const ArrayOfVector& level0_data,
38 const ArrayOfTime& level0_time,
39 const Vector& cold_temp,
40 const Vector& hot_temp,
41 const Index& first_c_index,
42 const Verbosity&)
43{
44 ARTS_USER_ERROR_IF(level0_data.nelem() not_eq cold_temp.nelem() or
45 level0_data.nelem() not_eq hot_temp.nelem(),
46 "Length of vectors must be correct");
47 ARTS_USER_ERROR_IF (level0_time.nelem() not_eq level0_data.nelem() and
48 level0_time.nelem() not_eq 0,
49 "Bad level0_time length, must be empty of same as level0_data");
50
51 ybatch = Raw::Calibration::caha(level0_data, cold_temp, hot_temp, first_c_index);
52
53 // Fix time using the same method as CAHA
54 if (level0_time.nelem()) {
55 sensor_time.resize(ybatch.nelem());
56 // First position as described by method is at H if this index is too large
57 const Index pos = first_c_index - ((first_c_index > 1) ? 2 : 0);
58 for (Index i=0; i<sensor_time.nelem(); i++) {
59 sensor_time[i] = level0_time[pos + 2*i];
60 }
61 }
62}
63
64
65void ybatchTimeAveraging(ArrayOfVector& ybatch,
66 ArrayOfTime& sensor_time,
67 const String& time_step,
68 const Index& disregard_first,
69 const Index& disregard_last,
70 const Verbosity&)
71{
72 // Size of problem
73 const Index n=sensor_time.nelem();
74 ARTS_USER_ERROR_IF (sensor_time.nelem() not_eq n,
75 "Time vector length must match input data length");
76
77 // Time is not decreasing
78 ARTS_USER_ERROR_IF (not std::is_sorted(sensor_time.cbegin(), sensor_time.cend()),
79 "Time vector cannot decrease");
80
81 // Find the limits of the range
82 auto lims = time_steps(sensor_time, time_stepper_selection(time_step));
83
84 // Output variables
85 ArrayOfVector ybatch_out;
86 ArrayOfTime sensor_time_out;
87
88 if (lims.front() == n) {
89 ybatch_out.resize(0);
90 sensor_time_out.resize(0);
91 } else {
92
93 // Frequency grids
94 const Index k = ybatch[0].nelem();
95 ARTS_USER_ERROR_IF (not std::all_of(ybatch.cbegin(), ybatch.cend(),
96 [k](auto& x){return x.nelem() == k;}),
97 "Bad frequency grid size in input data; expects all equal");
98
99 // Allocate output
100 const Index m = lims.nelem() - 1;
101 ARTS_USER_ERROR_IF (m < 0,
102 "Must include last if time step covers all of the range");
103 ybatch_out = ArrayOfVector(m, Vector(k));
104 sensor_time_out = ArrayOfTime(m);
105
106 // Fill output
107 #pragma omp parallel for if (not arts_omp_in_parallel()) schedule(guided)
108 for (Index i=0; i<m; i++) {
109 sensor_time_out[i] = mean_time(sensor_time, lims[i], lims[i+1]);
110 Raw::Average::nanavg(ybatch_out[i], ybatch, lims[i], lims[i+1]);
111 }
112 }
113
114 if (disregard_first) {
115 ybatch_out.erase(ybatch_out.begin());
116 sensor_time_out.erase(sensor_time_out.begin());
117 }
118
119 if (disregard_last) {
120 ybatch_out.pop_back();
121 sensor_time_out.pop_back();
122 }
123
124 ybatch = ybatch_out;
125 sensor_time = sensor_time_out;
126}
127
128
130 ArrayOfVector& ybatch,
131 const ArrayOfIndex& range,
132 const Vector& trop_temp,
133 const Numeric& targ_temp,
134 const Verbosity&)
135{
136 // Size of problem
137 const Index n=ybatch.nelem();
138 const Index m=n ? ybatch[0].nelem() : 0;
139
140 ARTS_USER_ERROR_IF (m == 0, "A frequency range is required")
141
142 ARTS_USER_ERROR_IF (std::any_of(ybatch.begin(), ybatch.end(),
143 [m](auto& y){return y.nelem() not_eq m;}),
144 "Bad input size, all of ybatch must match itself");
145
146 ARTS_USER_ERROR_IF (trop_temp.nelem() not_eq n and trop_temp.nelem() not_eq 1,
147 "Bad input size, trop_temp must match ybatch or be 1-long,\n"
148 "trop_temp: [", trop_temp, "]\ntrop_temp.nelem(): ",
149 trop_temp.nelem(), "\nybatch.nelem(): ", n);
150
151 // This algorithm stores partial transmission and median and tropospheric temperature in the correction terms
152 ybatch_corr = ArrayOfVector(n, Vector(3));
153
154 // Compute tropospheric correction
155 for (Index i=0; i<n; i++) {
156 ybatch_corr[i][2] = trop_temp.nelem() > 1 ? trop_temp[i] : trop_temp[0];
157 ybatch_corr[i][0] = Raw::Average::nanmedian(ybatch[i], range);
158 ybatch_corr[i][1] = std::exp(- std::log((ybatch_corr[i][2] - ybatch_corr[i][0]) / (ybatch_corr[i][2] - targ_temp)));
159 }
160
161 // Apply correction
162 for (Index i=0; i<n; i++) {
163 ybatch[i] *= ybatch_corr[i][1];
164 ybatch[i] += ybatch_corr[i][2] * (1 - ybatch_corr[i][1]);
165 }
166}
167
168
170 const ArrayOfVector& ybatch_corr,
171 const Verbosity&)
172{
173 // Size of problem
174 const Index n=ybatch.nelem();
175 ARTS_USER_ERROR_IF ((std::any_of(ybatch_corr.begin(), ybatch_corr.end(),
176 [](auto& corr){return corr.nelem() not_eq 3;})) or ybatch_corr.nelem() not_eq n,
177 "Bad input size, all of ybatch_corr must match ybatch and have three elements each");
178
179 // Apply inverse of correction
180 for (Index i=0; i<n; i++) {
181 ybatch[i] -= ybatch_corr[i][2] * (1 - ybatch_corr[i][1]);
182 ybatch[i] /= ybatch_corr[i][1];
183 }
184}
185
187 const Numeric& dx,
188 const Verbosity&) {
189 const Numeric median = Raw::Average::nanmedian(y);
190 auto mask = Raw::Mask::out_of_bounds(y, median-dx, median+dx);
191
192 for (Index i=0; i<y.nelem(); i++) {
193 y[i] = mask[i] ? std::numeric_limits<Numeric>::quiet_NaN() : y[i];
194 }
195}
196
197void ybatchMaskOutsideMedianRange(ArrayOfVector& ybatch,
198 const Numeric& dx,
199 const Verbosity& verbosity) {
200 for (Vector& y: ybatch) yMaskOutsideMedianRange(y, dx, verbosity);
201}
202
203void yDoublingMeanFocus(Vector& f_grid,
204 Vector& y,
205 const Numeric& f0,
206 const Numeric& df,
207 const Verbosity&) {
208 if (f_grid.nelem() not_eq y.nelem()) {
209 throw std::runtime_error("f_grid and y must have the same size");
210 }
211
212 if (not is_increasing(f_grid)) {
213 throw std::runtime_error("f_grid must be sorted and ever increasing");
214 }
215
216 if (f_grid.nelem() < 2) {
217 throw std::runtime_error("Must have at least 2 frequency grid points");
218 }
219
220 // Sets defaults by method description
221 const Numeric F0 = f0 > 0 ? f0 : mean(f_grid);
222 const Numeric DF = df > 0 ? df : 10 * (f_grid[1] - f_grid[0]);
223
224 if (f_grid[0] <= F0 or F0 >= last(f_grid)) {
225 throw std::runtime_error("Needs F0 in the range of f_grid");
226 }
227
228 auto red = Raw::Reduce::focus_doublescaling(f_grid, F0, DF);
229 y = Raw::Reduce::nanfocus(y, red);
230 f_grid = Raw::Reduce::nanfocus(f_grid, red);
231}
232
233void ybatchDoublingMeanFocus(Vector& f_grid,
234 ArrayOfVector& ybatch,
235 const Numeric& f0,
236 const Numeric& df,
237 const Verbosity&) {
238 for (Vector& y: ybatch) {
239 if (f_grid.nelem() not_eq y.nelem()) {
240 throw std::runtime_error("f_grid and all of ybatch must have the same size");
241 }
242 }
243
244 if (not is_increasing(f_grid)) {
245 throw std::runtime_error("f_grid must be sorted and ever increasing");
246 }
247
248 if (f_grid.nelem() < 2) {
249 throw std::runtime_error("Must have at least 2 frequency grid points");
250 }
251
252 // Sets defaults by method description
253 const Numeric F0 = f0 > 0 ? f0 : mean(f_grid);
254 const Numeric DF = df > 0 ? df : 10 * (f_grid[1] - f_grid[0]);
255
256 if (last(f_grid) <= F0 or f_grid[0] >= F0) {
257 throw std::runtime_error("Needs F0 in the range of f_grid");
258 }
259
260 auto red = Raw::Reduce::focus_doublescaling(f_grid, F0, DF);
261 for (Vector& y: ybatch) y = Raw::Reduce::nanfocus(y, red);
262 f_grid = Raw::Reduce::nanfocus(f_grid, red);
263}
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Definition artstime.cc:172
ArrayOfIndex time_steps(const ArrayOfTime &times, const TimeStep &DT)
Finds the index matching demands in a list of times.
Definition artstime.cc:59
TimeStep time_stepper_selection(const String &time_step)
Returns a time step from valid string.
Definition artstime.cc:19
Time mean_time(const ArrayOfTime &ts, Index s, Index E)
Computes the average time in a list.
Definition artstime.cc:133
TimeStep median(ArrayOfTimeStep dt)
Returns the median time step.
Definition artstime.cc:161
Stuff related to time in ARTS.
Array< Time > ArrayOfTime
List of times.
Definition artstime.h:122
This can be used to make arrays out of anything.
Definition array.h:31
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
void yDoublingMeanFocus(Vector &f_grid, Vector &y, const Numeric &f0, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: yDoublingMeanFocus.
Definition m_raw.cc:203
void yColdAtmHot(Vector &y, const Vector &cold, const Vector &atm, const Vector &hot, const Numeric &cold_temp, const Numeric &hot_temp, const Index &calib, const Verbosity &)
WORKSPACE METHOD: yColdAtmHot.
Definition m_raw.cc:14
void ybatchTimeAveraging(ArrayOfVector &ybatch, ArrayOfTime &sensor_time, const String &time_step, const Index &disregard_first, const Index &disregard_last, const Verbosity &)
WORKSPACE METHOD: ybatchTimeAveraging.
Definition m_raw.cc:65
void ybatchDoublingMeanFocus(Vector &f_grid, ArrayOfVector &ybatch, const Numeric &f0, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: ybatchDoublingMeanFocus.
Definition m_raw.cc:233
void ybatchColdAtmHotAtmCycle(ArrayOfVector &ybatch, ArrayOfTime &sensor_time, const ArrayOfVector &level0_data, const ArrayOfTime &level0_time, const Vector &cold_temp, const Vector &hot_temp, const Index &first_c_index, const Verbosity &)
WORKSPACE METHOD: ybatchColdAtmHotAtmCycle.
Definition m_raw.cc:35
void yMaskOutsideMedianRange(Vector &y, const Numeric &dx, const Verbosity &)
WORKSPACE METHOD: yMaskOutsideMedianRange.
Definition m_raw.cc:186
void ybatchTroposphericCorrectionNaiveMedianForward(ArrayOfVector &ybatch_corr, ArrayOfVector &ybatch, const ArrayOfIndex &range, const Vector &trop_temp, const Numeric &targ_temp, const Verbosity &)
WORKSPACE METHOD: ybatchTroposphericCorrectionNaiveMedianForward.
Definition m_raw.cc:129
void ybatchTroposphericCorrectionNaiveMedianInverse(ArrayOfVector &ybatch, const ArrayOfVector &ybatch_corr, const Verbosity &)
WORKSPACE METHOD: ybatchTroposphericCorrectionNaiveMedianInverse.
Definition m_raw.cc:169
void ybatchMaskOutsideMedianRange(ArrayOfVector &ybatch, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchMaskOutsideMedianRange.
Definition m_raw.cc:197
Numeric last(ConstVectorView x)
last
Numeric nanmedian(const ConstVectorView &v, const ArrayOfIndex &pos)
Get the median of the vector in the range ignoring all non-normal values.
Definition raw.cc:242
VectorView nanavg(VectorView y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the average of the ranged ys ignoring all non-normal values.
Definition raw.cc:112
ArrayOfVector caha(const ArrayOfVector &data, const Vector &tcvec, const Vector &thvec, const Index first_c_index)
Calibrate the data by CAHA.
Definition raw.cc:35
Vector systemtemp(const Vector &pc, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear receiver temperature formula elementwise.
Definition raw.cc:27
Vector calibration(const Vector &pc, const Vector &pa, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear calibration formula elementwise.
Definition raw.cc:19
std::vector< bool > out_of_bounds(const Vector &x, const Numeric xmin, const Numeric xmax)
Masks values that are out of bounds in x.
Definition raw.cc:383
Vector nanfocus(const Vector &x, const ArrayOfIndex &scaling)
Returns the re-averaged values of x according to scaling ignoring all non-normal values.
Definition raw.cc:341
ArrayOfIndex focus_doublescaling(const Vector &x, const Numeric x0, const Numeric dx)
Returns the relative position scale for each value in x.
Definition raw.cc:276
Stuff related to generating y-data from raw data.