← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to imagery/i.atcorr/altitude.cpp

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
imagery/i.atcorr/altitude.cpp
 
1
#include "common.h"
 
2
#include "altitude.h"
 
3
#include "atmosmodel.h"
 
4
#include "aerosolconcentration.h"
 
5
 
 
6
/* Update the atmospheric profile (P(z),T(z),H2O(z),O3(z)) in case the target 
 
7
   is not at sea level.
 
8
 
 
9
 
 
10
   Given the altitude of the target in kilometers as input, we transform the
 
11
   original atmospheric profile (Pressure, Temperature, Water Vapor, Ozone) 
 
12
   so that first level of the new profile is the one at the target altitude. 
 
13
   We also compute the new integrated content in water vapor and ozone, that
 
14
   are used as outputs or in computations when the user chooses to enter a
 
15
   specific amount of Ozone and Water Vapor.
 
16
*/
 
17
void Altitude::pressure(AtmosModel& atms, float& uw, float& uo3)
 
18
{
 
19
    /* log linear interpolation */
 
20
    if(xps >= 100) xps = 99.99f;
 
21
                
 
22
    int i;
 
23
    for(i = 0; atms.z[i] <= xps; i++);
 
24
                
 
25
    int isup = i;
 
26
    int iinf = i - 1;
 
27
 
 
28
    float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
 
29
    float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
 
30
    float ps = (float)exp((xps - xb) / xa);
 
31
 
 
32
    /* interpolating temperature wator vapor and ozone profile versus altitude */
 
33
    float xalt = xps;
 
34
    float xtemp = (atms.t[isup] - atms.t[iinf]) / (atms.z[isup] - atms.z[iinf]);
 
35
    xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
 
36
    float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
 
37
    xwo = xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
 
38
    float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
 
39
    xwh = xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];
 
40
 
 
41
    /* updating atmospheric profile
 
42
       1rst level: target     , complete to 34
 
43
       with interpolated layers */
 
44
    atms.z[0] = xalt;                                                          
 
45
    atms.p[0] = ps;
 
46
    atms.t[0] = xtemp;
 
47
    atms.wh[0] = xwh;
 
48
    atms.wo[0] = xwo;
 
49
 
 
50
    for (i = 1; i < 33 - iinf; ++i)
 
51
    {
 
52
        atms.z[i] = atms.z[i + iinf];
 
53
        atms.p[i] = atms.p[i + iinf];
 
54
        atms.t[i] = atms.t[i + iinf];
 
55
        atms.wh[i] = atms.wh[i + iinf];
 
56
        atms.wo[i] = atms.wo[i + iinf];
 
57
    }
 
58
    int l = 33 - iinf - 1;
 
59
    for (i = l; i < 34; ++i)
 
60
    {
 
61
        atms.z[i] = (atms.z[33] - atms.z[l]) * (i - l) / (33 - l) + atms.z[l];
 
62
        atms.p[i] = (atms.p[33] - atms.p[l]) * (i - l) / (33 - l) + atms.p[l];
 
63
        atms.t[i] = (atms.t[33] - atms.t[l]) * (i - l) / (33 - l) + atms.t[l];
 
64
        atms.wh[i] = (atms.wh[33] - atms.wh[l]) * (i - l) / (33 - l) + atms.wh[l];
 
65
        atms.wo[i] = (atms.wo[33] - atms.wo[l]) * (i - l) / (33 - l) + atms.wo[l];
 
66
    }
 
67
 
 
68
    /* compute modified h2o and o3 integrated content */
 
69
    uw = 0;
 
70
    uo3 = 0;
 
71
    const float g = 98.1f;
 
72
    const float air = 0.028964f/0.0224f;
 
73
    const float ro3 = 0.048f/0.0224f;
 
74
 
 
75
    float rmwh[34];
 
76
    float rmo3[34];
 
77
    int k;
 
78
    for (k = 0; k < 33; ++k)
 
79
    {
 
80
        float roair = air * 273.16f * atms.p[k] / (atms.t[k] * 1013.25f);
 
81
        rmwh[k] = atms.wh[k] / (roair * 1e3f);
 
82
        rmo3[k] = atms.wo[k] / (roair * 1e3f);
 
83
    }
 
84
 
 
85
    for (k = 1; k < 33; ++k)
 
86
    {
 
87
        float ds = (atms.p[k - 1] - atms.p[k]) / atms.p[0];
 
88
        uw += (rmwh[k] + rmwh[k - 1]) * ds / 2.f;
 
89
        uo3 += (rmo3[k] + rmo3[k - 1]) * ds / 2.f;
 
90
    }
 
