ARTS 2.5.0 (git: 9ee3ac6c)
m_raw.cc
Go to the documentation of this file.
1/* Copyright (C) 2020
2 * Richard Larsson <ric.larsson@gmail.com>
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2, or (at your option) any
7 * later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 * USA. */
18
28#include "artstime.h"
29#include "raw.h"
30
31
33 const Vector& cold,
34 const Vector& atm,
35 const Vector& hot,
36 const Numeric& cold_temp,
37 const Numeric& hot_temp,
38 const Index& calib,
39 const Verbosity&)
40{
41 ARTS_USER_ERROR_IF (cold.nelem() not_eq atm.nelem() or atm.nelem() not_eq hot.nelem(),
42 "Length of vectors must be correct");
43
44 y.resize(atm.nelem());
45 if (calib) {
46 y = Raw::Calibration::calibration(cold, atm, hot, cold_temp, hot_temp);
47 } else {
48 y = Raw::Calibration::systemtemp(cold, hot, cold_temp, hot_temp);
49 }
50}
51
52
54 ArrayOfTime& sensor_time,
55 const ArrayOfVector& level0_data,
56 const ArrayOfTime& level0_time,
57 const Vector& cold_temp,
58 const Vector& hot_temp,
59 const Index& first_c_index,
60 const Verbosity&)
61{
62 ARTS_USER_ERROR_IF(level0_data.nelem() not_eq cold_temp.nelem() or
63 level0_data.nelem() not_eq hot_temp.nelem(),
64 "Length of vectors must be correct");
65 ARTS_USER_ERROR_IF (level0_time.nelem() not_eq level0_data.nelem() and
66 level0_time.nelem() not_eq 0,
67 "Bad level0_time length, must be empty of same as level0_data");
68
69 ybatch = Raw::Calibration::caha(level0_data, cold_temp, hot_temp, first_c_index);
70
71 // Fix time using the same method as CAHA
72 if (level0_time.nelem()) {
73 sensor_time.resize(ybatch.nelem());
74 // First position as described by method is at H if this index is too large
75 const Index pos = first_c_index - ((first_c_index > 1) ? 2 : 0);
76 for (Index i=0; i<sensor_time.nelem(); i++) {
77 sensor_time[i] = level0_time[pos + 2*i];
78 }
79 }
80}
81
82
84 ArrayOfTime& sensor_time,
85 const String& time_step,
86 const Index& disregard_first,
87 const Index& disregard_last,
88 const Verbosity&)
89{
90 // Size of problem
91 const Index n=sensor_time.nelem();
92 ARTS_USER_ERROR_IF (sensor_time.nelem() not_eq n,
93 "Time vector length must match input data length");
94
95 // Time is not decreasing
96 ARTS_USER_ERROR_IF (not std::is_sorted(sensor_time.cbegin(), sensor_time.cend()),
97 "Time vector cannot decrease");
98
99 // Find the limits of the range
100 auto lims = time_steps(sensor_time, time_stepper_selection(time_step));
101
102 // Output variables
103 ArrayOfVector ybatch_out;
104 ArrayOfTime sensor_time_out;
105
106 if (lims.front() == n) {
107 ybatch_out.resize(0);
108 sensor_time_out.resize(0);
109 } else {
110
111 // Frequency grids
112 const Index k = ybatch[0].nelem();
113 ARTS_USER_ERROR_IF (not std::all_of(ybatch.cbegin(), ybatch.cend(),
114 [k](auto& x){return x.nelem() == k;}),
115 "Bad frequency grid size in input data; expects all equal");
116
117 // Allocate output
118 const Index m = lims.nelem() - 1;
119 ARTS_USER_ERROR_IF (m < 0,
120 "Must include last if time step covers all of the range");
121 ybatch_out = ArrayOfVector(m, Vector(k));
122 sensor_time_out = ArrayOfTime(m);
123
124 // Fill output
125 #pragma omp parallel for if (not arts_omp_in_parallel()) schedule(guided)
126 for (Index i=0; i<m; i++) {
127 sensor_time_out[i] = mean_time(sensor_time, lims[i], lims[i+1]);
128 Raw::Average::nanavg(ybatch_out[i], ybatch, lims[i], lims[i+1]);
129 }
130 }
131
132 if (disregard_first) {
133 ybatch_out.erase(ybatch_out.begin());
134 sensor_time_out.erase(sensor_time_out.begin());
135 }
136
137 if (disregard_last) {
138 ybatch_out.pop_back();
139 sensor_time_out.pop_back();
140 }
141
142 ybatch = ybatch_out;
143 sensor_time = sensor_time_out;
144}
145
146
148 ArrayOfVector& ybatch,
149 const ArrayOfIndex& range,
150 const Vector& trop_temp,
151 const Numeric& targ_temp,
152 const Verbosity&)
153{
154 // Size of problem
155 const Index n=ybatch.nelem();
156 const Index m=n ? ybatch[0].nelem() : 0;
157
158 ARTS_USER_ERROR_IF (m == 0, "A frequency range is required")
159
160 ARTS_USER_ERROR_IF (std::any_of(ybatch.begin(), ybatch.end(),
161 [m](auto& y){return y.nelem() not_eq m;}),
162 "Bad input size, all of ybatch must match itself");
163
164 ARTS_USER_ERROR_IF (trop_temp.nelem() not_eq n and trop_temp.nelem() not_eq 1,
165 "Bad input size, trop_temp must match ybatch or be 1-long,\n"
166 "trop_temp: [", trop_temp, "]\ntrop_temp.nelem(): ",
167 trop_temp.nelem(), "\nybatch.nelem(): ", n);
168
169 // This algorithm stores partial transmission and median and tropospheric temperature in the correction terms
170 ybatch_corr = ArrayOfVector(n, Vector(3));
171
172 // Compute tropospheric correction
173 for (Index i=0; i<n; i++) {
174 ybatch_corr[i][2] = trop_temp.nelem() > 1 ? trop_temp[i] : trop_temp[0];
175 ybatch_corr[i][0] = Raw::Average::nanmedian(ybatch[i], range);
176 ybatch_corr[i][1] = std::exp(- std::log((ybatch_corr[i][2] - ybatch_corr[i][0]) / (ybatch_corr[i][2] - targ_temp)));
177 }
178
179 // Apply correction
180 for (Index i=0; i<n; i++) {
181 ybatch[i] *= ybatch_corr[i][1];
182 ybatch[i] += ybatch_corr[i][2] * (1 - ybatch_corr[i][1]);
183 }
184}
185
186
188 const ArrayOfVector& ybatch_corr,
189 const Verbosity&)
190{
191 // Size of problem
192 const Index n=ybatch.nelem();
193 ARTS_USER_ERROR_IF ((std::any_of(ybatch_corr.begin(), ybatch_corr.end(),
194 [](auto& corr){return corr.nelem() not_eq 3;})) or ybatch_corr.nelem() not_eq n,
195 "Bad input size, all of ybatch_corr must match ybatch and have three elements each");
196
197 // Apply inverse of correction
198 for (Index i=0; i<n; i++) {
199 ybatch[i] -= ybatch_corr[i][2] * (1 - ybatch_corr[i][1]);
200 ybatch[i] /= ybatch_corr[i][1];
201 }
202}
203
205 const Numeric& dx,
206 const Verbosity&) {
209
210 for (Index i=0; i<y.nelem(); i++) {
211 y[i] = mask[i] ? std::numeric_limits<Numeric>::quiet_NaN() : y[i];
212 }
213}
214
216 const Numeric& dx,
217 const Verbosity& verbosity) {
218 for (Vector& y: ybatch) yMaskOutsideMedianRange(y, dx, verbosity);
219}
220
222 Vector& y,
223 const Numeric& f0,
224 const Numeric& df,
225 const Verbosity&) {
226 if (f_grid.nelem() not_eq y.nelem()) {
227 throw std::runtime_error("f_grid and y must have the same size");
228 }
229
230 if (not is_increasing(f_grid)) {
231 throw std::runtime_error("f_grid must be sorted and ever increasing");
232 }
233
234 if (f_grid.nelem() < 2) {
235 throw std::runtime_error("Must have at least 2 frequency grid points");
236 }
237
238 // Sets defaults by method description
239 const Numeric F0 = f0 > 0 ? f0 : mean(f_grid);
240 const Numeric DF = df > 0 ? df : 10 * (f_grid[1] - f_grid[0]);
241
242 if (f_grid[0] <= F0 or F0 >= last(f_grid)) {
243 throw std::runtime_error("Needs F0 in the range of f_grid");
244 }
245
246 auto red = Raw::Reduce::focus_doublescaling(f_grid, F0, DF);
247 y = Raw::Reduce::nanfocus(y, red);
248 f_grid = Raw::Reduce::nanfocus(f_grid, red);
249}
250
252 ArrayOfVector& ybatch,
253 const Numeric& f0,
254 const Numeric& df,
255 const Verbosity&) {
256 for (Vector& y: ybatch) {
257 if (f_grid.nelem() not_eq y.nelem()) {
258 throw std::runtime_error("f_grid and all of ybatch must have the same size");
259 }
260 }
261
262 if (not is_increasing(f_grid)) {
263 throw std::runtime_error("f_grid must be sorted and ever increasing");
264 }
265
266 if (f_grid.nelem() < 2) {
267 throw std::runtime_error("Must have at least 2 frequency grid points");
268 }
269
270 // Sets defaults by method description
271 const Numeric F0 = f0 > 0 ? f0 : mean(f_grid);
272 const Numeric DF = df > 0 ? df : 10 * (f_grid[1] - f_grid[0]);
273
274 if (last(f_grid) <= F0 or f_grid[0] >= F0) {
275 throw std::runtime_error("Needs F0 in the range of f_grid");
276 }
277
278 auto red = Raw::Reduce::focus_doublescaling(f_grid, F0, DF);
279 for (Vector& y: ybatch) y = Raw::Reduce::nanfocus(y, red);
280 f_grid = Raw::Reduce::nanfocus(f_grid, red);
281}
Array< Vector > ArrayOfVector
An array of vectors.
Definition: array.h:57
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Definition: artstime.cc:190
ArrayOfIndex time_steps(const ArrayOfTime &times, const TimeStep &DT)
Finds the index matching demands in a list of times.
Definition: artstime.cc:77
TimeStep time_stepper_selection(const String &time_step)
Returns a time step from valid string.
Definition: artstime.cc:37
Time mean_time(const ArrayOfTime &ts, Index s, Index E)
Computes the average time in a list.
Definition: artstime.cc:151
TimeStep median(ArrayOfTimeStep dt)
Returns the median time step.
Definition: artstime.cc:179
Stuff related to time in ARTS.
Array< Time > ArrayOfTime
List of times.
Definition: artstime.h:90
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define dx
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
bool is_sorted(ConstVectorView x)
Checks if a vector is sorted in ascending order.
Definition: logic.cc:199
void yDoublingMeanFocus(Vector &f_grid, Vector &y, const Numeric &f0, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: yDoublingMeanFocus.
Definition: m_raw.cc:221
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:32
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:83
void ybatchDoublingMeanFocus(Vector &f_grid, ArrayOfVector &ybatch, const Numeric &f0, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: ybatchDoublingMeanFocus.
Definition: m_raw.cc:251
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:53
void yMaskOutsideMedianRange(Vector &y, const Numeric &dx, const Verbosity &)
WORKSPACE METHOD: yMaskOutsideMedianRange.
Definition: m_raw.cc:204
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:147
void ybatchTroposphericCorrectionNaiveMedianInverse(ArrayOfVector &ybatch, const ArrayOfVector &ybatch_corr, const Verbosity &)
WORKSPACE METHOD: ybatchTroposphericCorrectionNaiveMedianInverse.
Definition: m_raw.cc:187
void ybatchMaskOutsideMedianRange(ArrayOfVector &ybatch, const Numeric &dx, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchMaskOutsideMedianRange.
Definition: m_raw.cc:215
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:159
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
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:260
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:130
ArrayOfVector caha(const ArrayOfVector &data, const Vector &tcvec, const Vector &thvec, const Index first_c_index)
Calibrate the data by CAHA.
Definition: raw.cc:53
Vector systemtemp(const Vector &pc, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear receiver temperature formula elementwise.
Definition: raw.cc:45
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:37
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:401
VectorView mask(VectorView x, const std::vector< bool > &masking)
Masks all true entries of masking in x by turning them into NaN.
Definition: raw.cc:408
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:359
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:294
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
Stuff related to generating y-data from raw data.