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