91
    uw = uw * atms.p[0] * 100.f / g;
 
92
    uo3 = uo3 * atms.p[0] * 100.f / g;
 
93
    uo3 = uo3 * 1e3f / ro3;
 
94
}
 
95
 
 
96
/*
 
97
  Function: Update the atmospheric profile (P(z),T(z),H2O(z),O3(z)) in case the observer is on
 
98
  board an aircraft.
 
99
 
 
100
  Description: Given the altitude or pressure at aircraft level as input, the first task is to
 
101
  compute the altitude (in case the pressure has been entered) or the pressure (in case the altitude has
 
102
  been entered) at plane level. Then, a new atmospheric profile is created (Pp,Tp,H2Op,O3p) for which
 
103
  the last level is located at the plane altitude. This profile is used in the gaseous absorption
 
104
  computation (ABSTRA.f) for the path from target to sensor (upward transmission). The ozone and
 
105
  water vapor integrated content of the "plane" atmospheric profile are also an output of this
 
106
  subroutine. The last output is the proportion of molecules below plane level which is useful in
 
107
  scattering computations (OS.f,ISO.f).
 
108
*/
 
109
void Altitude::presplane(AtmosModel& atms)
 
110
{
 
111
    /* log linear interpolation */
 
112
    xpp += atms.z[0];
 
113
    if(xpp >= 100) xpp = 1000;
 
114
 
 
115
    int i;
 
116
    for(i = 0; atms.z[i] <= xpp; i++);
 
117
 
 
118
    int isup = i;
 
119
    int iinf = i-1;
 
120
 
 
121
    float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
 
122
    float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
 
123
    float ps = (float)(exp((xpp - xb) / xa));
 
124
 
 
125
    /* interpolating temperature wator vapor and ozone profile versus altitud */
 
126
    float xalt = xpp;
 
127
    float xtemp  = (atms.t[isup] - atms.t[iinf])/ (atms.z[isup] - atms.z[iinf]);
 
128
    xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
 
129
    float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
 
130
    xwo =  xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
 
131
    float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
 
132
    xwh =  xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];
 
133
 
 
134
    /* updating atmospheric profile
 
135
       last level: plane     , complete to 34
 
136
       with interpolated layers */
 
137
    for(i = 0; i <= iinf; i++)
 
138
    {
 
139
        plane_sim.zpl[i] = atms.z[i];
 
140
        plane_sim.ppl[i] = atms.p[i];
 
141
        plane_sim.tpl[i] = atms.t[i];
 
142
        plane_sim.whpl[i] = atms.wh[i];
 
143
        plane_sim.wopl[i] = atms.wo[i];
 
144
    }
 
145
 
 
146
    for(i = iinf+1; i < 34; i++)
 
147
    {
 
148
        plane_sim.zpl[i] = xalt;
 
149
        plane_sim.ppl[i] = ps;
 
150
        plane_sim.tpl[i] = xtemp;
 
151
        plane_sim.whpl[i] = xwh;
 
152
        plane_sim.wopl[i] = xwo;
 
153
    }
 
154
 
 
155
    /* compute modified h2o and o3 integrated content
 
156
       compute conversion factor for rayleigh optical thickness computation
 
157
       ftray=rp/rt */
 
158
    atms.uw = 0;
 
159
    atms.uo3 = 0;
 
160
    const float g = 98.1f;
 
161
    const float air = 0.028964f/0.0224f;
 
162
    const float ro3 = 0.048f/0.0224f;
 
163
    float rt = 0;
 
164
    float rp = 0;
 
165
 
 
166
    float rmo3[34];
 
167
    float rmwh[34];
 
168
    int k;
 
169
    for(k = 0; k < 33; k++)
 
170
    {
 
171
        float roair = (float)(air * 273.16 * plane_sim.ppl[k] / (1013.25 * plane_sim.tpl[k]));
 
172
        rmwh[k] = atms.wh[k] / (roair * 1000);
 
173
        rmo3[k] = atms.wo[k] / (roair * 1000);
 
174
        rt += (atms.p[k+1] / atms.t[k+1] + atms.p[k] / atms.p[k]) * (atms.z[k+1] - atms.z[k]);
 
175
        rp += (plane_sim.ppl[k+1] / plane_sim.tpl[k+1] + plane_sim.ppl[k] / plane_sim.tpl[k]) 
 
176
            * (plane_sim.zpl[k+1] - plane_sim.zpl[k]);
 
177
    }
 
178
 
 
179
    ftray = rp / rt;
 
