source: grass/trunk/imagery/i.atcorr/computations.cpp

Last change on this file was 73852, checked in by neteler, 6 years ago

bulk fixing of typos in comments (using tools/fix_typos.sh)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-c++src
File size: 48.4 KB
Line 
1#include <cstring>
2
3extern "C" {
4#include <grass/gis.h>
5#include <grass/glocale.h>
6}
7
8#include "common.h"
9#include "geomcond.h"
10#include "atmosmodel.h"
11#include "aerosolmodel.h"
12#include "aerosolconcentration.h"
13#include "altitude.h"
14#include "iwave.h"
15#include "gauss.h"
16#include "transform.h"
17#ifdef WIN32
18#pragma warning (disable : 4305)
19#endif
20
21struct OpticalAtmosProperties
22{
23 double rorayl, romix, roaero;
24 double ddirtr, ddiftr;
25 double ddirtt, ddiftt;
26 double ddirta, ddifta;
27 double udirtr, udiftr;
28 double udirtt, udiftt;
29 double udirta, udifta;
30 double sphalbr, sphalbt, sphalba;
31};
32
33/* To compute the molecular optical depth as a function of wavelength for any
34 atmosphere defined by the pressure and temperature profiles. */
35double odrayl(const AtmosModel &atms, const double wl)
36{
37 /* air refraction index edlen 1966 / metrologia,2,71-80 putting pw=0 */
38 double ak=1/wl;
39 double awl= wl*wl*wl*wl * 0.0254743;
40 double a1 = 130 - ak * ak;
41 double a2 = 38.9 - ak * ak;
42 double a3 = 2406030 / a1;
43 double a4 = 15997 / a2;
44 double an = (8342.13 + a3 + a4) * 1.0e-08 + 1;
45 double a = (24 * M_PI * M_PI * M_PI) * ((an * an - 1) * (an * an - 1))
46 * (6 + 3 * delta) / (6 - 7 * delta) / ((an * an + 2) * (an * an + 2));
47
48 double tray = 0;
49 for(int k = 0; k < 33; k++)
50 {
51 double dppt = (288.15 / 1013.25) * (atms.p[k] / atms.t[k] + atms.p[k+1] / atms.t[k+1]) / 2;
52 double sr = a * dppt / awl;
53 tray += (double)((atms.z[k+1] - atms.z[k]) * sr);
54 }
55
56 return tray;
57}
58
59/*
60 decompose the aerosol phase function in series of Legendre polynomial used in
61 OS.f and ISO.f and compute truncation coefficient f to modify aerosol optical thickness t and single
62 scattering albedo w0 according to:
63 t' = (1-w0 f) t
64
65 w0' = w0 (1- f)
66 --------
67 (1-w0 f)
68*/
69double trunca()
70{
71 double ptemp[83];
72 double cosang[80];
73 double weight[80];
74 double rmu[83];
75 double ga[83];
76
77 int i;
78 for(i = 0; i < 83; i++) ptemp[i] = sixs_trunc.pha[i];
79
80 Gauss::gauss(-1,1,cosang,weight,80);
81
82 for(i = 0; i < 40; i++)
83 {
84 rmu[i+1] = cosang[i];
85 ga[i+1] = weight[i];
86 }
87
88 rmu[0] = -1;
89 ga[0] = 0;
90 rmu[41] = 0;
91 ga[41] = 0;
92
93 for(i = 40; i < 80; i++)
94 {
95 rmu[i+2] = cosang[i];
96 ga[i+2] = weight[i];
97 }
98
99 rmu[82] = 1;
100 ga[82] = 0;
101
102 int k = 0;
103 for(i = 0; i < 83; i++)
104 {
105 if(rmu[i] > 0.8) break;
106 k = i - 1;
107 }
108
109 int kk = 0;
110 for(i = 0; i < 83; i++)
111 {
112 if(rmu[i] > 0.94) break;
113 kk = i - 1;
114 }
115
116
117 double aa = (double)((log10(sixs_trunc.pha[kk]) - log10(sixs_trunc.pha[k])) /
118 (acos(rmu[kk]) - acos(rmu[k])));
119 double x1 = (double)(log10(sixs_trunc.pha[kk]));
120 double x2 = (double)acos(rmu[kk]);
121
122 for(i = kk + 1; i < 83; i++)
123 {
124 double a;
125 if(fabs(rmu[i] - 1) <= 1e-08) a = x1 - aa * x2;
126 else a = x1 + aa * (acos(rmu[i]) - x2);
127 ptemp[i] = (double)pow(10,a);
128 }
129
130
131 for(i = 0; i < 83; i++) sixs_trunc.pha[i] = ptemp[i];
132 for(i = 0; i < 80; i++) sixs_trunc.betal[i] = 0;
133
134 double pl[83];
135
136#define IPL(X) ((X)+1)
137
138
139 for(i = 0; i < 83; i++)
140 {
141 double x = sixs_trunc.pha[i] * ga[i];
142 double rm = rmu[i];
143 pl[IPL(-1)] = 0;
144 pl[IPL(0)] = 1;
145
146 for(k = 0; k <= 80; k++)
147 {
148 pl[IPL(k+1)] = ((2 * k + 1) * rm * pl[IPL(k)] - k * pl[IPL(k-1)]) / (k + 1);
149 sixs_trunc.betal[k] += x * pl[IPL(k)];
150 }
151 }
152
153 for(i = 0; i <= 80; i++) sixs_trunc.betal[i] *= (2 * i + 1) * 0.5f;
154
155 double z1 = sixs_trunc.betal[0];
156 for(i = 0; i <= 80; i++) sixs_trunc.betal[i] /= z1;
157 if(sixs_trunc.betal[80] < 0) sixs_trunc.betal[80] = 0;
158
159 return 1 - z1;
160}
161
162
163/*
164 Decompose the atmosphere in a finite number of layers. For each layer, DISCRE
165 provides the optical thickness, the proportion of molecules and aerosols assuming an exponential
166 distribution for each constituents. Figure 1 illustrate the way molecules and aerosols are mixed in a
167 realistic atmosphere. For molecules, the scale height is 8km. For aerosols it is assumed to be 2km
168 unless otherwise specified by the user (using aircraft measurements).
169*/
170
171double discre(const double ta, const double ha, const double tr, const double hr,
172 const int it, const int ntd, const double yy, const double dd,
173 const double ppp2, const double ppp1)
174{
175 if( ha >= 7 )
176 {
177 G_warning(_("Check aerosol measurements or plane altitude"));
178 return 0;
179 }
180
181 double dt;
182 if( it == 0 ) dt = 1e-17;
183 else dt = 2 * (ta + tr - yy) / (ntd - it + 1);
184
185 double zx; /* return value */
186 double ecart = 0;
187 do {
188 dt = dt / 2;
189 double ti = yy + dt;
190 double y1 = ppp2;
191 double y2;
192 double y3 = ppp1;
193
194 while(true)
195 {
196 y2 = (y1 + y3) * 0.5f;
197
198 double xx = -y2 / ha;
199 double x2;
200 if (xx < -18) x2 = tr * exp(-y2 / hr);
201 else x2 = ta * exp(xx) + tr * exp(-y2 / hr);
202
203 if(fabs(ti - x2) < 0.00001) break;
204
205 if(ti - x2 < 0) y3 = y2;
206 else y1 = y2;
207 }
208
209 zx = y2;
210 double cdelta = (double)(1. / (1 + ta * hr / tr / ha * exp((zx - ppp1) * (1. / hr - 1. / ha))));
211 if(dd != 0) ecart = (double)fabs((dd - cdelta) / dd);
212 } while((ecart > 0.75) && (it != 0));
213 return zx;
214
215}
216
217
218/* indexing macro for the psl variable */
219#define PSI(X) ((X)+1)
220/*
221 Compute the values of Legendre polynomials used in the successive order of
222 scattering method.
223*/
224void kernel(const int is, double (&xpl)[2*mu + 1], double (&bp)[26][2*mu + 1], Gauss &gauss)
225{
226 const double rac3 = 1.7320508075688772935274463415059;
227#define PSI(X) ((X)+1)
228 double psl[82][2*mu + 1];
229
230 if(is == 0)
231 {
232 for(int j = 0; j <= mu; j++)
233 {
234 psl[PSI(0)][STDI(-j)] = 1;
235 psl[PSI(0)][STDI(j)] = 1;
236 psl[PSI(1)][STDI(j)] = gauss.rm[STDI(j)];
237 psl[PSI(1)][STDI(-j)] = -gauss.rm[STDI(j)];
238
239 double xdb = (3 * gauss.rm[STDI(j)] * gauss.rm[STDI(j)] - 1) * 0.5;
240 if(fabs(xdb) < 1e-30) xdb = 0;
241 psl[PSI(2)][STDI(-j)] = (double)xdb;
242 psl[PSI(2)][STDI(j)] = (double)xdb;
243 }
244 psl[PSI(1)][STDI(0)] = gauss.rm[STDI(0)];
245 }
246 else if(is == 1)
247 {
248 for(int j = 0; j <= mu; j++)
249 {
250 double x = 1 - gauss.rm[STDI(j)] * gauss.rm[STDI(j)];
251 psl[PSI(0)][STDI(j)] = 0;
252 psl[PSI(0)][STDI(-j)] = 0;
253 psl[PSI(1)][STDI(-j)] = (double)sqrt(x * 0.5);
254 psl[PSI(1)][STDI(j)] = (double)sqrt(x * 0.5);
255 psl[PSI(2)][STDI(j)] = (double)(gauss.rm[STDI(j)] * psl[PSI(1)][STDI(j)] * rac3);
256 psl[PSI(2)][STDI(-j)] = -psl[PSI(2)][STDI(j)];
257
258 }
259 psl[PSI(2)][STDI(0)] = -psl[PSI(2)][STDI(0)];
260 }
261 else
262 {
263 double a = 1;
264 for(int i = 1; i <= is; i++) a *= sqrt((double)(i + is) / (double)i) * 0.5;
265/* double b = a * sqrt((double)is / (is + 1.)) * sqrt((is - 1.) / (is + 2.));*/
266
267 for(int j = 0; j <= mu; j++)
268 {
269 double xx = 1 - gauss.rm[STDI(j)] * gauss.rm[STDI(j)];
270 psl[PSI(is - 1)][STDI(j)] = 0;
271 double xdb = a * pow(xx, is * 0.5);
272 if(fabs(xdb) < 1e-30) xdb = 0;
273 psl[PSI(is)][STDI(-j)] = (double)xdb;
274 psl[PSI(is)][STDI(j)] = (double)xdb;
275 }
276 }
277
278 int k = 2;
279 int ip = 80;
280
281 if(is > 2) k = is;
282 if(k != ip)
283 {
284 int ig = -1;
285 if( is == 1 ) ig = 1;
286 for(int l = k; l < ip; l++)
287 {
288 double a = (2 * l + 1.) / sqrt((l + is + 1.) * (l - is + 1.));
289 double b = sqrt(double((l + is) * (l - is))) / (2. * l + 1.);
290
291 for(int j = 0; j <= mu; j++)
292 {
293 double xdb = a * (gauss.rm[STDI(j)] * psl[PSI(l)][STDI(j)] - b * psl[PSI(l-1)][STDI(j)]);
294 if (fabs(xdb) < 1e-30) xdb = 0;
295 psl[PSI(l+1)][STDI(j)] = (double)xdb;
296 if(j != 0) psl[PSI(l+1)][STDI(-j)] = ig * psl[PSI(l+1)][STDI(j)];
297 }
298 ig = -ig;
299 }
300 }
301
302 int j;
303 for(j = -mu; j <= mu; j++) xpl[STDI(j)] = psl[PSI(2)][STDI(j)];
304
305 for(j = 0; j <= mu; j++)
306 {
307 for(k = -mu; k <= mu; k++)
308 {
309 double sbp = 0;
310 if(is <= 80)
311 {
312 for(int l = is; l <= 80; l++)
313 sbp += psl[PSI(l)][STDI(j)] * psl[PSI(l)][STDI(k)] * sixs_trunc.betal[l];
314
315 if(fabs(sbp) < 1e-30) sbp = 0;
316 bp[j][STDI(k)] = (double)sbp;
317 }
318 }
319 }
320
321}
322
323
324
325#define accu 1e-20
326#define accu2 1e-3
327#define mum1 (mu - 1)
328
329void os(const double tamoy, const double trmoy, const double pizmoy,
330 const double tamoyp, const double trmoyp, double (&xl)[2*mu + 1][np],
331 Gauss &gauss, const Altitude &alt, const GeomCond &geom)
332{
333 double trp = trmoy - trmoyp;
334 double tap = tamoy - tamoyp;
335 int iplane = 0;
336
337 /* if plane observations recompute scale height for aerosol knowing:
338 the aerosol optical depth as measure from the plane = tamoyp
339 the rayleigh scale height = = hr (8km)
340 the rayleigh optical depth at plane level = trmoyp
341 the altitude of the plane = palt
342 the rayleigh optical depth for total atmos = trmoy
343 the aerosol optical depth for total atmos = tamoy
344 if not plane observations then ha is equal to 2.0km
345 ntp local variable: if ntp=nt no plane observation selected
346 ntp=nt-1 plane observation selected
347 it's a mixing rayleigh+aerosol */
348
349 double ha = 2;
350 int snt = nt;
351 int ntp = snt;
352 if(alt.palt <= 900 && alt.palt > 0)
353 {
354 if(tap > 1.e-03) ha = -alt.palt / (double)log(tap / tamoy);
355 ntp = snt - 1;
356 }
357
358 double xmus = -gauss.rm[STDI(0)];
359
360 /* compute mixing rayleigh, aerosol
361 case 1: pure rayleigh
362 case 2: pure aerosol
363 case 3: mixing rayleigh-aerosol */
364
365 double h[31];
366 memset(h, 0, sizeof(h));
367 double ch[31];
368 double ydel[31];
369
370 double xdel[31];
371 double altc[31];
372 if( (tamoy <= accu2) && (trmoy > tamoy) )
373 {
374 for(int j = 0; j <= ntp; j++)
375 {
376 h[j] = j * trmoy / ntp;
377 ch[j]= (double)exp(-h[j] / xmus) / 2;
378 ydel[j] = 1;
379 xdel[j] = 0;
380
381 if (j == 0) altc[j] = 300;
382 else altc[j] = -(double)log(h[j] / trmoy) * 8;
383 }
384 }
385
386
387 if( (trmoy <= accu2) && (tamoy > trmoy) )
388 {
389 for(int j = 0; j <= ntp; j++)
390 {
391 h[j] = j * tamoy / ntp;
392 ch[j]= (double)exp(-h[j] / xmus) / 2;
393 ydel[j] = 0;
394 xdel[j] = pizmoy;
395
396 if (j == 0) altc[j] = 300;
397 else altc[j] = -(double)log(h[j] / tamoy) * ha;
398 }
399 }
400
401 if(trmoy > accu2 && tamoy > accu2)
402 {
403 ydel[0] = 1;
404 xdel[0] = 0;
405 h[0] = 0;
406 ch[0] = 0.5;
407 altc[0] = 300;
408 double zx = 300;
409 iplane = 0;
410
411 for(int it = 0; it <= ntp; it++)
412 {
413 if(it == 0) zx = discre(tamoy, ha, trmoy, 8.0, it, ntp, 0, 0, 300, 0);
414 else zx = discre(tamoy, ha, trmoy, 8.0, it, ntp, h[it - 1], ydel[it - 1], 300, 0);
415
416 double xx = -zx / ha;
417 double ca;
418 if( xx <= -20 ) ca = 0;
419 else ca = tamoy * (double)exp(xx);
420
421 xx = -zx / 8;
422 double cr = trmoy * (double)exp(xx);
423 h[it] = cr + ca;
424
425
426 altc[it] = zx;
427 ch[it] = (double)exp(-h[it] / xmus) / 2;
428 cr = cr / 8;
429 ca = ca / ha;
430 double ratio = cr / (cr + ca);
431 xdel[it] = (1 - ratio) * pizmoy;
432 ydel[it] = ratio;
433 }
434 }
435
436 /* update plane layer if necessary */
437 if (ntp == (snt - 1))
438 {
439 /* compute position of the plane layer */
440 double taup = tap + trp;
441 iplane = -1;
442 for(int i = 0; i <= ntp; i++) if (taup >= h[i]) iplane = i;
443
444 /* update the layer from the end to the position to update if necessary */
445 double xt1 = (double)fabs(h[iplane] - taup);
446 double xt2 = (double)fabs(h[iplane+1] - taup);
447
448 if ((xt1 > 0.0005) && (xt2 > 0.0005))
449 {
450 for(int i = snt; i >= iplane+1; i--)
451 {
452 xdel[i] = xdel[i-1];
453 ydel[i] = ydel[i-1];
454 h[i] = h[i-1];
455 altc[i] = altc[i-1];
456 ch[i] = ch[i-1];
457 }
458 }
459 else
460 {
461 snt = ntp;
462 if (xt2 > xt1) iplane++;
463
464 }
465
466 h[iplane] = taup;
467 if ( trmoy > accu2 && tamoy > accu2 )
468 {
469
470 double ca = tamoy * (double)exp(-alt.palt / ha);
471 double cr = trmoy * (double)exp(-alt.palt / 8);
472 h[iplane] = ca + cr;
473 cr = cr / 8;
474 ca = ca / ha;
475 double ratio = cr / (cr + ca);
476 xdel[iplane] = (1 - ratio) * pizmoy;
477 ydel[iplane] = ratio;
478 altc[iplane] = alt.palt;
479 ch[iplane] = (double)exp(-h[iplane] / xmus) / 2;
480 }
481
482 if ( trmoy > accu2 && tamoy <= accu2 )
483 {
484 ydel[iplane] = 1;
485 xdel[iplane] = 0;
486 altc[iplane] = alt.palt;
487 }
488
489 if ( trmoy <= accu2 && tamoy > accu2 )
490 {
491 ydel[iplane] = 0;
492 xdel[iplane] = pizmoy;
493 altc[iplane] = alt.palt;
494 }
495 }
496
497
498 double phi = (double)geom.phirad;
499 int i, k;
500 for(i = 0; i < np; i++) for(int m = -mu; m <= mu; m++) xl[STDI(m)][i] = 0;
501
502 /* ************ incident angle mus ******* */
503
504 double aaaa = delta / (2 - delta);
505 double ron = (1 - aaaa) / (1 + 2 * aaaa);
506
507 /* rayleigh phase function */
508
509 double beta0 = 1;
510 double beta2 = 0.5f * ron;
511
512 /* fourier decomposition */
513 double i1[31][2*mu + 1];
514 double i2[31][2*mu + 1];
515 double i3[2*mu + 1];
516 double i4[2*mu + 1];
517 double in[2*mu + 1];
518 double inm1[2*mu + 1];
519 double inm2[2*mu + 1];
520 for(i = -mu; i <= mu; i++) i4[STDI(i)] = 0;
521
522 int iborm = 80;
523 if(fabs(xmus - 1.000000) < 1e-06) iborm = 0;
524
525 for(int is = 0; is <= iborm; is++)
526 {
527 /* primary scattering */
528 int ig = 1;
529 double roavion0 = 0;
530 double roavion1 = 0;
531 double roavion2 = 0;
532 double roavion = 0;
533
534 int j;
535 for(j = -mu; j <= mu; j++) i3[STDI(j)] = 0;
536
537 /* kernel computations */
538 double xpl[2*mu + 1];
539 double bp[26][2*mu + 1];
540 memset(xpl, 0, sizeof(double)*(2*mu+1));
541 memset(bp, 0, sizeof(double)*26*(2*mu+1));
542
543 kernel(is,xpl,bp,gauss);
544
545 if(is > 0) beta0 = 0;
546
547 for(j = -mu; j <= mu; j++)
548 {
549 double sa1;
550 double sa2;
551
552 if((is - 2) <= 0)
553 {
554 double spl = xpl[STDI(0)];
555 sa1 = beta0 + beta2 * xpl[STDI(j)] * spl;
556 sa2 = bp[0][STDI(j)];
557 }
558 else
559 {
560 sa2 = bp[0][STDI(j)];
561 sa1 = 0;
562 }
563 /* primary scattering source function at every level within the layer */
564
565 for(k = 0; k <= snt; k++)
566 {
567 double c = ch[k];
568 double a = ydel[k];
569 double b = xdel[k];
570 i2[k][STDI(j)] = (double)(c * (sa2 * b + sa1 * a));
571 }
572 }
573
574 /* vertical integration, primary upward radiation */
575 for(k = 1; k <= mu; k++)
576 {
577 i1[snt][STDI(k)] = 0;
578 double zi1 = i1[snt][STDI(k)];
579
580 for(i = snt - 1; i >= 0; i--)
581 {
582 double f = h[i + 1] - h[i];
583 double a = (i2[i + 1][STDI(k)] - i2[i][STDI(k)]) / f;
584 double b = i2[i][STDI(k)] - a * h[i];
585 double c = (double)exp(-f / gauss.rm[STDI(k)]);
586
587 double xx = h[i] - h[i + 1] * c;
588 zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
589 i1[i][STDI(k)] = zi1;
590 }
591 }
592
593 /* vertical integration, primary downward radiation */
594 for(k = -mu; k <= -1; k++)
595 {
596 i1[0][STDI(k)] = 0;
597 double zi1 = i1[0][STDI(k)];
598
599 for(i = 1; i <= snt; i++)
600 {
601 double f = h[i] - h[i - 1];
602 double c = (double)exp(f / gauss.rm[STDI(k)]);
603 double a = (i2[i][STDI(k)] -i2[i - 1][STDI(k)]) / f;
604 double b = i2[i][STDI(k)] - a * h[i];
605 double xx = h[i] - h[i - 1] * c;
606 zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx)/ 2);
607 i1[i][STDI(k)] = zi1;
608 }
609 }
610
611 /* inm2 is inialized with scattering computed at n-2
612 i3 is inialized with primary scattering */
613 for(k = -mu; k <= mu; k++)
614 {
615 if(k < 0)
616 {
617 inm1[STDI(k)] = i1[snt][STDI(k)];
618 inm2[STDI(k)] = i1[snt][STDI(k)];
619 i3[STDI(k)] = i1[snt][STDI(k)];
620 }
621 else if(k > 0)
622 {
623 inm1[STDI(k)] = i1[0][STDI(k)];
624 inm2[STDI(k)] = i1[0][STDI(k)];
625 i3[STDI(k)] = i1[0][STDI(k)];
626 }
627 }
628 roavion2 = i1[iplane][STDI(mu)];
629 roavion = i1[iplane][STDI(mu)];
630
631 do
632 {
633 /* loop on successive order */
634 ig++;
635
636 /* successive orders
637 multiple scattering source function at every level within the laye
638 if is < ou = 2 kernels are a mixing of aerosols and molecules kern
639 if is >2 aerosols kernels only */
640
641 if(is - 2 <= 0)
642 {
643 for(k = 1; k <= mu; k++)
644 {
645 for(i = 0; i <= snt; i++)
646 {
647 double ii1 = 0;
648 double ii2 = 0;
649
650 for(j = 1; j <= mu; j++)
651 {
652 double bpjk = bp[j][STDI(k)] * xdel[i] + ydel[i] * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(k)]);
653 double bpjmk = bp[j][STDI(-k)] * xdel[i] + ydel[i] * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(-k)]);
654 double xdb = gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjk + i1[i][STDI(-j)] * bpjmk);
655 ii2 += xdb;
656 xdb = gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjmk + i1[i][STDI(-j)] * bpjk);
657 ii1 += xdb;
658 }
659
660 if (ii2 < 1e-30) ii2 = 0;
661 if (ii1 < 1e-30) ii1 = 0;
662 i2[i][STDI(k)] = (double)ii2;
663 i2[i][STDI(-k)]= (double)ii1;
664 }
665 }
666 }
667 else
668 {
669 for(k = 1; k <= mu; k++)
670 {
671 double ii1;
672 double ii2;
673
674 for(i = 0; i <= snt; i++)
675 {
676 ii1 = 0;
677 ii2 = 0;
678
679 for(j = 1; j <= mu; j++)
680 {
681 double bpjk = bp[j][STDI(k)] * xdel[i];
682 double bpjmk = bp[j][STDI(-k)] * xdel[i];
683 double xdb = gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjk + i1[i][STDI(-j)] * bpjmk);
684 ii2 += xdb;
685 xdb = gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjmk + i1[i][STDI(-j)] * bpjk);
686 ii1 += xdb;
687 }
688
689 if (ii2 < 1e-30) ii2 = 0;
690 if (ii1 < 1e-30) ii1 = 0;
691 i2[i][STDI(k)] = (double)ii2;
692 i2[i][STDI(-k)] = (double)ii1;
693 }
694 }
695 }
696
697
698 /* vertical integration, upward radiation */
699 for(k = 1; k <= mu; k++)
700 {
701 i1[snt][STDI(k)] = 0;
702 double zi1 = i1[snt][STDI(k)];
703
704 for(i = snt-1; i >= 0; i--)
705 {
706 double f = h[i + 1] - h[i];
707 double a = (i2[i + 1][STDI(k)] - i2[i][STDI(k)]) / f;
708 double b = i2[i][STDI(k)] - a * h[i];
709 double c = (double)exp(-f / gauss.rm[STDI(k)]);
710 double xx = h[i] - h[i + 1] * c;
711 zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
712 if (fabs(zi1) <= 1e-20) zi1 = 0;
713 i1[i][STDI(k)] = zi1;
714 }
715 }
716
717 /* vertical integration, downward radiation */
718 for(k = -mu; k <= -1; k++)
719 {
720 i1[0][STDI(k)] = 0;
721 double zi1 = i1[0][STDI(k)];
722
723 for(i = 1; i <= snt; i++)
724 {
725 double f = h[i] - h[i - 1];
726 double c = (double)exp(f / gauss.rm[STDI(k)]);
727 double a = (i2[i][STDI(k)] - i2[i - 1][STDI(k)]) / f;
728 double b = i2[i][STDI(k)] - a * h[i];
729 double xx = h[i] - h[i - 1] * c;
730 zi1 = (double)(c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2);
731
732 if (fabs(zi1) <= 1e-20) zi1 = 0;
733 i1[i][STDI(k)] = zi1;
734 }
735 }
736
737 /* in is the nieme scattering order */
738 for(k = -mu; k <= mu; k++)
739 {
740 if(k < 0) in[STDI(k)] = i1[snt][STDI(k)];
741 else if(k > 0) in[STDI(k)] = i1[0][STDI(k)];
742 }
743 roavion0 = i1[iplane][STDI(mu)];
744
745 /* convergence test (geometrical serie) */
746 if(ig > 2)
747 {
748 double a1 = roavion2;
749 double d1 = roavion1;
750
751 double g1 = roavion0;
752
753
754 double z = 0;
755 if(a1 >= accu && d1 >= accu && roavion >= accu)
756 {
757 double y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / roavion));
758 y = fabs(y);
759 z = y > z ? y : z;
760 }
761
762 for(int l = -mu; l <= mu; l++)
763 {
764 if (l == 0) continue;
765 a1 = inm2[STDI(l)];
766 d1 = inm1[STDI(l)];
767 g1 = in[STDI(l)];
768 if(a1 <= accu) continue;
769 if(d1 <= accu) continue;
770 if(i3[STDI(l)] <= accu) continue;
771
772 double y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / i3[STDI(l)]));
773 y = fabs(y);
774 z = y > z ? y : z;
775 }
776
777 if(z < 0.0001)
778 {
779 /* successful test (geometrical serie) */
780 double y1;
781
782 for(int l = -mu; l <= mu; l++)
783 {
784 y1 = 1;
785 d1 = inm1[STDI(l)];
786 g1 = in[STDI(l)];
787 if(d1 <= accu) continue;
788 y1 = 1 - g1 / d1;
789 if(fabs(g1 - d1) <= accu) continue;
790 g1 /= y1;
791 i3[STDI(l)] += g1;
792 }
793
794
795 d1 = roavion1;
796 g1 = roavion0;
797 y1 = 1;
798 if(d1 >= accu)
799 {
800 if(fabs(g1 - d1) >= accu)
801 {
802 y1 = 1 - g1 / d1;
803 g1 /= y1;
804 }
805
806 roavion += g1;
807 }
808
809 break; /* break out of the while loop */
810 }
811
812 /* inm2 is the (n-2)ieme scattering order */
813 for(k = -mu; k <= mu; k++) inm2[STDI(k)] = inm1[STDI(k)];
814 roavion2 = roavion1;
815 }
816
817 /* inm1 is the (n-1)ieme scattering order */
818 for(k = -mu; k <= mu; k++) inm1[STDI(k)] = in[STDI(k)];
819 roavion1 = roavion0;
820
821 /* sum of the n-1 orders */
822 for(k = -mu; k <= mu; k++) i3[STDI(k)] += in[STDI(k)];
823 roavion += roavion0;
824
825 /* stop if order n is less than 1% of the sum */
826 double z = 0;
827 for(k = -mu; k <= mu; k++)
828 {
829 if (fabs(i3[STDI(k)]) >= accu)
830 {
831 double y = fabs(in[STDI(k)] / i3[STDI(k)]);
832 z = z >= y ? z : y;
833 }
834 }
835 if(z < 0.00001) break;
836
837 } while( ig <= 20 ); /* stop if order n is greater than 20 in any case */
838
839 /* sum of the fourier component s */
840 double delta0s = 1;
841 if(is != 0) delta0s = 2;
842 for(k = -mu; k <= mu; k++) i4[STDI(k)] += delta0s * i3[STDI(k)];
843
844 /* stop of the fourier decomposition */
845 int l;
846 for(l = 0; l < np; l++)
847 {
848 phi = gauss.rp[l];
849
850 for(int m = -mum1; m <= mum1; m++)
851 {
852 if(m > 0) xl[STDI(m)][l] += (double)(delta0s * i3[STDI(m)] * cos(is * (phi + M_PI)));
853 else xl[STDI(m)][l] += (double)(delta0s * i3[STDI(m)] * cos(is * phi));
854 }
855 }
856
857 if(is == 0)
858 for(k = 1; k <= mum1; k++) xl[STDI(0)][0] += gauss.rm[STDI(k)] * gauss.gb[STDI(k)] * i3[STDI(-k)];
859
860 xl[STDI(mu)][0] += (double)(delta0s * i3[STDI(mu)] * cos(is * (geom.phirad + M_PI)));
861 xl[STDI(-mu)][0] += (double)(delta0s * roavion * cos(is * (geom.phirad + M_PI)));
862
863 double z = 0;
864 for(l = -mu; l <= mu; l++)
865 {
866 if(l == 0) continue;
867 if (fabs(i4[STDI(l)]) > accu) continue;
868 double x = fabs(i3[STDI(l)] / i4[STDI(l)]);
869 z = z > x ? z : x;
870 }
871
872 if(z <= 0.001) break;
873 }
874}
875
876
877/*
878 Compute the atmospheric transmission for either a satellite or aircraft observation
879 as well as the spherical albedo of the atmosphere.
880*/
881void iso(const double tamoy, const double trmoy, const double pizmoy,
882 const double tamoyp, const double trmoyp, double (&xf)[3],
883 Gauss &gauss, const Altitude &alt)
884{
885 /* molecular ratio within the layer
886 computations are performed assuming a scale of 8km for
887 molecules and 2km for aerosols */
888
889 /* the optical thickness above plane are recomputed to give o.t above pla */
890 double trp = trmoy - trmoyp;
891 double tap = tamoy - tamoyp;
892
893 /* if plane observations recompute scale height for aerosol knowing:
894 the aerosol optical depth as measure from the plane = tamoyp
895 the rayleigh scale height = = hr (8km)
896 the rayleigh optical depth at plane level = trmoyp
897 the altitude of the plane = palt
898 the rayleigh optical depth for total atmos = trmoy
899 the aerosol optical depth for total atmos = tamoy
900 if not plane observations then ha is equal to 2.0km
901 sntp local variable: if sntp=snt no plane observation selected
902 sntp=snt-1 plane observation selected */
903
904 /* it's a mixing rayleigh+aerosol */
905 int snt = nt;
906 int iplane = 0;
907 int ntp = snt;
908 double ha = 2.0;
909 if(alt.palt <= 900. && alt.palt > 0.0)
910 {
911 if (tap > 1.e-03) ha = (double)(-alt.palt / log(tap / tamoy));
912 else ha = 2.;
913 ntp = snt - 1;
914 }
915
916 /* compute mixing rayleigh, aerosol
917 case 1: pure rayleigh
918 case 2: pure aerosol
919 case 3: mixing rayleigh-aerosol */
920
921 double h[31];
922 double ydel[31];
923 double xdel[31];
924 double altc[31];
925
926 if((tamoy <= accu2) && (trmoy > tamoy))
927 {
928 for(int j = 0; j <= ntp; j++)
929 {
930 h[j] = j * trmoy / ntp;
931 ydel[j] = 1.0;
932 xdel[j] = 0.0;
933 }
934 }
935
936 if((trmoy <= accu2) && (tamoy > trmoy))
937 {
938 for(int j = 0; j <= ntp; j++)
939 {
940 h[j] = j * tamoy / ntp;
941 ydel[j] = 0.0;
942 xdel[j] = pizmoy;
943 }
944 }
945
946 if(trmoy > accu2 && tamoy > accu2)
947 {
948 ydel[0] = 1.0;
949 xdel[0] = 0.0;
950 h[0] = 0;
951 altc[0] = 300;
952 double zx = 300;
953 iplane = 0;
954
955 for(int it = 0; it <= ntp; it++)
956
957 {
958 if (it == 0) zx = discre(tamoy,ha,trmoy,8.0,it,ntp,0,0,300.0,0.0);
959 else zx = discre(tamoy,ha,trmoy,8.0,it,ntp,h[it-1],ydel[it-1],300.0,0.0);
960
961 double ca;
962 if ((-zx / ha) < -18) ca = 0;
963 else ca = (double)(tamoy * exp(-zx / ha));
964
965 double cr = (double)(trmoy * exp(-zx / 8.0));
966 h[it] = cr + ca;
967 altc[it] = zx;
968
969 cr = cr / 8;
970 ca = ca / ha;
971 double ratio = cr / (cr + ca);
972 xdel[it] = (1 - ratio) * pizmoy;
973 ydel[it] = ratio;
974
975 }
976 }
977
978 /* update plane layer if necessary */
979 if (ntp == (snt-1))
980 {
981 /* compute position of the plane layer */
982 double taup = tap + trp;
983 iplane = -1;
984 for(int i = 0; i <= ntp; i++) if (taup >= h[i]) iplane = i;
985
986 /* update the layer from the end to the position to update if necessary */
987 double xt1 = (double)fabs(h[iplane] - taup);
988 double xt2 = (double)fabs(h[iplane + 1] - taup);
989 if ((xt1 > 0.005) && (xt2 > 0.005))
990 {
991 for(int i = snt; i >= iplane + 1; i--)
992 {
993 xdel[i] = xdel[i-1];
994
995
996 ydel[i] = ydel[i-1];
997 h[i] = h[i-1];
998 altc[i] = altc[i-1];
999 }
1000 }
1001 else
1002 {
1003 snt = ntp;
1004 if (xt2 < xt1) iplane = iplane + 1;
1005 }
1006
1007 h[iplane] = taup;
1008 if ( trmoy > accu2 && tamoy > accu2)
1009 {
1010 double ca = (double)(tamoy * exp(-alt.palt / ha));
1011 double cr = (double)(trmoy * exp(-alt.palt / 8.0));
1012 cr = cr / 8;
1013 ca = ca / ha;
1014 double ratio = cr / (cr + ca);
1015 xdel[iplane] = (1 - ratio) * pizmoy;
1016
1017 ydel[iplane] = ratio;
1018 altc[iplane] = alt.palt;
1019 }
1020
1021 if ( trmoy > accu2 && tamoy <= accu2)
1022 {
1023 ydel[iplane] = 1;
1024 xdel[iplane] = 0;
1025 altc[iplane] = alt.palt;
1026 }
1027
1028 if ( trmoy <= accu2 && tamoy > accu2)
1029 {
1030 ydel[iplane] = 0;
1031 xdel[iplane] = 1 * pizmoy;
1032 altc[iplane] = alt.palt;
1033 }
1034 }
1035
1036 double aaaa = delta / (2-delta);
1037 double ron = (1 - aaaa) / (1 + 2 * aaaa);
1038
1039 /* rayleigh phase function */
1040 double beta0 = 1;
1041 double beta2 = 0.5f * ron;
1042
1043 /* primary scattering */
1044 int ig = 1;
1045 double tavion0 = 0;
1046 double tavion1 = 0;
1047 double tavion2 = 0;
1048 double tavion = 0;
1049
1050 double i1[31][2*mu + 1];
1051 double i2[31][2*mu + 1];
1052 double i3[2*mu + 1];
1053 double in[2*mu + 1];
1054 double inm1[2*mu + 1];
1055 double inm2[2*mu + 1];
1056 int j;
1057 for(j = -mu; j <= mu; j++) i3[STDI(j)] = 0;
1058
1059 /* kernel computations */
1060 double xpl[2*mu + 1];
1061 double bp[26][2*mu + 1];
1062 kernel(0, xpl, bp, gauss);
1063
1064 for(j = -mu; j <= mu; j++)
1065 for(int k = 0; k <= snt; k++) i2[k][STDI(j)] = 0;
1066
1067 /* vertical integration, primary upward radiation */
1068 int k;
1069 for(k = 1; k <= mu; k++)
1070 {
1071 i1[snt][STDI(k)] = 1.0;
1072 for(int i = snt-1; i >= 0; i--)
1073 i1[i][STDI(k)] = (double)(exp(-(tamoy + trmoy - h[i]) / gauss.rm[STDI(k)]));
1074 }
1075
1076
1077 /* vertical integration, primary downward radiation */
1078 for(k = -mu; k <= -1; k++)
1079 for(int i = 0; i <= snt; i++) i1[i][STDI(k)] = 0.0;
1080
1081 /* inm2 is inialized with scattering computed at n-2
1082 i3 is inialized with primary scattering */
1083 for(k = -mu; k <= mu; k++)
1084 {
1085 if(k == 0) continue;
1086 if(k < 0)
1087 {
1088 inm1[STDI(k)] = i1[snt][STDI(k)];
1089 inm2[STDI(k)] = i1[snt][STDI(k)];
1090 i3[STDI(k)] = i1[snt][STDI(k)];
1091 }
1092 else
1093 {
1094 inm1[STDI(k)] = i1[0][STDI(k)];
1095 inm2[STDI(k)] = i1[0][STDI(k)];
1096 i3[STDI(k)] = i1[0][STDI(k)];
1097 }
1098 }
1099 tavion = i1[iplane][STDI(mu)];
1100 tavion2 = i1[iplane][STDI(mu)];
1101
1102 do {
1103 /* loop on successive order */
1104 ig++;
1105
1106 /* successive orders
1107 multiple scattering source function at every level within the layer */
1108 for(k = 1; k <= mu; k++)
1109 {
1110 for(int i = 0; i <= snt; i++)
1111 {
1112 double ii1 = 0;
1113 double ii2 = 0;
1114 double x = xdel[i];
1115 double y = ydel[i];
1116
1117 for(j = 1; j <= mu; j++)
1118 {
1119 double bpjk = bp[j][STDI(k)] * x + y * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(k)]);
1120 double bpjmk= bp[j][STDI(-k)] * x + y * (beta0 + beta2 * xpl[STDI(j)] * xpl[STDI(-k)]);
1121 ii2 += gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjk + i1[i][STDI(-j)] * bpjmk);
1122 ii1 += gauss.gb[STDI(j)] * (i1[i][STDI(j)] * bpjmk + i1[i][STDI(-j)] * bpjk);
1123 }
1124
1125 i2[i][STDI(k)] = (double)ii2;
1126 i2[i][STDI(-k)] = (double)ii1;
1127 }
1128 }
1129
1130 /* vertical integration, upward radiation */
1131 for(k = 1; k <= mu; k++)
1132 {
1133 i1[snt][STDI(k)] = 0.0;
1134 double zi1 = i1[snt][STDI(k)];
1135
1136 for(int i = snt-1; i >= 0; i--)
1137 {
1138 double f = h[i+1] - h[i];
1139 double a = (i2[i+1][STDI(k)] -i2[i][STDI(k)]) / f;
1140 double b = i2[i][STDI(k)] - a * h[i];
1141 double c = (double)exp(-f / gauss.rm[STDI(k)]);
1142 double xx = h[i] - h[i+1] * c;
1143
1144 zi1 = c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2;
1145 i1[i][STDI(k)] = zi1;
1146 }
1147 }
1148
1149 /* vertical integration, downward radiation */
1150 for(k = -mu; k <= -1; k++)
1151 {
1152 i1[0][STDI(k)] = 0;
1153 double zi1 = i1[0][STDI(k)];
1154
1155 for(int i = 1; i <= snt; i++)
1156 {
1157 double f = h[i] - h[i-1];
1158 double c = (double)exp(f / gauss.rm[STDI(k)]);
1159 double a = (i2[i][STDI(k)] - i2[i-1][STDI(k)]) / f;
1160 double b = i2[i][STDI(k)] - a * h[i];
1161 double xx = h[i] - h[i-1] * c;
1162 zi1 = c * zi1 + ((1 - c) * (b + a * gauss.rm[STDI(k)]) + a * xx) / 2;
1163 i1[i][STDI(k)] = zi1;
1164 }
1165 }
1166
1167 /* in is the nieme scattering order */
1168 for(k = -mu; k <= mu; k++)
1169 {
1170 if(k == 0) continue;
1171 if(k < 0) in[STDI(k)] = i1[snt][STDI(k)];
1172 else in[STDI(k)] = i1[0][STDI(k)];
1173 }
1174 tavion0 = i1[iplane][STDI(mu)];
1175
1176 /* convergence test (geometrical serie) */
1177 if(ig > 2)
1178 {
1179 double z = 0;
1180 double a1 = tavion2;
1181 double d1 = tavion1;
1182 double g1 = tavion0;
1183 if (a1 >= accu && d1 >= accu && tavion >= accu)
1184 {
1185 double y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / tavion));
1186 y = (double)fabs(y);
1187 z = z >= y ? z : y;
1188 }
1189
1190 for(int l = -mu; l <= mu; l++)
1191 {
1192 if (l == 0) continue;
1193 a1 = inm2[STDI(l)];
1194 d1 = inm1[STDI(l)];
1195 g1 = in[STDI(l)];
1196 if(a1 == 0) continue;
1197 if(d1 == 0) continue;
1198 if(i3[STDI(l)] == 0) continue;
1199
1200 double y = ((g1 / d1 - d1 / a1) / ((1 - g1 / d1) * (1 - g1 / d1)) * (g1 / i3[STDI(l)]));
1201 y = (double)fabs(y);
1202 z = z >= y ? z : y;
1203 }
1204
1205 if(z < 0.0001)
1206 {
1207 /* successful test (geometrical serie) */
1208
1209 for(int l = -mu; l <= mu; l++)
1210 {
1211 if (l == 0) continue;
1212 double y1 = 1;
1213 d1 = inm1[STDI(l)];
1214 g1 = in[STDI(l)];
1215 if(d1 == 0) continue;
1216 y1 = 1 - g1 / d1;
1217 g1 = g1 / y1;
1218 i3[STDI(l)] += g1;
1219 }
1220
1221 d1 = tavion1;
1222 g1 = tavion0;
1223 if (d1 >= accu)
1224 {
1225 if (fabs(g1 - d1) >= accu)
1226 {
1227 double y1 = 1 - g1 / d1;
1228 g1 = g1 / y1;
1229 }
1230 tavion = tavion + g1;
1231
1232 }
1233
1234 break; /* go to 505 */
1235 }
1236
1237 /* inm2 is the (n-2)ieme scattering order */
1238 for(k = -mu; k <= mu; k++) inm2[STDI(k)] = inm1[STDI(k)];
1239 tavion2 = tavion1;
1240 }
1241
1242 /* inm1 is the (n-1)ieme scattering order */
1243 for(k = -mu; k <= mu; k++) inm1[STDI(k)] = in[STDI(k)];
1244 tavion1 = tavion0;
1245
1246 /* sum of the n-1 orders */
1247 for(k = -mu; k <= mu; k++) i3[STDI(k)] += in[STDI(k)];
1248 tavion = tavion + tavion0;
1249
1250 /* stop if order n is less than 1% of the sum */
1251 double z = 0;
1252 for(k = -mu; k <= mu; k++)
1253 {
1254 if(i3[STDI(k)] != 0)
1255 {
1256 double y = (double)fabs(in[STDI(k)] / i3[STDI(k)]);
1257 z = z >= y ? z : y;
1258 }
1259 }
1260 if(z < 0.00001) break;
1261
1262 /* stop if order n is greater than 20 in any case */
1263 } while(ig <= 20);
1264
1265 /* dimension for os computation */
1266 xf[0] = tavion;
1267 xf[1] = 0;
1268 xf[2] = 0;
1269
1270 xf[2] += i3[STDI(mu)];
1271 for(k = 1; k <= mu; k++) xf[1] += gauss.rm[STDI(k)] * gauss.gb[STDI(k)] * i3[STDI(-k)];
1272
1273}
1274
1275/*
1276 To compute the atmospheric reflectance for the molecular atmosphere in
1277 case of satellite observation.
1278*/
1279double chand(const double xtau, const GeomCond &geom)
1280{
1281 /* input parameters: xphi,xmus,xmuv,xtau
1282 xphi: azimuthal difference between sun and observation (xphi=0,
1283 in backscattering) and expressed in degree (0.:360.)
1284 xmus: cosine of the sun zenith angle
1285 xmuv: cosine of the observation zenith angle
1286 xtau: molecular optical depth
1287 output parameter: xrray : molecular reflectance (0.:1.)
1288 constant : xdep: depolarization factor (0.0279) */
1289
1290 const double xdep = 0.0279;
1291
1292 static const double as0[10] = {
1293 .33243832,-6.777104e-02,.16285370,1.577425e-03,-.30924818,
1294 -1.240906e-02,-.10324388,3.241678e-02,.11493334,-3.503695e-02
1295 };
1296
1297 static const double as1[2] = { .19666292, -5.439061e-02 };
1298 static const double as2[2] = { .14545937,-2.910845e-02 };
1299
1300 double phios = (180 - geom.phi);
1301 double xcosf1 = 1;
1302 double xcosf2 = cos(phios * M_PI / 180);
1303 double xcosf3 = cos(2 * phios * M_PI / 180);
1304
1305
1306 double xfd = xdep / (2 - xdep);
1307 xfd = (1 - xfd) / (1 + 2 * xfd);
1308
1309 double xph1 = 1 + (3 * geom.xmus * geom.xmus - 1) * (3 * geom.xmuv * geom.xmuv - 1) * xfd / 8;
1310 double xph2 = -geom.xmus * geom.xmuv * sqrt(1 - geom.xmus * geom.xmus) * sqrt(1 - geom.xmuv * geom.xmuv);
1311 xph2 *= xfd * 0.75;
1312 double xph3 = (1 - geom.xmus * geom.xmus) * (1 - geom.xmuv * geom.xmuv);
1313 xph3 *= xfd * 0.1875;
1314
1315 double xitm = (1 - exp(-xtau * (1 / geom.xmus + 1 / geom.xmuv))) * geom.xmus / (4 * (geom.xmus + geom.xmuv));
1316 double xp1 = xph1 * xitm;
1317 double xp2 = xph2 * xitm;
1318 double xp3 = xph3 * xitm;
1319
1320 xitm = (1 - exp(-xtau / geom.xmus)) * (1 - exp(-xtau / geom.xmuv));
1321 double cfonc1 = xph1 * xitm;
1322 double cfonc2 = xph2 * xitm;
1323 double cfonc3 = xph3 * xitm;
1324 double xlntau = log(xtau);
1325
1326 double pl[10];
1327 pl[0] = 1;
1328 pl[1] = xlntau;
1329 pl[2] = geom.xmus + geom.xmuv;
1330 pl[3] = xlntau * pl[2];
1331 pl[4] = geom.xmus * geom.xmuv;
1332 pl[5] = xlntau * pl[4];
1333 pl[6] = geom.xmus * geom.xmus + geom.xmuv * geom.xmuv;
1334 pl[7] = xlntau * pl[6];
1335 pl[8] = geom.xmus * geom.xmus * geom.xmuv * geom.xmuv;
1336 pl[9] = xlntau * pl[8];
1337
1338 double fs0 = 0;
1339 for(int i = 0; i < 10; i++) fs0 += pl[i] * as0[i];
1340
1341 double fs1 = pl[0] * as1[0] + pl[1] * as1[1];
1342 double fs2 = pl[0] * as2[0] + pl[1] * as2[1];
1343 double xitot1 = xp1 + cfonc1 * fs0 * geom.xmus;
1344 double xitot2 = xp2 + cfonc2 * fs1 * geom.xmus;
1345 double xitot3 = xp3 + cfonc3 * fs2 * geom.xmus;
1346
1347 double xrray = (double)(xitot1 * xcosf1);
1348 xrray += (double)(xitot2 * xcosf2 * 2);
1349
1350 xrray += (double)(xitot3 * xcosf3 * 2);
1351 xrray /= (double)geom.xmus;
1352
1353 return xrray;
1354}
1355
1356/*
1357 To compute the atmospheric reflectance for the molecular and aerosol atmospheres
1358 and the mixed atmosphere. In 6S instead of an approximation as in 5S, we use the scalar Successive
1359 Order of Scattering method (subroutine OS.f). The polarization terms of aerosol or rayleigh phase
1360 are not accounted for in the computation of the aerosol reflectance and the mixed Rayleigh-aerosol
1361 reflectance. The polarization is addressed in computing the Rayleigh reflectance (Subroutine
1362 CHAND.f) by semi-empirical fitting of the vectorized Successive Orders of Scattering method
1363 (Deuzé et al, 1989).
1364*/
1365void atmref(const double tamoy, const double trmoy, const double pizmoy,
1366 const double tamoyp, const double trmoyp, OpticalAtmosProperties &oap,
1367 Gauss &gauss, const GeomCond &geom, const AerosolModel &aero,
1368 const Altitude &alt)
1369{
1370 double xlm1[2 * mu + 1][np];
1371 double xlm2[2 * mu + 1][np];
1372
1373 /* atmospheric reflectances */
1374 oap.rorayl = 0;
1375 oap.roaero = 0;
1376
1377 /* rayleigh reflectance 3 cases (satellite,plane,ground) */
1378 if(alt.palt < 900 && alt.palt > 0)
1379 {
1380 gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
1381 gauss.rm[STDI(mu)] = (double)geom.xmuv;
1382 gauss.rm[STDI(0)] = -(double)geom.xmus;
1383
1384 os(0, trmoy, pizmoy, 0, trmoyp, xlm1, gauss, alt, geom);
1385
1386 oap.rorayl = (double)(xlm1[STDI(-mu)][0] / geom.xmus);
1387 }
1388 else if(alt.palt <= 0) oap.rorayl = 0;
1389 else oap.rorayl = chand(trmoy, geom);
1390
1391 if (aero.iaer == 0)
1392 {
1393 oap.romix = oap.rorayl;
1394 return;
1395 }
1396
1397 /* rayleigh+aerosol=romix,aerosol=roaero reflectance computed
1398 using successive order of scattering method
1399 3 cases: satellite,plane,ground */
1400 if (alt.palt > 0)
1401 {
1402 gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
1403 gauss.rm[STDI(mu)] = (double)geom.xmuv;
1404 gauss.rm[STDI(0)] = -(double)geom.xmus;
1405
1406 os(tamoy, trmoy, pizmoy, tamoyp, trmoyp, xlm2, gauss, alt, geom);
1407 oap.romix = (double)(xlm2[STDI(-mu)][0] / geom.xmus);
1408
1409 os(tamoy, 0, pizmoy, tamoyp, 0, xlm2, gauss, alt, geom);
1410 oap.roaero = (double)(xlm2[STDI(-mu)][0] / geom.xmus);
1411 }
1412 else
1413 {
1414 oap.roaero = 0;
1415 oap.romix = 0;
1416 }
1417}
1418
1419
1420double fintexp1(const double xtau)
1421{
1422 /* accuracy 2e-07... for 0<xtau<1 */
1423 double a[6] = { -.57721566,0.99999193,-0.24991055,0.05519968,-0.00976004,0.00107857 };
1424 double xftau = 1;
1425 double xx = a[0];
1426 for(int i = 1; i <= 5; i++)
1427 {
1428 xftau *= xtau;
1429 xx += a[i] * xftau;
1430 }
1431 return (double)(xx - log(xtau));
1432}
1433
1434double fintexp3(const double xtau)
1435{
1436 return (double)((exp(-xtau) * (1. - xtau) + xtau * xtau * fintexp1(xtau)) / 2.);
1437}
1438
1439/* To compute the spherical albedo of the molecular layer. */
1440void csalbr(const double xtau, double& xalb)
1441{
1442 xalb = (double)((3. * xtau - fintexp3(xtau) * (4. + 2. * xtau) + 2. * exp(-xtau)));
1443 xalb = (double)(xalb / (4. + 3. * xtau));
1444}
1445
1446void scatra(const double taer, const double taerp,
1447 const double tray, const double trayp,
1448 const double piza, OpticalAtmosProperties& oap,
1449 Gauss &gauss, const GeomCond &geom, const Altitude &alt)
1450{
1451 /* computations of the direct and diffuse transmittances
1452 for downward and upward paths , and spherical albedo */
1453 double tamol,tamolp;
1454 double xtrans[3];
1455
1456 oap.ddirtt = 1; oap.ddiftt = 0;
1457 oap.udirtt = 1; oap.udiftt = 0;
1458 oap.ddirtr = 1; oap.ddiftr = 0;
1459 oap.udirtr = 1; oap.udiftr = 0;
1460 oap.ddirta = 1; oap.ddifta = 0;
1461 oap.udirta = 1; oap.udifta = 0;
1462 oap.sphalbt = 0;
1463 oap.sphalbr = 0;
1464 oap.sphalba = 0;
1465
1466 for(int it = 1; it <= 3; it++)
1467
1468 {
1469 /* it=1 rayleigh only, it=2 aerosol only, it=3 rayleigh+aerosol */
1470 if (it == 2 && taer <= 0) continue;
1471
1472 /* compute upward,downward diffuse transmittance for rayleigh,aerosol */
1473 if (it == 1)
1474 {
1475 if (alt.palt > 900)
1476 {
1477 oap.udiftt = (double)((2. / 3. + geom.xmuv) + (2. / 3. - geom.xmuv) * exp(-tray / geom.xmuv));
1478
1479 oap.udiftt = (double)(oap.udiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmuv));
1480 oap.ddiftt = (double)((2. / 3. + geom.xmus) + (2. / 3. - geom.xmus) * exp(-tray / geom.xmus));
1481 oap.ddiftt = (double)(oap.ddiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmus));
1482 oap.ddirtt = (double)exp(-tray / geom.xmus);
1483 oap.udirtt = (double)exp(-tray / geom.xmuv);
1484
1485 csalbr(tray, oap.sphalbt);
1486 }
1487 else if (alt.palt > 0 && alt.palt <= 900)
1488 {
1489 tamol = 0;
1490 tamolp= 0;
1491 gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
1492 gauss.rm[STDI(mu)] = (double)geom.xmuv;
1493 gauss.rm[STDI(0)] = (double)geom.xmus;
1494
1495 iso(tamol, tray, piza, tamolp, trayp, xtrans, gauss, alt);
1496
1497 oap.udiftt = (double)(xtrans[0] - exp(-trayp / geom.xmuv));
1498 oap.udirtt = (double)exp(-trayp / geom.xmuv);
1499 gauss.rm[STDI(-mu)] = -(double)geom.xmus;
1500 gauss.rm[STDI(mu)] = (double)geom.xmus;
1501 gauss.rm[STDI(0)] = (double)geom.xmus;
1502
1503 oap.ddiftt = (double)((2. / 3. + geom.xmus) + (2. / 3. - geom.xmus) * exp(-tray / geom.xmus));
1504 oap.ddiftt = (double)(oap.ddiftt / ((4. / 3.) + tray) - exp(-tray / geom.xmus));
1505 oap.ddirtt = (double)exp(-tray / geom.xmus);
1506 oap.udirtt = (double)exp(-tray / geom.xmuv);
1507
1508 csalbr(tray, oap.sphalbt);
1509 }
1510 else if (alt.palt <= 0.)
1511 {
1512 oap.udiftt = 0;
1513 oap.udirtt = 1;
1514 }
1515
1516 oap.ddirtr = oap.ddirtt;
1517 oap.ddiftr = oap.ddiftt;
1518 oap.udirtr = oap.udirtt;
1519 oap.udiftr = oap.udiftt;
1520 oap.sphalbr = oap.sphalbt;
1521 }
1522
1523 else if (it == 2)
1524 {
1525 tamol = 0;
1526 tamolp= 0;
1527 gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
1528 gauss.rm[STDI(mu)] = (double)geom.xmuv;
1529 gauss.rm[STDI(0)] = (double)geom.xmus;
1530
1531 iso(taer, tamol, piza, taerp, tamolp, xtrans, gauss, alt);
1532
1533 oap.udiftt = (double)(xtrans[0] - exp(-taerp / geom.xmuv));
1534 oap.udirtt = (double)exp(-taerp / geom.xmuv);
1535 gauss.rm[STDI(-mu)] = -(double)geom.xmus;
1536 gauss.rm[STDI(mu)] = (double)geom.xmus;
1537 gauss.rm[STDI(0)] = (double)geom.xmus;
1538
1539 double tmp_alt = alt.palt;
1540 alt.palt = 999;
1541 iso(taer, tamol, piza, taerp, tamolp, xtrans, gauss, alt);
1542 alt.palt = tmp_alt;
1543
1544 oap.ddirtt = (double)exp(-taer / geom.xmus);
1545 oap.ddiftt = (double)(xtrans[2] - exp(-taer / geom.xmus));
1546 oap.sphalbt= (double)(xtrans[1] * 2);
1547
1548 if(alt.palt <= 0)
1549 {
1550 oap.udiftt = 0;
1551 oap.udirtt = 1;
1552 }
1553
1554 oap.ddirta = oap.ddirtt;
1555 oap.ddifta = oap.ddiftt;
1556 oap.udirta = oap.udirtt;
1557 oap.udifta = oap.udiftt;
1558 oap.sphalba = oap.sphalbt;
1559 }
1560 else if(it == 3)
1561 {
1562 gauss.rm[STDI(-mu)] = -(double)geom.xmuv;
1563 gauss.rm[STDI(mu)] = (double)geom.xmuv;
1564 gauss.rm[STDI(0)] = (double)geom.xmus;
1565
1566 iso(taer, tray, piza, taerp, trayp, xtrans, gauss, alt);
1567
1568 oap.udirtt = (double)exp(-(taerp + trayp) / geom.xmuv);
1569 oap.udiftt = (double)(xtrans[0] - exp(-(taerp + trayp) / geom.xmuv));
1570 gauss.rm[STDI(-mu)] = -(double)geom.xmus;
1571 gauss.rm[STDI(mu)] = (double)geom.xmus;
1572 gauss.rm[STDI(0)] = (double)geom.xmus;
1573
1574 double tmp_alt = alt.palt;
1575 alt.palt = 999;
1576 iso(taer, tray, piza, taerp, trayp, xtrans, gauss, alt);
1577 alt.palt = tmp_alt;
1578
1579 oap.ddiftt = (double)(xtrans[2] - exp(-(taer + tray) / geom.xmus));
1580 oap.ddirtt = (double)exp(-(taer + tray) / geom.xmus);
1581
1582 oap.sphalbt= (double)(xtrans[1] * 2);
1583
1584 if (alt.palt <= 0)
1585 {
1586 oap.udiftt = 0;
1587 oap.udirtt = 1;
1588 }
1589 }
1590 }
1591}
1592
1593/* To compute the optical properties of the atmosphere at the 10 discrete
1594 wavelengths. */
1595void discom(const GeomCond &geom, const AtmosModel &atms,
1596 const AerosolModel &aero, const AerosolConcentration &aerocon,
1597 const Altitude &alt, const IWave &iwave)
1598{
1599 OpticalAtmosProperties oap;
1600 memset(&oap, 0, sizeof(oap));
1601
1602 Gauss gauss;
1603
1604 gauss.init(); /* discom is the only function that uses the gauss data */
1605
1606 memset(&sixs_trunc, 0, sizeof(sixs_trunc)); /* clear this to keep preconditions the same and output consistent */
1607
1608/* computation of all scattering parameters at wavelength
1609 discrete values,so we
1610 can interpolate at any wavelength */
1611 int i;
1612 for(i = 0; i < 10; i++)
1613 {
1614 if(!((((i < 2) && iwave.ffu.wlsup < sixs_disc.wldis[0])) || ((iwave.ffu.wlinf > sixs_disc.wldis[9]) && (i >= 8))))
1615 if (((i < 9) && (sixs_disc.wldis[i] < iwave.ffu.wlinf) && (sixs_disc.wldis[i+1] < iwave.ffu.wlinf)) ||
1616 ((i > 0) && (sixs_disc.wldis[i] > iwave.ffu.wlsup) && (sixs_disc.wldis[i-1] > iwave.ffu.wlsup))) continue;
1617
1618 double wl = sixs_disc.wldis[i];
1619 /* computation of rayleigh optical depth at wl */
1620 double tray = odrayl(atms, wl);
1621 double trayp;
1622
1623 /* plane case discussed here above */
1624 if (alt.idatmp == 0) trayp = 0;
1625 else if (alt.idatmp == 4) trayp = tray;
1626 else trayp = tray * alt.ftray;
1627
1628 sixs_disc.trayl[i] = tray;
1629 sixs_disc.traypl[i] = trayp;
1630
1631 /* computation of aerosol optical properties at wl */
1632
1633 double taer = aerocon.taer55 * sixs_aer.ext[i] / sixs_aer.ext[3];
1634 double taerp = alt.taer55p * sixs_aer.ext[i] / sixs_aer.ext[3];
1635 double piza = sixs_aer.ome[i];
1636
1637 /* computation of atmospheric reflectances
1638
1639 rorayl is rayleigh ref
1640 roaero is aerosol ref
1641 call plegen to decompose aerosol phase function in Betal */
1642
1643 double coeff = 0;
1644 if(aero.iaer != 0)
1645 {
1646 for(int k = 0; k < 83; k++) sixs_trunc.pha[k] = sixs_sos.phasel[i][k];
1647 coeff = trunca();
1648 }
1649
1650 double tamoy = taer * (1 - piza * coeff);
1651 double tamoyp = taerp * (1 - piza * coeff);
1652 double pizmoy = piza * (1 - coeff) / (1 - piza * coeff);
1653
1654 atmref(tamoy, tray, pizmoy, tamoyp, trayp, oap, gauss, geom, aero, alt);
1655
1656 /* computation of scattering transmitances (direct and diffuse)
1657 first time for rayleigh ,next total (rayleigh+aerosols) */
1658
1659 scatra(tamoy, tamoyp, tray, trayp, pizmoy, oap, gauss, geom, alt);
1660
1661 sixs_disc.roatm[0][i] = oap.rorayl;
1662 sixs_disc.roatm[1][i] = oap.romix;
1663 sixs_disc.roatm[2][i] = oap.roaero;
1664 sixs_disc.dtdir[0][i] = oap.ddirtr;
1665 sixs_disc.dtdif[0][i] = oap.ddiftr;
1666 sixs_disc.dtdir[1][i] = oap.ddirtt;
1667 sixs_disc.dtdif[1][i] = oap.ddiftt;
1668 sixs_disc.dtdir[2][i] = oap.ddirta;
1669 sixs_disc.dtdif[2][i] = oap.ddifta;
1670 sixs_disc.utdir[0][i] = oap.udirtr;
1671 sixs_disc.utdif[0][i] = oap.udiftr;
1672 sixs_disc.utdir[1][i] = oap.udirtt;
1673 sixs_disc.utdif[1][i] = oap.udiftt;
1674 sixs_disc.utdir[2][i] = oap.udirta;
1675 sixs_disc.utdif[2][i] = oap.udifta;
1676 sixs_disc.sphal[0][i] = oap.sphalbr;
1677 sixs_disc.sphal[1][i] = oap.sphalbt;
1678 sixs_disc.sphal[2][i] = oap.sphalba;
1679 }
1680
1681}
1682
1683/*
1684 To compute the atmospheric properties at the equivalent wavelength (see
1685 EQUIVWL.f) needed for the calculation of the downward radiation field used
1686 in the computation of the non lambertian target contribution (main.f).
1687*/
1688void specinterp(const double wl, double& tamoy, double& tamoyp, double& pizmoy, double& pizmoyp,
1689 const AerosolConcentration &aerocon, const Altitude &alt)
1690{
1691
1692 int linf = 0;
1693 int i = 8;
1694 while (i >= 0 && linf == 0) {
1695 if (wl >= sixs_disc.wldis[i] && wl <= sixs_disc.wldis[i+1]) linf = i;
1696 i--;
1697 }
1698 if(wl > sixs_disc.wldis[9]) linf = 8;
1699
1700 int lsup = linf + 1;
1701 double coef = (double)log(sixs_disc.wldis[lsup] / sixs_disc.wldis[linf]);
1702 double wlinf = sixs_disc.wldis[linf];
1703
1704 double alphaa = (double)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] /
1705 (sixs_aer.ext[linf] * sixs_aer.ome[linf])) / coef);
1706 double betaa = (double)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
1707 double tsca = (double)(aerocon.taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
1708 alphaa = (double)(log(sixs_aer.ext[lsup] / sixs_aer.ext[linf]) / coef);
1709 betaa = (double)(sixs_aer.ext[linf] / pow(wlinf,alphaa));
1710 tamoy = (double)(aerocon.taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
1711 tamoyp= (double)(alt.taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
1712 pizmoy= tsca / tamoy;
1713 pizmoyp = pizmoy;
1714
1715 for(int k = 0; k < 83; k++)
1716 {
1717 alphaa = (double)log(sixs_sos.phasel[lsup][k] / sixs_sos.phasel[linf][k]) / coef;
1718 betaa = (double)(sixs_sos.phasel[linf][k] / pow(wlinf,alphaa));
1719 sixs_trunc.pha[k] = (double)(betaa * pow(wl,alphaa));
1720 }
1721
1722 double coeff = trunca();
1723
1724 tamoy *= 1 - pizmoy * coeff;
1725 tamoyp *= 1 - pizmoyp * coeff;
1726 pizmoy *= (1 - coeff) / (1 - pizmoy * coeff);
1727}
1728
1729/**********************************************************************c
1730c c
1731c start of computations c
1732c c
1733c**********************************************************************/
1734
1735/*
1736 To compute the environment functions F(r) which allows us to account for an
1737 inhomogeneous ground.
1738*/
1739void enviro (const double difr, const double difa, const double r, const double palt,
1740 const double xmuv, double& fra, double& fae, double& fr)
1741{
1742 double fae0, fra0, xlnv;
1743 static const double alt[16] = {
1744 0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,
1745 10.0,12.0,14.0,16.0,18.0,20.0,60.0
1746 };
1747
1748 static const double cfr1[16] = {
1749 0.730,0.710,0.656,0.606,0.560,0.516,0.473,
1750 0.433,0.395,0.323,0.258,0.209,0.171,0.142,0.122,0.070
1751 };
1752
1753 static const double cfr2[16] = {
1754 2.8,1.51,0.845,0.634,0.524,0.465,0.429,
1755 0.405,0.390,0.386,0.409,0.445,0.488,0.545,0.608,0.868
1756 };
1757
1758 static const double cfa1[16] = {
1759 0.239,0.396,0.588,0.626,0.612,0.505,0.454,
1760 0.448,0.444,0.445,0.444,0.448,0.448,0.448,0.448,0.448
1761 };
1762
1763 static const double cfa2[16] = {
1764 1.40,1.20,1.02,0.86,0.74,0.56,0.46,0.42,
1765 0.38,0.34,0.3,0.28,0.27,0.27,0.27,0.27
1766 };
1767
1768 static const double cfa3[16] = {
1769 9.17,6.26,5.48,5.16,4.74,3.65,3.24,3.15,
1770 3.07,2.97,2.88,2.83,2.83,2.83,2.83,2.83
1771 };
1772
1773
1774/* calculation of the environmental function for
1775 rayleigh and aerosols contribution.
1776
1777 this calculation have been done for nadir observation
1778 and are corrected of the effect of the view zenith angle. */
1779
1780 const double a0 = 1.3347;
1781 const double b0 = 0.57757;
1782 const double a1 = -1.479;
1783 const double b1 = -1.5275;
1784
1785 if (palt >= 60)
1786 {
1787 fae0 = (double)(1. - 0.448 * exp( -r * 0.27) - 0.552 * exp( -r * 2.83));
1788 fra0 = (double)(1. - 0.930 * exp( -r * 0.080) - 0.070 * exp( -r * 1.100));
1789 }
1790 else
1791 {
1792 int i;
1793 for(i = 0; palt >= alt[i]; i++);
1794 double xcfr1 = 0, xcfr2 = 0, xcfa1 = 0, xcfa2 = 0, xcfa3 = 0;
1795
1796 if ((i > 0) && (i < 16))
1797 {
1798 double zmin = alt[i - 1];
1799 double zmax = alt[i];
1800 xcfr1 = cfr1[i - 1] + (cfr1[i] - cfr1[i - 1]) * (palt - zmin) / (zmax - zmin);
1801 xcfr2 = cfr2[i - 1] + (cfr2[i] - cfr2[i - 1]) * (palt - zmin) / (zmax - zmin);
1802 xcfa1 = cfa1[i - 1] + (cfa1[i] - cfa1[i - 1]) * (palt - zmin) / (zmax - zmin);
1803 xcfa2 = cfa2[i - 1] + (cfa2[i] - cfa2[i - 1]) * (palt - zmin) / (zmax - zmin);
1804 xcfa3 = cfa3[i - 1] + (cfa3[i] - cfa3[i - 1]) * (palt - zmin) / (zmax - zmin);
1805 }
1806
1807 if (i == 0)
1808 {
1809 xcfr1 = cfr1[0];
1810 xcfr2 = cfr2[0];
1811 xcfa1 = cfa1[0];
1812 xcfa2 = cfa2[0];
1813 xcfa3 = cfa3[0];
1814 }
1815
1816 fra0 = (double)(1. - xcfr1 * exp(-r * xcfr2) - (1. - xcfr1) * exp(-r * 0.08));
1817 fae0 = (double)(1. - xcfa1 * exp(-r * xcfa2) - (1. - xcfa1) * exp(-r * xcfa3));
1818 }
1819
1820 /* correction of the effect of the view zenith angle */
1821 xlnv = (double)log(xmuv);
1822 fra = (double)(fra0 * (xlnv * (1 - fra0) + 1));
1823 fae = (double)(fae0 * ((1 + a0 * xlnv + b0 * xlnv * xlnv) + fae0 * (a1 * xlnv + b1 * xlnv * xlnv) +
1824 fae0 * fae0 * ((-a1 - a0) * xlnv + ( - b1 - b0) * xlnv * xlnv)));
1825
1826 if ((difa + difr) > 1e-03) fr = (fae * difa + fra * difr) / (difa + difr);
1827 else fr = 1;
1828}
1829
1830
Note: See TracBrowser for help on using the repository browser.