Line data Source code
1 : /*******************************************************************************
2 : * ALMA - Atacama Large Millimiter Array
3 : * (c) Instituto de Estructura de la Materia, 2009
4 : *
5 : * This library is free software; you can redistribute it and/or
6 : * modify it under the terms of the GNU Lesser General Public
7 : * License as published by the Free Software Foundation; either
8 : * version 2.1 of the License, or (at your option) any later version.
9 : *
10 : * This library is distributed in the hope that it will be useful,
11 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 : * Lesser General Public License for more details.
14 : *
15 : * You should have received a copy of the GNU Lesser General Public
16 : * License along with this library; if not, write to the Free Software
17 : * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 : *
19 : * "@(#) $Id: ATMSkyStatus.cpp Exp $"
20 : *
21 : * who when what
22 : * -------- -------- ----------------------------------------------
23 : * pardo 24/03/09 created
24 : */
25 :
26 : #include "ATMSkyStatus.h"
27 :
28 : #include <iostream>
29 : #include <math.h>
30 :
31 : ATM_NAMESPACE_BEGIN
32 :
33 : // Constructors
34 :
35 18 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile) :
36 18 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(1.0),
37 18 : skyBackgroundTemperature_(2.73, Temperature::UnitKelvin)
38 : {
39 :
40 18 : iniSkyStatus();
41 :
42 18 : }
43 :
44 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
45 0 : double airMass) :
46 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
47 0 : skyBackgroundTemperature_(2.73, Temperature::UnitKelvin)
48 : {
49 :
50 0 : iniSkyStatus();
51 :
52 0 : }
53 :
54 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
55 0 : const Temperature &temperatureBackground) :
56 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(1.0),
57 0 : skyBackgroundTemperature_(temperatureBackground)
58 : {
59 :
60 0 : iniSkyStatus();
61 :
62 0 : }
63 :
64 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
65 0 : const Length &wh2o) :
66 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(1.0),
67 0 : skyBackgroundTemperature_(2.73, Temperature::UnitKelvin), wh2o_user_(wh2o)
68 : {
69 :
70 0 : iniSkyStatus();
71 :
72 0 : }
73 :
74 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
75 : const Temperature &temperatureBackground,
76 0 : double airMass) :
77 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
78 0 : skyBackgroundTemperature_(temperatureBackground)
79 : {
80 :
81 0 : iniSkyStatus();
82 :
83 0 : }
84 :
85 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
86 : double airMass,
87 0 : const Temperature &temperatureBackground) :
88 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
89 0 : skyBackgroundTemperature_(temperatureBackground)
90 : {
91 :
92 0 : iniSkyStatus();
93 :
94 0 : }
95 :
96 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
97 : const Length &wh2o,
98 0 : double airMass) :
99 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
100 0 : skyBackgroundTemperature_(2.73, Temperature::UnitKelvin), wh2o_user_(wh2o)
101 : {
102 :
103 0 : iniSkyStatus();
104 :
105 0 : }
106 :
107 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
108 : double airMass,
109 0 : const Length &wh2o) :
110 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
111 0 : skyBackgroundTemperature_(2.73, Temperature::UnitKelvin), wh2o_user_(wh2o)
112 : {
113 :
114 0 : iniSkyStatus();
115 :
116 0 : }
117 :
118 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
119 : const Length &wh2o,
120 0 : const Temperature &temperatureBackground) :
121 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(1.0),
122 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
123 : {
124 :
125 0 : iniSkyStatus();
126 :
127 0 : }
128 :
129 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
130 : const Temperature &temperatureBackground,
131 0 : const Length &wh2o) :
132 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(1.0),
133 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
134 : {
135 :
136 0 : iniSkyStatus();
137 :
138 0 : }
139 :
140 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
141 : const Temperature &temperatureBackground,
142 : double airMass,
143 0 : const Length &wh2o) :
144 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
145 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
146 : {
147 :
148 0 : iniSkyStatus();
149 :
150 0 : }
151 :
152 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
153 : const Temperature &temperatureBackground,
154 : const Length &wh2o,
155 0 : double airMass) :
156 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
157 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
158 : {
159 :
160 0 : iniSkyStatus();
161 :
162 0 : }
163 :
164 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
165 : double airMass,
166 : const Temperature &temperatureBackground,
167 0 : const Length &wh2o) :
168 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
169 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
170 : {
171 :
172 0 : iniSkyStatus();
173 :
174 0 : }
175 :
176 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
177 : double airMass,
178 : const Length &wh2o,
179 0 : const Temperature &temperatureBackground) :
180 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
181 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
182 : {
183 :
184 0 : iniSkyStatus();
185 :
186 0 : }
187 :
188 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
189 : const Length &wh2o,
190 : const Temperature &temperatureBackground,
191 0 : double airMass) :
192 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
193 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
194 : {
195 :
196 0 : iniSkyStatus();
197 :
198 0 : }
199 :
200 0 : SkyStatus::SkyStatus(const RefractiveIndexProfile &refractiveIndexProfile,
201 : const Length &wh2o,
202 : double airMass,
203 0 : const Temperature &temperatureBackground) :
204 0 : RefractiveIndexProfile(refractiveIndexProfile), airMass_(airMass),
205 0 : skyBackgroundTemperature_(temperatureBackground), wh2o_user_(wh2o)
206 : {
207 :
208 0 : iniSkyStatus();
209 :
210 0 : }
211 :
212 0 : SkyStatus::SkyStatus(const SkyStatus & a) : RefractiveIndexProfile(a)
213 : {
214 :
215 : // 2015-03-02 (DB) :
216 : // Use constructor of base class RefractiveIndexProfile()
217 : // and comment the following code because it is executed in RefractiveIndexProfile() constructor
218 :
219 : // groundTemperature_ = a.groundTemperature_;
220 : // tropoLapseRate_ = a.tropoLapseRate_;
221 : // groundPressure_ = a.groundPressure_;
222 : // relativeHumidity_ = a.relativeHumidity_;
223 : // wvScaleHeight_ = a.wvScaleHeight_;
224 : // pressureStep_ = a.pressureStep_;
225 : // pressureStepFactor_ = a.pressureStepFactor_;
226 : // altitude_ = a.altitude_;
227 : // topAtmProfile_ = a.topAtmProfile_;
228 : // numLayer_ = a.numLayer_;
229 : // newBasicParam_ = a.newBasicParam_;
230 : // v_layerThickness_.reserve(numLayer_);
231 : // v_layerTemperature_.reserve(numLayer_);
232 : // v_layerWaterVapor_.reserve(numLayer_);
233 : // v_layerCO_.reserve(numLayer_);
234 : // v_layerO3_.reserve(numLayer_);
235 : // v_layerN2O_.reserve(numLayer_);
236 : // v_layerNO2_.reserve(numLayer_);
237 : // v_layerSO2_.reserve(numLayer_);
238 : // for(unsigned int n = 0; n < numLayer_; n++) {
239 : // v_layerThickness_.push_back(a.v_layerThickness_[n]);
240 : // v_layerTemperature_.push_back(a.v_layerTemperature_[n]);
241 : // //v_layerDeltaT_.push_back(a.v_layerDeltaT_[n]);
242 : // //cout << "n=" << n << endl;
243 : // v_layerWaterVapor_.push_back(a.v_layerWaterVapor_[n]);
244 : // v_layerPressure_.push_back(a.v_layerPressure_[n]);
245 : // v_layerCO_.push_back(a.v_layerCO_[n]);
246 : // v_layerO3_.push_back(a.v_layerO3_[n]);
247 : // v_layerN2O_.push_back(a.v_layerN2O_[n]);
248 : // v_layerNO2_.push_back(a.v_layerNO2_[n]);
249 : // v_layerSO2_.push_back(a.v_layerSO2_[n]);
250 : // }
251 :
252 : // // level Spectral Grid
253 : // freqUnits_ = a.freqUnits_;
254 : // v_chanFreq_ = a.v_chanFreq_;
255 :
256 : // v_numChan_ = a.v_numChan_;
257 : // v_refChan_ = a.v_refChan_;
258 : // v_refFreq_ = a.v_refFreq_;
259 : // v_chanSep_ = a.v_chanSep_;
260 : // v_maxFreq_ = a.v_maxFreq_;
261 : // v_minFreq_ = a.v_minFreq_;
262 : // v_intermediateFrequency_ = a.v_intermediateFrequency_;
263 : // v_loFreq_ = a.v_loFreq_;
264 :
265 : // v_sidebandSide_ = a.v_sidebandSide_;
266 : // v_sidebandType_ = a.v_sidebandType_;
267 :
268 : // vv_assocSpwId_ = a.vv_assocSpwId_;
269 : // vv_assocNature_ = a.vv_assocNature_;
270 :
271 : // v_transfertId_ = a.v_transfertId_;
272 :
273 : // // level Absorption Profile
274 : // vv_N_H2OLinesPtr_.reserve(a.v_chanFreq_.size());
275 : // vv_N_H2OContPtr_.reserve(a.v_chanFreq_.size());
276 : // vv_N_O2LinesPtr_.reserve(a.v_chanFreq_.size());
277 : // vv_N_DryContPtr_.reserve(a.v_chanFreq_.size());
278 : // vv_N_O3LinesPtr_.reserve(a.v_chanFreq_.size());
279 : // vv_N_COLinesPtr_.reserve(a.v_chanFreq_.size());
280 : // vv_N_N2OLinesPtr_.reserve(a.v_chanFreq_.size());
281 : // vv_N_NO2LinesPtr_.reserve(a.v_chanFreq_.size());
282 : // vv_N_SO2LinesPtr_.reserve(a.v_chanFreq_.size());
283 :
284 : // vector<complex<double> >* v_N_H2OLinesPtr;
285 : // vector<complex<double> >* v_N_H2OContPtr;
286 : // vector<complex<double> >* v_N_O2LinesPtr;
287 : // vector<complex<double> >* v_N_DryContPtr;
288 : // vector<complex<double> >* v_N_O3LinesPtr;
289 : // vector<complex<double> >* v_N_COLinesPtr;
290 : // vector<complex<double> >* v_N_N2OLinesPtr;
291 : // vector<complex<double> >* v_N_NO2LinesPtr;
292 : // vector<complex<double> >* v_N_SO2LinesPtr;
293 :
294 : // for(unsigned int nc = 0; nc < v_chanFreq_.size(); nc++) {
295 :
296 : // v_N_H2OLinesPtr = new vector<complex<double> > ;
297 : // v_N_H2OLinesPtr->reserve(numLayer_);
298 : // v_N_H2OContPtr = new vector<complex<double> > ;
299 : // v_N_H2OContPtr->reserve(numLayer_);
300 : // v_N_O2LinesPtr = new vector<complex<double> > ;
301 : // v_N_O2LinesPtr->reserve(numLayer_);
302 : // v_N_DryContPtr = new vector<complex<double> > ;
303 : // v_N_DryContPtr->reserve(numLayer_);
304 : // v_N_O3LinesPtr = new vector<complex<double> > ;
305 : // v_N_O3LinesPtr->reserve(numLayer_);
306 : // v_N_COLinesPtr = new vector<complex<double> > ;
307 : // v_N_COLinesPtr->reserve(numLayer_);
308 : // v_N_N2OLinesPtr = new vector<complex<double> > ;
309 : // v_N_N2OLinesPtr->reserve(numLayer_);
310 : // v_N_NO2LinesPtr = new vector<complex<double> > ;
311 : // v_N_NO2LinesPtr->reserve(numLayer_);
312 : // v_N_SO2LinesPtr = new vector<complex<double> > ;
313 : // v_N_SO2LinesPtr->reserve(numLayer_);
314 :
315 : // for(unsigned int n = 0; n < numLayer_; n++) {
316 :
317 : // // cout << "numLayer_=" << nc << " " << n << endl; // COMMENTED OUT BY JUAN MAY/16/2005
318 :
319 : // v_N_H2OLinesPtr->push_back(a.vv_N_H2OLinesPtr_[nc]->at(n));
320 : // v_N_H2OContPtr->push_back(a.vv_N_H2OContPtr_[nc]->at(n));
321 : // v_N_O2LinesPtr->push_back(a.vv_N_O2LinesPtr_[nc]->at(n));
322 : // v_N_DryContPtr->push_back(a.vv_N_DryContPtr_[nc]->at(n));
323 : // v_N_O3LinesPtr->push_back(a.vv_N_O3LinesPtr_[nc]->at(n));
324 : // v_N_COLinesPtr->push_back(a.vv_N_COLinesPtr_[nc]->at(n));
325 : // v_N_N2OLinesPtr->push_back(a.vv_N_N2OLinesPtr_[nc]->at(n));
326 : // v_N_NO2LinesPtr->push_back(a.vv_N_NO2LinesPtr_[nc]->at(n));
327 : // v_N_SO2LinesPtr->push_back(a.vv_N_SO2LinesPtr_[nc]->at(n));
328 :
329 : // }
330 :
331 : // vv_N_H2OLinesPtr_.push_back(v_N_H2OLinesPtr);
332 : // vv_N_H2OContPtr_.push_back(v_N_H2OContPtr);
333 : // vv_N_O2LinesPtr_.push_back(v_N_O2LinesPtr);
334 : // vv_N_DryContPtr_.push_back(v_N_DryContPtr);
335 : // vv_N_O3LinesPtr_.push_back(v_N_O3LinesPtr);
336 : // vv_N_COLinesPtr_.push_back(v_N_COLinesPtr);
337 : // vv_N_N2OLinesPtr_.push_back(v_N_N2OLinesPtr);
338 : // vv_N_NO2LinesPtr_.push_back(v_N_NO2LinesPtr);
339 : // vv_N_SO2LinesPtr_.push_back(v_N_SO2LinesPtr);
340 :
341 : // }
342 :
343 : // level Atm Radiance
344 :
345 0 : airMass_ = a.airMass_;
346 0 : skyBackgroundTemperature_ = a.skyBackgroundTemperature_;
347 0 : wh2o_user_ = a.wh2o_user_;
348 :
349 0 : }
350 :
351 33 : SkyStatus::~SkyStatus()
352 : {
353 18 : rmSkyStatus();
354 33 : }
355 :
356 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude,
357 : const Pressure &groundPressure,
358 : const Temperature &groundTemperature,
359 : double tropoLapseRate,
360 : const Humidity &relativeHumidity,
361 : const Length &wvScaleHeight)
362 : {
363 0 : bool update = updateProfilesAndRadiance(altitude,
364 : groundPressure,
365 : groundTemperature,
366 : tropoLapseRate,
367 : relativeHumidity,
368 : wvScaleHeight);
369 0 : return update;
370 : }
371 :
372 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude,
373 : const Length &wvScaleHeight)
374 : {
375 0 : bool update = updateProfilesAndRadiance(altitude,
376 0 : groundPressure_,
377 0 : groundTemperature_,
378 : tropoLapseRate_,
379 0 : relativeHumidity_,
380 : wvScaleHeight);
381 0 : return update;
382 : }
383 :
384 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude)
385 : {
386 0 : bool update = updateProfilesAndRadiance(altitude,
387 0 : groundPressure_,
388 0 : groundTemperature_,
389 : tropoLapseRate_,
390 0 : relativeHumidity_,
391 0 : wvScaleHeight_);
392 0 : return update;
393 : }
394 :
395 0 : bool SkyStatus::setBasicAtmosphericParameters(const Temperature &groundTemperature)
396 : {
397 0 : bool update = updateProfilesAndRadiance(altitude_,
398 0 : groundPressure_,
399 : groundTemperature,
400 : tropoLapseRate_,
401 0 : relativeHumidity_,
402 0 : wvScaleHeight_);
403 0 : return update;
404 : }
405 :
406 120 : bool SkyStatus::setBasicAtmosphericParameters(const Pressure &groundPressure)
407 : {
408 :
409 240 : bool update = updateProfilesAndRadiance(altitude_,
410 : groundPressure,
411 120 : groundTemperature_,
412 : tropoLapseRate_,
413 120 : relativeHumidity_,
414 120 : wvScaleHeight_);
415 :
416 120 : return update;
417 : }
418 :
419 0 : bool SkyStatus::setBasicAtmosphericParameters(const Humidity &relativeHumidity)
420 : {
421 0 : bool update = updateProfilesAndRadiance(altitude_,
422 0 : groundPressure_,
423 0 : groundTemperature_,
424 : tropoLapseRate_,
425 : relativeHumidity,
426 0 : wvScaleHeight_);
427 0 : return update;
428 : }
429 :
430 0 : bool SkyStatus::setBasicAtmosphericParameters(double tropoLapseRate)
431 : {
432 0 : bool update = updateProfilesAndRadiance(altitude_,
433 0 : groundPressure_,
434 0 : groundTemperature_,
435 : tropoLapseRate,
436 0 : relativeHumidity_,
437 0 : wvScaleHeight_);
438 0 : return update;
439 : }
440 :
441 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude,
442 : const Temperature &groundTemperature)
443 : {
444 0 : bool update = updateProfilesAndRadiance(altitude,
445 0 : groundPressure_,
446 : groundTemperature,
447 : tropoLapseRate_,
448 0 : relativeHumidity_,
449 0 : wvScaleHeight_);
450 0 : return update;
451 : }
452 :
453 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude,
454 : const Pressure &groundPressure)
455 : {
456 0 : bool update = updateProfilesAndRadiance(altitude,
457 : groundPressure,
458 0 : groundTemperature_,
459 : tropoLapseRate_,
460 0 : relativeHumidity_,
461 0 : wvScaleHeight_);
462 0 : return update;
463 : }
464 :
465 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude,
466 : const Humidity &relativeHumidity)
467 : {
468 0 : bool update = updateProfilesAndRadiance(altitude,
469 0 : groundPressure_,
470 0 : groundTemperature_,
471 : tropoLapseRate_,
472 : relativeHumidity,
473 0 : wvScaleHeight_);
474 0 : return update;
475 : }
476 :
477 0 : bool SkyStatus::setBasicAtmosphericParameters(const Length &altitude,
478 : double tropoLapseRate)
479 : {
480 0 : bool update = updateProfilesAndRadiance(altitude,
481 0 : groundPressure_,
482 0 : groundTemperature_,
483 : tropoLapseRate,
484 0 : relativeHumidity_,
485 0 : wvScaleHeight_);
486 0 : return update;
487 : }
488 :
489 0 : bool SkyStatus::setBasicAtmosphericParameters(const Temperature &groundTemperature,
490 : const Pressure &groundPressure)
491 : {
492 0 : bool update = updateProfilesAndRadiance(altitude_,
493 : groundPressure,
494 : groundTemperature,
495 : tropoLapseRate_,
496 0 : relativeHumidity_,
497 0 : wvScaleHeight_);
498 0 : return update;
499 : }
500 0 : bool SkyStatus::setBasicAtmosphericParameters(const Pressure &groundPressure,
501 : const Temperature &groundTemperature)
502 : {
503 0 : bool update = updateProfilesAndRadiance(altitude_,
504 : groundPressure,
505 : groundTemperature,
506 : tropoLapseRate_,
507 0 : relativeHumidity_,
508 0 : wvScaleHeight_);
509 0 : return update;
510 : }
511 :
512 0 : bool SkyStatus::setBasicAtmosphericParameters(const Temperature &groundTemperature,
513 : const Humidity &relativeHumidity)
514 : {
515 0 : bool update = updateProfilesAndRadiance(altitude_,
516 0 : groundPressure_,
517 : groundTemperature,
518 : tropoLapseRate_,
519 : relativeHumidity,
520 0 : wvScaleHeight_);
521 0 : return update;
522 : }
523 0 : bool SkyStatus::setBasicAtmosphericParameters(const Humidity &relativeHumidity,
524 : const Temperature &groundTemperature)
525 : {
526 0 : bool update = updateProfilesAndRadiance(altitude_,
527 0 : groundPressure_,
528 : groundTemperature,
529 : tropoLapseRate_,
530 : relativeHumidity,
531 0 : wvScaleHeight_);
532 0 : return update;
533 : }
534 :
535 0 : bool SkyStatus::setBasicAtmosphericParameters(const Temperature &groundTemperature,
536 : double tropoLapseRate)
537 : {
538 0 : bool update = updateProfilesAndRadiance(altitude_,
539 0 : groundPressure_,
540 : groundTemperature,
541 : tropoLapseRate,
542 0 : relativeHumidity_,
543 0 : wvScaleHeight_);
544 0 : return update;
545 : }
546 0 : bool SkyStatus::setBasicAtmosphericParameters(double tropoLapseRate,
547 : const Temperature &groundTemperature)
548 : {
549 0 : bool update = updateProfilesAndRadiance(altitude_,
550 0 : groundPressure_,
551 : groundTemperature,
552 : tropoLapseRate,
553 0 : relativeHumidity_,
554 0 : wvScaleHeight_);
555 0 : return update;
556 : }
557 :
558 0 : bool SkyStatus::setBasicAtmosphericParameters(const Temperature &groundTemperature,
559 : const Length &wvScaleHeight)
560 : {
561 0 : bool update = updateProfilesAndRadiance(altitude_,
562 0 : groundPressure_,
563 : groundTemperature,
564 : tropoLapseRate_,
565 0 : relativeHumidity_,
566 : wvScaleHeight);
567 0 : return update;
568 : }
569 :
570 0 : bool SkyStatus::setBasicAtmosphericParameters(const Pressure &groundPressure,
571 : const Humidity &relativeHumidity)
572 : {
573 0 : bool update = updateProfilesAndRadiance(altitude_,
574 : groundPressure,
575 0 : groundTemperature_,
576 : tropoLapseRate_,
577 : relativeHumidity,
578 0 : wvScaleHeight_);
579 0 : return update;
580 : }
581 0 : bool SkyStatus::setBasicAtmosphericParameters(const Humidity &relativeHumidity,
582 : const Pressure &groundPressure)
583 : {
584 0 : bool update = updateProfilesAndRadiance(altitude_,
585 : groundPressure,
586 0 : groundTemperature_,
587 : tropoLapseRate_,
588 : relativeHumidity,
589 0 : wvScaleHeight_);
590 0 : return update;
591 : }
592 :
593 0 : bool SkyStatus::setBasicAtmosphericParameters(const Pressure &groundPressure,
594 : double tropoLapseRate)
595 : {
596 0 : bool update = updateProfilesAndRadiance(altitude_,
597 : groundPressure,
598 0 : groundTemperature_,
599 : tropoLapseRate,
600 0 : relativeHumidity_,
601 0 : wvScaleHeight_);
602 0 : return update;
603 : }
604 0 : bool SkyStatus::setBasicAtmosphericParameters(double tropoLapseRate,
605 : const Pressure &groundPressure)
606 : {
607 0 : bool update = updateProfilesAndRadiance(altitude_,
608 : groundPressure,
609 0 : groundTemperature_,
610 : tropoLapseRate,
611 0 : relativeHumidity_,
612 0 : wvScaleHeight_);
613 0 : return update;
614 : }
615 :
616 0 : bool SkyStatus::setBasicAtmosphericParameters(const Pressure &groundPressure,
617 : const Length &wvScaleHeight)
618 : {
619 0 : bool update = updateProfilesAndRadiance(altitude_,
620 : groundPressure,
621 0 : groundTemperature_,
622 : tropoLapseRate_,
623 0 : relativeHumidity_,
624 : wvScaleHeight);
625 0 : return update;
626 : }
627 :
628 0 : bool SkyStatus::setBasicAtmosphericParameters(const Humidity &relativeHumidity,
629 : double tropoLapseRate)
630 : {
631 0 : bool update = updateProfilesAndRadiance(altitude_,
632 0 : groundPressure_,
633 0 : groundTemperature_,
634 : tropoLapseRate,
635 : relativeHumidity,
636 0 : wvScaleHeight_);
637 0 : return update;
638 : }
639 0 : bool SkyStatus::setBasicAtmosphericParameters(double tropoLapseRate,
640 : const Humidity &relativeHumidity)
641 : {
642 0 : bool update = updateProfilesAndRadiance(altitude_,
643 0 : groundPressure_,
644 0 : groundTemperature_,
645 : tropoLapseRate,
646 : relativeHumidity,
647 0 : wvScaleHeight_);
648 0 : return update;
649 : }
650 :
651 0 : bool SkyStatus::setBasicAtmosphericParameters(const Humidity &relativeHumidity,
652 : const Length &wvScaleHeight)
653 : {
654 0 : bool update = updateProfilesAndRadiance(altitude_,
655 0 : groundPressure_,
656 0 : groundTemperature_,
657 : tropoLapseRate_,
658 : relativeHumidity,
659 : wvScaleHeight);
660 0 : return update;
661 : }
662 :
663 0 : bool SkyStatus::setBasicAtmosphericParameters(double tropoLapseRate,
664 : const Length &wvScaleHeight)
665 : {
666 0 : bool update = updateProfilesAndRadiance(altitude_,
667 0 : groundPressure_,
668 0 : groundTemperature_,
669 : tropoLapseRate,
670 0 : relativeHumidity_,
671 : wvScaleHeight);
672 0 : return update;
673 : }
674 :
675 120 : bool SkyStatus::updateProfilesAndRadiance(const Length &altitude,
676 : const Pressure &groundPressure,
677 : const Temperature &groundTemperature,
678 : double tropoLapseRate,
679 : const Humidity &relativeHumidity,
680 : const Length &wvScaleHeight)
681 : {
682 :
683 120 : bool updated = false;
684 :
685 120 : bool mkNewAbsPhaseProfile = updateRefractiveIndexProfile(altitude,
686 : groundPressure,
687 : groundTemperature,
688 : tropoLapseRate,
689 : relativeHumidity,
690 : wvScaleHeight);
691 :
692 120 : if(mkNewAbsPhaseProfile) {
693 120 : updated = true;
694 : }
695 120 : return updated;
696 : }
697 :
698 26 : Opacity SkyStatus::getH2OLinesOpacity(unsigned int nc)
699 : {
700 26 : if(!chanIndexIsValid(nc)) return (double) -999.0;
701 26 : double kv = 0;
702 656 : for(unsigned int j = 0; j < numLayer_; j++) {
703 630 : kv = kv + imag(vv_N_H2OLinesPtr_[nc]->at(j)) * v_layerThickness_[j];
704 : }
705 52 : return ((getUserWH2O().get()) / (getGroundWH2O().get())) * kv;
706 : }
707 :
708 0 : Opacity SkyStatus::getH2OLinesOpacityUpTo(unsigned int nc, Length refalti)
709 : {
710 0 : unsigned int ires; unsigned int numlayerold; Length alti;
711 0 : Opacity opacityout0; Opacity opacityout1; Opacity zeroOp(0.0,Opacity::UnitNeper);
712 : double fractionLast; double g1; double g2;
713 :
714 0 : if(refalti.get(Length::UnitKiloMeter) <= altitude_.get(Length::UnitKiloMeter)) {
715 0 : return zeroOp;
716 : }else{
717 0 : fractionLast = 1.0; numlayerold = numLayer_;
718 0 : g1=getGroundWH2O().get();
719 0 : opacityout0=getH2OLinesOpacity(nc); ires=numlayerold-1; alti=altitude_;
720 0 : for(unsigned int i=0; i<numLayer_; i++){
721 0 : if(alti.get(Length::UnitKiloMeter) < refalti.get(Length::UnitKiloMeter) && (alti.get(Length::UnitKiloMeter)+v_layerThickness_[i]/1000.0) >= refalti.get(Length::UnitKiloMeter))
722 0 : { ires=i; fractionLast = (refalti.get(Length::UnitMeter)-alti.get(Length::UnitMeter))/v_layerThickness_[i]; }
723 0 : alti = alti + Length(v_layerThickness_[i],Length::UnitMeter);
724 : }
725 0 : numLayer_ = ires; g2=getGroundWH2O().get();
726 0 : opacityout0=getH2OLinesOpacity(nc)*(g2/g1);
727 0 : numLayer_ = ires+1; g2=getGroundWH2O().get();
728 0 : opacityout1=getH2OLinesOpacity(nc)*(g2/g1);
729 0 : numLayer_ = numlayerold;
730 0 : return opacityout0+(opacityout1-opacityout0)*fractionLast;
731 : }
732 0 : }
733 :
734 : /*
735 : Opacity SkyStatus::getTotalOpacityUpTo(unsigned int nc, Length refalti) //15NOV2017
736 : {
737 : unsigned int ires; unsigned int numlayerold; Length alti;
738 : Opacity opacityout; Opacity opacityout0;
739 : Opacity opacityout1; Opacity zeroOp(0.0,Opacity::UnitNeper);
740 : double fractionLast; double g1; double g2;
741 :
742 : if(refalti.get(Length::UnitKiloMeter) <= altitude_.get(Length::UnitKiloMeter)) {
743 : return zeroOp;
744 : }else{
745 : fractionLast = 1.0; numlayerold = numLayer_;
746 : g1=getGroundWH2O().get(); ires = numlayerold-1; alti=altitude_;
747 : for(unsigned int i=0; i<numLayer_; i++){
748 : if(alti.get(Length::UnitKiloMeter) < refalti.get(Length::UnitKiloMeter) && (alti.get(Length::UnitKiloMeter)+v_layerThickness_[i]/1000.0) >= refalti.get(Length::UnitKiloMeter))
749 : { ires=i; fractionLast = (refalti.get(Length::UnitMeter)-alti.get(Length::UnitMeter))/v_layerThickness_[i]; }
750 : alti = alti + Length(v_layerThickness_[i],Length::UnitMeter);
751 : }
752 : numLayer_ = ires+1; g2=getGroundWH2O().get();
753 : // std::cout << "En getTotalOpacityUpTo numlayerold numLayer_ fractionLast " << numlayerold << " " << numLayer_ << " " << fractionLast << std::endl;
754 : // std::cout << "calculando opacityout1" << std::endl;
755 : opacityout1 = getWetOpacity(nc)*(g2/g1) + getDryOpacity(nc);
756 : // std::cout << "opacityout1= " << opacityout1.get("np") << std::endl;
757 : // std::cout << "calculando opacityout0" << std::endl;
758 : // std::cout << "ires=" << ires << std::endl;
759 : if(ires == 0){
760 : opacityout0 = zeroOp;
761 : }else{
762 : numLayer_ = ires; g2=getGroundWH2O().get();
763 : opacityout0 = getWetOpacity(nc)*(g2/g1) + getDryOpacity(nc);
764 : }
765 : // std::cout << "opacityout0= " << opacityout0.get("np") << std::endl;
766 : opacityout = opacityout0 + (opacityout1-opacityout0)*fractionLast;
767 : numLayer_ = numlayerold;
768 : return opacityout;
769 : }
770 : }
771 : */
772 :
773 26 : Opacity SkyStatus::getH2OContOpacity(unsigned int nc)
774 : {
775 26 : if(!chanIndexIsValid(nc)) return (double) -999.0;
776 26 : double kv = 0;
777 656 : for(unsigned int j = 0; j < numLayer_; j++) {
778 630 : kv = kv + imag(vv_N_H2OContPtr_[nc]->at(j)) * v_layerThickness_[j];
779 : }
780 52 : return ((getUserWH2O().get()) / (getGroundWH2O().get())) * kv;
781 : }
782 :
783 0 : Opacity SkyStatus::getH2OContOpacityUpTo(unsigned int nc, Length refalti)
784 : {
785 0 : unsigned int ires; unsigned int numlayerold; Length alti;
786 0 : Opacity opacityout0; Opacity opacityout1; Opacity zeroOp(0.0,Opacity::UnitNeper);
787 : double fractionLast; double g1; double g2;
788 :
789 :
790 0 : if(refalti.get(Length::UnitKiloMeter) <= altitude_.get(Length::UnitKiloMeter)) {
791 0 : return zeroOp;
792 : }else{
793 0 : fractionLast = 1.0; numlayerold = numLayer_;
794 0 : g1=getGroundWH2O().get();
795 0 : opacityout0=getH2OContOpacity(nc); ires=numlayerold-1; alti=altitude_;
796 0 : for(unsigned int i=0; i<numLayer_; i++){
797 0 : if(alti.get(Length::UnitKiloMeter) < refalti.get(Length::UnitKiloMeter) && (alti.get(Length::UnitKiloMeter)+v_layerThickness_[i]/1000.0) >= refalti.get(Length::UnitKiloMeter))
798 0 : { ires=i; fractionLast = (refalti.get(Length::UnitMeter)-alti.get(Length::UnitMeter))/v_layerThickness_[i]; }
799 0 : alti = alti + Length(v_layerThickness_[i],Length::UnitMeter);
800 : }
801 0 : numLayer_ = ires; g2=getGroundWH2O().get();
802 0 : opacityout0=getH2OContOpacity(nc)*(g2/g1);
803 0 : numLayer_ = ires+1; g2=getGroundWH2O().get();
804 0 : opacityout1=getH2OContOpacity(nc)*(g2/g1);
805 0 : numLayer_ = numlayerold;
806 0 : return opacityout0+(opacityout1-opacityout0)*fractionLast;
807 : }
808 0 : }
809 :
810 4800 : Angle SkyStatus::getDispersiveH2OPhaseDelay(unsigned int nc)
811 : {
812 4800 : if(!chanIndexIsValid(nc)) {
813 0 : Angle aa(0.0, Angle::UnitDegree);
814 0 : return aa;
815 0 : }
816 4800 : double kv = 0;
817 100400 : for(unsigned int j = 0; j < numLayer_; j++) {
818 95600 : kv = kv + real(vv_N_H2OLinesPtr_[nc]->at(j)) * v_layerThickness_[j];
819 : }
820 9600 : Angle aa(((getUserWH2O().get()) / (getGroundWH2O().get())) * kv * 57.29578,
821 4800 : Angle::UnitDegree);
822 4800 : return aa;
823 4800 : }
824 :
825 4800 : Length SkyStatus::getDispersiveH2OPathLength(unsigned int nc)
826 : {
827 4800 : if(!chanIndexIsValid(nc)) {
828 0 : Length ll(0.0, Length::UnitMilliMeter);
829 0 : return ll;
830 0 : }
831 4800 : double wavelength = 299792458.0 / v_chanFreq_[nc]; // in m
832 4800 : Length ll((wavelength / 360.0) * getDispersiveH2OPhaseDelay(nc).get(Angle::UnitDegree),
833 4800 : Length::UnitMeter);
834 4800 : return ll;
835 4800 : }
836 :
837 4800 : Angle SkyStatus::getNonDispersiveH2OPhaseDelay(unsigned int nc)
838 : {
839 4800 : double kv = 0;
840 4800 : if(!chanIndexIsValid(nc)) {
841 0 : Angle aa(0.0, Angle::UnitDegree);
842 0 : return aa;
843 0 : }
844 100400 : for(unsigned int j = 0; j < numLayer_; j++) {
845 95600 : kv = kv + real(vv_N_H2OContPtr_[nc]->at(j)) * v_layerThickness_[j];
846 : }
847 9600 : Angle aa(((getUserWH2O().get()) / (getGroundWH2O().get())) * kv * 57.29578,
848 4800 : Angle::UnitDegree);
849 4800 : return aa;
850 4800 : }
851 :
852 4800 : Length SkyStatus::getNonDispersiveH2OPathLength(unsigned int nc)
853 : {
854 4800 : if(!chanIndexIsValid(nc)) {
855 0 : Length ll(0.0, Length::UnitMilliMeter);
856 0 : return ll;
857 0 : }
858 4800 : double wavelength = 299792458.0 / v_chanFreq_[nc]; // in m
859 : Length
860 4800 : ll((wavelength / 360.0) * getNonDispersiveH2OPhaseDelay(nc).get(Angle::UnitDegree),
861 4800 : Length::UnitMeter);
862 4800 : return ll;
863 4800 : }
864 :
865 0 : Angle SkyStatus::getAverageDispersiveH2OPhaseDelay(unsigned int spwid)
866 : {
867 0 : if(!spwidAndIndexAreValid(spwid, 0)) {
868 0 : Angle aa(-999.0, Angle::UnitDegree);
869 0 : return aa;
870 0 : }
871 0 : double av = 0.0;
872 0 : for(unsigned int i = 0; i < getNumChan(spwid); i++) {
873 0 : av = av + getDispersiveH2OPhaseDelay(v_transfertId_[spwid] + i).get(Angle::UnitDegree);
874 : }
875 0 : av = av / getNumChan(spwid);
876 0 : Angle average(av, Angle::UnitDegree);
877 0 : return average;
878 0 : }
879 :
880 240 : Length SkyStatus::getAverageDispersiveH2OPathLength(unsigned int spwid)
881 : {
882 240 : if(!spwidAndIndexAreValid(spwid, 0)) {
883 0 : Length ll(0.0, Length::UnitMilliMeter);
884 0 : return ll;
885 0 : }
886 240 : double av = 0.0;
887 5040 : for(unsigned int i = 0; i < getNumChan(spwid); i++) {
888 4800 : av = av + getDispersiveH2OPathLength(v_transfertId_[spwid] + i).get(Length::UnitMilliMeter);
889 : }
890 240 : av = av / getNumChan(spwid);
891 240 : Length average(av, Length::UnitMilliMeter);
892 240 : return average;
893 240 : }
894 :
895 0 : Angle SkyStatus::getAverageNonDispersiveH2OPhaseDelay(unsigned int spwid)
896 : {
897 0 : if(!spwidAndIndexAreValid(spwid, 0)) {
898 0 : Angle aa(0.0, Angle::UnitDegree);
899 0 : return aa;
900 0 : }
901 0 : double av = 0.0;
902 0 : for(unsigned int i = 0; i < getNumChan(spwid); i++) {
903 0 : av = av
904 0 : + getNonDispersiveH2OPhaseDelay(v_transfertId_[spwid] + i).get(Angle::UnitDegree);
905 : }
906 0 : av = av / getNumChan(spwid);
907 0 : Angle average(av, Angle::UnitDegree);
908 0 : return average;
909 0 : }
910 :
911 240 : Length SkyStatus::getAverageNonDispersiveH2OPathLength(unsigned int spwid)
912 : {
913 240 : if(!spwidAndIndexAreValid(spwid, 0)) {
914 0 : Length ll(0.0, Length::UnitMilliMeter);
915 0 : return ll;
916 0 : }
917 240 : double av = 0.0;
918 5040 : for(unsigned int i = 0; i < getNumChan(spwid); i++) {
919 4800 : av = av
920 4800 : + getNonDispersiveH2OPathLength(v_transfertId_[spwid] + i).get(Length::UnitMilliMeter);
921 : }
922 240 : av = av / getNumChan(spwid);
923 240 : Length average(av, Length::UnitMilliMeter);
924 240 : return average;
925 240 : }
926 :
927 4 : Temperature SkyStatus::getAverageTebbSky(unsigned int spwid,
928 : const Length &wh2o,
929 : double airmass,
930 : double skycoupling,
931 : const Temperature &Tspill)
932 : {
933 4 : Temperature tt(-999, Temperature::UnitKelvin);
934 4 : if(!spwidAndIndexAreValid(spwid, 0)) {
935 0 : return tt;
936 : }
937 4 : if(wh2o.get() < 0.0) {
938 0 : return tt;
939 : }
940 : // if(skycoupling<0.0 || skycoupling>1.0){return tt;}
941 4 : if(airmass < 1.0) {
942 0 : return tt;
943 : }
944 4 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
945 0 : return tt;
946 : }
947 8 : return RT(((wh2o.get()) / (getGroundWH2O().get())),
948 : skycoupling,
949 : Tspill.get(Temperature::UnitKelvin),
950 : airmass,
951 8 : spwid);
952 4 : }
953 :
954 0 : Temperature SkyStatus::getAverageTebbSky(unsigned int spwid,
955 : const Length &wh2o,
956 : double airmass,
957 : double skycoupling,
958 : double signalgain, // adition
959 : const Temperature &Tspill)
960 : {
961 0 : Temperature tt(-999, Temperature::UnitKelvin);
962 0 : if(!spwidAndIndexAreValid(spwid, 0)) {
963 0 : return tt;
964 : }
965 0 : if(wh2o.get() < 0.0) {
966 0 : return tt;
967 : }
968 : // if(skycoupling<0.0 || skycoupling>1.0){return tt;}
969 0 : if(airmass < 1.0) {
970 0 : return tt;
971 : }
972 0 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
973 0 : return tt;
974 : }
975 0 : return signalgain*RT(((wh2o.get()) / (getGroundWH2O().get())),
976 : skycoupling,
977 : Tspill.get(Temperature::UnitKelvin),
978 : airmass,
979 : spwid)+
980 0 : +(1-signalgain)*RT(((wh2o.get()) / (getGroundWH2O().get())),
981 : skycoupling,
982 : Tspill.get(Temperature::UnitKelvin),
983 : airmass,
984 0 : getAssocSpwId(spwid)[0]);
985 0 : }
986 :
987 11238 : Temperature SkyStatus::getTebbSky(unsigned int spwid,
988 : unsigned int nc,
989 : const Length &wh2o,
990 : double airmass,
991 : double skycoupling,
992 : const Temperature &Tspill)
993 : {
994 11238 : Temperature tt(-999, Temperature::UnitKelvin);
995 11238 : if(!spwidAndIndexAreValid(spwid, nc)) {
996 0 : return tt;
997 : }
998 11238 : if(wh2o.get() < 0.0) {
999 0 : return tt;
1000 : }
1001 11238 : if(skycoupling < 0.0 || skycoupling > 1.0) {
1002 0 : return tt;
1003 : }
1004 11238 : if(airmass < 1.0) {
1005 0 : return tt;
1006 : }
1007 11238 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
1008 0 : return tt;
1009 : }
1010 22476 : return RT(((wh2o.get()) / (getGroundWH2O().get())),
1011 : skycoupling,
1012 : Tspill.get(Temperature::UnitKelvin),
1013 : airmass,
1014 : spwid,
1015 22476 : nc);
1016 11238 : }
1017 :
1018 :
1019 0 : Temperature SkyStatus::getAverageTrjSky(unsigned int spwid,
1020 : const Length &wh2o,
1021 : double airmass,
1022 : double skycoupling,
1023 : const Temperature &Tspill)
1024 : {
1025 0 : Temperature tt(-999, Temperature::UnitKelvin);
1026 0 : if(!spwidAndIndexAreValid(spwid, 0)) {
1027 0 : return tt;
1028 : }
1029 0 : if(wh2o.get() < 0.0) {
1030 0 : return tt;
1031 : }
1032 : // if(skycoupling<0.0 || skycoupling>1.0){return tt;}
1033 0 : if(airmass < 1.0) {
1034 0 : return tt;
1035 : }
1036 0 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
1037 0 : return tt;
1038 : }
1039 0 : return RTRJ(((wh2o.get()) / (getGroundWH2O().get())),
1040 : skycoupling,
1041 : Tspill.get(Temperature::UnitKelvin),
1042 : airmass,
1043 0 : spwid);
1044 0 : }
1045 :
1046 0 : Temperature SkyStatus::getAverageTrjSky(unsigned int spwid,
1047 : const Length &wh2o,
1048 : double airmass,
1049 : double skycoupling,
1050 : double signalgain, // adition
1051 : const Temperature &Tspill)
1052 : {
1053 0 : Temperature tt(-999, Temperature::UnitKelvin);
1054 0 : if(!spwidAndIndexAreValid(spwid, 0)) {
1055 0 : return tt;
1056 : }
1057 0 : if(wh2o.get() < 0.0) {
1058 0 : return tt;
1059 : }
1060 : // if(skycoupling<0.0 || skycoupling>1.0){return tt;}
1061 0 : if(airmass < 1.0) {
1062 0 : return tt;
1063 : }
1064 0 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
1065 0 : return tt;
1066 : }
1067 0 : return signalgain*RTRJ(((wh2o.get()) / (getGroundWH2O().get())),
1068 : skycoupling,
1069 : Tspill.get(Temperature::UnitKelvin),
1070 : airmass,
1071 : spwid)+
1072 0 : +(1-signalgain)*RTRJ(((wh2o.get()) / (getGroundWH2O().get())),
1073 : skycoupling,
1074 : Tspill.get(Temperature::UnitKelvin),
1075 : airmass,
1076 0 : getAssocSpwId(spwid)[0]);
1077 0 : }
1078 :
1079 0 : Temperature SkyStatus::getTrjSky(unsigned int spwid,
1080 : unsigned int nc,
1081 : const Length &wh2o,
1082 : double airmass,
1083 : double skycoupling,
1084 : const Temperature &Tspill)
1085 : {
1086 0 : Temperature tt(-999, Temperature::UnitKelvin);
1087 0 : if(!spwidAndIndexAreValid(spwid, nc)) {
1088 0 : return tt;
1089 : }
1090 0 : if(wh2o.get() < 0.0) {
1091 0 : return tt;
1092 : }
1093 0 : if(skycoupling < 0.0 || skycoupling > 1.0) {
1094 0 : return tt;
1095 : }
1096 0 : if(airmass < 1.0) {
1097 0 : return tt;
1098 : }
1099 0 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
1100 0 : return tt;
1101 : }
1102 0 : return RTRJ(((wh2o.get()) / (getGroundWH2O().get())),
1103 : skycoupling,
1104 : Tspill.get(Temperature::UnitKelvin),
1105 : airmass,
1106 : spwid,
1107 0 : nc);
1108 0 : }
1109 :
1110 :
1111 :
1112 0 : Angle SkyStatus::getDispersiveH2OPhaseDelay(unsigned int spwid, unsigned int nc)
1113 : {
1114 0 : if(!spwidAndIndexAreValid(spwid, nc)) {
1115 0 : Angle aa(0.0, Angle::UnitDegree);
1116 0 : return aa;
1117 0 : }
1118 0 : return getDispersiveH2OPhaseDelay(v_transfertId_[spwid] + nc);
1119 : }
1120 :
1121 0 : Length SkyStatus::getDispersiveH2OPathLength(unsigned int spwid,
1122 : unsigned int nc)
1123 : {
1124 0 : if(!spwidAndIndexAreValid(spwid, nc)) {
1125 0 : Length ll(0.0, Length::UnitMilliMeter);
1126 0 : return ll;
1127 0 : }
1128 0 : return getDispersiveH2OPathLength(v_transfertId_[spwid] + nc);
1129 : }
1130 :
1131 0 : Angle SkyStatus::getNonDispersiveH2OPhaseDelay(unsigned int spwid,
1132 : unsigned int nc)
1133 : {
1134 0 : if(!spwidAndIndexAreValid(spwid, nc)) {
1135 0 : Angle aa(0.0, Angle::UnitDegree);
1136 0 : return aa;
1137 0 : }
1138 0 : return getNonDispersiveH2OPhaseDelay(v_transfertId_[spwid] + nc);
1139 : }
1140 :
1141 0 : Length SkyStatus::getNonDispersiveH2OPathLength(unsigned int spwid,
1142 : unsigned int nc)
1143 : {
1144 0 : if(!spwidAndIndexAreValid(spwid, nc)) return (double) 0.0;
1145 0 : return getNonDispersiveH2OPathLength(v_transfertId_[spwid] + nc);
1146 : }
1147 :
1148 0 : Length SkyStatus::WaterVaporRetrieval_fromFTS(unsigned int spwId,
1149 : const vector<double> &v_transmission,
1150 : const Frequency &f1,
1151 : const Frequency &f2)
1152 : {
1153 0 : if(f1.get() > f2.get()) {
1154 0 : return Length(-999, Length::UnitMilliMeter);
1155 : }
1156 0 : if(v_transmission.size() == getSpectralWindow(spwId).size()) {
1157 : return mkWaterVaporRetrieval_fromFTS(spwId,
1158 : v_transmission,
1159 : // getAirMass(), // unused parameter
1160 : f1,
1161 0 : f2);
1162 : } else {
1163 0 : return Length(-999.0, Length::UnitMilliMeter);
1164 : }
1165 : }
1166 :
1167 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1168 : const vector<Temperature> &v_tebb,
1169 : const vector<double> &skycoupling,
1170 : const vector<Temperature> &tspill)
1171 : {
1172 0 : vector<vector<double> > spwId_filters;
1173 0 : vector<double> spwId_filter;
1174 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1175 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1176 0 : spwId_filter.push_back(1.0);
1177 : }
1178 0 : spwId_filters.push_back(spwId_filter);
1179 0 : spwId_filter.clear();
1180 : }
1181 : return WaterVaporRetrieval_fromTEBB(spwId,
1182 : v_tebb,
1183 : spwId_filters,
1184 : skycoupling,
1185 0 : tspill);
1186 0 : }
1187 :
1188 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1189 : const vector<Temperature> &v_tebb,
1190 : double skycoupling,
1191 : const Temperature &tspill)
1192 : {
1193 0 : vector<double> spwId_filter;
1194 0 : for(unsigned int n = 0; n < v_numChan_[spwId]; n++) {
1195 0 : spwId_filter.push_back(1.0);
1196 : }
1197 : return WaterVaporRetrieval_fromTEBB(spwId,
1198 : v_tebb,
1199 : spwId_filter,
1200 : skycoupling,
1201 0 : tspill);
1202 0 : }
1203 :
1204 :
1205 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1206 : const vector< vector<Temperature> > &vv_tebb,
1207 : const vector<double> &skycoupling,
1208 : const vector<Temperature> &tspill)
1209 : {
1210 0 : vector<vector<double> > spwId_filters;
1211 0 : vector<double> spwId_filter;
1212 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1213 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1214 0 : spwId_filter.push_back(1.0);
1215 : }
1216 0 : spwId_filters.push_back(spwId_filter);
1217 0 : spwId_filter.clear();
1218 : }
1219 : return WaterVaporRetrieval_fromTEBB(spwId,
1220 : vv_tebb,
1221 : spwId_filters,
1222 : skycoupling,
1223 0 : tspill);
1224 0 : }
1225 :
1226 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1227 : const Percent &signalGain,
1228 : const vector<Temperature> &v_tebb,
1229 : double skycoupling,
1230 : const Temperature &tspill)
1231 : {
1232 0 : vector<double> spwId_filter;
1233 0 : for(unsigned int n = 0; n < v_numChan_[spwId]; n++) {
1234 0 : spwId_filter.push_back(1.0);
1235 : }
1236 : return WaterVaporRetrieval_fromTEBB(spwId,
1237 : signalGain,
1238 : v_tebb,
1239 : spwId_filter,
1240 : skycoupling,
1241 0 : tspill);
1242 0 : }
1243 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1244 : const vector<Percent> &signalGain,
1245 : const vector<vector<Temperature> > &vv_tebb,
1246 : const vector<double> &skycoupling,
1247 : const vector<Temperature> &tspill)
1248 : {
1249 0 : vector<vector<double> > spwId_filters;
1250 0 : vector<double> spwId_filter;
1251 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1252 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1253 0 : spwId_filter.push_back(1.0);
1254 : }
1255 0 : spwId_filters.push_back(spwId_filter);
1256 0 : spwId_filter.clear();
1257 : }
1258 : return WaterVaporRetrieval_fromTEBB(spwId,
1259 : signalGain,
1260 : vv_tebb,
1261 : spwId_filters,
1262 : skycoupling,
1263 0 : tspill);
1264 0 : }
1265 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1266 : const vector<Percent> &signalGain,
1267 : const vector<Temperature> &v_tebb,
1268 : const vector<double> &skycoupling,
1269 : const vector<Temperature> &tspill)
1270 : {
1271 0 : vector<vector<double> > spwId_filters;
1272 0 : vector<double> spwId_filter;
1273 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1274 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1275 0 : spwId_filter.push_back(1.0);
1276 : }
1277 0 : spwId_filters.push_back(spwId_filter);
1278 0 : spwId_filter.clear();
1279 : }
1280 : return WaterVaporRetrieval_fromTEBB(spwId,
1281 : signalGain,
1282 : v_tebb,
1283 : spwId_filters,
1284 : skycoupling,
1285 0 : tspill);
1286 0 : }
1287 :
1288 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1289 : const vector<Temperature> &v_tebb,
1290 : const vector<double> &spwId_filter,
1291 : double skycoupling,
1292 : const Temperature &tspill)
1293 : {
1294 0 : if(v_tebb.size() == getSpectralWindow(spwId).size()) {
1295 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1296 0 : Percent(100.0, Percent::UnitPercent),
1297 : v_tebb,
1298 : getAirMass(),
1299 : spwId_filter,
1300 : skycoupling,
1301 0 : tspill);
1302 : } else {
1303 0 : return Length(-999.0, Length::UnitMilliMeter);
1304 : }
1305 : }
1306 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1307 : const vector<vector<Temperature> > &vv_tebb,
1308 : const vector<vector<double> > &spwId_filters,
1309 : const vector<double> &skycoupling,
1310 : const vector<Temperature> &tspill)
1311 : {
1312 0 : vector<Percent> signalGain;
1313 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1314 0 : signalGain.push_back(Percent(100.0, Percent::UnitPercent));
1315 : }
1316 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1317 : signalGain,
1318 : vv_tebb,
1319 : getAirMass(),
1320 : spwId_filters,
1321 : skycoupling,
1322 0 : tspill);
1323 0 : }
1324 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1325 : const vector<Temperature> &v_tebb,
1326 : const vector<vector<double> > &spwId_filters,
1327 : const vector<double> &skycoupling,
1328 : const vector<Temperature> &tspill)
1329 : {
1330 0 : vector<Percent> signalGain;
1331 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1332 0 : signalGain.push_back(Percent(100.0, Percent::UnitPercent));
1333 : }
1334 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1335 : signalGain,
1336 : v_tebb,
1337 : getAirMass(),
1338 : spwId_filters,
1339 : skycoupling,
1340 0 : tspill);
1341 0 : }
1342 :
1343 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1344 : const Percent &signalGain,
1345 : const vector<Temperature> &v_tebb,
1346 : const vector<double> &spwId_filter,
1347 : double skycoupling,
1348 : const Temperature &tspill)
1349 : {
1350 0 : if(v_tebb.size() == getSpectralWindow(spwId).size()) {
1351 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1352 : signalGain,
1353 : v_tebb,
1354 : getAirMass(),
1355 : spwId_filter,
1356 : skycoupling,
1357 0 : tspill);
1358 : } else {
1359 0 : return Length(-999.0, Length::UnitMilliMeter);
1360 : }
1361 : }
1362 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1363 : const vector<Percent> &signalGain,
1364 : const vector<vector<Temperature> > &vv_tebb,
1365 : const vector<vector<double> > &spwId_filters,
1366 : const vector<double> &skycoupling,
1367 : const vector<Temperature> &tspill)
1368 : {
1369 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1370 : signalGain,
1371 : vv_tebb,
1372 : getAirMass(),
1373 : spwId_filters,
1374 : skycoupling,
1375 0 : tspill);
1376 : }
1377 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1378 : const vector<Percent> &signalGain,
1379 : const vector<Temperature> &v_tebb,
1380 : const vector<vector<double> > &spwId_filters,
1381 : const vector<double> &skycoupling,
1382 : const vector<Temperature> &tspill)
1383 : {
1384 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1385 : signalGain,
1386 : v_tebb,
1387 : getAirMass(),
1388 : spwId_filters,
1389 : skycoupling,
1390 0 : tspill);
1391 : }
1392 :
1393 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1394 : const vector<Temperature> &v_tebb,
1395 : double airmass,
1396 : double skycoupling,
1397 : const Temperature &tspill)
1398 : {
1399 0 : vector<double> spwId_filter;
1400 0 : for(unsigned int n = 0; n < v_numChan_[spwId]; n++) {
1401 0 : spwId_filter.push_back(1.0);
1402 : }
1403 : return WaterVaporRetrieval_fromTEBB(spwId,
1404 : v_tebb,
1405 : spwId_filter,
1406 : airmass,
1407 : skycoupling,
1408 0 : tspill);
1409 0 : }
1410 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1411 : const vector<vector<Temperature> > &vv_tebb,
1412 : double airmass,
1413 : const vector<double> &skycoupling,
1414 : const vector<Temperature> &tspill)
1415 : {
1416 0 : vector<vector<double> > spwId_filters;
1417 0 : vector<double> spwId_filter;
1418 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1419 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1420 0 : spwId_filter.push_back(1.0);
1421 : }
1422 0 : spwId_filters.push_back(spwId_filter);
1423 0 : spwId_filter.clear();
1424 : }
1425 : return WaterVaporRetrieval_fromTEBB(spwId,
1426 : vv_tebb,
1427 : spwId_filters,
1428 : airmass,
1429 : skycoupling,
1430 0 : tspill);
1431 0 : }
1432 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1433 : const vector<Temperature> &v_tebb,
1434 : double airmass,
1435 : const vector<double> &skycoupling,
1436 : const vector<Temperature> &tspill)
1437 : {
1438 0 : vector<vector<double> > spwId_filters;
1439 0 : vector<double> spwId_filter;
1440 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1441 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1442 0 : spwId_filter.push_back(1.0);
1443 : }
1444 0 : spwId_filters.push_back(spwId_filter);
1445 0 : spwId_filter.clear();
1446 : }
1447 : return WaterVaporRetrieval_fromTEBB(spwId,
1448 : v_tebb,
1449 : spwId_filters,
1450 : airmass,
1451 : skycoupling,
1452 0 : tspill);
1453 0 : }
1454 :
1455 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1456 : const Percent &signalGain,
1457 : const vector<Temperature> &v_tebb,
1458 : double airmass,
1459 : double skycoupling,
1460 : const Temperature &tspill)
1461 : {
1462 0 : vector<double> spwId_filter;
1463 0 : for(unsigned int n = 0; n < v_numChan_[spwId]; n++) {
1464 0 : spwId_filter.push_back(1.0);
1465 : }
1466 : return WaterVaporRetrieval_fromTEBB(spwId,
1467 : signalGain,
1468 : v_tebb,
1469 : spwId_filter,
1470 : airmass,
1471 : skycoupling,
1472 0 : tspill);
1473 0 : }
1474 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1475 : const vector<Percent> &signalGain,
1476 : const vector<vector<Temperature> > &vv_tebb,
1477 : double airmass,
1478 : const vector<double> &skycoupling,
1479 : const vector<Temperature> &tspill)
1480 : {
1481 0 : vector<vector<double> > spwId_filters;
1482 0 : vector<double> spwId_filter;
1483 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1484 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1485 0 : spwId_filter.push_back(1.0);
1486 : }
1487 0 : spwId_filters.push_back(spwId_filter);
1488 0 : spwId_filter.clear();
1489 : }
1490 : return WaterVaporRetrieval_fromTEBB(spwId,
1491 : signalGain,
1492 : vv_tebb,
1493 : spwId_filters,
1494 : airmass,
1495 : skycoupling,
1496 0 : tspill);
1497 0 : }
1498 :
1499 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1500 : const Percent &signalGain,
1501 : const Temperature &tebb,
1502 : double airmass,
1503 : double skycoupling,
1504 : const Temperature &tspill)
1505 : {
1506 0 : vector<unsigned int> spwIdv;
1507 0 : vector<Percent> signalGainv;
1508 0 : vector<Temperature> tebbv;
1509 0 : vector<double> skycouplingv;
1510 0 : vector<Temperature> tspillv;
1511 0 : spwIdv.push_back(spwId);
1512 0 : signalGainv.push_back(signalGain);
1513 0 : tebbv.push_back(tebb);
1514 0 : skycouplingv.push_back(skycoupling);
1515 0 : tspillv.push_back(tspill);
1516 0 : return WaterVaporRetrieval_fromTEBB(spwIdv,signalGainv,tebbv,airmass,skycouplingv,tspillv);
1517 0 : }
1518 :
1519 :
1520 :
1521 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1522 : const vector<Percent> &signalGain,
1523 : const vector<Temperature> &v_tebb,
1524 : double airmass,
1525 : const vector<double> &skycoupling,
1526 : const vector<Temperature> &tspill)
1527 : {
1528 0 : vector<vector<double> > spwId_filters;
1529 0 : vector<double> spwId_filter;
1530 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1531 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1532 0 : spwId_filter.push_back(1.0);
1533 : }
1534 0 : spwId_filters.push_back(spwId_filter);
1535 0 : spwId_filter.clear();
1536 : }
1537 : return WaterVaporRetrieval_fromTEBB(spwId,
1538 : signalGain,
1539 : v_tebb,
1540 : spwId_filters,
1541 : airmass,
1542 : skycoupling,
1543 0 : tspill);
1544 0 : }
1545 :
1546 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1547 : const vector<Temperature> &v_tebb,
1548 : const vector<double> &spwId_filter,
1549 : double airmass,
1550 : double skycoupling,
1551 : const Temperature &tspill)
1552 : {
1553 0 : if(v_tebb.size() == getSpectralWindow(spwId).size()) {
1554 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1555 0 : Percent(100.0, Percent::UnitPercent),
1556 : v_tebb,
1557 : airmass,
1558 : spwId_filter,
1559 : skycoupling,
1560 0 : tspill);
1561 : } else {
1562 0 : return Length(-999.0, Length::UnitMilliMeter);
1563 : }
1564 : }
1565 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1566 : const vector<vector<Temperature> > &vv_tebb,
1567 : const vector<vector<double> > &spwId_filters,
1568 : double airmass,
1569 : const vector<double> &skycoupling,
1570 : const vector<Temperature> &tspill)
1571 : {
1572 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
1573 0 : if(vv_tebb[j].size() != getSpectralWindow(spwId[j]).size()) {
1574 0 : return Length(-999.0, Length::UnitMilliMeter);
1575 : }
1576 : }
1577 0 : vector<Percent> signalGain;
1578 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1579 0 : signalGain.push_back(Percent(100.0, Percent::UnitPercent));
1580 : }
1581 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1582 : signalGain,
1583 : vv_tebb,
1584 : airmass,
1585 : spwId_filters,
1586 : skycoupling,
1587 0 : tspill);
1588 0 : }
1589 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1590 : const vector<Temperature> &v_tebb,
1591 : const vector<vector<double> > &spwId_filters,
1592 : double airmass,
1593 : const vector<double> &skycoupling,
1594 : const vector<Temperature> &tspill)
1595 : {
1596 0 : vector<Percent> signalGain;
1597 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1598 0 : signalGain.push_back(Percent(100.0, Percent::UnitPercent));
1599 : }
1600 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1601 : signalGain,
1602 : v_tebb,
1603 : airmass,
1604 : spwId_filters,
1605 : skycoupling,
1606 0 : tspill);
1607 0 : }
1608 :
1609 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(unsigned int spwId,
1610 : const Percent &signalGain,
1611 : const vector<Temperature> &v_tebb,
1612 : const vector<double> &spwId_filter,
1613 : double airmass,
1614 : double skycoupling,
1615 : const Temperature &tspill)
1616 : {
1617 0 : if(v_tebb.size() == getSpectralWindow(spwId).size()) {
1618 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1619 : signalGain,
1620 : v_tebb,
1621 : airmass,
1622 : spwId_filter,
1623 : skycoupling,
1624 0 : tspill);
1625 : } else {
1626 0 : return Length(-999.0, Length::UnitMilliMeter);
1627 : }
1628 : }
1629 :
1630 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1631 : const vector<Percent> &signalGain,
1632 : const vector<vector<Temperature> > &vv_tebb,
1633 : const vector<vector<double> > &spwId_filter,
1634 : double airmass,
1635 : const vector<double> &skycoupling,
1636 : const vector<Temperature> &tspill)
1637 : {
1638 :
1639 0 : if(spwId.size() != signalGain.size()) {
1640 0 : return Length(-999.0, Length::UnitMilliMeter);
1641 : }
1642 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
1643 0 : if(vv_tebb[j].size() != getSpectralWindow(spwId[j]).size()) {
1644 0 : return Length(-999.0, Length::UnitMilliMeter);
1645 : }
1646 : }
1647 0 : if(spwId.size() != spwId_filter.size()) {
1648 0 : return Length(-999.0, Length::UnitMilliMeter);
1649 : }
1650 0 : if(spwId.size() != skycoupling.size()) {
1651 0 : return Length(-999.0, Length::UnitMilliMeter);
1652 : }
1653 0 : if(spwId.size() != tspill.size()) {
1654 0 : return Length(-999.0, Length::UnitMilliMeter);
1655 : }
1656 :
1657 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1658 : signalGain,
1659 : vv_tebb,
1660 : airmass,
1661 : spwId_filter,
1662 : skycoupling,
1663 0 : tspill);
1664 : }
1665 :
1666 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1667 : const vector<Percent> &signalGain,
1668 : const vector<Temperature> &v_tebb,
1669 : const vector<vector<double> > &spwId_filter,
1670 : double airmass,
1671 : const vector<double> &skycoupling,
1672 : const vector<Temperature> &tspill)
1673 : {
1674 0 : if(spwId.size() != signalGain.size()) {
1675 0 : return Length(-999.0, Length::UnitMilliMeter);
1676 : }
1677 0 : if(spwId.size() != v_tebb.size()) {
1678 0 : return Length(-999.0, Length::UnitMilliMeter);
1679 : }
1680 0 : if(spwId.size() != spwId_filter.size()) {
1681 0 : return Length(-999.0, Length::UnitMilliMeter);
1682 : }
1683 0 : if(spwId.size() != skycoupling.size()) {
1684 0 : return Length(-999.0, Length::UnitMilliMeter);
1685 : }
1686 0 : if(spwId.size() != tspill.size()) {
1687 0 : return Length(-999.0, Length::UnitMilliMeter);
1688 : }
1689 :
1690 : return mkWaterVaporRetrieval_fromTEBB(spwId,
1691 : signalGain,
1692 : v_tebb,
1693 : airmass,
1694 : spwId_filter,
1695 : skycoupling,
1696 0 : tspill);
1697 : }
1698 :
1699 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1700 : const Percent &signalGain,
1701 : const vector<Temperature> &v_tebb,
1702 : const vector<vector<double> > &spwId_filters,
1703 : double airmass,
1704 : double skycoupling,
1705 : const Temperature &tspill)
1706 : {
1707 0 : vector<Percent> v_signalGain;
1708 0 : vector<double> v_skycoupling;
1709 0 : vector<Temperature> v_tspill;
1710 0 : v_signalGain.reserve(spwId.size());
1711 0 : v_skycoupling.reserve(spwId.size());
1712 0 : v_tspill.reserve(spwId.size());
1713 :
1714 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
1715 0 : v_signalGain.push_back(signalGain);
1716 0 : v_skycoupling.push_back(skycoupling);
1717 0 : v_tspill.push_back(tspill);
1718 : }
1719 : return WaterVaporRetrieval_fromTEBB(spwId,
1720 : v_signalGain,
1721 : v_tebb,
1722 : spwId_filters,
1723 : airmass,
1724 : v_skycoupling,
1725 0 : v_tspill);
1726 0 : }
1727 0 : Length SkyStatus::WaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
1728 : const Percent &signalGain,
1729 : const vector<Temperature> &v_tebb,
1730 : double airmass,
1731 : double skycoupling,
1732 : const Temperature &tspill)
1733 : {
1734 0 : vector<vector<double> > spwId_filters;
1735 0 : vector<double> spwId_filter;
1736 0 : for(unsigned int i = 0; i < spwId.size(); i++) {
1737 0 : for(unsigned int n = 0; n < v_numChan_[spwId[i]]; n++) {
1738 0 : spwId_filter.push_back(1.0);
1739 : }
1740 0 : spwId_filters.push_back(spwId_filter);
1741 0 : spwId_filter.clear();
1742 : }
1743 : return WaterVaporRetrieval_fromTEBB(spwId,
1744 : signalGain,
1745 : v_tebb,
1746 : spwId_filters,
1747 : airmass,
1748 : skycoupling,
1749 0 : tspill);
1750 0 : }
1751 :
1752 0 : double SkyStatus::SkyCouplingRetrieval_fromTEBB(unsigned int spwId,
1753 : const vector<Temperature> &v_tebb,
1754 : double skycoupling,
1755 : const Temperature &tspill)
1756 : {
1757 0 : vector<double> spwId_filter;
1758 0 : for(unsigned int n = 0; n < v_numChan_[spwId]; n++) {
1759 0 : spwId_filter.push_back(1.0);
1760 : }
1761 0 : return SkyCouplingRetrieval_fromTEBB(spwId,
1762 : v_tebb,
1763 : spwId_filter,
1764 : skycoupling,
1765 0 : tspill);
1766 0 : }
1767 :
1768 0 : double SkyStatus::SkyCouplingRetrieval_fromTEBB(unsigned int spwId,
1769 : const vector<Temperature> &v_tebb,
1770 : const vector<double> &spwId_filter,
1771 : double skycoupling,
1772 : const Temperature &tspill)
1773 : {
1774 0 : if(v_tebb.size() == getSpectralWindow(spwId).size()) {
1775 0 : return mkSkyCouplingRetrieval_fromTEBB(spwId,
1776 0 : Percent(100, Percent::UnitPercent),
1777 : v_tebb,
1778 : getAirMass(),
1779 : spwId_filter,
1780 : skycoupling,
1781 0 : tspill);
1782 : } else {
1783 0 : return -999.0;
1784 : }
1785 : }
1786 :
1787 0 : double SkyStatus::SkyCouplingRetrieval_fromTEBB(unsigned int spwId,
1788 : const vector<Temperature> &v_tebb,
1789 : double airmass,
1790 : double skycoupling,
1791 : const Temperature &tspill)
1792 : {
1793 0 : vector<double> spwId_filter;
1794 0 : for(unsigned int n = 0; n < v_numChan_[spwId]; n++) {
1795 0 : spwId_filter.push_back(1.0);
1796 : }
1797 0 : return SkyCouplingRetrieval_fromTEBB(spwId,
1798 : v_tebb,
1799 : spwId_filter,
1800 : airmass,
1801 : skycoupling,
1802 0 : tspill);
1803 0 : }
1804 :
1805 0 : double SkyStatus::SkyCouplingRetrieval_fromTEBB(unsigned int spwId,
1806 : const vector<Temperature> &v_tebb,
1807 : const vector<double> &spwId_filter,
1808 : double airmass,
1809 : double skycoupling,
1810 : const Temperature &tspill)
1811 : {
1812 0 : if(v_tebb.size() == getSpectralWindow(spwId).size()) {
1813 0 : return mkSkyCouplingRetrieval_fromTEBB(spwId,
1814 0 : Percent(100, Percent::UnitPercent),
1815 : v_tebb,
1816 : airmass,
1817 : spwId_filter,
1818 : skycoupling,
1819 0 : tspill);
1820 : } else {
1821 0 : return -999.0;
1822 : }
1823 : }
1824 :
1825 0 : double SkyStatus::getSigmaTransmissionFit(unsigned int spwId,
1826 : const vector<double> &v_transmission,
1827 : double airm,
1828 : const Frequency &f1,
1829 : const Frequency &f2)
1830 : {
1831 0 : if(f1.get() > f2.get()) {
1832 0 : return -999.0;
1833 : }
1834 0 : if(v_transmission.size() == getSpectralWindow(spwId).size()) {
1835 0 : double rms = 0.0;
1836 0 : unsigned int num = 0;
1837 :
1838 0 : for(unsigned int i = 0; i < v_transmission.size(); i++) {
1839 0 : if((getSpectralWindow(spwId)[i] * 1E-09 >= f1.get(Frequency::UnitGigaHertz)
1840 0 : && getSpectralWindow(spwId)[i] * 1E-09 <= f2.get(Frequency::UnitGigaHertz))) {
1841 0 : num++;
1842 0 : rms = rms + pow((v_transmission[i] - exp(-airm
1843 0 : * ( (getDryContOpacity(spwId, i).get()+getO2LinesOpacity(spwId, i).get()+0.65*getO3LinesOpacity(spwId, i).get() ) // getDryOpacity(spwId, i).get()
1844 0 : + getWetOpacity(spwId, i).get()))),
1845 : 2);
1846 : }
1847 : }
1848 0 : rms = sqrt(rms / num);
1849 0 : return rms;
1850 : } else {
1851 0 : return -999.0;
1852 : }
1853 : }
1854 :
1855 0 : Temperature SkyStatus::getSigmaFit(unsigned int spwId,
1856 : const vector<Temperature> &v_tebbspec,
1857 : const Length &wh2o,
1858 : double airmass,
1859 : double skyCoupling,
1860 : const Temperature &Tspill)
1861 : {
1862 0 : Temperature ttt(-999, Temperature::UnitKelvin);
1863 0 : if(!spwidAndIndexAreValid(spwId, 0)) {
1864 0 : return ttt;
1865 : }
1866 0 : if(v_tebbspec.size() != getSpectralWindow(spwId).size()) {
1867 0 : return ttt;
1868 : }
1869 0 : if(wh2o.get(Length::UnitMilliMeter) < 0.0) {
1870 0 : return ttt;
1871 : }
1872 0 : if(skyCoupling < 0.0 || skyCoupling > 1.0) {
1873 0 : return ttt;
1874 : }
1875 0 : if(airmass < 1.0) {
1876 0 : return ttt;
1877 : }
1878 0 : if(Tspill.get(Temperature::UnitKelvin) < 0.0 || Tspill.get(Temperature::UnitKelvin) > 350.0) {
1879 0 : return ttt;
1880 : }
1881 0 : double rms = 0.0;
1882 0 : unsigned int num = 0;
1883 0 : for(unsigned int i = 0; i < v_tebbspec.size(); i++) {
1884 0 : if(v_tebbspec[i].get() > 0.0) {
1885 0 : num++;
1886 0 : rms = rms + pow((v_tebbspec[i].get(Temperature::UnitKelvin) - getTebbSky(spwId,
1887 : i,
1888 : wh2o,
1889 : airmass,
1890 : skyCoupling,
1891 0 : Tspill).get(Temperature::UnitKelvin)),
1892 : 2);
1893 : }
1894 : }
1895 0 : rms = sqrt(rms / num);
1896 0 : return Temperature(rms, Temperature::UnitKelvin);
1897 0 : }
1898 :
1899 0 : Length SkyStatus::mkWaterVaporRetrieval_fromFTS(unsigned int spwId,
1900 : const vector<double> &measuredSkyTransmission,
1901 : //double airm, // unused parameter
1902 : const Frequency &fre1,
1903 : const Frequency &fre2)
1904 : {
1905 : double pfit_wh2o;
1906 0 : double deltaa = 0.02;
1907 0 : double sig_fit = -999.0;
1908 0 : double eps = 0.01;
1909 0 : vector<double> transmission_fit;
1910 0 : transmission_fit.reserve(measuredSkyTransmission.size()); // allows resizing
1911 :
1912 : unsigned int num;
1913 : double flamda;
1914 0 : unsigned int niter = 20;
1915 : double alpha;
1916 : double beta;
1917 : double array;
1918 : double f1;
1919 : double psave;
1920 : double f2;
1921 : double deriv;
1922 : double chisq1;
1923 : double chisqr;
1924 : double pfit_wh2o_b;
1925 : double res;
1926 0 : Length wh2o_retrieved(-999.0, Length::UnitMilliMeter);
1927 0 : Length werr(-888, Length::UnitMilliMeter);
1928 : double sigma_fit_transm0;
1929 0 : Length sigma_wh2o;
1930 :
1931 0 : num = 0;
1932 :
1933 0 : flamda = 0.001;
1934 0 : pfit_wh2o = 1.0; // (getUserWH2O().get(Length::UnitMilliMeter))/(getGroundWH2O().get(Length::UnitMilliMeter));
1935 :
1936 0 : unsigned int nl = 0;
1937 :
1938 0 : if(fre1.get(Frequency::UnitGigaHertz) < 0) {
1939 0 : nl = getSpectralWindow(spwId).size();
1940 : } else {
1941 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId).size(); i++) {
1942 0 : if(getSpectralWindow(spwId)[i] * 1E-09 >= fre1.get(Frequency::UnitGigaHertz)
1943 0 : && getSpectralWindow(spwId)[i] * 1E-09 <= fre2.get(Frequency::UnitGigaHertz)) {
1944 0 : nl = nl + 1;
1945 : }
1946 : }
1947 : }
1948 :
1949 0 : for(unsigned int kite = 0; kite < niter; kite++) {
1950 :
1951 0 : num = num + 1;
1952 :
1953 0 : beta = 0.0;
1954 0 : alpha = 0.0;
1955 :
1956 : // for(unsigned int i=0; i<getSpectralWindow(spwId).size(); i++){
1957 :
1958 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId).size(); i++) {
1959 :
1960 0 : if(nl == getSpectralWindow(spwId).size() || (getSpectralWindow(spwId)[i]
1961 0 : * 1E-09 >= fre1.get(Frequency::UnitGigaHertz) && getSpectralWindow(spwId)[i] * 1E-09
1962 0 : <= fre2.get(Frequency::UnitGigaHertz))) {
1963 : //
1964 0 : transmission_fit[i] = exp(-( (getDryContOpacity(spwId, i).get()+getO2LinesOpacity(spwId, i).get()+0.65*getO3LinesOpacity(spwId, i).get() ) //getDryOpacity(spwId, i).get()
1965 0 : + pfit_wh2o * getWetOpacity(spwId, i).get() ) );
1966 : //
1967 0 : f1 = transmission_fit[i];
1968 0 : psave = pfit_wh2o;
1969 0 : pfit_wh2o = pfit_wh2o + deltaa;
1970 0 : f2 = exp(-( (getDryContOpacity(spwId, i).get()+getO2LinesOpacity(spwId, i).get()+0.65*getO3LinesOpacity(spwId, i).get() ) //getDryOpacity(spwId, i).get()
1971 0 : + pfit_wh2o* getWetOpacity(spwId, i).get()));
1972 0 : deriv = (f2 - f1) / deltaa;
1973 0 : pfit_wh2o = psave;
1974 0 : beta = beta + (measuredSkyTransmission[i] - transmission_fit[i])
1975 0 : * deriv;
1976 0 : alpha = alpha + deriv * deriv;
1977 : }
1978 : }
1979 :
1980 0 : chisq1 = 0;
1981 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId).size(); i++) {
1982 0 : if(nl == getSpectralWindow(spwId).size() || (getSpectralWindow(spwId)[i]
1983 0 : * 1E-09 >= fre1.get(Frequency::UnitGigaHertz) && getSpectralWindow(spwId)[i] * 1E-09
1984 0 : <= fre2.get(Frequency::UnitGigaHertz))) {
1985 0 : res = -transmission_fit[i] + measuredSkyTransmission[i];
1986 0 : chisq1 = chisq1 + res * res;
1987 : }
1988 : }
1989 0 : if(nl > 1) {
1990 0 : chisq1 = chisq1 / (nl - 1);
1991 : }
1992 :
1993 0 : adjust: array = 1.0 / (1.0 + flamda);
1994 0 : pfit_wh2o_b = pfit_wh2o;
1995 0 : pfit_wh2o_b = pfit_wh2o_b + beta * array / alpha;
1996 0 : if(pfit_wh2o_b < 0.0) pfit_wh2o_b = 0.9 * pfit_wh2o;
1997 :
1998 0 : chisqr = 0;
1999 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId).size(); i++) {
2000 0 : if(nl == getSpectralWindow(spwId).size() || (getSpectralWindow(spwId)[i]
2001 0 : * 1E-09 >= fre1.get(Frequency::UnitGigaHertz) && getSpectralWindow(spwId)[i] * 1E-09
2002 0 : <= fre2.get(Frequency::UnitGigaHertz))) {
2003 0 : transmission_fit[i] = exp(-( ( getDryContOpacity(spwId, i).get()+getO2LinesOpacity(spwId, i).get()+0.65*getO3LinesOpacity(spwId, i).get() ) // (getDryOpacity(spwId, i).get())
2004 0 : + pfit_wh2o_b * getWetOpacity(spwId, i).get()));
2005 0 : res = -transmission_fit[i] + measuredSkyTransmission[i];
2006 0 : chisqr = chisqr + res * res;
2007 : }
2008 : }
2009 :
2010 0 : if(nl > 1) {
2011 0 : chisqr = chisqr / (nl - 1);
2012 : }
2013 :
2014 0 : if(fabs(chisq1 - chisqr) > 0.001) {
2015 0 : if(chisq1 < chisqr) {
2016 0 : flamda = flamda * 10.0;
2017 0 : goto adjust;
2018 : }
2019 : }
2020 :
2021 0 : flamda = flamda / 10.0;
2022 0 : sig_fit = sqrt(chisqr);
2023 0 : pfit_wh2o = pfit_wh2o_b;
2024 0 : sigma_wh2o = Length(sqrt(array / alpha) * sig_fit * pfit_wh2o
2025 0 : * (getUserWH2O().get()), Length::UnitMilliMeter);
2026 :
2027 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
2028 :
2029 : /* for(unsigned int i=0; i<getSpectralWindow(spwId).size(); i++){
2030 : cout << getSpectralWindow(spwId)[i]*1E-09 << " " << measuredSkyTransmission[i] << " " << transmission_fit[i] << endl;
2031 : } */
2032 :
2033 0 : sigma_fit_transm0 = sig_fit;
2034 :
2035 0 : wh2o_retrieved = Length(pfit_wh2o * getUserWH2O().get(Length::UnitMilliMeter), Length::UnitMilliMeter);
2036 :
2037 0 : goto salir;
2038 :
2039 : }
2040 :
2041 : }
2042 :
2043 0 : wh2o_retrieved = werr; // Extra error code, fit not reached after 20 iterations
2044 0 : sigma_fit_transm0 = -888.0; // Extra error code, fit not reached after 20 iterations
2045 0 : sigma_wh2o = werr; // Extra error code, fit not reached after 20 iterations
2046 :
2047 0 : salir:
2048 :
2049 0 : sigma_transmission_FTSfit_ = sigma_fit_transm0;
2050 :
2051 0 : if(wh2o_retrieved.get() > 0.0) {
2052 0 : wh2o_user_ = wh2o_retrieved;
2053 : }
2054 0 : return wh2o_retrieved;
2055 :
2056 0 : }
2057 :
2058 0 : Length SkyStatus::mkWaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
2059 : const vector<Percent> &signalGain,
2060 : const vector<Temperature> &measuredAverageSkyTEBB,
2061 : double airm,
2062 : const vector<vector<double> > &spwId_filter,
2063 : const vector<double> &skycoupling,
2064 : const vector<Temperature> &tspill)
2065 : {
2066 :
2067 : double pfit_wh2o;
2068 0 : double deltaa = 0.02;
2069 0 : double sig_fit = -999.0;
2070 0 : double eps = 0.001;
2071 0 : vector<Temperature> v_tebb_fit;
2072 :
2073 : unsigned int num;
2074 : double flamda;
2075 0 : unsigned int niter = 20;
2076 : double alpha;
2077 : double beta;
2078 : double array;
2079 : double f1;
2080 : double psave;
2081 : double f2;
2082 : double deriv;
2083 : double chisq1;
2084 : double chisqr;
2085 : double pfit_wh2o_b;
2086 : double res;
2087 0 : Length wh2o_retrieved(-999.0, Length::UnitMilliMeter);
2088 0 : Length werr(-888, Length::UnitMilliMeter);
2089 : double sigma_fit_transm0;
2090 0 : Length sigma_wh2o;
2091 0 : v_tebb_fit.reserve(measuredAverageSkyTEBB.size());
2092 :
2093 0 : num = 0;
2094 0 : flamda = 0.001;
2095 :
2096 0 : pfit_wh2o = (getUserWH2O().get(Length::UnitMilliMeter)) / (getGroundWH2O().get(Length::UnitMilliMeter));
2097 :
2098 0 : for(unsigned int kite = 0; kite < niter; kite++) {
2099 :
2100 0 : num = num + 1;
2101 :
2102 0 : beta = 0.0;
2103 0 : alpha = 0.0;
2104 :
2105 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
2106 :
2107 :
2108 0 : f1 = RT(pfit_wh2o,
2109 0 : skycoupling[j],
2110 0 : tspill[j].get(Temperature::UnitKelvin),
2111 : airm,
2112 0 : spwId[j],
2113 0 : spwId_filter[j],
2114 0 : signalGain[j]); // if signalGain[j] < 1.0 then there is an image side
2115 : // band and it is correctly taken into account in RT
2116 :
2117 0 : v_tebb_fit[j] = Temperature(f1, Temperature::UnitKelvin);
2118 :
2119 0 : psave = pfit_wh2o;
2120 0 : pfit_wh2o = pfit_wh2o + deltaa;
2121 :
2122 0 : f2 = RT(pfit_wh2o,
2123 0 : skycoupling[j],
2124 0 : tspill[j].get(Temperature::UnitKelvin),
2125 : airm,
2126 0 : spwId[j],
2127 0 : spwId_filter[j],
2128 0 : signalGain[j]);
2129 :
2130 0 : deriv = (f2 - f1) / deltaa;
2131 :
2132 0 : pfit_wh2o = psave;
2133 0 : beta = beta + ((measuredAverageSkyTEBB[j]).get(Temperature::UnitKelvin) - f1) * deriv;
2134 0 : alpha = alpha + deriv * deriv;
2135 :
2136 : }
2137 :
2138 0 : chisq1 = 0;
2139 :
2140 0 : for(unsigned int j = 0; j < measuredAverageSkyTEBB.size(); j++) {
2141 :
2142 0 : res = -v_tebb_fit[j].get(Temperature::UnitKelvin) + (measuredAverageSkyTEBB[j]).get(Temperature::UnitKelvin);
2143 0 : chisq1 = chisq1 + res * res;
2144 :
2145 : }
2146 :
2147 0 : if(measuredAverageSkyTEBB.size() > 1) {
2148 0 : chisq1 = chisq1 / (measuredAverageSkyTEBB.size() - 1);
2149 : }
2150 :
2151 0 : adjust: array = 1.0 / (1.0 + flamda);
2152 0 : pfit_wh2o_b = pfit_wh2o;
2153 0 : pfit_wh2o_b = pfit_wh2o_b + beta * array / alpha;
2154 0 : if(pfit_wh2o_b < 0.0) pfit_wh2o_b = 0.9 * pfit_wh2o;
2155 :
2156 0 : chisqr = 0;
2157 :
2158 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
2159 :
2160 0 : v_tebb_fit[j] = Temperature(RT(pfit_wh2o_b,
2161 0 : skycoupling[j],
2162 0 : tspill[j].get(Temperature::UnitKelvin),
2163 : airm,
2164 0 : spwId[j],
2165 0 : spwId_filter[j],
2166 0 : signalGain[j]), Temperature::UnitKelvin);
2167 :
2168 0 : res = -v_tebb_fit[j].get(Temperature::UnitKelvin) + (measuredAverageSkyTEBB[j]).get(Temperature::UnitKelvin);
2169 0 : chisqr = chisqr + res * res;
2170 :
2171 : }
2172 :
2173 0 : if(spwId.size() > 1) {
2174 0 : chisqr = chisqr / (spwId.size() - 1);
2175 : }
2176 :
2177 0 : if(fabs(chisq1 - chisqr) > 0.001) {
2178 0 : if(chisq1 < chisqr) {
2179 0 : flamda = flamda * 10.0;
2180 0 : goto adjust;
2181 : }
2182 : }
2183 :
2184 0 : flamda = flamda / 10.0;
2185 0 : sig_fit = sqrt(chisqr);
2186 0 : pfit_wh2o = pfit_wh2o_b;
2187 0 : sigma_wh2o = Length(sqrt(array / alpha) * sig_fit * pfit_wh2o
2188 0 : * (getGroundWH2O().get(Length::UnitMilliMeter)), Length::UnitMilliMeter);
2189 :
2190 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
2191 :
2192 0 : sigma_fit_transm0 = sig_fit;
2193 :
2194 0 : wh2o_retrieved = Length(pfit_wh2o * getGroundWH2O().get(Length::UnitMilliMeter), Length::UnitMilliMeter);
2195 :
2196 0 : goto salir;
2197 :
2198 : }
2199 :
2200 : }
2201 :
2202 0 : wh2o_retrieved = werr; // Extra error code, fit not reached after 20 iterations
2203 0 : sigma_fit_transm0 = -888.0; // Extra error code, fit not reached after 20 iterations
2204 0 : sigma_wh2o = werr; // Extra error code, fit not reached after 20 iterations
2205 :
2206 0 : salir:
2207 :
2208 0 : sigma_TEBBfit_ = Temperature(sigma_fit_transm0, Temperature::UnitKelvin);
2209 0 : if(wh2o_retrieved.get() > 0.0) {
2210 0 : wh2o_user_ = wh2o_retrieved;
2211 : }
2212 0 : return wh2o_retrieved;
2213 :
2214 0 : }
2215 :
2216 0 : Length SkyStatus::mkWaterVaporRetrieval_fromTEBB(const vector<unsigned int> &spwId,
2217 : const vector<Percent> &signalGain,
2218 : const vector<vector<Temperature> > &measuredSkyTEBB,
2219 : double airm,
2220 : const vector<vector<double> > &spwId_filter,
2221 : const vector<double> &skycoupling,
2222 : const vector<Temperature> &tspill)
2223 : {
2224 :
2225 : double pfit_wh2o;
2226 0 : double deltaa = 0.02;
2227 0 : double sig_fit = -999.0;
2228 0 : double eps = 0.01;
2229 0 : vector<vector<Temperature> > vv_tebb_fit;
2230 0 : vector<Temperature> v_tebb_fit;
2231 :
2232 : unsigned int num;
2233 : double flamda;
2234 0 : unsigned int niter = 20;
2235 : double alpha;
2236 : double beta;
2237 : double array;
2238 : double f1;
2239 : double psave;
2240 : double f2;
2241 : double deriv;
2242 : double chisq1;
2243 : double chisqr;
2244 : double pfit_wh2o_b;
2245 : double res;
2246 0 : Length wh2o_retrieved(-999.0, Length::UnitMilliMeter);
2247 0 : Length werr(-888, Length::UnitMilliMeter);
2248 : double sigma_fit_transm0;
2249 0 : Length sigma_wh2o;
2250 0 : vector<double> spwIdNorm;
2251 0 : vector<double> validchannels;
2252 0 : spwIdNorm.reserve(spwId.size());
2253 0 : validchannels.reserve(spwId.size());
2254 0 : vv_tebb_fit.reserve(spwId.size());
2255 :
2256 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
2257 0 : spwIdNorm[j] = 0.0;
2258 0 : validchannels[j] = 0.0;
2259 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId[j]).size(); i++) {
2260 0 : if(spwId_filter[j][i] > 0) {
2261 0 : spwIdNorm[j] = spwIdNorm[j] + spwId_filter[j][i];
2262 0 : validchannels[j] = validchannels[j] + 1.0;
2263 : }
2264 : }
2265 : }
2266 :
2267 0 : num = 0;
2268 0 : flamda = 0.001;
2269 0 : pfit_wh2o = (getUserWH2O().get(Length::UnitMilliMeter)) / (getGroundWH2O().get(Length::UnitMilliMeter));
2270 :
2271 0 : for(unsigned int kite = 0; kite < niter; kite++) {
2272 :
2273 0 : num = num + 1;
2274 :
2275 0 : beta = 0.0;
2276 0 : alpha = 0.0;
2277 :
2278 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
2279 :
2280 0 : v_tebb_fit.clear();
2281 :
2282 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId[j]).size(); i++) {
2283 :
2284 0 : if(spwId_filter[j][i] > 0) {
2285 :
2286 0 : if(signalGain[j].get() < 1.0) {
2287 0 : f1 = RT(pfit_wh2o,
2288 0 : skycoupling[j],
2289 0 : tspill[j].get(Temperature::UnitKelvin),
2290 : airm,
2291 0 : spwId[j],
2292 0 : i) * signalGain[j].get() + RT(pfit_wh2o,
2293 0 : skycoupling[j],
2294 0 : tspill[j].get(Temperature::UnitKelvin),
2295 : airm,
2296 0 : getAssocSpwId(spwId[j])[0],
2297 : i)
2298 0 : * (1 - signalGain[j].get());
2299 : } else {
2300 0 : f1 = RT(pfit_wh2o,
2301 0 : skycoupling[j],
2302 0 : tspill[j].get(Temperature::UnitKelvin),
2303 : airm,
2304 0 : spwId[j],
2305 : i);
2306 : }
2307 :
2308 0 : psave = pfit_wh2o;
2309 0 : pfit_wh2o = pfit_wh2o + deltaa;
2310 :
2311 0 : v_tebb_fit.push_back(Temperature(f1, Temperature::UnitKelvin));
2312 :
2313 0 : if(signalGain[j].get() < 1.0) {
2314 0 : f2 = (RT(pfit_wh2o,
2315 0 : skycoupling[j],
2316 0 : tspill[j].get(Temperature::UnitKelvin),
2317 : airm,
2318 0 : spwId[j],
2319 0 : i) * signalGain[j].get() + RT(pfit_wh2o,
2320 0 : skycoupling[j],
2321 0 : tspill[j].get(Temperature::UnitKelvin),
2322 : airm,
2323 0 : getAssocSpwId(spwId[j])[0],
2324 0 : i) * (1
2325 0 : - signalGain[j].get()));
2326 : } else {
2327 0 : f2 = (RT(pfit_wh2o,
2328 0 : skycoupling[j],
2329 0 : tspill[j].get(Temperature::UnitKelvin),
2330 : airm,
2331 0 : spwId[j],
2332 0 : i) * signalGain[j].get());
2333 : }
2334 :
2335 0 : f2 = f2 * spwId_filter[j][i] * validchannels[j] / spwIdNorm[j];
2336 0 : f1 = f1 * spwId_filter[j][i] * validchannels[j] / spwIdNorm[j];
2337 0 : deriv = (f2 - f1) / deltaa;
2338 :
2339 0 : pfit_wh2o = psave;
2340 0 : beta = beta + ((measuredSkyTEBB[j][i]).get(Temperature::UnitKelvin) - f1) * deriv; //*spwId_filter[j][i]/spwIdNorm[j];
2341 0 : alpha = alpha + deriv * deriv;
2342 :
2343 : } else {
2344 :
2345 0 : v_tebb_fit.push_back(Temperature(-999, Temperature::UnitKelvin));
2346 :
2347 : }
2348 :
2349 : }
2350 :
2351 0 : vv_tebb_fit.push_back(v_tebb_fit);
2352 :
2353 : }
2354 :
2355 0 : chisq1 = 0;
2356 :
2357 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
2358 :
2359 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId[j]).size(); i++) {
2360 0 : if(spwId_filter[j][i] > 0) {
2361 : res
2362 0 : = (-vv_tebb_fit[j][i].get(Temperature::UnitKelvin) + (measuredSkyTEBB[j][i]).get(Temperature::UnitKelvin))
2363 0 : * spwId_filter[j][i] * validchannels[j] / spwIdNorm[j];
2364 0 : chisq1 = chisq1 + res * res;
2365 : }
2366 : }
2367 :
2368 : }
2369 :
2370 0 : if(spwId.size() > 1) {
2371 0 : chisq1 = chisq1 / (spwId.size() - 1);
2372 : }
2373 :
2374 0 : adjust: array = 1.0 / (1.0 + flamda);
2375 0 : pfit_wh2o_b = pfit_wh2o;
2376 0 : pfit_wh2o_b = pfit_wh2o_b + beta * array / alpha;
2377 0 : if(pfit_wh2o_b < 0.0) pfit_wh2o_b = 0.9 * pfit_wh2o;
2378 :
2379 0 : chisqr = 0;
2380 :
2381 0 : for(unsigned int j = 0; j < spwId.size(); j++) {
2382 :
2383 0 : for(unsigned int i = 0; i < getSpectralWindow(spwId[j]).size(); i++) {
2384 :
2385 0 : if(spwId_filter[j][i] > 0) {
2386 :
2387 0 : if(signalGain[j].get() < 1.0) {
2388 0 : vv_tebb_fit[j][i] = Temperature(RT(pfit_wh2o_b,
2389 0 : skycoupling[j],
2390 0 : tspill[j].get(Temperature::UnitKelvin),
2391 : airm,
2392 0 : spwId[j],
2393 0 : i) * signalGain[j].get()
2394 0 : + RT(pfit_wh2o_b,
2395 0 : skycoupling[j],
2396 0 : tspill[j].get(Temperature::UnitKelvin),
2397 : airm,
2398 0 : getAssocSpwId(spwId[j])[0],
2399 0 : i) * (1 - signalGain[j].get()), Temperature::UnitKelvin);
2400 : } else {
2401 0 : vv_tebb_fit[j][i] = Temperature(RT(pfit_wh2o_b,
2402 0 : skycoupling[j],
2403 0 : tspill[j].get(Temperature::UnitKelvin),
2404 : airm,
2405 0 : spwId[j],
2406 0 : i) * signalGain[j].get(), Temperature::UnitKelvin);
2407 : }
2408 :
2409 : res
2410 0 : = (-vv_tebb_fit[j][i].get(Temperature::UnitKelvin) + (measuredSkyTEBB[j][i]).get(Temperature::UnitKelvin))
2411 0 : * spwId_filter[j][i] * validchannels[j] / spwIdNorm[j];
2412 :
2413 0 : chisqr = chisqr + res * res;
2414 : }
2415 : }
2416 :
2417 : }
2418 :
2419 0 : if(spwId.size() > 1) {
2420 0 : chisqr = chisqr / (spwId.size() - 1);
2421 : }
2422 :
2423 0 : if(fabs(chisq1 - chisqr) > 0.001) {
2424 0 : if(chisq1 < chisqr) {
2425 0 : flamda = flamda * 10.0;
2426 0 : goto adjust;
2427 : }
2428 : }
2429 :
2430 0 : flamda = flamda / 10.0;
2431 0 : sig_fit = sqrt(chisqr);
2432 0 : pfit_wh2o = pfit_wh2o_b;
2433 0 : sigma_wh2o = Length(sqrt(array / alpha) * sig_fit * pfit_wh2o
2434 0 : * (getGroundWH2O().get(Length::UnitMilliMeter)), Length::UnitMilliMeter);
2435 :
2436 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
2437 :
2438 0 : sigma_fit_transm0 = sig_fit;
2439 :
2440 0 : wh2o_retrieved = Length(pfit_wh2o * getGroundWH2O().get(Length::UnitMilliMeter), Length::UnitMilliMeter);
2441 :
2442 0 : goto salir;
2443 :
2444 : }
2445 :
2446 : }
2447 :
2448 0 : wh2o_retrieved = werr; // Extra error code, fit not reached after 20 iterations
2449 0 : sigma_fit_transm0 = -888.0; // Extra error code, fit not reached after 20 iterations
2450 0 : sigma_wh2o = werr; // Extra error code, fit not reached after 20 iterations
2451 :
2452 0 : salir:
2453 :
2454 0 : sigma_TEBBfit_ = Temperature(sigma_fit_transm0, Temperature::UnitKelvin);
2455 0 : if(wh2o_retrieved.get() > 0.0) {
2456 0 : wh2o_user_ = wh2o_retrieved;
2457 : }
2458 0 : return wh2o_retrieved;
2459 :
2460 0 : }
2461 :
2462 0 : double SkyStatus::mkSkyCouplingRetrieval_fromTEBB(unsigned int spwId,
2463 : const Percent &signalGain,
2464 : const vector<Temperature> &measuredSkyTEBB,
2465 : double airm,
2466 : const vector<double> &spwId_filter,
2467 : double skycoupling,
2468 : const Temperature &tspill)
2469 : {
2470 :
2471 : // double pfit_wh2o; // [-Wunused_but_set_variable]
2472 : double pfit_skycoupling;
2473 : double pfit_skycoupling_b;
2474 0 : double deltaa = 0.02;
2475 : // double sig_fit = -999.0; // [-Wunused_but_set_variable]
2476 0 : double eps = 0.01;
2477 0 : vector<Temperature> tebb_fit;
2478 0 : tebb_fit.reserve(measuredSkyTEBB.size()); // allows resizing
2479 :
2480 : unsigned int num;
2481 : double flamda;
2482 0 : unsigned int niter = 20;
2483 : double alpha;
2484 : double beta;
2485 : double array;
2486 : double f1;
2487 : double psave;
2488 : double f2;
2489 : double deriv;
2490 : double chisq1;
2491 : double chisqr;
2492 : double res;
2493 0 : Length wh2o_retrieved(-999.0, Length::UnitMilliMeter);
2494 0 : Length werr(-888, Length::UnitMilliMeter);
2495 : // double sigma_fit_transm0; // [-Wunused_but_set_variable]
2496 0 : Length sigma_wh2o;
2497 :
2498 0 : num = 0;
2499 :
2500 0 : flamda = 0.001;
2501 : // pfit_wh2o = (getUserWH2O().get(Length::UnitMilliMeter)) / (getGroundWH2O().get(Length::UnitMilliMeter)); // [-Wunused_but_set_variable]
2502 0 : pfit_skycoupling = 1.0;
2503 :
2504 0 : for(unsigned int kite = 0; kite < niter; kite++) {
2505 :
2506 0 : num = num + 1;
2507 :
2508 0 : beta = 0.0;
2509 0 : alpha = 0.0;
2510 :
2511 : // for(unsigned int i=0; i<getSpectralWindow(spwId).size(); i++){
2512 :
2513 0 : mkWaterVaporRetrieval_fromTEBB(spwId,
2514 : signalGain,
2515 : measuredSkyTEBB,
2516 : airm,
2517 : spwId_filter,
2518 : pfit_skycoupling * skycoupling,
2519 : tspill);
2520 0 : f1 = sigma_TEBBfit_.get(Temperature::UnitKelvin);
2521 : // cout << "pfit_skycoupling =" << pfit_skycoupling << " f1=" << f1 << " wh2o_user_=" << wh2o_user_.get(Length::UnitMilliMeter) << " mm" << endl;
2522 0 : psave = pfit_skycoupling;
2523 0 : pfit_skycoupling = pfit_skycoupling + deltaa;
2524 0 : mkWaterVaporRetrieval_fromTEBB(spwId,
2525 : signalGain,
2526 : measuredSkyTEBB,
2527 : airm,
2528 : spwId_filter,
2529 : pfit_skycoupling * skycoupling,
2530 : tspill);
2531 0 : f2 = sigma_TEBBfit_.get(Temperature::UnitKelvin);
2532 : // cout << "pfit_skycoupling =" << pfit_skycoupling << " f2=" << f2 << " wh2o_user_=" << wh2o_user_.get(Length::UnitMilliMeter) << " mm" << endl;
2533 0 : deriv = (f2 - f1) / deltaa;
2534 0 : pfit_skycoupling = psave;
2535 0 : beta = beta - f1 * deriv;
2536 0 : alpha = alpha + deriv * deriv;
2537 :
2538 0 : chisq1 = 0;
2539 0 : res = f1;
2540 0 : chisq1 = chisq1 + res * res;
2541 :
2542 0 : adjust: array = 1.0 / (1.0 + flamda);
2543 0 : pfit_skycoupling_b = pfit_skycoupling;
2544 0 : pfit_skycoupling_b = pfit_skycoupling_b + beta * array / alpha;
2545 0 : if(pfit_skycoupling_b < 0.0) pfit_skycoupling_b = 0.9 * pfit_skycoupling;
2546 0 : chisqr = 0;
2547 0 : mkWaterVaporRetrieval_fromTEBB(spwId,
2548 : signalGain,
2549 : measuredSkyTEBB,
2550 : airm,
2551 : spwId_filter,
2552 : pfit_skycoupling_b * skycoupling,
2553 : tspill);
2554 0 : res = sigma_TEBBfit_.get(Temperature::UnitKelvin);
2555 : // cout << "pfit_skycoupling =" << pfit_skycoupling_b << " res=" << res << " wh2o_user_=" << wh2o_user_.get(Length::UnitMilliMeter) << " mm" << endl;
2556 0 : chisqr = chisqr + res * res;
2557 :
2558 0 : if(fabs(chisq1 - chisqr) > 0.001) {
2559 0 : if(chisq1 < chisqr) {
2560 0 : flamda = flamda * 10.0;
2561 0 : goto adjust;
2562 : }
2563 : }
2564 :
2565 0 : flamda = flamda / 10.0;
2566 : // sig_fit = sqrt(chisqr); // [-Wunused_but_set_variable]
2567 0 : pfit_skycoupling = pfit_skycoupling_b;
2568 :
2569 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
2570 0 : goto salir;
2571 : }
2572 :
2573 : }
2574 :
2575 0 : wh2o_retrieved = werr; // Extra error code, fit not reached after 20 iterations
2576 : //sigma_fit_transm0 = -888.0; // Extra error code, fit not reached after 20 iterations [-Wunused_but_set_variable]
2577 0 : sigma_wh2o = werr; // Extra error code, fit not reached after 20 iterations
2578 :
2579 0 : salir: return pfit_skycoupling * skycoupling;
2580 :
2581 0 : }
2582 :
2583 :
2584 0 : double SkyStatus::RT(double pfit_wh2o,
2585 : double skycoupling,
2586 : double tspill,
2587 : double airmass,
2588 : unsigned int spwid,
2589 : const vector<double> &spwId_filter,
2590 : const Percent &signalgain)
2591 : {
2592 :
2593 0 : double tebb_channel = 0.0;
2594 : double rtr;
2595 0 : double norm = 0.0;
2596 :
2597 0 : for(unsigned int n = 0; n < v_numChan_[spwid]; n++) {
2598 0 : if(spwId_filter[n] > 0) {
2599 0 : norm = norm + spwId_filter[n];
2600 : }
2601 : }
2602 :
2603 0 : if(norm == 0.0) {
2604 0 : return norm;
2605 : }
2606 :
2607 0 : for(unsigned int n = 0; n < v_numChan_[spwid]; n++) {
2608 :
2609 0 : if(spwId_filter[n] > 0) {
2610 :
2611 0 : if(signalgain.get() < 1.0) {
2612 0 : rtr = RT(pfit_wh2o, skycoupling, tspill, airmass, spwid, n)
2613 0 : * signalgain.get() + RT(pfit_wh2o,
2614 : skycoupling,
2615 : tspill,
2616 : airmass,
2617 0 : getAssocSpwId(spwid)[0],
2618 0 : n) * (1.0 - signalgain.get());
2619 : } else {
2620 0 : rtr = RT(pfit_wh2o, skycoupling, tspill, airmass, spwid, n);
2621 : }
2622 0 : tebb_channel = tebb_channel + rtr * spwId_filter[n] / norm;
2623 : }
2624 : }
2625 :
2626 0 : return tebb_channel;
2627 : }
2628 :
2629 :
2630 :
2631 11242 : double SkyStatus::RT(double pfit_wh2o,
2632 : double skycoupling,
2633 : double tspill,
2634 : double airm,
2635 : unsigned int spwid,
2636 : unsigned int nc)
2637 : {
2638 :
2639 : double radiance;
2640 : double singlefreq;
2641 : // double chanwidth; // [-Wunused_but_set_variable]
2642 : double tebb;
2643 11242 : double h_div_k = 0.04799274551; /* plank=6.6262e-34,boltz=1.3806E-23 */
2644 : double kv;
2645 : double tau_layer;
2646 11242 : double tbgr = skyBackgroundTemperature_.get(Temperature::UnitKelvin);
2647 11242 : double ratioWater = pfit_wh2o;
2648 :
2649 11242 : singlefreq = getChanFreq(spwid, nc).get(Frequency::UnitGigaHertz);
2650 : // chanwidth = getChanWidth(spwid, nc).get(Frequency::UnitGigaHertz); // [-Wunused_but_set_variable]
2651 : // cout << "Chan freq. =" << singlefreq << " GHz" << endl;
2652 : // cout << "Chan width =" << chanwidth << " GHz" << endl;
2653 :
2654 11242 : tebb=0.0;
2655 11242 : kv = 0.0;
2656 11242 : radiance = 0.0;
2657 :
2658 236222 : for(unsigned int i = 0; i < getNumLayer(); i++) {
2659 :
2660 224980 : tau_layer = ((getAbsTotalWet(spwid, nc, i).get()) * ratioWater+ getAbsTotalDry(spwid, nc, i).get()) * getLayerThickness(i).get();
2661 :
2662 : // cout << i << " " << getAbsTotalWet(spwid, nc, i).get() << " " << getAbsTotalDry(spwid, nc, i).get() << endl;
2663 :
2664 224980 : radiance = radiance + (1.0 / (exp(h_div_k * singlefreq/ getLayerTemperature(i).get()) - 1.0)) * exp(-kv * airm) * (1.0- exp(-airm * tau_layer));
2665 :
2666 224980 : kv = kv + tau_layer;
2667 :
2668 : }
2669 :
2670 11242 : radiance = skycoupling * (radiance + (1.0 / (exp(h_div_k * singlefreq / tbgr)- 1.0)) * exp(-kv * airm)) + (1.0 / (exp(h_div_k * singlefreq / tspill)- 1.0)) * (1 - skycoupling);
2671 :
2672 11242 : tebb = h_div_k * singlefreq / log(1 + (1 / radiance));
2673 : // tebb = tebb+ h_div_k * singlefreq * radiance;
2674 : // cout << "singlefreq = " << singlefreq << " total opacity = " << kv << " tebb = " << tebb << endl;
2675 :
2676 11242 : return tebb;
2677 :
2678 : }
2679 :
2680 :
2681 0 : double SkyStatus::RTRJ(double pfit_wh2o,
2682 : double skycoupling,
2683 : double tspill,
2684 : double airm,
2685 : unsigned int spwid,
2686 : unsigned int nc)
2687 : {
2688 :
2689 : double radiance;
2690 : double singlefreq;
2691 : // double chanwidth; // [-Wunused_but_set_variable]
2692 : double trj;
2693 0 : double h_div_k = 0.04799274551; /* plank=6.6262e-34,boltz=1.3806E-23 */
2694 : double kv;
2695 : double tau_layer;
2696 0 : double tbgr = skyBackgroundTemperature_.get(Temperature::UnitKelvin);
2697 0 : double ratioWater = pfit_wh2o;
2698 :
2699 0 : singlefreq = getChanFreq(spwid, nc).get(Frequency::UnitGigaHertz);
2700 : // chanwidth = getChanWidth(spwid, nc).get(Frequency::UnitGigaHertz); // [-Wunused_but_set_variable]
2701 0 : trj=0.0;
2702 0 : kv = 0.0;
2703 0 : radiance = 0.0;
2704 :
2705 0 : for(unsigned int i = 0; i < getNumLayer(); i++) {
2706 :
2707 0 : tau_layer = ((getAbsTotalWet(spwid, nc, i).get()) * ratioWater+ getAbsTotalDry(spwid, nc, i).get()) * getLayerThickness(i).get();
2708 :
2709 0 : radiance = radiance + (1.0 / (exp(h_div_k * singlefreq/ getLayerTemperature(i).get()) - 1.0)) * exp(-kv * airm) * (1.0- exp(-airm * tau_layer));
2710 :
2711 0 : kv = kv + tau_layer;
2712 :
2713 : }
2714 :
2715 0 : radiance = skycoupling * (radiance + (1.0 / (exp(h_div_k * singlefreq / tbgr)- 1.0)) * exp(-kv * airm)) + (1.0 / (exp(h_div_k * singlefreq / tspill)- 1.0)) * (1 - skycoupling);
2716 :
2717 0 : trj = h_div_k * singlefreq * radiance;
2718 :
2719 0 : return trj;
2720 :
2721 : }
2722 :
2723 :
2724 :
2725 18 : void SkyStatus::iniSkyStatus()
2726 : {
2727 :
2728 18 : Length wh2o_default(1, Length::UnitMilliMeter);
2729 18 : Length wh2o_default_neg(-999, Length::UnitMilliMeter);
2730 18 : Temperature temp_default_neg(-999, Temperature::UnitKelvin);
2731 :
2732 36 : if(wh2o_user_.get() <= 0.0 || wh2o_user_.get() > (getGroundWH2O().get())
2733 18 : * (200 / (getRelativeHumidity().get(Percent::UnitPercent)))) {
2734 18 : wh2o_user_ = wh2o_default;
2735 : }
2736 :
2737 18 : }
2738 :
2739 0 : double SkyStatus::getAverageNonDispersiveDryPathLength_GroundPressureDerivative(unsigned int spwid)
2740 : {
2741 0 : Pressure ref = getGroundPressure();
2742 : Length a =
2743 0 : RefractiveIndexProfile::getAverageNonDispersiveDryPathLength(spwid);
2744 0 : setBasicAtmosphericParameters(ref + Pressure(1.0, Pressure::UnitMilliBar));
2745 : Length b =
2746 0 : RefractiveIndexProfile::getAverageNonDispersiveDryPathLength(spwid);
2747 0 : setBasicAtmosphericParameters(ref);
2748 0 : return b.get(Length::UnitMicrons) - a.get(Length::UnitMicrons);
2749 0 : }
2750 :
2751 0 : double SkyStatus::getAverageNonDispersiveDryPathLength_GroundTemperatureDerivative(unsigned int spwid)
2752 : {
2753 0 : Temperature ref = getGroundTemperature();
2754 0 : double oldLapseRate = getTropoLapseRate();
2755 : Length a =
2756 0 : RefractiveIndexProfile::getAverageNonDispersiveDryPathLength(spwid);
2757 : double newLapseRate =
2758 0 : ((getTropopauseTemperature() - ref - Temperature(1.0, Temperature::UnitKelvin)).get(Temperature::UnitKelvin))
2759 0 : / ((getTropopauseAltitude() - getAltitude()).get(Length::UnitKiloMeter));
2760 0 : setBasicAtmosphericParameters(ref + Temperature(1.0, Temperature::UnitKelvin), newLapseRate);
2761 : Length b =
2762 0 : RefractiveIndexProfile::getAverageNonDispersiveDryPathLength(spwid);
2763 0 : setBasicAtmosphericParameters(ref, oldLapseRate);
2764 0 : return b.get(Length::UnitMicrons) - a.get(Length::UnitMicrons);
2765 0 : }
2766 :
2767 0 : double SkyStatus::getAverageDispersiveDryPathLength_GroundPressureDerivative(unsigned int spwid)
2768 : {
2769 0 : Pressure ref = getGroundPressure();
2770 0 : Length a = RefractiveIndexProfile::getAverageDispersiveDryPathLength(spwid);
2771 : // scanf("%d",&e);
2772 0 : setBasicAtmosphericParameters(ref + Pressure(1.0, Pressure::UnitMilliBar));
2773 0 : Length b = RefractiveIndexProfile::getAverageDispersiveDryPathLength(spwid);
2774 0 : setBasicAtmosphericParameters(ref);
2775 0 : return b.get(Length::UnitMicrons) - a.get(Length::UnitMicrons);
2776 0 : }
2777 :
2778 0 : double SkyStatus::getAverageDispersiveDryPathLength_GroundTemperatureDerivative(unsigned int spwid)
2779 : {
2780 0 : Temperature ref = getGroundTemperature();
2781 0 : double oldLapseRate = getTropoLapseRate();
2782 0 : Length a = RefractiveIndexProfile::getAverageDispersiveDryPathLength(spwid);
2783 : double newLapseRate =
2784 0 : ((getTropopauseTemperature() - ref - Temperature(1.0, Temperature::UnitKelvin)).get(Temperature::UnitKelvin))
2785 0 : / ((getTropopauseAltitude() - getAltitude()).get(Length::UnitKiloMeter));
2786 0 : setBasicAtmosphericParameters(ref + Temperature(1.0, Temperature::UnitKelvin), newLapseRate);
2787 0 : Length b = RefractiveIndexProfile::getAverageDispersiveDryPathLength(spwid);
2788 0 : setBasicAtmosphericParameters(ref, oldLapseRate);
2789 0 : return b.get(Length::UnitMicrons) - a.get(Length::UnitMicrons);
2790 0 : }
2791 :
2792 0 : Temperature SkyStatus::getWVRAverageSigmaTskyFit(const vector<WVRMeasurement> &RadiometerData,
2793 : unsigned int n,
2794 : unsigned int m)
2795 : {
2796 0 : double sigma = 0.0;
2797 : double tr;
2798 0 : Temperature sigmaT;
2799 0 : if(m < n) {
2800 0 : return Temperature(-999, Temperature::UnitKelvin);
2801 : }
2802 0 : for(unsigned int i = n; i < m; i++) {
2803 0 : tr = RadiometerData[i].getSigmaFit().get(Temperature::UnitKelvin);
2804 0 : if(tr < 0) {
2805 0 : return Temperature(-999, Temperature::UnitKelvin);
2806 : }
2807 0 : sigma = sigma + tr * tr;
2808 : }
2809 0 : if(m == n) {
2810 0 : sigmaT = RadiometerData[n].getSigmaFit();
2811 : } else {
2812 0 : sigma = sqrt(sigma / (m - n));
2813 0 : sigmaT = Temperature(sigma, Temperature::UnitKelvin);
2814 : }
2815 0 : return sigmaT;
2816 0 : }
2817 0 : Temperature SkyStatus::getWVRSigmaChannelTskyFit(const vector<WVRMeasurement> &RadiometerData,
2818 : unsigned int ichan,
2819 : unsigned int n,
2820 : unsigned int m)
2821 : {
2822 0 : double sigma = 0.0;
2823 : double dtr;
2824 0 : Temperature sigmaT;
2825 0 : if(m <= n) {
2826 0 : return Temperature(-999, Temperature::UnitKelvin);
2827 : }
2828 0 : for(unsigned int i = n; i < m; i++) {
2829 0 : dtr = RadiometerData[i].getmeasuredSkyBrightness()[ichan].get(Temperature::UnitKelvin)
2830 0 : -RadiometerData[i].getfittedSkyBrightness()[ichan].get(Temperature::UnitKelvin);
2831 0 : sigma = sigma + dtr * dtr;
2832 : }
2833 0 : sigma = sqrt(sigma / (m - n));
2834 0 : sigmaT = Temperature(sigma, Temperature::UnitKelvin);
2835 0 : return sigmaT;
2836 0 : }
2837 :
2838 0 : WVRMeasurement SkyStatus::mkWaterVaporRetrieval_fromWVR(const vector<Temperature> &measuredSkyBrightnessVector,
2839 : const vector<unsigned int> &IdChannels,
2840 : const vector<double> &skyCoupling,
2841 : const vector<Percent> &signalGain,
2842 : const Temperature &spilloverTemperature,
2843 : const Angle &elevation)
2844 : {
2845 0 : double tspill = spilloverTemperature.get(Temperature::UnitKelvin);
2846 : double pfit_wh2o;
2847 0 : double deltaa = 0.02;
2848 0 : double sig_fit = -999.0;
2849 0 : double eps = 0.01;
2850 0 : vector<double> tebb_fit;
2851 0 : tebb_fit.reserve(measuredSkyBrightnessVector.size());
2852 0 : double airm = 1.0 / sin((3.1415926 * elevation.get(Angle::UnitDegree)) / 180.0);
2853 : unsigned int num;
2854 : double flamda;
2855 0 : unsigned int niter = 20;
2856 : double alpha;
2857 : double beta;
2858 : double array;
2859 : double f1;
2860 : double psave;
2861 : double f2;
2862 : double deriv;
2863 : double chisq1;
2864 : double chisqr;
2865 : double pfit_wh2o_b;
2866 : double res;
2867 0 : Length wh2o_retrieved(-999.0, Length::UnitMilliMeter);
2868 0 : Length werr(-888, Length::UnitMilliMeter);
2869 0 : Temperature sigma_fit_temp0;
2870 0 : Temperature t_astro;
2871 0 : Length sigma_wh2o;
2872 :
2873 0 : num = 0;
2874 :
2875 0 : flamda = 0.001;
2876 :
2877 0 : pfit_wh2o = (getUserWH2O().get(Length::UnitMilliMeter)) / (getGroundWH2O().get(Length::UnitMilliMeter));
2878 :
2879 : // cout << "pfit_wh2o=" << pfit_wh2o << endl;
2880 :
2881 0 : for(unsigned int kite = 0; kite < niter; kite++) {
2882 :
2883 0 : num = num + 1;
2884 :
2885 0 : beta = 0.0;
2886 0 : alpha = 0.0;
2887 :
2888 0 : for(unsigned int i = 0; i < IdChannels.size(); i++) {
2889 :
2890 0 : tebb_fit[i] = RT(pfit_wh2o,
2891 0 : skyCoupling[i],
2892 : tspill,
2893 : airm,
2894 0 : IdChannels[i],
2895 0 : signalGain[i]);
2896 : // cout << i << " " << tebb_fit[i] << endl;
2897 :
2898 :
2899 0 : f1 = tebb_fit[i];
2900 0 : psave = pfit_wh2o;
2901 0 : pfit_wh2o = pfit_wh2o + deltaa;
2902 0 : f2 = RT(pfit_wh2o,
2903 0 : skyCoupling[i],
2904 : tspill,
2905 : airm,
2906 0 : IdChannels[i],
2907 0 : signalGain[i]);
2908 0 : deriv = (f2 - f1) / deltaa;
2909 0 : pfit_wh2o = psave;
2910 0 : beta = beta + (measuredSkyBrightnessVector[i].get(Temperature::UnitKelvin) - tebb_fit[i])
2911 0 : * deriv;
2912 :
2913 0 : alpha = alpha + deriv * deriv;
2914 :
2915 : }
2916 :
2917 0 : chisq1 = 0;
2918 0 : for(unsigned int i = 0; i < measuredSkyBrightnessVector.size(); i++) {
2919 0 : res = -tebb_fit[i] + measuredSkyBrightnessVector[i].get(Temperature::UnitKelvin);
2920 0 : chisq1 = chisq1 + res * res;
2921 : }
2922 0 : if(measuredSkyBrightnessVector.size() > 1) {
2923 0 : chisq1 = chisq1 / (measuredSkyBrightnessVector.size() - 1);
2924 : }
2925 :
2926 0 : adjust: array = 1.0 / (1.0 + flamda);
2927 0 : pfit_wh2o_b = pfit_wh2o;
2928 0 : pfit_wh2o_b = pfit_wh2o_b + beta * array / alpha;
2929 0 : if(pfit_wh2o_b < 0.0) pfit_wh2o_b = 0.9 * pfit_wh2o;
2930 :
2931 0 : for(unsigned int i = 0; i < IdChannels.size(); i++) {
2932 :
2933 0 : tebb_fit[i] = RT(pfit_wh2o_b,
2934 0 : skyCoupling[i],
2935 : tspill,
2936 : airm,
2937 0 : IdChannels[i],
2938 0 : signalGain[i]);
2939 :
2940 : }
2941 :
2942 0 : chisqr = 0;
2943 0 : for(unsigned int i = 0; i < IdChannels.size(); i++) {
2944 0 : res = -tebb_fit[i] + measuredSkyBrightnessVector[i].get(Temperature::UnitKelvin);
2945 0 : chisqr = chisqr + res * res;
2946 : }
2947 0 : if(IdChannels.size() > 1) {
2948 0 : chisqr = chisqr / (IdChannels.size() - 1);
2949 : }
2950 :
2951 0 : if(fabs(chisq1 - chisqr) > 0.001) {
2952 0 : if(chisq1 < chisqr) {
2953 0 : flamda = flamda * 10.0;
2954 0 : goto adjust;
2955 : }
2956 : }
2957 :
2958 0 : flamda = flamda / 10.0;
2959 0 : sig_fit = sqrt(chisqr);
2960 0 : pfit_wh2o = pfit_wh2o_b;
2961 0 : sigma_wh2o = Length(sqrt(array / alpha) * sig_fit * pfit_wh2o
2962 0 : * (getGroundWH2O().get(Length::UnitMilliMeter)), Length::UnitMilliMeter);
2963 :
2964 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
2965 :
2966 0 : sigma_fit_temp0 = Temperature(sig_fit, Temperature::UnitKelvin);
2967 :
2968 0 : wh2o_retrieved = Length(pfit_wh2o * getGroundWH2O().get(Length::UnitMilliMeter), Length::UnitMilliMeter);
2969 :
2970 0 : goto salir;
2971 :
2972 : }
2973 :
2974 : }
2975 :
2976 0 : wh2o_retrieved = werr; // Extra error code, fit not reached after 20 iterations
2977 0 : sigma_fit_temp0 = Temperature(sig_fit, Temperature::UnitKelvin); // Extra error code, fit not reached after 20 iterations
2978 0 : sigma_wh2o = werr; // Extra error code, fit not reached after 20 iterations
2979 :
2980 0 : salir:
2981 :
2982 0 : vector<Temperature> ttt;
2983 :
2984 0 : for(unsigned int i = 0; i < IdChannels.size(); i++) {
2985 0 : ttt.push_back(Temperature(tebb_fit[i], Temperature::UnitKelvin));
2986 : }
2987 :
2988 0 : if(wh2o_retrieved.get() > 0.0) {
2989 0 : wh2o_user_ = wh2o_retrieved;
2990 : }
2991 : return WVRMeasurement(elevation,
2992 : measuredSkyBrightnessVector,
2993 : ttt,
2994 : wh2o_retrieved,
2995 0 : sigma_fit_temp0);
2996 :
2997 0 : }
2998 :
2999 0 : void SkyStatus::WaterVaporRetrieval_fromWVR(vector<WVRMeasurement> &RadiometerData,
3000 : unsigned int n,
3001 : unsigned int m)
3002 : {
3003 :
3004 0 : for(unsigned int i = n; i < m; i++) {
3005 :
3006 0 : WaterVaporRetrieval_fromWVR(RadiometerData[i]);
3007 :
3008 : }
3009 :
3010 0 : }
3011 :
3012 0 : void SkyStatus::WaterVaporRetrieval_fromWVR(WVRMeasurement &RadiometerData)
3013 : {
3014 :
3015 0 : WVRMeasurement RadiometerData_withRetrieval;
3016 :
3017 : // cout << waterVaporRadiometer_.getIdChannels().size() << endl;
3018 : // cout << RadiometerData.getmeasuredSkyBrightness()[0].get(Temperature::UnitKelvin) << " K" << endl;
3019 : // cout << waterVaporRadiometer_.getIdChannels()[1] << endl;
3020 : // cout << getAssocSpwIds(waterVaporRadiometer_.getIdChannels())[1] << endl;
3021 :
3022 :
3023 : // cout << "zz=" << waterVaporRadiometer_.getIdChannels().size() << endl;
3024 :
3025 : RadiometerData_withRetrieval
3026 0 : = mkWaterVaporRetrieval_fromWVR(RadiometerData.getmeasuredSkyBrightness(),
3027 0 : waterVaporRadiometer_.getIdChannels(),
3028 0 : waterVaporRadiometer_.getSkyCoupling(),
3029 0 : waterVaporRadiometer_.getsignalGain(),
3030 0 : waterVaporRadiometer_.getSpilloverTemperature(),
3031 0 : RadiometerData.getElevation());
3032 :
3033 : // cout << "_fromWVR Sky Coupling = " << waterVaporRadiometer_.getSkyCoupling()[0] << endl;
3034 : // cout << "Signal Gain = " << waterVaporRadiometer_.getsignalGain()[0].get(Percent::UnitPercent) << " %" << endl;
3035 : // cout << "Spillover Temp. = " << waterVaporRadiometer_.getSpilloverTemperature().get(Temperature::UnitKelvin) << " K" << endl;
3036 : // cout << "Elevation = " << RadiometerData.getElevation().get(Angle::UnitDegree) << endl;
3037 : // cout << "PWV=" << RadiometerData_withRetrieval.getretrievedWaterVaporColumn().get(Length::UnitMilliMeter) << " mm" << endl;
3038 :
3039 0 : RadiometerData.setretrievedWaterVaporColumn(RadiometerData_withRetrieval.getretrievedWaterVaporColumn());
3040 0 : RadiometerData.setfittedSkyBrightness(RadiometerData_withRetrieval.getfittedSkyBrightness());
3041 0 : RadiometerData.setSigmaFit(RadiometerData_withRetrieval.getSigmaFit());
3042 :
3043 0 : }
3044 :
3045 0 : void SkyStatus::updateSkyCoupling_fromWVR(vector<WVRMeasurement> &RadiometerData,
3046 : unsigned int n,
3047 : unsigned int m)
3048 : {
3049 : double pfit;
3050 0 : double deltaa = 0.02;
3051 : // double sig_fit = -999.0; // [-Wunused_but_set_variable]
3052 0 : double eps = 0.01;
3053 :
3054 : unsigned int num;
3055 : double flamda;
3056 0 : unsigned int niter = 20;
3057 : double alpha;
3058 : double beta;
3059 : double array;
3060 : double f1;
3061 : double psave;
3062 : double f2;
3063 : double deriv;
3064 : double chisq1;
3065 : double chisqr;
3066 : double pfit_b;
3067 : double res;
3068 :
3069 0 : num = 0;
3070 :
3071 0 : flamda = 0.001;
3072 0 : pfit = 0.5;
3073 :
3074 : // Find the maximum of skycoupling on the WVR channels
3075 : // This value is used to assure that the found skycoupling does not exceed 1
3076 0 : double maxCoupling =0.;
3077 0 : for (unsigned int i=0; i<waterVaporRadiometer_.getSkyCoupling().size(); i++)
3078 0 : if (waterVaporRadiometer_.getSkyCoupling()[i]>maxCoupling)
3079 0 : maxCoupling = waterVaporRadiometer_.getSkyCoupling()[i];
3080 :
3081 :
3082 0 : for(unsigned int kite = 0; kite < niter; kite++) {
3083 :
3084 0 : num = num + 1;
3085 0 : beta = 0.0;
3086 0 : alpha = 0.0;
3087 :
3088 0 : if (pfit*maxCoupling>1.5)
3089 0 : pfit = 1.-deltaa;
3090 0 : f1 = sigmaSkyCouplingRetrieval_fromWVR(pfit,
3091 0 : waterVaporRadiometer_,
3092 : RadiometerData,
3093 : n,
3094 : m);
3095 0 : psave = pfit;
3096 0 : pfit = pfit + deltaa;
3097 0 : f2 = sigmaSkyCouplingRetrieval_fromWVR(pfit,
3098 0 : waterVaporRadiometer_,
3099 : RadiometerData,
3100 : n,
3101 : m);
3102 0 : deriv = (f2 - f1) / deltaa;
3103 0 : pfit = psave;
3104 0 : beta = beta - f1 * deriv;
3105 0 : alpha = alpha + deriv * deriv;
3106 0 : chisq1 = f1 * f1;
3107 :
3108 0 : adjust: array = 1.0 / (1.0 + flamda);
3109 0 : pfit_b = pfit;
3110 0 : pfit_b = pfit_b + beta * array / alpha;
3111 :
3112 0 : chisqr = 0.;
3113 :
3114 0 : if(pfit_b < 0.0) {
3115 0 : pfit_b = 0.9 * pfit;
3116 : }
3117 0 : if (pfit_b*maxCoupling>1.5)
3118 0 : pfit_b = 1;
3119 :
3120 0 : res = sigmaSkyCouplingRetrieval_fromWVR(pfit_b,
3121 0 : waterVaporRadiometer_,
3122 : RadiometerData,
3123 : n,
3124 : m);
3125 0 : chisqr = chisqr + res * res;
3126 :
3127 0 : if(fabs(chisq1 - chisqr) > 0.001) {
3128 0 : if(chisq1 < chisqr) {
3129 0 : flamda = flamda * 10.0;
3130 0 : goto adjust;
3131 : }
3132 : }
3133 :
3134 0 : flamda = flamda / 10.0;
3135 : // sig_fit = sqrt(chisqr); // [-Wunused_but_set_variable]
3136 0 : pfit = pfit_b;
3137 :
3138 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
3139 0 : goto salir;
3140 : }
3141 :
3142 : }
3143 :
3144 0 : salir: waterVaporRadiometer_.multiplySkyCoupling(pfit);
3145 :
3146 : // cout << "pfit=" << pfit << " sky_coupling: " << waterVaporRadiometer_.getSkyCoupling()[0] << endl;
3147 :
3148 :
3149 0 : }
3150 0 : void SkyStatus::updateSkyCouplingChannel_fromWVR(vector<WVRMeasurement> &RadiometerData,
3151 : unsigned int ichan,
3152 : unsigned int n,
3153 : unsigned int m)
3154 : {
3155 : double pfit;
3156 0 : double deltaa = 0.02;
3157 : // double sig_fit = -999.0; // [-Wunused_but_set_variable]
3158 0 : double eps = 0.01;
3159 :
3160 : unsigned int num;
3161 : double flamda;
3162 0 : unsigned int niter = 20;
3163 : double alpha;
3164 : double beta;
3165 : double array;
3166 : double f1;
3167 : double psave;
3168 : double f2;
3169 : double deriv;
3170 : double chisq1;
3171 : double chisqr;
3172 : double pfit_b;
3173 : double res;
3174 :
3175 0 : num = 0;
3176 :
3177 0 : flamda = 0.001;
3178 0 : pfit = 1.00;
3179 :
3180 : // This value is used to assure that the found skycoupling does not exceed 1
3181 0 : double maxCoupling = waterVaporRadiometer_.getSkyCoupling()[ichan];
3182 :
3183 :
3184 0 : for(unsigned int kite = 0; kite < niter; kite++) {
3185 :
3186 0 : num = num + 1;
3187 0 : beta = 0.0;
3188 0 : alpha = 0.0;
3189 :
3190 0 : if (pfit*maxCoupling>1)
3191 0 : pfit = 1.-deltaa;
3192 0 : f1 = sigmaSkyCouplingChannelRetrieval_fromWVR(pfit,
3193 0 : waterVaporRadiometer_,
3194 : RadiometerData,
3195 : ichan,
3196 : n,
3197 : m);
3198 0 : psave = pfit;
3199 0 : pfit = pfit + deltaa;
3200 0 : f2 = sigmaSkyCouplingChannelRetrieval_fromWVR(pfit,
3201 0 : waterVaporRadiometer_,
3202 : RadiometerData,
3203 : ichan,
3204 : n,
3205 : m);
3206 0 : deriv = (f2 - f1) / deltaa;
3207 0 : pfit = psave;
3208 0 : beta = beta - f1 * deriv;
3209 0 : alpha = alpha + deriv * deriv;
3210 0 : chisq1 = f1 * f1;
3211 :
3212 0 : adjust: array = 1.0 / (1.0 + flamda);
3213 0 : pfit_b = pfit;
3214 0 : pfit_b = pfit_b + beta * array / alpha;
3215 :
3216 0 : chisqr = 0.;
3217 :
3218 0 : if(pfit_b < 0.0) {
3219 0 : pfit_b = 0.9 * pfit;
3220 : }
3221 0 : if (pfit_b*maxCoupling>1)
3222 0 : pfit_b = 1/maxCoupling;
3223 :
3224 0 : res = sigmaSkyCouplingChannelRetrieval_fromWVR(pfit_b,
3225 0 : waterVaporRadiometer_,
3226 : RadiometerData,
3227 : ichan,
3228 : n,
3229 : m);
3230 0 : chisqr = chisqr + res * res;
3231 :
3232 0 : if(fabs(chisq1 - chisqr) > 0.001) {
3233 0 : if(chisq1 < chisqr) {
3234 0 : flamda = flamda * 10.0;
3235 0 : goto adjust;
3236 : }
3237 : }
3238 :
3239 0 : flamda = flamda / 10.0;
3240 : // sig_fit = sqrt(chisqr); // [-Wunused_but_set_variable]
3241 0 : pfit = pfit_b;
3242 :
3243 0 : if(fabs(sqrt(chisq1) - sqrt(chisqr)) < eps) {
3244 0 : goto salir;
3245 : }
3246 :
3247 : }
3248 :
3249 0 : salir: waterVaporRadiometer_.multiplySkyCouplingChannel(ichan, pfit);
3250 :
3251 : // cout << "pfit=" << pfit << " sky_coupling: " << waterVaporRadiometer_.getSkyCoupling()[0] << endl;
3252 :
3253 :
3254 0 : }
3255 :
3256 0 : double SkyStatus::sigmaSkyCouplingRetrieval_fromWVR(double par_fit,
3257 : const WaterVaporRadiometer &wvr,
3258 : vector<WVRMeasurement> &RadiometerData,
3259 : unsigned int n,
3260 : unsigned int m)
3261 : {
3262 :
3263 0 : vector<double> skyCoupling = wvr.getSkyCoupling();
3264 :
3265 0 : for(unsigned int i = 0; i < skyCoupling.size(); i++) {
3266 0 : skyCoupling[i] = skyCoupling[i] * par_fit;
3267 : }
3268 :
3269 0 : WVRMeasurement RadiometerData_withRetrieval;
3270 :
3271 0 : for(unsigned int i = n; i < m; i++) {
3272 :
3273 : RadiometerData_withRetrieval
3274 0 : = mkWaterVaporRetrieval_fromWVR(RadiometerData[i].getmeasuredSkyBrightness(),
3275 0 : wvr.getIdChannels(),
3276 : skyCoupling,
3277 0 : wvr.getsignalGain(),
3278 0 : wvr.getSpilloverTemperature(),
3279 0 : RadiometerData[i].getElevation());
3280 :
3281 0 : RadiometerData[i].setretrievedWaterVaporColumn(RadiometerData_withRetrieval.getretrievedWaterVaporColumn());
3282 0 : RadiometerData[i].setfittedSkyBrightness(RadiometerData_withRetrieval.getfittedSkyBrightness());
3283 0 : RadiometerData[i].setSigmaFit(RadiometerData_withRetrieval.getSigmaFit());
3284 :
3285 : }
3286 :
3287 0 : return getWVRAverageSigmaTskyFit(RadiometerData, n, m).get(Temperature::UnitKelvin);
3288 :
3289 0 : }
3290 :
3291 0 : double SkyStatus::sigmaSkyCouplingChannelRetrieval_fromWVR(double par_fit,
3292 : const WaterVaporRadiometer &wvr,
3293 : vector<WVRMeasurement> &RadiometerData,
3294 : unsigned int ichan,
3295 : unsigned int n,
3296 : unsigned int m)
3297 : {
3298 :
3299 0 : vector<double> skyCoupling = wvr.getSkyCoupling();
3300 :
3301 : //for(unsigned int i = 0; i < skyCoupling.size(); i++) {
3302 0 : skyCoupling[ichan] = skyCoupling[ichan] * par_fit;
3303 : //}
3304 :
3305 0 : WVRMeasurement RadiometerData_withRetrieval;
3306 :
3307 0 : for(unsigned int i = n; i < m; i++) {
3308 :
3309 : RadiometerData_withRetrieval
3310 0 : = mkWaterVaporRetrieval_fromWVR(RadiometerData[i].getmeasuredSkyBrightness(),
3311 0 : wvr.getIdChannels(),
3312 : skyCoupling,
3313 0 : wvr.getsignalGain(),
3314 0 : wvr.getSpilloverTemperature(),
3315 0 : RadiometerData[i].getElevation());
3316 :
3317 0 : RadiometerData[i].setretrievedWaterVaporColumn(RadiometerData_withRetrieval.getretrievedWaterVaporColumn());
3318 0 : RadiometerData[i].setfittedSkyBrightness(RadiometerData_withRetrieval.getfittedSkyBrightness());
3319 0 : RadiometerData[i].setSigmaFit(RadiometerData_withRetrieval.getSigmaFit());
3320 :
3321 : }
3322 :
3323 0 : return getWVRAverageSigmaTskyFit(RadiometerData, n, m).get(Temperature::UnitKelvin);
3324 : // return getWVRSigmaChannelTskyFit(RadiometerData, ichan, n, m).get(Temperature::UnitKelvin);
3325 :
3326 0 : }
3327 :
3328 18 : void SkyStatus::rmSkyStatus()
3329 : {
3330 :
3331 : /* CODE À ÉCRIRE */
3332 :
3333 18 : }
3334 : ATM_NAMESPACE_END
3335 :
|