180
    for(k = 1; k < 33; k++)
 
181
    {
 
182
        float ds = (plane_sim.ppl[k-1] - plane_sim.ppl[k]) / plane_sim.ppl[0];
 
183
        atms.uw += (rmwh[k] + rmwh[k-1])*ds/2;
 
184
        atms.uo3+= (rmo3[k] + rmo3[k-1])*ds/2;
 
185
    }
 
186
 
 
187
    atms.uw *= plane_sim.ppl[0] * 100 / g;
 
188
    atms.uo3*= plane_sim.ppl[0] * 100 / g;
 
189
    atms.uo3*= 1000 / ro3;
 
190
}
 
191
 
 
192
void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
 
193
{
 
194
    xps = original_xps;
 
195
    xpp = original_xpp;
 
196
 
 
197
    float uwus;
 
198
    float uo3us;
 
199
    if(xps <= 0)
 
200
    {
 
201
        xps = 0;
 
202
        uwus = 1.424f;
 
203
        uo3us = 0.344f;
 
204
    }
 
205
    else if(atms.idatm != 8) pressure(atms, atms.uw, atms.uo3);
 
206
    else pressure(atms, uwus, uo3us);
 
207
 
 
208
    if(xpp <= 0)
 
209
    {
 
210
        /* ground measurement option */
 
211
        palt = 0;
 
212
        pps = atms.p[0];
 
213
        idatmp = 0;
 
214
        original_taer55p = taer55p = 0;
 
215
        puw = 0;
 
216
    }
 
217
    else if(xpp >= 100)
 
218
    {
 
219
        /* satellite case of equivalent */
 
220
        palt = 1000;
 
221
        pps = 0;
 
222
        original_taer55p = taer55p = aerocon.taer55;
 
223
        puw = 0;
 
224
        ftray = 1;
 
225
        idatmp = 4;
 
226
    }
 
227
    else
 
228
    {
 
229
        /* "real" plane case */
 
230
        cin >> original_puw;
 
231
        cin >> original_puo3;
 
232
        cin.ignore(numeric_limits < int >::max(), '\n');        /* ignore comments */
 
233
 
 
234
        puw = original_puw;
 
235
        puo3 = original_puo3;
 
236
        if ( puw < 0 )
 
237
        {
 
238
            presplane(atms);
 
239
            idatmp = 2;
 
240
 
 
241
            if (atms.idatm == 8)
 
242
            {
 
243
                puwus = puw;
 
244
                puo3us = puo3;
 
245
                puw *= atms.uw / uwus;
 
246
                puo3 *= atms.uo3 / uo3us;
 
247
                idatmp = 8;
 
248
            }
 
249
        }
 
250
        else
 
251
        {
 
252
            presplane(atms);
 
253
            idatmp = 8;
 
254
        }
 
255
 
 
256
        palt = plane_sim.zpl[33] - atms.z[0];
 
257
        pps = plane_sim.ppl[33];
 
258
        cin >> original_taer55p;
 
259
        taer55p = original_taer55p;
 
260
 
 
261
        if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03))
 
262
        {
 
263
            /* a scale heigh of 2km is assumed in case no value is given for taer55p */
 
264
            taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 2)));
 
265
        }
 
266
        else
 
267
        {
 
268
            /* compute effective scale heigh */
 
269
            double sham = exp(-palt / 4);
 
270
            double sha = 1 - (taer55p / aerocon.taer55);
 
271
 
 
272
            if( sha >= sham) taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
 
273
            else {
 
274
                sha = -palt/log(sha);
 
275
                taer55p = (float)(aerocon.taer55 * (1 - exp(-palt/sha)));
 
276
            }
 
277
        }
 
278
    }
 
279
}
 
280
 
 
281
void Altitude::update_hv(AtmosModel & atms, const AerosolConcentration & aerocon)
 
