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: ATMProfile.cpp Exp $"
20 : *
21 : * who when what
22 : * -------- -------- ----------------------------------------------
23 : * pardo 24/03/09 created
24 : */
25 :
26 : #include "ATMProfile.h"
27 : #include "ATMException.h"
28 : #include <iostream>
29 : #include <math.h>
30 : #include <sstream>
31 :
32 :
33 :
34 : ATM_NAMESPACE_BEGIN
35 :
36 0 : AtmProfile::AtmProfile(unsigned int n)
37 : {
38 0 : numLayer_ = n;
39 0 : initBasicAtmosphericParameterThresholds();
40 0 : for(unsigned int i = 0; i < numLayer_; ++i) {
41 0 : v_layerO3_.push_back(0.0);
42 0 : v_layerCO_.push_back(0.0);
43 0 : v_layerN2O_.push_back(0.0);
44 0 : v_layerNO2_.push_back(0.0);
45 0 : v_layerSO2_.push_back(0.0);
46 0 : v_layerThickness_.push_back(0.0);
47 0 : v_layerTemperature_.push_back(0.0);
48 0 : v_layerTemperature0_.push_back(0.0);
49 0 : v_layerTemperature1_.push_back(0.0);
50 0 : v_layerPressure_.push_back(0.0);
51 0 : v_layerPressure0_.push_back(0.0);
52 0 : v_layerPressure1_.push_back(0.0);
53 0 : v_layerWaterVapor_.push_back(0.0);
54 0 : v_layerWaterVapor0_.push_back(0.0);
55 0 : v_layerWaterVapor1_.push_back(0.0);
56 : }
57 0 : }
58 :
59 :
60 :
61 825 : AtmProfile::AtmProfile(const Length &altitude,
62 : const Pressure &groundPressure,
63 : const Temperature &groundTemperature,
64 : double tropoLapseRate,
65 : const Humidity &relativeHumidity,
66 : const Length &wvScaleHeight,
67 : const Pressure &pressureStep,
68 : double pressureStepFactor,
69 : const Length &topAtmProfile,
70 825 : unsigned int atmType) :
71 : // Atmospheretype atmType):
72 825 : typeAtm_(atmType), groundTemperature_(groundTemperature),
73 825 : tropoLapseRate_(tropoLapseRate), groundPressure_(groundPressure),
74 825 : relativeHumidity_(relativeHumidity), wvScaleHeight_(wvScaleHeight),
75 825 : pressureStep_(pressureStep), pressureStepFactor_(pressureStepFactor),
76 1650 : altitude_(altitude), topAtmProfile_(topAtmProfile)
77 : {
78 825 : numLayer_ = 0;
79 825 : numLayer_ = mkAtmProfile();
80 : // std::cout << "ad numLayer_=" << numLayer_ << std::endl;
81 825 : initBasicAtmosphericParameterThresholds();
82 825 : newBasicParam_ = true;
83 825 : }
84 :
85 27 : AtmProfile::AtmProfile(const Length &altitude,
86 : const Pressure &groundPressure,
87 : const Temperature &groundTemperature,
88 : double tropoLapseRate,
89 : const Humidity &relativeHumidity,
90 : const Length &wvScaleHeight,
91 : const Pressure &pressureStep,
92 : double pressureStepFactor,
93 : const Length &topAtmProfile,
94 : unsigned int atmType,
95 : const vector<Length> &v_layerBoundaries,
96 27 : const vector<Temperature> &v_layerTemperature):
97 27 : typeAtm_(atmType), groundTemperature_(groundTemperature),
98 27 : tropoLapseRate_(tropoLapseRate), groundPressure_(groundPressure),
99 27 : relativeHumidity_(relativeHumidity), wvScaleHeight_(wvScaleHeight),
100 27 : pressureStep_(pressureStep), pressureStepFactor_(pressureStepFactor),
101 54 : altitude_(altitude), topAtmProfile_(topAtmProfile)
102 : {
103 27 : numLayer_ = 0;
104 27 : numLayer_ = mkAtmProfile();
105 27 : initBasicAtmosphericParameterThresholds();
106 27 : newBasicParam_ = true;
107 27 : unsigned int nL1 = v_layerBoundaries.size();
108 27 : unsigned int nL2 = v_layerTemperature.size();
109 27 : if(nL1 == nL2 ) {
110 27 : double h=altitude_.get(Length::UnitMeter);
111 : double h0;
112 : double h1;
113 : double counter;
114 : double avT;
115 758 : for(unsigned int n = 0; n < numLayer_; n++) {
116 731 : h0=h;
117 731 : h1 = h0 + v_layerThickness_[n];
118 731 : counter = 0.0;
119 731 : avT = 0.0;
120 843 : for(unsigned int m = 0; m < nL1; m++) {
121 112 : if( h0 <= v_layerBoundaries[m].get(Length::UnitMeter) && h1 >= v_layerBoundaries[m].get(Length::UnitMeter) ){
122 : // std::cout << "n=" << n << " h0=" << h0 << " h1=" << h1 << " v_layerBoundaries[" << m << "]=" << v_layerBoundaries[m].get(Length::UnitMeter) << std::endl;
123 0 : avT = avT + v_layerTemperature[m].get(Temperature::UnitKelvin);
124 0 : counter = counter + 1.0;
125 : }
126 : }
127 731 : if(avT > 0.0){
128 : // std::cout << "layer" << n << "old average temperature:" << v_layerTemperature_[n] << std::endl;
129 0 : v_layerTemperature_[n] = avT/counter;
130 : // std::cout << "layer" << n << "new average temperature:" << v_layerTemperature_[n] << std::endl;
131 : }
132 731 : h=h1;
133 : }
134 : }
135 27 : }
136 :
137 0 : AtmProfile::AtmProfile(const Length &altitude,
138 : const Pressure &groundPressure,
139 : const Temperature &groundTemperature,
140 : double tropoLapseRate,
141 : const Humidity &relativeHumidity,
142 : const Length &wvScaleHeight,
143 0 : unsigned int atmType) :
144 : // Atmospheretype atmType):
145 0 : typeAtm_(atmType), groundTemperature_(groundTemperature),
146 0 : tropoLapseRate_(tropoLapseRate), groundPressure_(groundPressure),
147 0 : relativeHumidity_(relativeHumidity), wvScaleHeight_(wvScaleHeight),
148 0 : pressureStep_(10.0, Pressure::UnitMilliBar), pressureStepFactor_(1.2), altitude_(altitude),
149 0 : topAtmProfile_(48.0, Length::UnitKiloMeter)
150 : {
151 0 : numLayer_ = 0;
152 0 : numLayer_ = mkAtmProfile();
153 0 : initBasicAtmosphericParameterThresholds();
154 0 : newBasicParam_ = true;
155 0 : }
156 :
157 0 : AtmProfile::AtmProfile(const Length &altitude,
158 : const vector<Length> &v_layerThickness,
159 : const vector<Pressure> &v_layerPressure,
160 : const vector<Temperature> &v_layerTemperature,
161 0 : const vector<MassDensity> &v_layerWaterVapor)
162 : {
163 0 : newBasicParam_ = true;
164 0 : unsigned int nL1 = v_layerThickness.size();
165 0 : unsigned int nL2 = v_layerPressure.size();
166 0 : unsigned int nL3 = v_layerTemperature.size();
167 0 : unsigned int nL4 = v_layerWaterVapor.size();
168 :
169 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4) {
170 0 : numLayer_ = nL1;
171 0 : for(unsigned int n = 0; n < numLayer_; n++) {
172 0 : v_layerO3_.push_back(0.0);
173 0 : v_layerCO_.push_back(0.0);
174 0 : v_layerN2O_.push_back(0.0);
175 0 : v_layerNO2_.push_back(0.0);
176 0 : v_layerSO2_.push_back(0.0);
177 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
178 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
179 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
180 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
181 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
182 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
183 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
184 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
185 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
186 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
187 : }
188 0 : } else {
189 0 : numLayer_ = 0;
190 : }
191 0 : initBasicAtmosphericParameterThresholds();
192 0 : }
193 :
194 0 : AtmProfile::AtmProfile(const vector<Length> &v_layerBoundaries,
195 : const vector<Pressure> &v_layerPressure,
196 : const vector<Temperature> &v_layerTemperature,
197 0 : const vector<MassDensity> &v_layerWaterVapor)
198 : {
199 0 : newBasicParam_ = true;
200 0 : unsigned int nL1 = v_layerBoundaries.size();
201 0 : unsigned int nL2 = v_layerPressure.size();
202 0 : unsigned int nL3 = v_layerTemperature.size();
203 0 : unsigned int nL4 = v_layerWaterVapor.size();
204 :
205 0 : if(nL1 == nL2 && nL2 == nL3 && nL3 == nL4) {
206 0 : numLayer_ = nL1 - 1;
207 0 : altitude_ = v_layerBoundaries[0];
208 0 : for(unsigned int n = 0; n < numLayer_; n++) {
209 0 : v_layerO3_.push_back(0.0);
210 0 : v_layerCO_.push_back(0.0);
211 0 : v_layerN2O_.push_back(0.0);
212 0 : v_layerNO2_.push_back(0.0);
213 0 : v_layerSO2_.push_back(0.0);
214 0 : v_layerThickness_.push_back(v_layerBoundaries[n+1].get(Length::UnitMeter) - v_layerBoundaries[n].get(Length::UnitMeter));
215 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
216 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
217 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
218 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
219 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
220 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
221 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
222 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
223 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
224 : }
225 0 : } else {
226 0 : numLayer_ = 0;
227 : }
228 0 : initBasicAtmosphericParameterThresholds();
229 0 : }
230 :
231 0 : AtmProfile::AtmProfile(const Length &altitude,
232 : const vector<Length> &v_layerThickness,
233 : const vector<Pressure> &v_layerPressure,
234 : const vector<Temperature> &v_layerTemperature,
235 0 : const vector<NumberDensity> &v_layerWaterVapor)
236 : {
237 0 : newBasicParam_ = true;
238 0 : unsigned int nL1 = v_layerThickness.size();
239 0 : unsigned int nL2 = v_layerPressure.size();
240 0 : unsigned int nL3 = v_layerTemperature.size();
241 0 : unsigned int nL4 = v_layerWaterVapor.size();
242 :
243 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4) {
244 0 : numLayer_ = nL1;
245 0 : for(unsigned int n = 0; n < numLayer_; n++) {
246 0 : v_layerO3_.push_back(0.0);
247 0 : v_layerCO_.push_back(0.0);
248 0 : v_layerN2O_.push_back(0.0);
249 0 : v_layerNO2_.push_back(0.0);
250 0 : v_layerSO2_.push_back(0.0);
251 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
252 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
253 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
254 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
255 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
256 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
257 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
258 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
259 0 : / (1000.0 * 6.023e23));
260 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
261 0 : / (1000.0 * 6.023e23));
262 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
263 0 : / (1000.0 * 6.023e23));
264 : }
265 0 : } else {
266 0 : numLayer_ = 0;
267 : }
268 0 : initBasicAtmosphericParameterThresholds();
269 0 : }
270 :
271 0 : AtmProfile::AtmProfile(const vector<Length> &v_layerBoundaries,
272 : const vector<Pressure> &v_layerPressure,
273 : const vector<Temperature> &v_layerTemperature,
274 0 : const vector<NumberDensity> &v_layerWaterVapor)
275 : {
276 0 : newBasicParam_ = true;
277 0 : unsigned int nL1 = v_layerBoundaries.size();
278 0 : unsigned int nL2 = v_layerPressure.size();
279 0 : unsigned int nL3 = v_layerTemperature.size();
280 0 : unsigned int nL4 = v_layerWaterVapor.size();
281 :
282 0 : if(nL1 == nL2 && nL2 == nL3 && nL3 == nL4) {
283 0 : numLayer_ = nL1 - 1;
284 0 : altitude_ = v_layerBoundaries[0];
285 0 : for(unsigned int n = 0; n < numLayer_; n++) {
286 0 : v_layerO3_.push_back(0.0);
287 0 : v_layerCO_.push_back(0.0);
288 0 : v_layerN2O_.push_back(0.0);
289 0 : v_layerNO2_.push_back(0.0);
290 0 : v_layerSO2_.push_back(0.0);
291 0 : v_layerThickness_.push_back(v_layerBoundaries[n+1].get(Length::UnitMeter) - v_layerBoundaries[n].get(Length::UnitMeter));
292 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
293 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
294 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
295 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
296 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
297 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
298 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
299 0 : / (1000.0 * 6.023e23));
300 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
301 0 : / (1000.0 * 6.023e23));
302 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
303 0 : / (1000.0 * 6.023e23));
304 : }
305 0 : } else {
306 0 : numLayer_ = 0;
307 : }
308 0 : initBasicAtmosphericParameterThresholds();
309 0 : }
310 :
311 0 : AtmProfile::AtmProfile(const Length &altitude,
312 : const vector<Length> &v_layerThickness,
313 : const vector<Pressure> &v_layerPressure,
314 : const vector<Temperature> &v_layerTemperature,
315 : const vector<MassDensity> &v_layerWaterVapor,
316 0 : const vector<NumberDensity> &v_layerO3)
317 : {
318 0 : newBasicParam_ = true;
319 0 : unsigned int nL1 = v_layerThickness.size();
320 0 : unsigned int nL2 = v_layerPressure.size();
321 0 : unsigned int nL3 = v_layerTemperature.size();
322 0 : unsigned int nL4 = v_layerWaterVapor.size();
323 0 : unsigned int nL5 = v_layerO3.size();
324 :
325 0 : if(nL1 * 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5) {
326 0 : numLayer_ = nL1;
327 0 : for(unsigned int n = 0; n < numLayer_; n++) {
328 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
329 0 : v_layerCO_.push_back(0.0);
330 0 : v_layerN2O_.push_back(0.0);
331 0 : v_layerNO2_.push_back(0.0);
332 0 : v_layerSO2_.push_back(0.0);
333 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
334 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
335 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
336 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
337 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
338 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
339 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
340 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
341 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
342 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
343 : }
344 0 : } else {
345 0 : numLayer_ = 0;
346 : }
347 0 : initBasicAtmosphericParameterThresholds();
348 0 : }
349 :
350 0 : AtmProfile::AtmProfile(const Length &altitude,
351 : const vector<Length> &v_layerThickness,
352 : const vector<Pressure> &v_layerPressure,
353 : const vector<Temperature> &v_layerTemperature,
354 : const vector<NumberDensity> &v_layerWaterVapor,
355 0 : const vector<NumberDensity> &v_layerO3)
356 : {
357 0 : newBasicParam_ = true;
358 0 : unsigned int nL1 = v_layerThickness.size();
359 0 : unsigned int nL2 = v_layerPressure.size();
360 0 : unsigned int nL3 = v_layerTemperature.size();
361 0 : unsigned int nL4 = v_layerWaterVapor.size();
362 0 : unsigned int nL5 = v_layerO3.size();
363 :
364 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5) {
365 0 : numLayer_ = nL1;
366 0 : for(unsigned int n = 0; n < numLayer_; n++) {
367 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
368 0 : v_layerCO_.push_back(0.0);
369 0 : v_layerN2O_.push_back(0.0);
370 0 : v_layerNO2_.push_back(0.0);
371 0 : v_layerSO2_.push_back(0.0);
372 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
373 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
374 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
375 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
376 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
377 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
378 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
379 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
380 0 : / (1000.0 * 6.023e23));
381 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
382 0 : / (1000.0 * 6.023e23));
383 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
384 0 : / (1000.0 * 6.023e23));
385 : }
386 0 : } else {
387 0 : numLayer_ = 0;
388 : }
389 0 : initBasicAtmosphericParameterThresholds();
390 0 : }
391 :
392 0 : AtmProfile::AtmProfile(const Length &altitude,
393 : const vector<Length> &v_layerThickness,
394 : const vector<Pressure> &v_layerPressure,
395 : const vector<Temperature> &v_layerTemperature,
396 : const vector<MassDensity> &v_layerWaterVapor,
397 : const vector<NumberDensity> &v_layerO3,
398 : const vector<NumberDensity> &v_layerCO,
399 : const vector<NumberDensity> &v_layerN2O,
400 : const vector<NumberDensity> &v_layerNO2,
401 0 : const vector<NumberDensity> &v_layerSO2)
402 : {
403 0 : newBasicParam_ = true;
404 0 : unsigned int nL1 = v_layerThickness.size();
405 0 : unsigned int nL2 = v_layerPressure.size();
406 0 : unsigned int nL3 = v_layerTemperature.size();
407 0 : unsigned int nL4 = v_layerWaterVapor.size();
408 0 : unsigned int nL5 = v_layerO3.size();
409 0 : unsigned int nL6 = v_layerCO.size();
410 0 : unsigned int nL7 = v_layerN2O.size();
411 0 : unsigned int nL8 = v_layerNO2.size();
412 0 : unsigned int nL9 = v_layerSO2.size();
413 :
414 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5 && nL5 == nL6 && nL6
415 0 : == nL7 && nL7 == nL8 && nL8 == nL9) {
416 0 : numLayer_ = nL1;
417 0 : for(unsigned int n = 0; n < numLayer_; n++) {
418 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
419 0 : v_layerCO_.push_back(v_layerCO[n].get(NumberDensity::UnitInverseCubicMeter));
420 0 : v_layerN2O_.push_back(v_layerN2O[n].get(NumberDensity::UnitInverseCubicMeter));
421 0 : v_layerNO2_.push_back(v_layerNO2[n].get(NumberDensity::UnitInverseCubicMeter));
422 0 : v_layerSO2_.push_back(v_layerSO2[n].get(NumberDensity::UnitInverseCubicMeter));
423 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
424 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
425 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
426 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
427 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
428 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
429 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
430 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
431 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
432 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
433 : }
434 0 : } else {
435 0 : numLayer_ = 0;
436 : }
437 0 : initBasicAtmosphericParameterThresholds();
438 0 : }
439 :
440 0 : AtmProfile::AtmProfile(const Length &altitude,
441 : const vector<Length> &v_layerThickness,
442 : const vector<Pressure> &v_layerPressure,
443 : const vector<Temperature> &v_layerTemperature,
444 : const vector<NumberDensity> &v_layerWaterVapor,
445 : const vector<NumberDensity> &v_layerO3,
446 : const vector<NumberDensity> &v_layerCO,
447 : const vector<NumberDensity> &v_layerN2O,
448 : const vector<NumberDensity> &v_layerNO2,
449 0 : const vector<NumberDensity> &v_layerSO2)
450 : {
451 0 : newBasicParam_ = true;
452 0 : unsigned int nL1 = v_layerThickness.size();
453 0 : unsigned int nL2 = v_layerPressure.size();
454 0 : unsigned int nL3 = v_layerTemperature.size();
455 0 : unsigned int nL4 = v_layerWaterVapor.size();
456 0 : unsigned int nL5 = v_layerO3.size();
457 0 : unsigned int nL6 = v_layerCO.size();
458 0 : unsigned int nL7 = v_layerN2O.size();
459 0 : unsigned int nL8 = v_layerNO2.size();
460 0 : unsigned int nL9 = v_layerSO2.size();
461 :
462 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5 && nL5 == nL6 && nL6
463 0 : == nL7 && nL7 == nL8 && nL8 == nL9) {
464 0 : numLayer_ = nL1;
465 0 : for(unsigned int n = 0; n < numLayer_; n++) {
466 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
467 0 : v_layerCO_.push_back(v_layerCO[n].get(NumberDensity::UnitInverseCubicMeter));
468 0 : v_layerN2O_.push_back(v_layerN2O[n].get(NumberDensity::UnitInverseCubicMeter));
469 0 : v_layerNO2_.push_back(v_layerNO2[n].get(NumberDensity::UnitInverseCubicMeter));
470 0 : v_layerSO2_.push_back(v_layerSO2[n].get(NumberDensity::UnitInverseCubicMeter));
471 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
472 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
473 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
474 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
475 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
476 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
477 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
478 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
479 0 : / (1000.0 * 6.023e23));
480 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
481 0 : / (1000.0 * 6.023e23));
482 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
483 0 : / (1000.0 * 6.023e23));
484 : }
485 0 : } else {
486 0 : numLayer_ = 0;
487 : }
488 0 : initBasicAtmosphericParameterThresholds();
489 0 : }
490 :
491 2600 : AtmProfile::AtmProfile(const AtmProfile &a)
492 : { //:AtmType(a.type_){
493 : // std::cout<<"AtmProfile copy constructor"<<endl; COMMENTED OUT BY JUAN MAY/16/2005
494 2600 : typeAtm_ = a.typeAtm_;
495 2600 : groundTemperature_ = a.groundTemperature_;
496 2600 : tropoLapseRate_ = a.tropoLapseRate_;
497 2600 : groundPressure_ = a.groundPressure_;
498 2600 : relativeHumidity_ = a.relativeHumidity_;
499 2600 : wvScaleHeight_ = a.wvScaleHeight_;
500 2600 : pressureStep_ = a.pressureStep_;
501 2600 : pressureStepFactor_ = a.pressureStepFactor_;
502 2600 : altitude_ = a.altitude_;
503 2600 : topAtmProfile_ = a.topAtmProfile_;
504 2600 : numLayer_ = a.numLayer_;
505 2600 : newBasicParam_ = a.newBasicParam_;
506 2600 : v_layerThickness_.reserve(numLayer_);
507 2600 : v_layerPressure_.reserve(numLayer_);
508 2600 : v_layerPressure0_.reserve(numLayer_);
509 2600 : v_layerPressure1_.reserve(numLayer_);
510 2600 : v_layerTemperature_.reserve(numLayer_);
511 2600 : v_layerTemperature0_.reserve(numLayer_);
512 2600 : v_layerTemperature1_.reserve(numLayer_);
513 2600 : v_layerWaterVapor_.reserve(numLayer_);
514 2600 : v_layerWaterVapor0_.reserve(numLayer_);
515 2600 : v_layerWaterVapor1_.reserve(numLayer_);
516 2600 : v_layerCO_.reserve(numLayer_);
517 2600 : v_layerO3_.reserve(numLayer_);
518 2600 : v_layerN2O_.reserve(numLayer_);
519 2600 : v_layerNO2_.reserve(numLayer_);
520 2600 : v_layerSO2_.reserve(numLayer_);
521 : // std::cout << "numLayer_=" << numLayer_ << std::endl; COMMENTED OUT BY JUAN MAY/16/2005
522 79594 : for(unsigned int n = 0; n < numLayer_; n++) {
523 76994 : v_layerThickness_.push_back(a.v_layerThickness_[n]);
524 76994 : v_layerTemperature_.push_back(a.v_layerTemperature_[n]);
525 76994 : v_layerTemperature0_.push_back(a.v_layerTemperature0_[n]);
526 76994 : v_layerTemperature1_.push_back(a.v_layerTemperature1_[n]);
527 : //cout << "n=" << n << std::endl;
528 76994 : v_layerWaterVapor_.push_back(a.v_layerWaterVapor_[n]);
529 76994 : v_layerWaterVapor0_.push_back(a.v_layerWaterVapor0_[n]);
530 76994 : v_layerWaterVapor1_.push_back(a.v_layerWaterVapor1_[n]);
531 76994 : v_layerPressure_.push_back(a.v_layerPressure_[n]);
532 76994 : v_layerPressure0_.push_back(a.v_layerPressure1_[n]);
533 76994 : v_layerPressure1_.push_back(a.v_layerPressure1_[n]);
534 76994 : v_layerCO_.push_back(a.v_layerCO_[n]);
535 76994 : v_layerO3_.push_back(a.v_layerO3_[n]);
536 76994 : v_layerN2O_.push_back(a.v_layerN2O_[n]);
537 76994 : v_layerNO2_.push_back(a.v_layerNO2_[n]);
538 76994 : v_layerSO2_.push_back(a.v_layerSO2_[n]);
539 : }
540 2600 : altitudeThreshold_ = a.altitudeThreshold_;
541 2600 : groundPressureThreshold_ = a.groundPressureThreshold_;
542 2600 : groundTemperatureThreshold_ = a.groundTemperatureThreshold_;
543 2600 : tropoLapseRateThreshold_ = a.tropoLapseRateThreshold_;
544 2600 : relativeHumidityThreshold_ = a.relativeHumidityThreshold_;
545 2600 : wvScaleHeightThreshold_ = a.wvScaleHeightThreshold_;
546 2600 : }
547 27 : void AtmProfile::setBasicAtmosphericParameterThresholds(const Length &altitudeThreshold,
548 : const Pressure &groundPressureThreshold,
549 : const Temperature &groundTemperatureThreshold,
550 : double tropoLapseRateThreshold,
551 : const Humidity &relativeHumidityThreshold,
552 : const Length &wvScaleHeightThreshold)
553 : {
554 27 : altitudeThreshold_ = altitudeThreshold;
555 27 : groundPressureThreshold_ = groundPressureThreshold;
556 27 : groundTemperatureThreshold_ = groundTemperatureThreshold;
557 27 : tropoLapseRateThreshold_ = tropoLapseRateThreshold;
558 27 : relativeHumidityThreshold_ = relativeHumidityThreshold;
559 27 : wvScaleHeightThreshold_ =wvScaleHeightThreshold;
560 27 : }
561 :
562 852 : void AtmProfile::initBasicAtmosphericParameterThresholds()
563 : {
564 852 : altitudeThreshold_ = Length(1.0,Length::UnitMeter);
565 852 : groundPressureThreshold_ = Pressure(99.,Pressure::UnitPascal); // DB 2014-03-06 : Choose 99 Pascal instead of 100 Pascal because the threshold must be lower than the value of delta_pressure used in getAverageDispersiveDryPathLength_GroundPressureDerivative()
566 852 : groundTemperatureThreshold_ = Temperature(0.3,Temperature::UnitKelvin);
567 852 : tropoLapseRateThreshold_ = 0.01;
568 852 : relativeHumidityThreshold_ = Humidity(100.,Percent::UnitPercent);
569 852 : wvScaleHeightThreshold_ = Length(20.,Length::UnitMeter);
570 852 : }
571 :
572 163 : bool AtmProfile::updateAtmProfile(const Length &altitude,
573 : const Pressure &groundPressure,
574 : const Temperature &groundTemperature,
575 : double tropoLapseRate,
576 : const Humidity &relativeHumidity,
577 : const Length &wvScaleHeight)
578 : {
579 :
580 : /* TODO A faire: pour decider s'il faut recalculer le profile on devrait plutot donner des seuils, eg
581 : if(fabs(altitude_.get()-altitude.get())>0.1)mkNewProfile=true;
582 : */
583 :
584 : unsigned int numLayer;
585 163 : bool mkNewProfile = false;
586 :
587 : // if(altitude_.get() != altitude.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"altitude has changed" <<endl;
588 : // if(groundPressure_.get() != groundPressure.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"ground pressure has changed" <<endl;
589 : // if(groundTemperature_.get() != groundTemperature.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"ground temperature has changed"<<endl;
590 : // if(wvScaleHeight_.get() != wvScaleHeight.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"wv scale height has changed" <<endl;
591 : // if(tropoLapseRate_ != tropoLapseRate) mkNewProfile = true; //if(mkNewProfile)cout<<"tropo lapse rate has changed" <<endl;
592 : // if(relativeHumidity_.get() != relativeHumidity.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"relative humidity has changed" <<endl;
593 : //! all these thresholds shoudl be parametrized.
594 163 : if(fabs(altitude_.get()-altitude.get())> altitudeThreshold_.get()) { // Length(0.1,Length::UnitMeter).get()) {
595 0 : mkNewProfile = true;
596 : //cout<<"altitude has changed" <<endl;
597 : }
598 163 : if(fabs(groundPressure_.get()-groundPressure.get())> groundPressureThreshold_.get()) { // Pressure(100.,Pressure::UnitPascal).get()) {
599 163 : mkNewProfile = true;
600 : //cout<<"ground pressure has changed" <<endl;
601 : }
602 163 : if(fabs(groundTemperature_.get()-groundTemperature.get()) > groundTemperatureThreshold_.get()) { // Temperature(0.3,Temperature::UnitKelvin).get()) {
603 43 : mkNewProfile = true;
604 : //cout<<"ground temperature has changed"<<endl;
605 : }
606 163 : if(fabs(wvScaleHeight_.get()-wvScaleHeight.get()) > wvScaleHeightThreshold_.get()) { // Length(20.,Length::UnitMeter).get() ) {
607 0 : mkNewProfile = true;
608 : //cout<<"wv scale height has changed" <<endl;
609 : }
610 163 : if(fabs(tropoLapseRate_-tropoLapseRate) > tropoLapseRateThreshold_) { // 0.01) {
611 0 : mkNewProfile = true;
612 : //cout<<"tropo lapse rate has changed" <<endl;
613 : }
614 : // we use WVR...
615 163 : if(fabs(relativeHumidity_.get()-relativeHumidity.get())> relativeHumidityThreshold_.get()) { // Humidity(100.,Percent::UnitPercent).get()) {
616 43 : mkNewProfile = true;
617 : //cout<<"relative humidity has changed" <<endl;
618 : }
619 163 : if(mkNewProfile) {
620 163 : newBasicParam_ = true;
621 163 : altitude_ = altitude;
622 163 : groundPressure_ = groundPressure;
623 163 : groundTemperature_ = groundTemperature;
624 163 : tropoLapseRate_ = tropoLapseRate;
625 163 : relativeHumidity_ = relativeHumidity;
626 163 : wvScaleHeight_ = wvScaleHeight;
627 163 : numLayer = mkAtmProfile();
628 163 : numLayer_ = numLayer;
629 : // std::cout << "There are new basic parameters, with " << numLayer_ << " layers " << std::endl;
630 : } else {
631 0 : numLayer = getNumLayer();
632 0 : numLayer_ = numLayer;
633 : }
634 :
635 163 : return mkNewProfile;
636 : }
637 :
638 : // Note that this setBasicAtmosphericParameters will be overrided by the subclasses.
639 0 : bool AtmProfile::setBasicAtmosphericParameters(const Length &altitude,
640 : const Pressure &groundPressure,
641 : const Temperature &groundTemperature,
642 : double tropoLapseRate,
643 : const Humidity &relativeHumidity,
644 : const Length &wvScaleHeight)
645 : {
646 0 : return updateAtmProfile(altitude,
647 : groundPressure,
648 : groundTemperature,
649 : tropoLapseRate,
650 : relativeHumidity,
651 0 : wvScaleHeight);
652 : }
653 :
654 813 : string AtmProfile::getAtmosphereType() const
655 : {
656 813 : return getAtmosphereType(typeAtm_);
657 : }
658 :
659 813 : string AtmProfile::getAtmosphereType(unsigned int typeAtm)
660 : {
661 813 : string type;
662 : string typeNames[] = { "TROPICAL", "MIDLATSUMMER", "MIDLATWINTER",
663 6504 : "SUBARTSUMMER", "SUBARTWINTER","US_ST76"};
664 :
665 813 : if(typeAtm < typeATM_end) {
666 813 : type = typeNames[typeAtm-1];
667 : } else {
668 0 : type = "DEFAULT";
669 : }
670 1626 : return type;
671 5691 : }
672 :
673 1 : vector<Temperature> AtmProfile::getTemperatureProfile() const
674 : {
675 1 : vector<Temperature> t;
676 1 : t.reserve(v_layerTemperature_.size());
677 31 : for(unsigned int i = 0; i < v_layerTemperature_.size(); i++) {
678 30 : Temperature tt(v_layerTemperature_[i], Temperature::UnitKelvin);
679 30 : t.push_back(tt);
680 30 : }
681 1 : return t;
682 0 : }
683 :
684 427457548 : Temperature AtmProfile::getLayerTemperature(unsigned int i) const
685 : {
686 : /*if(i > v_layerTemperature_.size() - 1) {
687 : Temperature t(-999.0, Temperature::UnitKelvin);
688 : return t;
689 : } else {
690 : Temperature t(v_layerTemperature_[i], Temperature::UnitKelvin);
691 : return t;
692 : }*/
693 427457548 : if(i > v_layerTemperature_.size() - 1) {
694 0 : std::ostringstream oss;
695 0 : oss << "Not a valid layer: " << i;
696 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
697 0 : }
698 427457548 : return Temperature(v_layerTemperature_[i], Temperature::UnitKelvin);
699 : }
700 :
701 0 : Temperature AtmProfile::getLayerBottomTemperature(unsigned int i) const
702 : {
703 0 : if(i > v_layerTemperature0_.size() - 1) {
704 0 : std::ostringstream oss;
705 0 : oss << "Not a valid layer: " << i;
706 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
707 0 : }
708 0 : return Temperature(v_layerTemperature0_[i], Temperature::UnitKelvin);
709 : }
710 :
711 0 : Temperature AtmProfile::getLayerTopTemperature(unsigned int i) const
712 : {
713 0 : if(i > v_layerTemperature1_.size() - 1) {
714 0 : std::ostringstream oss;
715 0 : oss << "Not a valid layer: " << i;
716 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
717 0 : }
718 0 : return Temperature(v_layerTemperature1_[i], Temperature::UnitKelvin);
719 : }
720 :
721 :
722 0 : void AtmProfile::setAltitude(const Length &groundaltitude)
723 : {
724 :
725 0 : if (groundaltitude <= altitude_){
726 :
727 : // std::cout << "extrapolar para abajo" << std::endl;
728 :
729 0 : int nextralayers = int( 0.50001+(altitude_-groundaltitude).get(Length::UnitMeter)/v_layerThickness_[0]);
730 0 : if (nextralayers == 0) {nextralayers=1;}
731 : // std::cout << "aaa=" << ( 0.50001+(altitude_-groundaltitude).get(Length::UnitMeter)/v_layerThickness_[0]) << std::endl;
732 0 : double newThickness = (altitude_-groundaltitude).get(Length::UnitMeter) / nextralayers;
733 :
734 : // std::cout << "Number of extra layers: " << nextralayers << " of thickness = " << newThickness << " m" << std::endl;
735 : // std::cout << "nextralayers-1 " << nextralayers-1 << std::endl;
736 :
737 0 : for(int i=nextralayers; i>0; i--){
738 0 : v_layerThickness_.push_back(newThickness);
739 0 : v_layerTemperature_.push_back(v_layerTemperature_[0]);
740 0 : v_layerTemperature0_.push_back(v_layerTemperature0_[0]);
741 0 : v_layerTemperature1_.push_back(v_layerTemperature1_[0]);
742 0 : v_layerWaterVapor_.push_back(v_layerWaterVapor_[0]);
743 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor0_[0]);
744 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor1_[0]);
745 0 : v_layerPressure_.push_back(v_layerPressure_[0]);
746 0 : v_layerPressure0_.push_back(v_layerPressure0_[0]);
747 0 : v_layerPressure1_.push_back(v_layerPressure1_[0]);
748 0 : v_layerCO_.push_back(v_layerCO_[0]);
749 0 : v_layerO3_.push_back(v_layerO3_[0]);
750 0 : v_layerN2O_.push_back(v_layerN2O_[0]);
751 0 : v_layerNO2_.push_back(v_layerNO2_[0]);
752 0 : v_layerSO2_.push_back(v_layerSO2_[0]);
753 : }
754 :
755 0 : for(int i=v_layerThickness_.size()-1; i>nextralayers-1; i--){
756 0 : v_layerThickness_[i] = v_layerThickness_[i-nextralayers];
757 0 : v_layerTemperature_[i] = v_layerTemperature_[i-nextralayers];
758 0 : v_layerTemperature0_[i] = v_layerTemperature0_[i-nextralayers];
759 0 : v_layerTemperature1_[i] = v_layerTemperature1_[i-nextralayers];
760 0 : v_layerWaterVapor_[i] = v_layerWaterVapor_[i-nextralayers];
761 0 : v_layerWaterVapor0_[i] = v_layerWaterVapor0_[i-nextralayers];
762 0 : v_layerWaterVapor1_[i] = v_layerWaterVapor1_[i-nextralayers];
763 0 : v_layerPressure_[i] = v_layerPressure_[i-nextralayers];
764 0 : v_layerPressure0_[i] = v_layerPressure0_[i-nextralayers];
765 0 : v_layerPressure1_[i] = v_layerPressure1_[i-nextralayers];
766 0 : v_layerCO_[i] = v_layerCO_[i-nextralayers];
767 0 : v_layerO3_[i] = v_layerO3_[i-nextralayers];
768 0 : v_layerN2O_[i] = v_layerN2O_[i-nextralayers];
769 0 : v_layerNO2_[i] = v_layerNO2_[i-nextralayers];
770 0 : v_layerSO2_[i] = v_layerSO2_[i-nextralayers];
771 : }
772 :
773 : // std::cout << "nextralayers=" << nextralayers << std::endl;
774 :
775 0 : for(int i=nextralayers; i>0 ; i--){
776 :
777 0 : v_layerThickness_[i-1] = newThickness;
778 0 : v_layerTemperature1_[i-1] = v_layerTemperature0_[i];
779 0 : v_layerTemperature0_[i-1] = v_layerTemperature0_[i]-tropoLapseRate_*0.001*v_layerThickness_[i-1];
780 0 : v_layerTemperature_[i-1] = (v_layerTemperature0_[i]+v_layerTemperature1_[i])/2.0;
781 :
782 0 : v_layerPressure1_[i-1] = v_layerPressure0_[i];
783 0 : v_layerPressure0_[i-1] = v_layerPressure0_[i]* exp(-0.0341695 *
784 0 : (-newThickness)/ v_layerTemperature_[i-1] );
785 0 : v_layerPressure_[i-1] = exp((log(v_layerPressure0_[i-1])+ log(v_layerPressure1_[i-1]))/2.0);
786 :
787 0 : v_layerWaterVapor1_[i-1] = v_layerWaterVapor0_[i];
788 0 : v_layerWaterVapor0_[i-1] = v_layerWaterVapor0_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
789 0 : v_layerWaterVapor_[i-1] = exp((log(v_layerWaterVapor0_[i-1])+ log(v_layerWaterVapor1_[i-1]))/2.0);
790 :
791 0 : v_layerCO_[i-1] = v_layerCO_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
792 0 : v_layerO3_[i-1] = v_layerO3_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
793 0 : v_layerN2O_[i-1] = v_layerN2O_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
794 0 : v_layerNO2_[i-1] = v_layerNO2_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
795 0 : v_layerSO2_[i-1] = v_layerSO2_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
796 :
797 : }
798 :
799 0 : tropoLayer_ = tropoLayer_ + (v_layerThickness_.size()-numLayer_);
800 0 : numLayer_=v_layerThickness_.size();
801 :
802 : }else{
803 : // std::cout << "extrapolar para arriba" << std::endl;
804 0 : double extraheight = fabs(groundaltitude.get(Length::UnitMeter)-altitude_.get(Length::UnitMeter));
805 0 : double cumulheight = 0.0;
806 :
807 : // std::cout << "before v_layerTemperature_.size()=" << v_layerTemperature_.size() << std::endl;
808 :
809 0 : for(unsigned int i=0; i<numLayer_; i++){
810 0 : cumulheight = cumulheight + v_layerThickness_[i];
811 : // std::cout << i << " cumulheight = " << cumulheight << " extraheight= " << extraheight << std::endl;
812 0 : if (cumulheight>=extraheight) {
813 :
814 0 : double remainingheight = fabs(cumulheight-extraheight);
815 : // std::cout << "remainingheight = " << remainingheight << " m" << std::endl;
816 0 : v_layerThickness_.erase(v_layerThickness_.begin(),v_layerThickness_.begin()+i);
817 0 : v_layerThickness_[0] = remainingheight;
818 :
819 :
820 0 : v_layerTemperature0_.erase(v_layerTemperature0_.begin(),v_layerTemperature0_.begin()+i);
821 0 : v_layerTemperature1_.erase(v_layerTemperature1_.begin(),v_layerTemperature1_.begin()+i);
822 0 : v_layerTemperature_.erase(v_layerTemperature_.begin(),v_layerTemperature_.begin()+i);
823 0 : v_layerTemperature0_[0] = v_layerTemperature0_[1]-tropoLapseRate_*0.001*v_layerThickness_[0];
824 0 : v_layerTemperature_[0] = (v_layerTemperature0_[0]+v_layerTemperature1_[0])/2.0;
825 :
826 0 : v_layerPressure0_.erase(v_layerPressure0_.begin(),v_layerPressure0_.begin()+i);
827 : // std::cout << "v_layerPressure1_[0]=" << v_layerPressure1_[0] << " remainingheight=" << remainingheight << std::endl;
828 0 : v_layerPressure1_.erase(v_layerPressure1_.begin(),v_layerPressure1_.begin()+i);
829 0 : v_layerPressure_.erase(v_layerPressure_.begin(),v_layerPressure_.begin()+i);
830 0 : v_layerPressure0_[0] = v_layerPressure1_[0]* exp(-0.0341695 *
831 0 : (-remainingheight)/ v_layerTemperature_[0] );
832 0 : v_layerPressure_[0] = exp((log(v_layerPressure0_[0])+ log(v_layerPressure1_[0]))/2.0);
833 :
834 0 : v_layerWaterVapor0_.erase(v_layerWaterVapor0_.begin(),v_layerWaterVapor0_.begin()+i);
835 0 : v_layerWaterVapor1_.erase(v_layerWaterVapor1_.begin(),v_layerWaterVapor1_.begin()+i);
836 0 : v_layerWaterVapor_.erase(v_layerWaterVapor_.begin(),v_layerWaterVapor_.begin()+i);
837 0 : v_layerWaterVapor0_[0] = v_layerWaterVapor1_[0]*(v_layerPressure0_[0]/v_layerPressure1_[0]);
838 0 : v_layerWaterVapor_[0] = exp((log(v_layerWaterVapor0_[0])+ log(v_layerWaterVapor1_[0]))/2.0);
839 :
840 0 : v_layerCO_.erase(v_layerCO_.begin(),v_layerCO_.begin()+i);
841 0 : v_layerO3_.erase(v_layerO3_.begin(),v_layerO3_.begin()+i);
842 0 : v_layerN2O_.erase(v_layerN2O_.begin(),v_layerN2O_.begin()+i);
843 0 : v_layerNO2_.erase(v_layerNO2_.begin(),v_layerNO2_.begin()+i);
844 0 : v_layerSO2_.erase(v_layerSO2_.begin(),v_layerSO2_.begin()+i);
845 :
846 0 : tropoLayer_ = tropoLayer_ - i;
847 :
848 : // std::cout << "after v_layerTemperature_.size()=" << v_layerTemperature_.size() << std::endl;
849 0 : break;
850 : }
851 : }
852 :
853 0 : numLayer_ = v_layerThickness_.size();
854 :
855 : }
856 :
857 0 : altitude_ = groundaltitude;
858 0 : groundTemperature_ = v_layerTemperature0_[0];
859 :
860 : // std::cout << "v_layerPressure0_[0]=" << v_layerPressure0_[0] << std::endl;
861 :
862 0 : groundPressure_ = v_layerPressure0_[0]*100.0; // default units for Pressure are Pascals
863 :
864 : /*
865 : const Temperature Tdif(((true_antenna_altitude.get(Length::UnitKiloMeter)-altitude.get(Length::UnitKiloMeter))*tropoLapseRate),Temperature::UnitKelvin);
866 : groundTemperature_ = groundTemperature+Tdif;
867 : groundPressure_ = groundPressure * exp(-0.0341695 * pow((6.371/(6.371+altitude.get(Length::UnitKiloMeter))),2.0)
868 : * (true_antenna_altitude.get(Length::UnitMeter)-altitude.get(Length::UnitMeter))
869 : / ((groundTemperature.get(Temperature::UnitKelvin)+Tdif.get(Temperature::UnitKelvin))/2.0) );
870 : */
871 :
872 :
873 :
874 0 : }
875 :
876 0 : void AtmProfile::setLayerTemperature(unsigned int i, const Temperature &layerTemperature)
877 : {
878 0 : if(i < v_layerTemperature_.size()) {
879 0 : v_layerTemperature_[i] = layerTemperature.get(Temperature::UnitKelvin);
880 : }
881 0 : }
882 :
883 0 : vector<Length> AtmProfile::getThicknessProfile() const
884 : {
885 0 : vector<Length> l;
886 0 : l.reserve(v_layerThickness_.size());
887 0 : for(unsigned int i = 0; i < v_layerThickness_.size(); i++) {
888 0 : Length ll(v_layerThickness_[i], Length::UnitMeter);
889 0 : l.push_back(ll);
890 0 : }
891 0 : return l;
892 0 : }
893 :
894 427455158 : Length AtmProfile::getLayerThickness(unsigned int i) const
895 : {
896 : /*if(i > v_layerThickness_.size() - 1) {
897 : Length l(-999.0, Length::UnitMeter);
898 : return l;
899 : } else {
900 : Length l(v_layerThickness_[i], Length::UnitMeter);
901 : return l;
902 : }*/
903 427455158 : if(i > v_layerThickness_.size() - 1) {
904 0 : std::ostringstream oss;
905 0 : oss << "Not a valid layer: " << i;
906 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
907 0 : }
908 427455158 : return Length(v_layerThickness_[i], Length::UnitMeter);
909 : }
910 :
911 0 : Length AtmProfile::getLayerBottomHeightAboveGround(unsigned int i) const
912 : {
913 0 : if(i > v_layerThickness_.size() - 1) {
914 0 : std::ostringstream oss;
915 0 : oss << "Not a valid layer: " << i;
916 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
917 0 : }
918 0 : double h=0.0;
919 0 : for(unsigned int j = 0; j < i; j++) {
920 0 : h = h + v_layerThickness_[j];
921 : }
922 0 : return Length(h, Length::UnitMeter);
923 : }
924 :
925 0 : Length AtmProfile::getLayerTopHeightAboveGround(unsigned int i) const
926 : {
927 0 : if(i > v_layerThickness_.size() - 1) {
928 0 : std::ostringstream oss;
929 0 : oss << "Not a valid layer: " << i;
930 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
931 0 : }
932 0 : double h=0.0;
933 0 : for(unsigned int j = 0; j < i+1; j++) {
934 0 : h = h + v_layerThickness_[j];
935 : }
936 0 : return Length(h, Length::UnitMeter);
937 : }
938 :
939 0 : Length AtmProfile::getLayerBottomHeightAboveSeaLevel(unsigned int i) const
940 : {
941 0 : if(i > v_layerThickness_.size() - 1) {
942 0 : std::ostringstream oss;
943 0 : oss << "Not a valid layer: " << i;
944 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
945 0 : }
946 0 : double h=altitude_.get(Length::UnitMeter);
947 0 : for(unsigned int j = 0; j < i; j++) {
948 0 : h = h + v_layerThickness_[j];
949 : }
950 0 : return Length(h, Length::UnitMeter);
951 : }
952 :
953 0 : Length AtmProfile::getLayerTopHeightAboveSeaLevel(unsigned int i) const
954 : {
955 0 : if(i > v_layerThickness_.size() - 1) {
956 0 : std::ostringstream oss;
957 0 : oss << "Not a valid layer: " << i;
958 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
959 0 : }
960 0 : double h=altitude_.get(Length::UnitMeter);
961 0 : for(unsigned int j = 0; j < i+1; j++) {
962 0 : h = h + v_layerThickness_[j];
963 : }
964 0 : return Length(h, Length::UnitMeter);
965 : }
966 :
967 :
968 0 : void AtmProfile::setLayerThickness(unsigned int i, const Length &layerThickness)
969 : {
970 0 : if(i < v_layerThickness_.size()) {
971 0 : v_layerThickness_[i] = layerThickness.get(Length::UnitMeter);
972 : }
973 0 : }
974 :
975 290 : MassDensity AtmProfile::getLayerWaterVaporMassDensity(unsigned int i) const
976 : {
977 : /*if(i > v_layerWaterVapor_.size() - 1) {
978 : MassDensity m(-999.0, MassDensity::UnitKiloGramPerCubicMeter);
979 : return m;
980 : } else {
981 : MassDensity m(v_layerWaterVapor_[i], MassDensity::UnitKiloGramPerCubicMeter);
982 : return m;
983 : }*/
984 290 : if(i > v_layerWaterVapor_.size() - 1) {
985 0 : std::ostringstream oss;
986 0 : oss << "Not a valid layer: " << i;
987 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
988 0 : }
989 290 : return MassDensity(v_layerWaterVapor_[i], MassDensity::UnitKiloGramPerCubicMeter);
990 : }
991 0 : MassDensity AtmProfile::getLayerBottomWaterVaporMassDensity(unsigned int i) const
992 : {
993 0 : if(i > v_layerWaterVapor0_.size() - 1) {
994 0 : std::ostringstream oss;
995 0 : oss << "Not a valid layer: " << i;
996 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
997 0 : }
998 0 : return MassDensity(v_layerWaterVapor0_[i], MassDensity::UnitKiloGramPerCubicMeter);
999 : }
1000 0 : MassDensity AtmProfile::getLayerTopWaterVaporMassDensity(unsigned int i) const
1001 : {
1002 0 : if(i > v_layerWaterVapor1_.size() - 1) {
1003 0 : std::ostringstream oss;
1004 0 : oss << "Not a valid layer: " << i;
1005 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1006 0 : }
1007 0 : return MassDensity(v_layerWaterVapor1_[i], MassDensity::UnitKiloGramPerCubicMeter);
1008 : }
1009 :
1010 290 : NumberDensity AtmProfile::getLayerWaterVaporNumberDensity(unsigned int i) const
1011 : {
1012 : /*if(i > v_layerWaterVapor_.size() - 1) {
1013 : NumberDensity m(-999.0, NumberDensity::UnitInverseCubicMeter);
1014 : return m;
1015 : } else {
1016 : NumberDensity
1017 : m(v_layerWaterVapor_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1018 : return m;
1019 : }*/
1020 290 : if(i > v_layerWaterVapor_.size() - 1) {
1021 0 : std::ostringstream oss;
1022 0 : oss << "Not a valid layer: " << i;
1023 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1024 0 : }
1025 290 : return NumberDensity(v_layerWaterVapor_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1026 : }
1027 0 : NumberDensity AtmProfile::getLayerBottomWaterVaporNumberDensity(unsigned int i) const
1028 : {
1029 0 : if(i > v_layerWaterVapor0_.size() - 1) {
1030 0 : std::ostringstream oss;
1031 0 : oss << "Not a valid layer: " << i;
1032 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1033 0 : }
1034 0 : return NumberDensity(v_layerWaterVapor0_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1035 : }
1036 0 : NumberDensity AtmProfile::getLayerTopWaterVaporNumberDensity(unsigned int i) const
1037 : {
1038 0 : if(i > v_layerWaterVapor1_.size() - 1) {
1039 0 : std::ostringstream oss;
1040 0 : oss << "Not a valid layer: " << i;
1041 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1042 0 : }
1043 0 : return NumberDensity(v_layerWaterVapor1_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1044 : }
1045 :
1046 0 : void AtmProfile::setLayerWaterVaporMassDensity(unsigned int i, const MassDensity &layerWaterVapor)
1047 : {
1048 0 : if(i <= v_layerWaterVapor_.size() - 1) {
1049 0 : v_layerWaterVapor_[i] = layerWaterVapor.get(MassDensity::UnitKiloGramPerCubicMeter);
1050 : }
1051 0 : }
1052 :
1053 0 : void AtmProfile::setLayerWaterVaporNumberDensity(unsigned int i, const NumberDensity &layerWaterVapor)
1054 : {
1055 0 : if(i <= v_layerWaterVapor_.size() - 1) {
1056 0 : v_layerWaterVapor_[i] = layerWaterVapor.get(NumberDensity::UnitInverseCubicMeter) * 18.0 / (6.023e23 * 1000.0);
1057 : }
1058 0 : }
1059 :
1060 0 : vector<Pressure> AtmProfile::getPressureProfile() const
1061 : {
1062 0 : vector<Pressure> p;
1063 0 : p.reserve(v_layerPressure_.size());
1064 0 : for(unsigned int i = 0; i < v_layerPressure_.size(); i++) {
1065 0 : Pressure pp(v_layerPressure_[i], Pressure::UnitMilliBar);
1066 0 : p.push_back(pp);
1067 0 : }
1068 0 : return p;
1069 0 : }
1070 :
1071 2680 : Pressure AtmProfile::getLayerPressure(unsigned int i) const
1072 : {
1073 : /*if(i > v_layerPressure_.size() - 1) {
1074 : Pressure p(-999.0, Pressure::UnitMilliBar);
1075 : return p;
1076 : } else {
1077 : Pressure p(v_layerPressure_[i], Pressure::UnitMilliBar);
1078 : return p;
1079 : }*/
1080 2680 : if(i > v_layerPressure_.size() - 1) {
1081 0 : std::ostringstream oss;
1082 0 : oss << "Not a valid layer: " << i;
1083 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1084 0 : }
1085 2680 : return Pressure(v_layerPressure_[i], Pressure::UnitMilliBar);
1086 : }
1087 :
1088 0 : Pressure AtmProfile::getLayerBottomPressure(unsigned int i) const
1089 : {
1090 0 : if(i > v_layerPressure0_.size() - 1) {
1091 0 : std::ostringstream oss;
1092 0 : oss << "Not a valid layer: " << i;
1093 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1094 0 : }
1095 0 : return Pressure(v_layerPressure0_[i], Pressure::UnitMilliBar);
1096 : }
1097 :
1098 0 : Pressure AtmProfile::getLayerTopPressure(unsigned int i) const
1099 : {
1100 0 : if(i > v_layerPressure1_.size() - 1) {
1101 0 : std::ostringstream oss;
1102 0 : oss << "Not a valid layer: " << i;
1103 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1104 0 : }
1105 0 : return Pressure(v_layerPressure1_[i], Pressure::UnitMilliBar);
1106 : }
1107 :
1108 :
1109 32060498 : Length AtmProfile::getGroundWH2O() const
1110 : {
1111 32060498 : double wm = 0;
1112 900700477 : for(unsigned int j = 0; j < numLayer_; j++) {
1113 868639979 : wm = wm + v_layerWaterVapor_[j] * v_layerThickness_[j]; // kg/m**2 or mm (from m*kg/m**3 IS units)
1114 : }
1115 32060498 : wm = wm * 1e-3; // (pasar de mm a m)
1116 32060498 : Length wh2o(wm);
1117 32060498 : return wh2o;
1118 : }
1119 :
1120 1015 : MassDensity AtmProfile::rwat(const Temperature &tt, const Humidity &rh, const Pressure &pp) const
1121 : {
1122 1015 : double t = tt.get(Temperature::UnitKelvin);
1123 1015 : double p = pp.get(Pressure::UnitMilliBar);
1124 1015 : double u = rh.get(Percent::UnitPercent);
1125 : double e, es, rwat0;
1126 :
1127 1015 : if(p <= 0 || t <= 0 || u <= 0) {
1128 0 : return MassDensity(0.0, MassDensity::UnitGramPerCubicMeter);
1129 : } else {
1130 1015 : es = 6.105 * exp(25.22 / t * (t - 273.0) - 5.31 * log(t / 273.0));
1131 1015 : e = 1.0 - (1.0 - u / 100.0) * es / p;
1132 1015 : e = es * u / 100.0 / e;
1133 1015 : rwat0 = e * 216.502 / t; //(en g/m*3)
1134 : }
1135 1015 : return MassDensity(rwat0, MassDensity::UnitGramPerCubicMeter);
1136 : }
1137 :
1138 0 : Humidity AtmProfile::rwat_inv(const Temperature &tt, const MassDensity &dd, const Pressure &pp) const
1139 : {
1140 0 : double p = pp.get(Pressure::UnitMilliBar);
1141 0 : double t = tt.get(Temperature::UnitKelvin);
1142 0 : double r = dd.get(MassDensity::UnitGramPerCubicMeter);
1143 : double es, e, rinv;
1144 :
1145 0 : if(p <= 0 || t <= 0 || r <= 0) {
1146 0 : rinv = 0.0;
1147 : } else {
1148 0 : es = 6.105 * exp(25.22 / t * (t - 273.0) - 5.31 * log(t / 273.0));
1149 0 : e = r * t / 216.502;
1150 0 : rinv = 100 * (e * (p - es) / (es * (p - e)));
1151 0 : if(rinv < 0 && p < 3) rinv = 0.0;
1152 : }
1153 0 : return Humidity(rinv, Percent::UnitPercent);
1154 : }
1155 :
1156 29917 : vector<NumberDensity> AtmProfile::st76(const Length &h, unsigned int tip) const
1157 : {
1158 : unsigned int i1, i2, i3, i_layer;
1159 : double x1, x2, x3, d;
1160 29917 : vector<NumberDensity> minorden;
1161 29917 : NumberDensity o3den, n2oden, coden, no2den, so2den;
1162 : static const double avogad = 6.022045E+23;
1163 : static const double airmwt = 28.964;
1164 : // static const double h2omwt=18.015;
1165 29917 : double ha = h.get(Length::UnitKiloMeter);
1166 :
1167 :
1168 : static const double
1169 : alt[50] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1170 : 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0,
1171 : 24.0, 25.0, 27.5, 30.0, 32.5, 35.0, 37.5, 40.0, 42.5, 45.0, 47.5,
1172 : 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0,
1173 : 105.0, 110.0, 115.0, 120.0 }; // REFERENCE LEVELS IN KM
1174 :
1175 : static const double
1176 : ozone[6][50] = { { 2.869E-02, 3.150E-02, 3.342E-02, 3.504E-02, 3.561E-02, 3.767E-02, 3.989E-02, 4.223E-02,
1177 : 4.471E-02, 5.000E-02, 5.595E-02, 6.613E-02, 7.815E-02, 9.289E-02, 1.050E-01, 1.256E-01,
1178 : 1.444E-01, 2.500E-01, 5.000E-01, 9.500E-01, 1.400E+00, 1.800E+00, 2.400E+00, 3.400E+00,
1179 : 4.300E+00, 5.400E+00, 7.800E+00, 9.300E+00, 9.850E+00, 9.700E+00, 8.800E+00, 7.500E+00,
1180 : 5.900E+00, 4.500E+00, 3.450E+00, 2.800E+00, 1.800E+00, 1.100E+00, 6.500E-01, 3.000E-01,
1181 : 1.800E-01, 3.300E-01, 5.000E-01, 5.200E-01, 5.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1182 : 5.000E-03, 5.000E-04 },
1183 : { 3.017E-02, 3.337E-02, 3.694E-02, 4.222E-02, 4.821E-02, 5.512E-02, 6.408E-02, 7.764E-02,
1184 : 9.126E-02, 1.111E-01, 1.304E-01, 1.793E-01, 2.230E-01, 3.000E-01, 4.400E-01, 5.000E-01,
1185 : 6.000E-01, 7.000E-01, 1.000E+00, 1.500E+00, 2.000E+00, 2.400E+00, 2.900E+00, 3.400E+00,
1186 : 4.000E+00, 4.800E+00, 6.000E+00, 7.000E+00, 8.100E+00, 8.900E+00, 8.700E+00, 7.550E+00,
1187 : 5.900E+00, 4.500E+00, 3.500E+00, 2.800E+00, 1.800E+00, 1.300E+00, 8.000E-01, 4.000E-01,
1188 : 1.900E-01, 2.000E-01, 5.700E-01, 7.500E-01, 7.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1189 : 5.000E-03, 5.000E-04 },
1190 : { 2.778E-02, 2.800E-02, 2.849E-02, 3.200E-02, 3.567E-02, 4.720E-02, 5.837E-02, 7.891E-02,
1191 : 1.039E-01, 1.567E-01, 2.370E-01, 3.624E-01, 5.232E-01, 7.036E-01, 8.000E-01, 9.000E-01,
1192 : 1.100E+00, 1.400E+00, 1.800E+00, 2.300E+00, 2.900E+00, 3.500E+00, 3.900E+00, 4.300E+00,
1193 : 4.700E+00, 5.100E+00, 5.600E+00, 6.100E+00, 6.800E+00, 7.100E+00, 7.200E+00, 6.900E+00,
1194 : 5.900E+00, 4.600E+00, 3.700E+00, 2.750E+00, 1.700E+00, 1.000E-00, 5.500E-01, 3.200E-01,
1195 : 2.500E-01, 2.300E-01, 5.500E-01, 8.000E-01, 8.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1196 : 5.000E-03, 5.000E-04 },
1197 : { 2.412E-02, 2.940E-02, 3.379E-02, 3.887E-02, 4.478E-02, 5.328E-02, 6.564E-02, 7.738E-02,
1198 : 9.114E-02, 1.420E-01, 1.890E-01, 3.050E-01, 4.100E-01, 5.000E-01, 6.000E-01, 7.000E-01,
1199 : 8.500E-01, 1.000E+00, 1.300E+00, 1.700E+00, 2.100E+00, 2.700E+00, 3.300E+00, 3.700E+00,
1200 : 4.200E+00, 4.500E+00, 5.300E+00, 5.700E+00, 6.900E+00, 7.700E+00, 7.800E+00, 7.000E+00,
1201 : 5.400E+00, 4.200E+00, 3.200E+00, 2.500E+00, 1.700E+00, 1.200E+00, 8.000E-01, 4.000E-01,
1202 : 2.000E-01, 1.800E-01, 6.500E-01, 9.000E-01, 8.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1203 : 5.000E-03, 5.000E-04 },
1204 : { 1.802E-02, 2.072E-02, 2.336E-02, 2.767E-02, 3.253E-02, 3.801E-02, 4.446E-02, 7.252E-02,
1205 : 1.040E-01, 2.100E-01, 3.000E-01, 3.500E-01, 4.000E-01, 6.500E-01, 9.000E-01, 1.200E+00,
1206 : 1.500E+00, 1.900E+00, 2.450E+00, 3.100E+00, 3.700E+00, 4.000E+00, 4.200E+00, 4.500E+00,
1207 : 4.600E+00, 4.700E+00, 4.900E+00, 5.400E+00, 5.900E+00, 6.200E+00, 6.250E+00, 5.900E+00,
1208 : 5.100E+00, 4.100E+00, 3.000E+00, 2.600E+00, 1.600E+00, 9.500E-01, 6.500E-01, 5.000E-01,
1209 : 3.300E-01, 1.300E-01, 7.500E-01, 8.000E-01, 8.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1210 : 5.000E-03, 5.000E-04 },
1211 : { 2.660E-02, 2.931E-02, 3.237E-02, 3.318E-02, 3.387E-02, 3.768E-02, 4.112E-02, 5.009E-02,
1212 : 5.966E-02, 9.168E-02, 1.313E-01, 2.149E-01, 3.095E-01, 3.846E-01, 5.030E-01, 6.505E-01,
1213 : 8.701E-01, 1.187E+00, 1.587E+00, 2.030E+00, 2.579E+00, 3.028E+00, 3.647E+00, 4.168E+00,
1214 : 4.627E+00, 5.118E+00, 5.803E+00, 6.553E+00, 7.373E+00, 7.837E+00, 7.800E+00, 7.300E+00,
1215 : 6.200E+00, 5.250E+00, 4.100E+00, 3.100E+00, 1.800E+00, 1.100E+00, 7.000E-01, 3.000E-01,
1216 : 2.500E-01, 3.000E-01, 5.000E-01, 7.000E-01, 7.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1217 : 5.000E-03, 5.000E-04 } };
1218 :
1219 : static const double
1220 : n2o[6][50] = { { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1221 : 3.200E-01, 3.195E-01, 3.179E-01, 3.140E-01, 3.095E-01, 3.048E-01, 2.999E-01, 2.944E-01,
1222 : 2.877E-01, 2.783E-01, 2.671E-01, 2.527E-01, 2.365E-01, 2.194E-01, 2.051E-01, 1.967E-01,
1223 : 1.875E-01, 1.756E-01, 1.588E-01, 1.416E-01, 1.165E-01, 9.275E-02, 6.693E-02, 4.513E-02,
1224 : 2.751E-02, 1.591E-02, 9.378E-03, 4.752E-03, 3.000E-03, 2.065E-03, 1.507E-03, 1.149E-03,
1225 : 8.890E-04, 7.056E-04, 5.716E-04, 4.708E-04, 3.932E-04, 3.323E-04, 2.837E-04, 2.443E-04,
1226 : 2.120E-04, 1.851E-04 },
1227 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1228 : 3.195E-01, 3.163E-01, 3.096E-01, 2.989E-01, 2.936E-01, 2.860E-01, 2.800E-01, 2.724E-01,
1229 : 2.611E-01, 2.421E-01, 2.174E-01, 1.843E-01, 1.607E-01, 1.323E-01, 1.146E-01, 1.035E-01,
1230 : 9.622E-02, 8.958E-02, 8.006E-02, 6.698E-02, 4.958E-02, 3.695E-02, 2.519E-02, 1.736E-02,
1231 : 1.158E-02, 7.665E-03, 5.321E-03, 3.215E-03, 2.030E-03, 1.397E-03, 1.020E-03, 7.772E-04,
1232 : 6.257E-04, 5.166E-04, 4.352E-04, 3.727E-04, 3.237E-04, 2.844E-04, 2.524E-04, 2.260E-04,
1233 : 2.039E-04, 1.851E-04 },
1234 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1235 : 3.195E-01, 3.163E-01, 3.096E-01, 2.989E-01, 2.936E-01, 2.860E-01, 2.800E-01, 2.724E-01,
1236 : 2.611E-01, 2.421E-01, 2.174E-01, 1.843E-01, 1.621E-01, 1.362E-01, 1.230E-01, 1.124E-01,
1237 : 1.048E-01, 9.661E-02, 8.693E-02, 7.524E-02, 6.126E-02, 5.116E-02, 3.968E-02, 2.995E-02,
1238 : 2.080E-02, 1.311E-02, 8.071E-03, 4.164E-03, 2.629E-03, 1.809E-03, 1.321E-03, 1.007E-03,
1239 : 7.883E-04, 6.333E-04, 5.194E-04, 4.333E-04, 3.666E-04, 3.140E-04, 2.717E-04, 2.373E-04,
1240 : 2.089E-04, 1.851E-04 },
1241 : { 3.100E-01, 3.100E-01, 3.100E-01, 3.100E-01, 3.079E-01, 3.024E-01, 2.906E-01, 2.822E-01,
1242 : 2.759E-01, 2.703E-01, 2.651E-01, 2.600E-01, 2.549E-01, 2.494E-01, 2.433E-01, 2.355E-01,
1243 : 2.282E-01, 2.179E-01, 2.035E-01, 1.817E-01, 1.567E-01, 1.350E-01, 1.218E-01, 1.102E-01,
1244 : 9.893E-02, 8.775E-02, 7.327E-02, 5.941E-02, 4.154E-02, 3.032E-02, 1.949E-02, 1.274E-02,
1245 : 9.001E-03, 6.286E-03, 4.558E-03, 2.795E-03, 1.765E-03, 1.214E-03, 8.866E-04, 6.756E-04,
1246 : 5.538E-04, 4.649E-04, 3.979E-04, 3.459E-04, 3.047E-04, 2.713E-04, 2.439E-04, 2.210E-04,
1247 : 2.017E-04, 1.851E-04 },
1248 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1249 : 3.195E-01, 3.163E-01, 3.096E-01, 2.989E-01, 2.936E-01, 2.860E-01, 2.800E-01, 2.724E-01,
1250 : 2.611E-01, 2.421E-01, 2.174E-01, 1.843E-01, 1.621E-01, 1.362E-01, 1.230E-01, 1.122E-01,
1251 : 1.043E-01, 9.570E-02, 8.598E-02, 7.314E-02, 5.710E-02, 4.670E-02, 3.439E-02, 2.471E-02,
1252 : 1.631E-02, 1.066E-02, 7.064E-03, 3.972E-03, 2.508E-03, 1.726E-03, 1.260E-03, 9.602E-04,
1253 : 7.554E-04, 6.097E-04, 5.024E-04, 4.210E-04, 3.579E-04, 3.080E-04, 2.678E-04, 2.350E-04,
1254 : 2.079E-04, 1.851E-04 },
1255 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1256 : 3.200E-01, 3.195E-01, 3.179E-01, 3.140E-01, 3.095E-01, 3.048E-01, 2.999E-01, 2.944E-01,
1257 : 2.877E-01, 2.783E-01, 2.671E-01, 2.527E-01, 2.365E-01, 2.194E-01, 2.051E-01, 1.967E-01,
1258 : 1.875E-01, 1.756E-01, 1.588E-01, 1.416E-01, 1.165E-01, 9.275E-02, 6.693E-02, 4.513E-02,
1259 : 2.751E-02, 1.591E-02, 9.378E-03, 4.752E-03, 3.000E-03, 2.065E-03, 1.507E-03, 1.149E-03,
1260 : 8.890E-04, 7.056E-04, 5.716E-04, 4.708E-04, 3.932E-04, 3.323E-04, 2.837E-04, 2.443E-04,
1261 : 2.120E-04, 1.851E-04 } };
1262 :
1263 : static const double
1264 : co[6][50] = { { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1265 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1266 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1267 : 1.400E-02, 1.521E-02, 1.722E-02, 1.995E-02, 2.266E-02, 2.487E-02, 2.738E-02, 3.098E-02,
1268 : 3.510E-02, 3.987E-02, 4.482E-02, 5.092E-02, 5.985E-02, 6.960E-02, 9.188E-02, 1.938E-01,
1269 : 5.688E-01, 1.549E+00, 3.849E+00, 6.590E+00, 1.044E+01, 1.705E+01, 2.471E+01, 3.358E+01,
1270 : 4.148E+01, 5.000E+01 },
1271 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1272 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1273 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1274 : 1.400E-02, 1.521E-02, 1.722E-02, 1.995E-02, 2.266E-02, 2.487E-02, 2.716E-02, 2.962E-02,
1275 : 3.138E-02, 3.307E-02, 3.487E-02, 3.645E-02, 3.923E-02, 4.673E-02, 6.404E-02, 1.177E-01,
1276 : 2.935E-01, 6.815E-01, 1.465E+00, 2.849E+00, 5.166E+00, 1.008E+01, 1.865E+01, 2.863E+01,
1277 : 3.890E+01, 5.000E+01 },
1278 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1279 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1280 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1281 : 1.400E-02, 1.498E-02, 1.598E-02, 1.710E-02, 1.850E-02, 1.997E-02, 2.147E-02, 2.331E-02,
1282 : 2.622E-02, 3.057E-02, 3.803E-02, 6.245E-02, 1.480E-01, 2.926E-01, 5.586E-01, 1.078E+00,
1283 : 1.897E+00, 2.960E+00, 4.526E+00, 6.862E+00, 1.054E+01, 1.709E+01, 2.473E+01, 3.359E+01,
1284 : 4.149E+01, 5.000E+01 },
1285 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1286 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1287 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1288 : 1.400E-02, 1.510E-02, 1.649E-02, 1.808E-02, 1.997E-02, 2.183E-02, 2.343E-02, 2.496E-02,
1289 : 2.647E-02, 2.809E-02, 2.999E-02, 3.220E-02, 3.650E-02, 4.589E-02, 6.375E-02, 1.176E-01,
1290 : 3.033E-01, 7.894E-01, 1.823E+00, 3.402E+00, 5.916E+00, 1.043E+01, 1.881E+01, 2.869E+01,
1291 : 3.892E+01, 5.000E+01 },
1292 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1293 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1294 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1295 : 1.400E-02, 1.521E-02, 1.722E-02, 2.037E-02, 2.486E-02, 3.168E-02, 4.429E-02, 6.472E-02,
1296 : 1.041E-01, 1.507E-01, 2.163E-01, 3.141E-01, 4.842E-01, 7.147E-01, 1.067E+00, 1.516E+00,
1297 : 2.166E+00, 3.060E+00, 4.564E+00, 6.877E+00, 1.055E+01, 1.710E+01, 2.473E+01, 3.359E+01,
1298 : 4.149E+01, 5.000E+01 },
1299 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1300 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1301 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1302 : 1.400E-02, 1.498E-02, 1.598E-02, 1.710E-02, 1.850E-02, 2.009E-02, 2.220E-02, 2.497E-02,
1303 : 2.824E-02, 3.241E-02, 3.717E-02, 4.597E-02, 6.639E-02, 1.073E-01, 1.862E-01, 3.059E-01,
1304 : 6.375E-01, 1.497E+00, 3.239E+00, 5.843E+00, 1.013E+01, 1.692E+01, 2.467E+01, 3.356E+01,
1305 : 4.148E+01, 5.000E+01 } };
1306 :
1307 : static const double
1308 : den[6][50] = { { 2.450E+19, 2.231E+19, 2.028E+19, 1.827E+19, 1.656E+19, 1.499E+19, 1.353E+19, 1.218E+19,
1309 : 1.095E+19, 9.789E+18, 8.747E+18, 7.780E+18, 6.904E+18, 6.079E+18, 5.377E+18, 4.697E+18,
1310 : 4.084E+18, 3.486E+18, 2.877E+18, 2.381E+18, 1.981E+18, 1.651E+18, 1.381E+18, 1.169E+18,
1311 : 9.920E+17, 8.413E+17, 5.629E+17, 3.807E+17, 2.598E+17, 1.789E+17, 1.243E+17, 8.703E+16,
1312 : 6.147E+16, 4.352E+16, 3.119E+16, 2.291E+16, 1.255E+16, 6.844E+15, 3.716E+15, 1.920E+15,
1313 : 9.338E+14, 4.314E+14, 1.801E+14, 7.043E+13, 2.706E+13, 1.098E+13, 4.445E+12, 1.941E+12,
1314 : 8.706E+11, 4.225E+11 },
1315 : { .496E+19, 2.257E+19, 2.038E+19, 1.843E+19, 1.666E+19, 1.503E+19, 1.351E+19, 1.212E+19,
1316 : 1.086E+19, 9.716E+18, 8.656E+18, 7.698E+18, 6.814E+18, 6.012E+18, 5.141E+18, 4.368E+18,
1317 : 3.730E+18, 3.192E+18, 2.715E+18, 2.312E+18, 1.967E+18, 1.677E+18, 1.429E+18, 1.223E+18,
1318 : 1.042E+18, 8.919E+17, 6.050E+17, 4.094E+17, 2.820E+17, 1.927E+17, 1.338E+17, 9.373E+16,
1319 : 6.624E+16, 4.726E+16, 3.398E+16, 2.500E+16, 1.386E+16, 7.668E+15, 4.196E+15, 2.227E+15,
1320 : 1.109E+15, 4.996E+14, 1.967E+14, 7.204E+13, 2.541E+13, 9.816E+12, 3.816E+12, 1.688E+12,
1321 : 8.145E+11, 4.330E+11 },
1322 : { .711E+19, 2.420E+19, 2.158E+19, 1.922E+19, 1.724E+19, 1.542E+19, 1.376E+19, 1.225E+19,
1323 : 1.086E+19, 9.612E+18, 8.472E+18, 7.271E+18, 6.237E+18, 5.351E+18, 4.588E+18, 3.931E+18,
1324 : 3.368E+18, 2.886E+18, 2.473E+18, 2.115E+18, 1.809E+18, 1.543E+18, 1.317E+18, 1.125E+18,
1325 : 9.633E+17, 8.218E+17, 5.536E+17, 3.701E+17, 2.486E+17, 1.647E+17, 1.108E+17, 7.540E+16,
1326 : 5.202E+16, 3.617E+16, 2.570E+16, 1.863E+16, 1.007E+16, 5.433E+15, 2.858E+15, 1.477E+15,
1327 : 7.301E+14, 3.553E+14, 1.654E+14, 7.194E+13, 3.052E+13, 1.351E+13, 6.114E+12, 2.952E+12,
1328 : 1.479E+12, 7.836E+11 },
1329 : { .549E+19, 2.305E+19, 2.080E+19, 1.873E+19, 1.682E+19, 1.508E+19, 1.357E+19, 1.216E+19,
1330 : 1.088E+19, 9.701E+18, 8.616E+18, 7.402E+18, 6.363E+18, 5.471E+18, 4.699E+18, 4.055E+18,
1331 : 3.476E+18, 2.987E+18, 2.568E+18, 2.208E+18, 1.899E+18, 1.632E+18, 1.403E+18, 1.207E+18,
1332 : 1.033E+18, 8.834E+17, 6.034E+17, 4.131E+17, 2.839E+17, 1.938E+17, 1.344E+17, 9.402E+16,
1333 : 6.670E+16, 4.821E+16, 3.516E+16, 2.581E+16, 1.421E+16, 7.946E+15, 4.445E+15, 2.376E+15,
1334 : 1.198E+15, 5.311E+14, 2.022E+14, 7.221E+13, 2.484E+13, 9.441E+12, 3.624E+12, 1.610E+12,
1335 : 7.951E+11, 4.311E+11 },
1336 : { .855E+19, 2.484E+19, 2.202E+19, 1.950E+19, 1.736E+19, 1.552E+19, 1.383E+19, 1.229E+19,
1337 : 1.087E+19, 9.440E+18, 8.069E+18, 6.898E+18, 5.893E+18, 5.039E+18, 4.308E+18, 3.681E+18,
1338 : 3.156E+18, 2.704E+18, 2.316E+18, 1.982E+18, 1.697E+18, 1.451E+18, 1.241E+18, 1.061E+18,
1339 : 9.065E+17, 7.742E+17, 5.134E+17, 3.423E+17, 2.292E+17, 1.533E+17, 1.025E+17, 6.927E+16,
1340 : 4.726E+16, 3.266E+16, 2.261E+16, 1.599E+16, 8.364E+15, 4.478E+15, 2.305E+15, 1.181E+15,
1341 : 6.176E+14, 3.127E+14, 1.531E+14, 7.244E+13, 3.116E+13, 1.403E+13, 6.412E+12, 3.099E+12,
1342 : 1.507E+12, 7.814E+11 },
1343 : { .548E+19, 2.313E+19, 2.094E+19, 1.891E+19, 1.704E+19, 1.532E+19, 1.373E+19, 1.228E+19,
1344 : 1.094E+19, 9.719E+18, 8.602E+18, 7.589E+18, 6.489E+18, 5.546E+18, 4.739E+18, 4.050E+18,
1345 : 3.462E+18, 2.960E+18, 2.530E+18, 2.163E+18, 1.849E+18, 1.575E+18, 1.342E+18, 1.144E+18,
1346 : 9.765E+17, 8.337E+17, 5.640E+17, 3.830E+17, 2.524E+17, 1.761E+17, 1.238E+17, 8.310E+16,
1347 : 5.803E+16, 4.090E+16, 2.920E+16, 2.136E+16, 1.181E+16, 6.426E+15, 3.386E+15, 1.723E+15,
1348 : 8.347E+14, 3.832E+14, 1.711E+14, 7.136E+13, 2.924E+13, 1.189E+13, 5.033E+12, 2.144E+12,
1349 : 9.688E+11, 5.114E+11 } };
1350 :
1351 : static const double so2[50] = { 3.00E-04, 2.74E-04, 2.36E-04, 1.90E-04, 1.46E-04,
1352 : 1.18E-04, 9.71E-05, 8.30E-05, 7.21E-05, 6.56E-05,
1353 : 6.08E-05, 5.79E-05, 5.60E-05, 5.59E-05, 5.64E-05,
1354 : 5.75E-05, 5.75E-05, 5.37E-05, 4.78E-05, 3.97E-05,
1355 : 3.19E-05, 2.67E-05, 2.28E-05, 2.07E-05, 1.90E-05,
1356 : 1.75E-05, 1.54E-05, 1.34E-05, 1.21E-05, 1.16E-05,
1357 : 1.21E-05, 1.36E-05, 1.65E-05, 2.10E-05, 2.77E-05,
1358 : 3.56E-05, 4.59E-05, 5.15E-05, 5.11E-05, 4.32E-05,
1359 : 2.83E-05, 1.33E-05, 5.56E-06, 2.24E-06, 8.96E-07,
1360 : 3.58E-07, 1.43E-07, 5.73E-08, 2.29E-08, 9.17E-09};
1361 :
1362 : static const double no2[50] = { 2.30E-05, 2.30E-05, 2.30E-05, 2.30E-05, 2.30E-05,
1363 : 2.30E-05, 2.30E-05, 2.30E-05, 2.30E-05, 2.32E-05,
1364 : 2.38E-05, 2.62E-05, 3.15E-05, 4.45E-05, 7.48E-05,
1365 : 1.71E-04, 3.19E-04, 5.19E-04, 7.71E-04, 1.06E-03,
1366 : 1.39E-03, 1.76E-03, 2.16E-03, 2.58E-03, 3.06E-03,
1367 : 3.74E-03, 4.81E-03, 6.16E-03, 7.21E-03, 7.28E-03,
1368 : 6.26E-03, 4.03E-03, 2.17E-03, 1.15E-03, 6.66E-04,
1369 : 4.43E-04, 3.39E-04, 2.85E-04, 2.53E-04, 2.31E-04,
1370 : 2.15E-04, 2.02E-04, 1.92E-04, 1.83E-04, 1.76E-04,
1371 : 1.70E-04, 1.64E-04, 1.59E-04, 1.55E-04, 1.51E-04};
1372 :
1373 : // --------------------------------------------[ altitude for interpolation ]
1374 :
1375 29917 : if((ha < 0.0) || (ha > 120.0)) {
1376 :
1377 0 : o3den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1378 0 : n2oden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1379 0 : coden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1380 0 : no2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1381 0 : so2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1382 :
1383 : } else {
1384 :
1385 29917 : i1 = 0; i2 = 0; i3 = 0;
1386 29917 : x1 = 0.0; x2 = 0.0; x3 = 0.0;
1387 :
1388 451026 : for(i_layer = 0; i_layer < 50; i_layer++) {
1389 :
1390 451026 : if(ha < alt[i_layer]) {
1391 :
1392 29917 : if(i_layer == 49) {
1393 :
1394 0 : i1 = i_layer - 2; i2 = i_layer - 1; i3 = i_layer;
1395 :
1396 : } else {
1397 :
1398 29917 : if(i_layer == 0) {
1399 :
1400 0 : i1 = i_layer; i2 = i_layer + 1; i3 = i_layer + 2;
1401 :
1402 : } else {
1403 :
1404 29917 : i1 = i_layer - 1; i2 = i_layer; i3 = i_layer + 1;
1405 :
1406 : }
1407 :
1408 : }
1409 :
1410 29917 : x1 = alt[i1]; x2 = alt[i2]; x3 = alt[i3];
1411 29917 : goto calc;
1412 : }
1413 : }
1414 :
1415 0 : calc:
1416 :
1417 29917 : if(x1 == 0.0 && x2 == 0.0 && x3 == 0.0) {
1418 :
1419 0 : o3den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1420 0 : n2oden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1421 0 : coden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1422 0 : no2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1423 0 : so2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1424 :
1425 : } else {
1426 :
1427 : // ----- ---------------------------------------[ pressure ]
1428 :
1429 : // p=poli2(ha,x1,x2,x3,p[tip][i1],p[tip][i2],p[tip][i3]);
1430 :
1431 : // --------------------------------------------[ TEMPERATURE ]
1432 :
1433 : // t=poli2(ha,x1,x2,x3,t[tip][i1],t[tip][i2],t[tip][i3]);
1434 :
1435 : // --------------------------------------------[ DENSITY ]
1436 :
1437 29917 : d = poli2(ha,x1,x2,x3,den[tip - 1][i1],den[tip - 1][i2],den[tip - 1][i3])
1438 29917 : * airmwt * 1e6 / avogad; // en g/m**3
1439 :
1440 : // D=A+B*HA+C*HA2 en cm**-3, pero dividido por 1E19 D=D*airmwt*1e6/avogad // en g/m**3
1441 :
1442 : // --------------------------------------------[ H2O ]
1443 :
1444 : // rr=poli2(ha,x1,x2,x3,agua[tip][i1],agua[tip][i2],agua[tip][i3]); // RR=A+B*HA+C*HA2 en ppmv
1445 : // rr=rr*1e-6*(h2omwt/airmwt)*d; // en g/m**3
1446 :
1447 : // --------------------------------------------[ OZONE ]
1448 :
1449 59834 : o3den= NumberDensity(poli2(ha, x1, x2, x3, ozone[tip - 1][i1], ozone[tip- 1][i2], ozone[tip - 1][i3])
1450 29917 : * 1e-12 * d * avogad / airmwt,NumberDensity::UnitInverseCubicCentiMeter);
1451 29917 : if(i2 < 30){ // FIX 15/11/2018
1452 25997 : o3den=o3den*0.82; // FIX 15/11/2018
1453 : }else{ // FIX 15/11/2018
1454 3920 : o3den=o3den*1.65; // FIX 15/11/2018
1455 : } // FIX 15/11/2018
1456 :
1457 : // std::cout << "ha o3den " << ha << " " << o3den.get() << std::endl;
1458 :
1459 : // OZONO=A+B*HA+C*HA2 // en ppmv
1460 :
1461 : // --------------------------------------------[ N2O ]
1462 :
1463 59834 : n2oden = NumberDensity(poli2(ha, x1, x2, x3, n2o[tip - 1][i1], n2o[tip - 1][i2], n2o[tip - 1][i3])
1464 29917 : * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1465 : // N2O=A+B*HA+C*HA2 // en ppmv
1466 :
1467 : // --------------------------------------------[ NO2 ]
1468 :
1469 29917 : no2den = NumberDensity(poli2(ha, x1, x2, x3, no2[i1], no2[i2], no2[i3]) * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1470 :
1471 : // --------------------------------------------[ SO2 ]
1472 :
1473 29917 : so2den = NumberDensity(poli2(ha, x1, x2, x3, no2[i1], so2[i2], so2[i3]) * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1474 :
1475 : // --------------------------------------------[ CO ]
1476 :
1477 59834 : coden = NumberDensity(poli2(ha,x1,x2,x3,co[tip-1][i1],co[tip-1][i2],co[tip-1][i3])
1478 29917 : * 1e-12 * d * avogad / airmwt,NumberDensity::UnitInverseCubicCentiMeter);
1479 :
1480 : }
1481 :
1482 : }
1483 :
1484 29917 : minorden.push_back(o3den);
1485 29917 : minorden.push_back(n2oden);
1486 29917 : minorden.push_back(coden);
1487 29917 : minorden.push_back(no2den);
1488 29917 : minorden.push_back(so2den);
1489 :
1490 59834 : return minorden;
1491 29917 : }
1492 :
1493 179502 : double AtmProfile::poli2(double ha,
1494 : double x1,
1495 : double x2,
1496 : double x3,
1497 : double y1,
1498 : double y2,
1499 : double y3) const
1500 : {
1501 : double a, b, c;
1502 :
1503 179502 : c = (y3 - y2) * (x2 - x1) - (y2 - y1) * (x3 - x2);
1504 179502 : b = (x2 - x1) * (x3 * x3 - x2 * x2) - (x2 * x2 - x1 * x1) * (x3 - x2);
1505 179502 : c = c / b;
1506 179502 : b = (y2 - y1) - c * (x2 * x2 - x1 * x1);
1507 179502 : b = b / (x2 - x1);
1508 179502 : a = y1 - c * x1 * x1 - b * x1;
1509 :
1510 179502 : return a + b * ha + c * pow(ha, 2);
1511 : }
1512 :
1513 1015 : unsigned int AtmProfile::mkAtmProfile()
1514 : {
1515 1015 : const int MAXNHX = 20;
1516 : static const double
1517 : hx[MAXNHX] = { 9.225, 10.225, 11.225, 12.850, 14.850, 16.850, 18.850, 22.600, 26.600, 30.600,
1518 : 34.850, 40.850, 46.850, 52.850, 58.850, 65.100, 73.100, 81.100, 89.100, 95.600 };
1519 :
1520 : static const double
1521 : px[6][MAXNHX] = { { 0.3190E+03, 0.2768E+03, 0.2391E+03, 0.1864E+03, 0.1354E+03, 0.9613E+02, 0.6833E+02,
1522 : 0.3726E+02, 0.2023E+02, 0.1121E+02, 0.6142E+01, 0.2732E+01, 0.1260E+01, 0.6042E+00,
1523 : 0.2798E+00, 0.1202E+00, 0.3600E-01, 0.9162E-02, 0.2076E-02, 0.6374E-03 },
1524 : { 0.3139E+03, 0.2721E+03, 0.2350E+03, 0.1833E+03, 0.1332E+03, 0.9726E+02, 0.7115E+02,
1525 : 0.3992E+02, 0.2185E+02, 0.1216E+02, 0.6680E+01, 0.2985E+01, 0.1400E+01, 0.6780E+00,
1526 : 0.3178E+00, 0.1380E+00, 0.4163E-01, 0.9881E-02, 0.2010E-02, 0.5804E-03 },
1527 : { 0.2892E+03, 0.2480E+03, 0.2124E+03, 0.1649E+03, 0.1206E+03, 0.8816E+02, 0.6433E+02,
1528 : 0.3558E+02, 0.1901E+02, 0.1014E+02, 0.5316E+01, 0.2255E+01, 0.1022E+01, 0.4814E+00,
1529 : 0.2206E+00, 0.9455E-01, 0.3000E-01, 0.8729E-02, 0.2332E-02, 0.8164E-03 },
1530 : { 0.3006E+03, 0.2587E+03, 0.2223E+03, 0.1739E+03, 0.1288E+03, 0.9495E+02, 0.7018E+02,
1531 : 0.3983E+02, 0.2199E+02, 0.1232E+02, 0.6771E+01, 0.3056E+01, 0.1452E+01, 0.7051E+00,
1532 : 0.3353E+00, 0.1459E+00, 0.4431E-01, 0.1024E-01, 0.1985E-02, 0.5627E-03 },
1533 : { 0.2731E+03, 0.2335E+03, 0.1995E+03, 0.1546E+03, 0.1130E+03, 0.8252E+02, 0.6017E+02,
1534 : 0.3314E+02, 0.1751E+02, 0.9306E+01, 0.4826E+01, 0.1988E+01, 0.8645E+00, 0.4000E+00,
1535 : 0.1819E+00, 0.7866E-01, 0.2639E-01, 0.8264E-02, 0.2364E-02, 0.8439E-03 },
1536 : { 0.2978E+03, 0.2559E+03, 0.2191E+03, 0.1698E+03, 0.1240E+03, 0.9060E+02, 0.6621E+02,
1537 : 0.3688E+02, 0.1999E+02, 0.1087E+02, 0.5862E+01, 0.2565E+01, 0.1182E+01, 0.5572E+00,
1538 : 0.2551E+00, 0.1074E+00, 0.3224E-01, 0.8697E-02, 0.2158E-02, 0.6851E-03 } };
1539 :
1540 : static const double
1541 : tx[6][MAXNHX] = { { 0.2421E+03, 0.2354E+03, 0.2286E+03, 0.2180E+03, 0.2046E+03, 0.1951E+03, 0.2021E+03,
1542 : 0.2160E+03, 0.2250E+03, 0.2336E+03, 0.2428E+03, 0.2558E+03, 0.2686E+03, 0.2667E+03,
1543 : 0.2560E+03, 0.2357E+03, 0.2083E+03, 0.1826E+03, 0.1767E+03, 0.1841E+03 },
1544 : { 0.2403E+03, 0.2338E+03, 0.2273E+03, 0.2167E+03, 0.2157E+03, 0.2157E+03, 0.2177E+03,
1545 : 0.2223E+03, 0.2271E+03, 0.2349E+03, 0.2448E+03, 0.2596E+03, 0.2740E+03, 0.2727E+03,
1546 : 0.2603E+03, 0.2394E+03, 0.2045E+03, 0.1715E+03, 0.1644E+03, 0.1785E+03 },
1547 : { 0.2239E+03, 0.2196E+03, 0.2191E+03, 0.2183E+03, 0.2173E+03, 0.2163E+03, 0.2153E+03,
1548 : 0.2152E+03, 0.2152E+03, 0.2178E+03, 0.2275E+03, 0.2458E+03, 0.2637E+03, 0.2633E+03,
1549 : 0.2531E+03, 0.2407E+03, 0.2243E+03, 0.2072E+03, 0.1991E+03, 0.2091E+03 },
1550 : { 0.2301E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03,
1551 : 0.2252E+03, 0.2298E+03, 0.2361E+03, 0.2468E+03, 0.2648E+03, 0.2756E+03, 0.2763E+03,
1552 : 0.2662E+03, 0.2392E+03, 0.2023E+03, 0.1681E+03, 0.1609E+03, 0.1770E+03 },
1553 : { 0.2172E+03, 0.2172E+03, 0.2172E+03, 0.2172E+03, 0.2172E+03, 0.2161E+03, 0.2149E+03,
1554 : 0.2126E+03, 0.2127E+03, 0.2165E+03, 0.2222E+03, 0.2368E+03, 0.2516E+03, 0.2601E+03,
1555 : 0.2523E+03, 0.2486E+03, 0.2388E+03, 0.2215E+03, 0.2033E+03, 0.2113E+03 },
1556 : { 228.260, 221.838, 216.778, 216.700, 216.700, 216.700, 216.700, 219.200, 223.136,
1557 : 227.340, 236.110, 252.746, 268.936, 265.057, 250.174, 233.026, 212.656, 196.466,
1558 : 187.260, 189.204} };
1559 :
1560 1015 : double T_ground = groundTemperature_.get(Temperature::UnitKelvin); // std::cout<<"T_ground: " << T_ground <<"K"<<endl;
1561 1015 : double P_ground = groundPressure_.get(Pressure::UnitMilliBar); // std::cout<<"P_ground: " << P_ground <<"mb"<<endl;
1562 1015 : double rh = relativeHumidity_.get(Percent::UnitPercent); // std::cout<<"rh: " << rh <<"%"<<endl;
1563 1015 : double h0 = wvScaleHeight_.get(Length::UnitKiloMeter); // std::cout<<"h0: " << h0 <<"km"<<endl;
1564 1015 : double dp = pressureStep_.get(Pressure::UnitMilliBar); // std::cout<<"dp: " << dp <<"mb"<<endl;
1565 1015 : double alti = altitude_.get(Length::UnitKiloMeter); // std::cout<<"alti: " << alti <<"km"<<endl;
1566 1015 : double atmh = topAtmProfile_.get(Length::UnitKiloMeter); // std::cout<<"atmh: " << atmh <<"km"<<endl;
1567 1015 : double dp1 = pressureStepFactor_; // std::cout<<"dp1: " << dp1 <<" "<<endl;
1568 1015 : double dt = tropoLapseRate_; // TODO implementer des unites (K/km) ici localement
1569 : double prLimit;
1570 : double minmin;
1571 :
1572 1015 : if(typeAtm_ == 1) {
1573 246 : prLimit = 100.0; // 100.0 230.2;
1574 769 : } else if(typeAtm_ == 2) {
1575 69 : prLimit = 100.0; // 198.0;
1576 700 : } else if(typeAtm_ == 3) {
1577 700 : prLimit = 100.0; // 300.0
1578 0 : } else if(typeAtm_ == 4) {
1579 0 : prLimit = 100.0; // 311.0
1580 0 : } else if(typeAtm_ == 5) {
1581 0 : prLimit = 100.0; // 332.0
1582 0 : } else if(typeAtm_ == 6) {
1583 0 : prLimit = 100.0; // 332.0
1584 : } else {
1585 0 : typeAtm_ = 6; // Default: US St76 Atm.
1586 0 : prLimit = 100.0; // 250.0;
1587 : }
1588 :
1589 1015 : unsigned int npp = 0; // number of layers initialized
1590 :
1591 1015 : double rt = 6371.2E+0; // Earth radius in km
1592 1015 : double g0 = 9.80665E+0; // Earth gravity at the surface (m/s**2)
1593 :
1594 : // static bool first = true; // [-Wunused_but_set_variable]
1595 :
1596 : unsigned int i;
1597 : unsigned int i0;
1598 : unsigned int j;
1599 : unsigned int k;
1600 : // double wh2o, wh2o0; // [-Wunused_but_set_variable]
1601 : double wgr0;
1602 1015 : double g = g0;
1603 : double www, altura;
1604 : // double abun_ozono, abun_n2o, abun_co;
1605 1015 : NumberDensity ozono, n2o, co;
1606 : double wgr, dh;
1607 : // double humrel; // [-Wunused_but_set_variable]
1608 : // bool errcode;
1609 :
1610 : //int nmaxLayers=40; // FV peut etre devrions nos avoir un garde-fou au cas ou le nb de couches serait stupidement trop grand
1611 :
1612 1015 : vector<double> v_layerPressure;
1613 1015 : vector<double> v_layerPressure0;
1614 1015 : vector<double> v_layerPressure1;
1615 1015 : vector<double> v_layerTemperature;
1616 1015 : vector<double> v_layerTemperature0;
1617 1015 : vector<double> v_layerTemperature1;
1618 1015 : vector<double> v_layerThickness;
1619 1015 : vector<double> v_layerWaterVapor;
1620 1015 : vector<double> v_layerWaterVapor0;
1621 1015 : vector<double> v_layerWaterVapor1;
1622 1015 : vector<double> v_layerO3;
1623 1015 : vector<double> v_layerCO;
1624 1015 : vector<double> v_layerN2O;
1625 1015 : vector<double> v_layerNO2;
1626 1015 : vector<double> v_layerSO2;
1627 :
1628 1015 : vector<double> v_layerO3_aux;
1629 1015 : vector<double> v_layerCO_aux;
1630 1015 : vector<double> v_layerN2O_aux;
1631 1015 : vector<double> v_layerNO2_aux;
1632 1015 : vector<double> v_layerSO2_aux;
1633 1015 : vector<double> v_layerPressure_aux;
1634 1015 : vector<double> v_layerPressure0_aux;
1635 1015 : vector<double> v_layerPressure1_aux;
1636 1015 : vector<double> v_layerTemperature_aux;
1637 1015 : vector<double> v_layerTemperature0_aux;
1638 1015 : vector<double> v_layerTemperature1_aux;
1639 1015 : vector<double> v_layerThickness_aux;
1640 1015 : vector<double> v_layerWaterVapor_aux;
1641 1015 : vector<double> v_layerWaterVapor0_aux;
1642 1015 : vector<double> v_layerWaterVapor1_aux;
1643 1015 : vector<NumberDensity> minorden;
1644 :
1645 1015 : v_layerPressure.push_back(P_ground);
1646 1015 : v_layerThickness.push_back(alti * 1000);
1647 1015 : v_layerTemperature.push_back(T_ground);
1648 :
1649 2030 : wgr = rwat(Temperature(v_layerTemperature[0], Temperature::UnitKelvin),
1650 2030 : Humidity(rh, Percent::UnitPercent),
1651 3045 : Pressure(v_layerPressure[0], Pressure::UnitMilliBar)).get(MassDensity::UnitGramPerCubicMeter);
1652 :
1653 : // Absolute Humidity in gr/m**3 at altitude alti
1654 1015 : wgr0 = wgr * exp(alti / h0); // h0 in km ==> wgr0 would be the absolute humidity (gr/m**3) at sea level
1655 : // wh2o0 = wgr0 * h0; // wh2o0 is the zenith column of water above sea level (kgr/m**2 or mm) // [-Wunused_but_set_variable]
1656 : // wh2o = wh2o0 * exp(-alti / h0); // wh2o is the zenith column of water above alti (kgr/m**2 or mm) // [-Wunused_but_set_variable]
1657 :
1658 1015 : v_layerWaterVapor.push_back(wgr); // in gr/m**3
1659 :
1660 1015 : i0 = 0;
1661 1015 : i = 0;
1662 1015 : j = 0;
1663 1015 : npp = 0;
1664 :
1665 : // std::cout << "layer " << i << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1666 :
1667 1015 : bool control = true;
1668 : // bool control2 = true;
1669 : while(true) {
1670 30932 : i++;
1671 :
1672 : // std::cout << "layer: " << i << " " << v_layerPressure[i - 1] - dp * pow(dp1, i - 1) << " " << prLimit << std::endl;
1673 :
1674 30932 : if( (v_layerPressure[i - 1] - dp * pow(dp1, i - 1) <= prLimit ) ) { //&& (v_layerTemperature[i - 1] <= tx[typeAtm_ - 1][0]) ) { // &&(fabs(prLimit-v_layerPressure[i - 1]-dp*pow(dp1,i-1))>=20.0) ) {
1675 :
1676 11335 : if(control) {
1677 :
1678 : //std::cout << "capa " << i-1 << " temperature=" << v_layerTemperature[i - 1] << " K " << " tx[0]=" << tx[typeAtm_ - 1][0] << std::endl;
1679 :
1680 : // std::cout << "prLimit=" << prLimit << std::endl;
1681 : // std::cout << "aaa " << v_layerPressure[i - 1] - dp * pow(dp1, i - 1) << std::endl;
1682 :
1683 : /* for(k = 0; k < 20; k++) {
1684 : if(control2) {
1685 : if(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) >= px[typeAtm_ - 1][k]) {
1686 : j = k;
1687 : std::cout << "P=" << v_layerPressure[i - 1] - dp * pow(dp1, i - 1)
1688 : << " prLimit=" << prLimit << " px[" << typeAtm_ - 1 << "][" << k << "]="
1689 : << px[typeAtm_ - 1][k] << std::endl;
1690 : control2 = false;
1691 : }
1692 : }
1693 : } */
1694 :
1695 1015 : minmin = 20000.0;
1696 :
1697 21315 : for(k = 0; k < MAXNHX; k++) {
1698 :
1699 35835 : if( (fabs(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) > 1.05*px[typeAtm_ - 1][k])) &&
1700 15535 : (fabs(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) - px[typeAtm_ - 1][k])) <= minmin ) {
1701 :
1702 1015 : j = k;
1703 :
1704 : /* std::cout << "P=" << v_layerPressure[i - 1] - dp * pow(dp1, i - 1)
1705 : << " prLimit=" << prLimit << " px[" << typeAtm_ - 1 << "][" << k << "]="
1706 : << px[typeAtm_ - 1][k] << std::endl; */
1707 1015 : minmin = fabs(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) - px[typeAtm_ - 1][k]);
1708 : /* std::cout << " minmin=" << minmin << std::endl; */
1709 :
1710 : }
1711 :
1712 : }
1713 :
1714 :
1715 : }
1716 :
1717 : // std::cout << "i,j,v_layerPressure.size()-1=" << i << "," << j << "," << v_layerPressure.size() - 1 << std::endl;
1718 :
1719 11335 : if(i0 == 0) i0 = i - 1;
1720 :
1721 : // std::cout << "i0,j " << i0 << "," << j << std::endl;
1722 :
1723 :
1724 11335 : if(i < v_layerPressure.size() - 1) {
1725 0 : if(control) {
1726 0 : int j0 = (j>0) ? j-1:0;
1727 0 : v_layerPressure[i] = px[typeAtm_ - 1][j0];
1728 0 : v_layerTemperature[i] = tx[typeAtm_ - 1][j0];
1729 :
1730 0 : std::cout << "v_layerTemperature[i]=" << v_layerTemperature[i] << std::endl;
1731 : // - tx[typeAtm_ - 1][0] + v_layerTemperature[i0]; COMMENTED OUT 31/5/2017
1732 :
1733 0 : www = v_layerWaterVapor[i - 1] / 1000.0;
1734 0 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1735 0 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1736 0 : v_layerThickness[i] = v_layerThickness[i - 1] + dh;
1737 : } else {
1738 0 : j++;
1739 0 : v_layerPressure[i] = px[typeAtm_ - 1][j - 1];
1740 0 : v_layerTemperature[i] = tx[typeAtm_ - 1][j - 1];
1741 : // - tx[typeAtm_ - 1][0] + v_layerTemperature[i0]; COMMENTED OUT 31/5/2017
1742 0 : if(j < MAXNHX) {
1743 0 : v_layerThickness[i] = (hx[j] - hx[j - 1]) * 1000.0 + v_layerThickness[i - 1];
1744 : }
1745 : }
1746 : // std::cout << "layer " << i << " j=" << j << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1747 0 : v_layerWaterVapor[i] = wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0));
1748 : } else {
1749 11335 : if(control) {
1750 1015 : int j0 = (j>0) ? j-1:0;
1751 1015 : v_layerPressure.push_back(px[typeAtm_ - 1][j0]);
1752 1015 : v_layerTemperature.push_back(tx[typeAtm_ - 1][j0]);
1753 : //- tx[typeAtm_ - 1][0] + v_layerTemperature[i0]); COMMENTED OUT 31/5/2017
1754 :
1755 1015 : www = v_layerWaterVapor[i - 1] / 1000.0;
1756 1015 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1757 1015 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1758 1015 : v_layerThickness.push_back(v_layerThickness[i - 1] + dh);
1759 : } else {
1760 10320 : j++;
1761 10320 : v_layerPressure.push_back(px[typeAtm_ - 1][j - 1]);
1762 10320 : v_layerTemperature.push_back(tx[typeAtm_ - 1][j - 1]);
1763 : // - tx[typeAtm_ - 1][0] + v_layerTemperature[i0]); COMMENTED OUT 31/5/2017
1764 10320 : if(j < MAXNHX) {
1765 10250 : v_layerThickness.push_back((hx[j] - hx[j - 1]) * 1000.0 + v_layerThickness[i - 1]);
1766 : }
1767 :
1768 : }
1769 11335 : v_layerWaterVapor.push_back(wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0)));
1770 : }
1771 : // std::cout << "layer " << i << " j=" << j << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1772 11335 : if(control) {
1773 1015 : tropoLayer_ = i - 1;
1774 : // std::cout << "tropoLayer_=" << tropoLayer_ << std::endl;
1775 1015 : control = false;
1776 : }
1777 : // std::cout << "2 layer " << i << " alt at end of layer = " << v_layerThickness[i] << " atmh=" << atmh*1000 << std::endl;
1778 :
1779 : // besides the thickness threshold, if j goes >=MAXNHX this loop cannot iterate further
1780 11335 : if((j >= MAXNHX) || (v_layerThickness[i] > (atmh * 1000.0)))
1781 1015 : break;
1782 : } else {
1783 :
1784 : // std::cout << "i,j,v_layerPressure.size()-1=" << i << "," << j << "," << v_layerPressure.size() - 1 << std::endl;
1785 :
1786 19597 : if(i > v_layerPressure.size() - 1) {
1787 19597 : v_layerPressure.push_back(v_layerPressure[i - 1] - dp * pow(dp1, i - 1));
1788 :
1789 19597 : www = v_layerWaterVapor[i - 1] / 1000.0; // in kg/m**3
1790 19597 : g = g0 * pow(1. + ((v_layerThickness[i - 1] / 1000.0)) / rt, -2); // gravity corrected for the height
1791 19597 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1792 19597 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1793 19597 : v_layerThickness.push_back(v_layerThickness[i - 1] + dh);
1794 : //std::cout << "layer " << i << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1795 19597 : v_layerTemperature.push_back(v_layerTemperature[i - 1] + dt * dh / 1000.0);
1796 19597 : v_layerWaterVapor.push_back(wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0))); //r[i] in kgr/(m**2*1000m) [gr/m**3]
1797 : } else {
1798 0 : v_layerPressure[i] = v_layerPressure[i - 1] - dp * pow(dp1, i - 1);
1799 0 : www = v_layerWaterVapor[i - 1] / 1000.0; // in kg/m**3
1800 0 : g = g0 * pow(1. + ((v_layerThickness[i - 1] / 1000.0)) / rt, -2); // gravity corrected for the height
1801 0 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1802 0 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1803 0 : v_layerThickness[i] = v_layerThickness[i - 1] + dh;
1804 : //std::cout << "layer " << i << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1805 0 : v_layerTemperature[i] = v_layerTemperature[i - 1] + dt * dh / 1000.0;
1806 0 : v_layerWaterVapor[i] = wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0)); //r[i] in kgr/(m**2*1000m) [gr/m**3]
1807 : }
1808 :
1809 : //std::cout << "dh=" << dh << std::endl;
1810 19597 : if( v_layerTemperature[i] <= tx[typeAtm_ - 1][0] && v_layerPressure[i] <= px[typeAtm_ - 1][0] )
1811 : {
1812 : // std::cout << "v_layerThickness[i]=" << v_layerThickness[i] << " " << v_layerThickness[i-1] <<std::endl;
1813 : // std::cout << "dh=" << dh << std::endl;
1814 : // std::cout << "old prLimit =" << prLimit << std::endl;
1815 1015 : prLimit=v_layerPressure[i]; //14NOV2017
1816 : // std::cout << "new prLimit =" << prLimit << std::endl;
1817 : }
1818 :
1819 : // humrel = rwat_inv(Temperature(v_layerTemperature[i], Temperature::UnitKelvin),
1820 : // MassDensity(v_layerWaterVapor[i], MassDensity::UnitGramPerCubicMeter),
1821 : // Pressure(v_layerPressure[i], Pressure::UnitMilliBar)).get(Percent::UnitPercent);
1822 :
1823 :
1824 : /* cout << "layer " << i
1825 : << " P " << v_layerPressure[i] << " prLimit " << prLimit
1826 : << " T " << v_layerTemperature[i]
1827 : << " Alti (layer top in m) " << v_layerThickness[i]
1828 : << " WaterVapor " << v_layerWaterVapor[i] << std::endl; */
1829 : // std::cout << "1 layer " << i << " alt at end of layer = " << v_layerThickness[i] << " atmh=" << atmh*1000 << std::endl;
1830 19597 : if(v_layerThickness[i] > (atmh * 1000.0)) break;
1831 : }
1832 : // std::cout << "v_layerPressure[" << i << "]" << v_layerPressure[i] << " v_layerTemperature[" << i << "] " << v_layerTemperature[i] << std::endl;
1833 29917 : }
1834 :
1835 1015 : npp = i - 1;
1836 :
1837 : // std::cout << "npp=" << npp << std::endl;
1838 : // std::cout << "tropoLayer=" << tropoLayer_ << std::endl;
1839 1015 : tropoTemperature_ = Temperature(v_layerTemperature[tropoLayer_], Temperature::UnitKelvin);
1840 1015 : tropoAltitude_ = Length(v_layerThickness[tropoLayer_], Length::UnitMeter);
1841 : // std::cout << "tropoAltitude=" << tropoAltitude_.get(Length::UnitKiloMeter) << " km" << std::endl;
1842 : // std::cout << "tropoTemperature=" << tropoTemperature_.get(Temperature::UnitKelvin) << " K" << std::endl;
1843 : // std::cout << "ground Altitude=" << altitude_.get(Length::UnitKiloMeter) << " km" << std::endl;
1844 : // std::cout << "ground Temperature=" << groundTemperature_.get(Temperature::UnitKelvin) << " K" << std::endl;
1845 : // std::cout << "Calculated Lapse Rate=" << (tropoTemperature_.get(Temperature::UnitKelvin)-groundTemperature_.get(Temperature::UnitKelvin))/(tropoAltitude_.get(Length::UnitKiloMeter)-altitude_.get(Length::UnitKiloMeter)) << " K/km" << std::endl;
1846 :
1847 :
1848 1015 : altura = alti;
1849 :
1850 : /*
1851 : if(first){
1852 : v_layerO3_.reserve(npp);
1853 : v_layerCO_.reserve(npp);
1854 : v_layerN2O_.reserve(npp);
1855 : v_layerNO2_.reserve(npp);
1856 : v_layerSO2_.reserve(npp);
1857 : }
1858 : */
1859 :
1860 :
1861 30932 : for(unsigned int jj = 0; jj < npp; jj++) {
1862 :
1863 29917 : v_layerTemperature0.push_back(v_layerTemperature[jj]);
1864 29917 : v_layerTemperature1.push_back(v_layerTemperature[jj + 1]);
1865 29917 : v_layerPressure0.push_back(v_layerPressure[jj]);
1866 29917 : v_layerPressure1.push_back(v_layerPressure[jj + 1]);
1867 29917 : v_layerWaterVapor0.push_back(1.0E-3*v_layerWaterVapor[jj]);
1868 29917 : v_layerWaterVapor1.push_back(1.0E-3*v_layerWaterVapor[jj + 1]);
1869 :
1870 : }
1871 :
1872 30932 : for(j = 0; j < npp; j++) {
1873 :
1874 29917 : v_layerThickness[j] = (v_layerThickness[j + 1] - v_layerThickness[j]); // in m
1875 29917 : altura = altura + v_layerThickness[j] / 2.0E3; // in km
1876 29917 : v_layerTemperature[j] = (v_layerTemperature[j + 1] + v_layerTemperature[j]) // AVERAGING MAY/30/2017
1877 29917 : / 2.; // in K
1878 29917 : v_layerPressure[j] = exp((log(v_layerPressure[j + 1])
1879 29917 : + log(v_layerPressure[j])) / 2.0); // in mb
1880 29917 : v_layerWaterVapor[j] = 1.0E-3 * exp((log(v_layerWaterVapor[j + 1]) + // in kg/m**3
1881 29917 : log(v_layerWaterVapor[j])) / 2.0); // 1.0E-3 ?
1882 :
1883 :
1884 : // std::cout << "type_=" << type_ << std::endl;
1885 : // std::cout << "typeAtm_=" << typeAtm_ << std::endl;
1886 :
1887 29917 : unsigned int atmType = typeAtm_; // conversion in int
1888 :
1889 : // std::cout << "going to minorden with atmType=" << atmType << std::endl;
1890 :
1891 29917 : minorden = st76(Length(altura, Length::UnitKiloMeter), atmType);
1892 :
1893 : // std::cout << "Ozone: " << abun_ozono << " " << ozono.get(NumberDensity::UnitInverseCubicCentiMeter) << std::endl;
1894 : // std::cout << "N2O : " << abun_n2o << " " << n2o.get(NumberDensity::UnitInverseCubicCentiMeter) << std::endl;
1895 : // std::cout << "CO : " << abun_co << " " << co.get(NumberDensity::UnitInverseCubicCentiMeter) << std::endl;
1896 :
1897 :
1898 : /* std::cout << j << " " << v_layerO3_.size()-1 << std::endl; */
1899 :
1900 : /* if(j>v_layerO3_.size()-1){
1901 : v_layerO3_.push_back( 1.E6*abun_ozono); // in m**-3
1902 : v_layerCO_.push_back( 1.E6*abun_co ); // in m**-3
1903 : v_layerN2O_.push_back( 1.E6*abun_n2o ); // in m**-3
1904 : v_layerNO2_.push_back( 1.E6*abun_no2 ); // in m**-3
1905 : v_layerSO2_.push_back( 1.E6*abun_so2 ); // in m**-3
1906 : std::cout << "uno" << std::endl;
1907 : }else{ */
1908 :
1909 : /* v_layerO3_[j] = 1.E6*abun_ozono; // in m**-3
1910 : v_layerCO_[j] = 1.E6*abun_co; // in m**-3
1911 : v_layerN2O_[j] = 1.E6*abun_n2o; // in m**-3
1912 : v_layerNO2_[j] = 1.E6*abun_no2; // in m**-3
1913 : v_layerSO2_[j] = 1.E6*abun_so2; // in m**-3 */
1914 :
1915 : // v_layerO3.push_back(1.E6*abun_ozono); // in m**-3
1916 : // v_layerCO.push_back(1.E6*abun_co); // in m**-3
1917 : // v_layerN2O.push_back(1.E6*abun_n2o); // in m**-3
1918 : // v_layerNO2.push_back(1.E6*abun_no2); // in m**-3
1919 : // v_layerSO2.push_back(1.E6*abun_so2); // in m**-3
1920 :
1921 29917 : v_layerO3.push_back(1.E6 * minorden[0].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
1922 29917 : v_layerCO.push_back(1.E6 * minorden[2].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
1923 29917 : v_layerN2O.push_back(1.E6 * minorden[1].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
1924 29917 : v_layerNO2.push_back(1.E6 * minorden[3].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
1925 29917 : v_layerSO2.push_back(1.E6 * minorden[4].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
1926 :
1927 :
1928 :
1929 : /* } */
1930 :
1931 29917 : altura = altura + v_layerThickness[j] / 2.0E3;
1932 :
1933 : /* std::cout << "j=" << j << "v_layerThickness.size()=" << v_layerThickness.size() << std::endl; */
1934 :
1935 : }
1936 :
1937 : /* std::cout << "j=" << j << " v_layerThickness_aux.size()= " << v_layerThickness_aux.size() << std::endl; */
1938 :
1939 : /* if(v_layerThickness.size()>npp){
1940 : for(j=npp; j<v_layerThickness.size(); j++){
1941 : delete v_layerThickness[j];
1942 : }
1943 : } */
1944 :
1945 30932 : for(j = 0; j < npp; j++) {
1946 :
1947 29917 : v_layerPressure_aux.push_back(v_layerPressure[j]);
1948 29917 : v_layerPressure0_aux.push_back(v_layerPressure0[j]);
1949 29917 : v_layerPressure1_aux.push_back(v_layerPressure1[j]);
1950 29917 : v_layerTemperature_aux.push_back(v_layerTemperature[j]);
1951 29917 : v_layerTemperature0_aux.push_back(v_layerTemperature0[j]);
1952 29917 : v_layerTemperature1_aux.push_back(v_layerTemperature1[j]);
1953 29917 : v_layerThickness_aux.push_back(v_layerThickness[j]);
1954 29917 : v_layerWaterVapor_aux.push_back(v_layerWaterVapor[j]);
1955 29917 : v_layerWaterVapor0_aux.push_back(v_layerWaterVapor0[j]);
1956 29917 : v_layerWaterVapor1_aux.push_back(v_layerWaterVapor1[j]);
1957 29917 : v_layerO3_aux.push_back(v_layerO3[j]);
1958 29917 : v_layerCO_aux.push_back(v_layerCO[j]);
1959 29917 : v_layerN2O_aux.push_back(v_layerN2O[j]);
1960 29917 : v_layerNO2_aux.push_back(v_layerNO2[j]);
1961 29917 : v_layerSO2_aux.push_back(v_layerSO2[j]);
1962 :
1963 : }
1964 :
1965 1015 : if(j > v_layerPressure_.size() - 1) { // ?????
1966 123 : v_layerPressure_.reserve(npp);
1967 123 : v_layerPressure0_.reserve(npp);
1968 123 : v_layerPressure1_.reserve(npp);
1969 123 : v_layerTemperature_.reserve(npp);
1970 123 : v_layerTemperature0_.reserve(npp);
1971 123 : v_layerTemperature1_.reserve(npp);
1972 123 : v_layerThickness_.reserve(npp);
1973 123 : v_layerWaterVapor_.reserve(npp);
1974 123 : v_layerWaterVapor0_.reserve(npp);
1975 123 : v_layerWaterVapor1_.reserve(npp);
1976 123 : v_layerO3_.reserve(npp);
1977 123 : v_layerCO_.reserve(npp);
1978 123 : v_layerN2O_.reserve(npp);
1979 123 : v_layerNO2_.reserve(npp);
1980 123 : v_layerSO2_.reserve(npp);
1981 : }
1982 :
1983 1015 : v_layerPressure_ = v_layerPressure_aux;
1984 1015 : v_layerPressure0_ = v_layerPressure0_aux;
1985 1015 : v_layerPressure1_ = v_layerPressure1_aux;
1986 1015 : v_layerTemperature_ = v_layerTemperature_aux;
1987 1015 : v_layerTemperature0_ = v_layerTemperature0_aux;
1988 1015 : v_layerTemperature1_ = v_layerTemperature1_aux;
1989 1015 : v_layerThickness_ = v_layerThickness_aux;
1990 1015 : v_layerWaterVapor_ = v_layerWaterVapor_aux;
1991 1015 : v_layerWaterVapor0_ = v_layerWaterVapor0_aux;
1992 1015 : v_layerWaterVapor1_ = v_layerWaterVapor1_aux;
1993 1015 : v_layerO3_ = v_layerO3_aux;
1994 1015 : v_layerCO_ = v_layerCO_aux;
1995 1015 : v_layerN2O_ = v_layerN2O_aux;
1996 1015 : v_layerNO2_ = v_layerNO2_aux;
1997 1015 : v_layerSO2_ = v_layerSO2_aux;
1998 :
1999 : // first = false; // ????? [-Wunused_but_set_variable]
2000 :
2001 1015 : return npp;
2002 1015 : }
2003 :
2004 : ATM_NAMESPACE_END
2005 :
|