282
{
 
283
    xps = original_xps;
 
284
    xpp = original_xpp;
 
285
 
 
286
    float uwus;
 
287
    float uo3us;
 
288
 
 
289
    if (xps <= 0) {
 
290
        xps = 0;
 
291
        uwus = 1.424f;
 
292
        uo3us = 0.344f;
 
293
    }
 
294
    else if (atms.idatm != 8)
 
295
        pressure(atms, atms.uw, atms.uo3);
 
296
    else
 
297
        pressure(atms, uwus, uo3us);
 
298
 
 
299
    if (xpp <= 0) {
 
300
        /* ground measurement option */
 
301
        palt = 0;
 
302
        pps = atms.p[0];
 
303
        idatmp = 0;
 
304
        taer55p = 0;
 
305
        puw = 0;
 
306
    }
 
307
    else if (xpp >= 100) {
 
308
        /* satellite case of equivalent */
 
309
        palt = 1000;
 
310
        pps = 0;
 
311
        taer55p = aerocon.taer55;
 
312
        puw = 0;
 
313
        ftray = 1;
 
314
        idatmp = 4;
 
315
    }
 
316
    else {
 
317
        /* "real" plane case */
 
318
 
 
319
        puw = original_puw;
 
320
        puo3 = original_puo3;
 
321
 
 
322
        if (puw < 0) {
 
323
            presplane(atms);
 
324
            idatmp = 2;
 
325
 
 
326
            if (atms.idatm == 8) {
 
327
                puwus = puw;
 
328
                puo3us = puo3;
 
329
                puw *= atms.uw / uwus;
 
330
                puo3 *= atms.uo3 / uo3us;
 
331
                idatmp = 8;
 
332
            }
 
333
        }
 
334
        else {
 
335
            presplane(atms);
 
336
            idatmp = 8;
 
337
        }
 
338
 
 
339
        palt = plane_sim.zpl[33] - atms.z[0];
 
340
        pps = plane_sim.ppl[33];
 
341
        taer55p = original_taer55p;
 
342
 
 
343
        if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03)) {
 
344
            /* a scale heigh of 2km is assumed in case no value is given for taer55p */
 
345
            taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 2)));
 
346
        }
 
347
        else {
 
348
            /* compute effective scale heigh */
 
349
            double sham = exp(-palt / 4);
 
350
            double sha = 1 - (taer55p / aerocon.taer55);
 
351
 
 
352
            if (sha >= sham)
 
353
                taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
 
354
            else {
 
355
                sha = -palt / log(sha);
 
356
                taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / sha)));
 
357
            }
 
358
        }
 
359
    }
 
360
}
 
361
 
 
362
void Altitude::parse()
 
363
{
 
364
    cin >> original_xps;
 
365
    cin.ignore(numeric_limits<int>::max(),'\n');        /* ignore comments */
 
366
    original_xps = -original_xps;
 
367
    
 
368
    cin >> original_xpp;
 
369
    cin.ignore(numeric_limits<int>::max(),'\n');        /* ignore comments */
 
370
    original_xpp = -original_xpp;
 
371
}
 
372
 
 
373
/* --- plane simulation output if selected ---- */
 
374
void Altitude::print()
 
375
{
 
376
    if(palt < 1000)
 
377
    {
 
378
        Output::Ln();
 
379
        Output::WriteLn(22," plane simulation description ");
 
380
        Output::WriteLn(22," ---------------------------- ");
 
381
                
 
382
        ostringstream s1;
 
383
        s1.setf(ios::fixed, ios::floatfield);
 
384
        s1.precision(2);
 
385
        s1 << " plane  pressure          [mb] " << setw(9) << pps << ends;
 
386
        Output::WriteLn(10,s1.str());
 
387
 
 
388
        ostringstream s2;
 
389
        s2.setf(ios::fixed, ios::floatfield);
 
390
        s2.precision(3);
 
391
        s2 << " plane  altitude absolute [km] " << setw(9) << plane_sim.zpl[33] << ends;
 
392
        Output::WriteLn(10,s2.str());
 
393
 
 
394
                
 
395
        Output::WriteLn(15," atmosphere under plane description: ");
 
396
 
 
397
        ostringstream s3;
 
398
        s3.setf(ios::fixed, ios::floatfield);
 
399
        s3.precision(3);
 
400
        s3 << " ozone content            " << setw(9) << puo3 << ends;
 
401
        Output::WriteLn(15,s3.str());
 
402
 
 
403
 
 
404
        ostringstream s4;
 
405
        s4.setf(ios::fixed, ios::floatfield);
 
406
        s4.precision(3);
 
407
        s4 << " h2o   content            " << setw(9) << puw << ends;
 
408
        Output::WriteLn(15,s4.str());
 
409
 
 
410
        ostringstream s5;
 
411
        s5.setf(ios::fixed, ios::floatfield);
 
412
        s5.precision(3);
 
413
        s5 << "aerosol opt. thick. 550nm " << setw(9) << taer55p << ends;
 
414
        Output::WriteLn(15,s5.str());
 
415
    }
 
416
}
 
417
 
 
418
Altitude Altitude::Parse()
 
419
{
 
420
    Altitude alt;
 
421
    alt.parse();
 
422
    return alt;
 
423
}