1 extern "C" {
2 #include <grass/gis.h>
3 #include <grass/glocale.h>
4 }
5
6 #include "common.h"
7 #include "aerosolmodel.h"
8 #include "atmosmodel.h"
9 #ifdef WIN32
10 #pragma warning(disable:4305) /* disable warning about initialization of a double by a double */
11 #endif /* WIN32 */
12
13 /* (background desert model...) */
bdm()14 void AerosolModel::bdm()
15 {
16 sixs_aerbas.ph = &sixs_aerbas.bdm_ph;
17 }
18
19 /* (biomass burning model...) */
bbm()20 void AerosolModel::bbm()
21 {
22 sixs_aerbas.ph = &sixs_aerbas.bbm_ph;
23 }
24
25 /* (stratospherique aerosol model...) */
stm()26 void AerosolModel::stm()
27 {
28 sixs_aerbas.ph = &sixs_aerbas.stm_ph;
29 }
30
31 /* dust model */
dust()32 void AerosolModel::dust()
33 {
34 sixs_aerbas.ph = &(sixs_aerbas.dust_ph);
35 }
36
37 /* water model */
wate()38 void AerosolModel::wate()
39 {
40 sixs_aerbas.ph = &(sixs_aerbas.wate_ph);
41 }
42
43 /* ocean model */
ocea()44 void AerosolModel::ocea()
45 {
46 sixs_aerbas.ph = &sixs_aerbas.ocea_ph;
47 }
48
49 /* soot model */
soot()50 void AerosolModel::soot()
51 {
52 sixs_aerbas.ph = &sixs_aerbas.soot_ph;
53 }
54
55 /* (user defined model from size distribution) */
56 /* To compute, using the scattering of electromagnetic waves by a homogeneous
57 isotropic sphere, the physical properties of particles whose sizes are
58 comparable to or larger than the wavelength, and to generate mixture of
59 dry particles. */
mie(double (& ex)[4][10],double (& sc)[4][10],double (& asy)[4][10])60 void AerosolModel::mie(double (&ex)[4][10], double (&sc)[4][10], double (&asy)[4][10])
61 {
62 double np[4];
63 double ext[10][4];
64 double sca2[10][4];
65 double p1[10][4][83];
66
67 const double rmul = 0.99526231496887960135245539673954; /*rlogpas = 0.030; (10**rlogpas-1.D+00)*/
68
69 int i;
70 for(i = 0; i < mie_in.icp; i++)
71 {
72
73
74
75 np[i] = 0;
76 for(int j = 0; j < 10; j++)
77 {
78 ex[i][j] = 0;
79 sc[i][j] = 0;
80 ext[j][i] = 0;
81 sca2[j][i] = 0;
82 for(int k = 0; k < 83; k++) p1[j][i][k] = 0;
83 }
84 }
85
86 double r;
87 double dr;
88 double nr = 0;
89 /* LOOPS ON THE NUMBER OF PARTICLE TYPE (4 max) */
90 for(i = 0; i < mie_in.icp; i++)
91 {
92 r = mie_in.rmin;
93 dr = r*rmul;
94
95 /* LOOPS ON THE RADIUS OF THE PARTICLE */
96 do {
97
98 /* call of the size distribution nr. For our computation, we need dn/dr for */
99 /* all functions except for sun-photometer inputs for which we need dV/dlog(r) */
100
101 switch(iaer-7)
102 {
103 case 1:
104 {
105 /* --- Mixture of particles (Log-Normal distribution functions, up to 5)*/
106 const double sqrt2PI = 2.506628274631000502415765284811;
107 const double ln10 = 2.3025850929940456840179914546844;
108 double log10_x2 = log10(mie_in.x2[i]);
109 double sq = log10(r/mie_in.x1[i])/log10_x2;
110 nr = exp(-0.5*sq*sq);
111 nr /= sqrt2PI * log10_x2 * ln10 * r;
112 break;
113 }
114 case 2:
115 {
116 /* --- Modified Gamma distribution function */
117 const double ldexp = -300.;
118 double arg=-mie_in.x2[i]*pow(r,(double)mie_in.x3[i]);
119 if(arg > ldexp) nr = pow(r,(double)mie_in.x1[i])*exp(arg);
120 else nr = 0;
121 break;
122 }
123 case 3:
124 {
125 /* --- Junge power-law function */
126 nr = pow(0.1,-(double)mie_in.x1[i]);
127 if(r > 0.1) nr = pow(r,-(double)mie_in.x1[i]);
128 break;
129 }
130 case 4:
131 {
132 /* --- from sun photometer */
133 nr = 0;
134 for(int j = 1; j < mie_in.irsunph; j++)
135 if((r-mie_in.rsunph[j]) < 0.000001)
136 {
137 nr = (r - mie_in.rsunph[j-1])/(mie_in.rsunph[j]-mie_in.rsunph[j-1]);
138 nr = mie_in.nrsunph[j-1] + nr*(mie_in.nrsunph[j] - mie_in.nrsunph[j-1]);
139 break;
140 }
141 }
142 }
143
144 /* The Mie's calculations have to be called several times (min=2, max=10 for
145 each type of particle): at wavelengths bounding the range of the selected
146 wavelengths,and at 0.550 microns to normalized the extinction coefficient
147 (if it's not in the selected range of wavelengths). */
148
149 double xndpr2 = nr * dr * M_PI * (r * r);
150 /* relatif number of particle for each type of particle (has to be equal to 1) */
151 np[i]+= nr * dr;
152
153 for(int j = 0; j < 10; j++)
154 {
155 if( (xndpr2*mie_in.cij[i]) < (1e-8 / sqrt(sixs_disc.wldis[j])) ) break;
156
157 double alpha = 2. * M_PI * r / sixs_disc.wldis[j];
158 double Qext, Qsca;
159 double p11[83];
160 exscphase(alpha, mie_in.rn[j][i], mie_in.ri[j][i], Qext, Qsca, p11);
161 ext[j][i] += xndpr2 * Qext;
162 sca2[j][i] += xndpr2 * Qsca;
163
164 /* phase function for each type of particle */
165 for(int k = 0; k < 83; k++) p1[j][i][k] += p11[k]*xndpr2;
166 }
167
168 r += dr;
169 dr = r * rmul;
170 } while (r < mie_in.rmax);
171 }
172
173
174 /* NOW WE MIXTE THE DIFFERENT TYPES OF PARTICLE
175 computation of the scattering and extinction coefficients. We first start
176 at 0.550 micron (the extinction coefficient is normalized at 0.550 micron) */
177 int j;
178 for(j = 0; j < 10; j++)
179 for(i = 0; i < mie_in.icp; i++)
180 {
181 ext[j][i] /= np[i] * 1000;
182 sca2[j][i] /= np[i] * 1000;
183 ex[0][j] += (double)(mie_in.cij[i] * ext[j][i]);
184 sc[0][j] += (double)(mie_in.cij[i] * sca2[j][i]);
185 }
186
187 /* computation of the phase function and the asymetry coefficient
188 of the mixture of particles */
189
190 for(j = 0; j < 10; j++)
191 {
192 double asy_n = 0;
193 double asy_d = 0;
194
195 for(int k = 0; k < 83; k++)
196 {
197 sixs_aerbas.usr_ph[j][k] = 0;
198 for(i = 0; i < mie_in.icp; i++)
199 sixs_aerbas.usr_ph[j][k] += (double)(mie_in.cij[i] * p1[j][i][k] / np[i] / 1000);
200
201 sixs_aerbas.usr_ph[j][k] += (double)sc[0][j];
202 asy_n += sixs_sos.cgaus[k] * sixs_aerbas.usr_ph[j][k] * sixs_sos.pdgs[k] / 10.;
203 asy_d += sixs_aerbas.usr_ph[j][k] * sixs_sos.pdgs[k] / 10.;
204 }
205
206 asy[0][j] = (double)(asy_n / asy_d);
207 }
208
209 sixs_aerbas.ph = &sixs_aerbas.usr_ph;
210 }
211
212 /***************************************************************************
213 Using the Mie's theory, this subroutine compute the scattering and
214 extinction efficiency factors (usually written Qsca and Qext) and it also
215 compute the scattering intensity efficiency
216 ***************************************************************************/
exscphase(const double X,const double nr,const double ni,double & Qext,double & Qsca,double (& p11)[83])217 void AerosolModel::exscphase(const double X, const double nr,
218 const double ni, double& Qext,
219 double& Qsca, double (&p11)[83])
220 {
221 double Ren=nr/(nr*nr+ni*ni);
222 double Imn=ni/(nr*nr+ni*ni);
223
224 /* ---Identification of the greater order of computation (=mu)
225 as defined by F.J. Corbato, J. Assoc. Computing Machinery, 1959,
226 6, 366-375 */
227 int N=int(0.5*(-1+sqrt(1+4*X*X)))+1;
228 if (N == 1) N = 2;
229
230 int mu2 = 1000000;
231 double Up = 2 * X / (2 * N + 1);
232 int mu1 = int(N + 30. * (0.1 + 0.35 * Up * (2. - Up * Up) / 2 / (1 - Up)));
233 int Np = int(X - 0.5 * sqrt(30. * 0.35 * X));
234
235 if(Np > N)
236 {
237 Up = 2 * X / (2 * Np + 1);
238 mu2 = int(Np + 30. * (0.1 + 0.35 * Up * (2 - Up * Up)/ 2 / (1 - Up)));
239 }
240
241 int mu = (mu1 < mu2) ? mu1 : mu2; /* min(mu1, mu2) */
242
243 /* --- Identification of the transition line. Below this line the Bessel
244 function j behaves as oscillating functions. Above the behavior
245 becomes monotonic. We start at a order greater than this transition
246 line (order max=mu) because a downward recursion is called for. */
247
248 double Rn[10001], xj[10001];
249 int k = mu + 1;
250
251 Rn[mu] = 0;
252 int mub;
253 while(true)
254 {
255 k--;
256 xj[k] = 0;
257 Rn[k-1] = X / (2 * k + 1 - X * Rn[k]);
258
259 if ( k == 2)
260 {
261 xj[mu + 1] = 0;
262 xj[mu] = 1;
263 mub = mu;
264 break;
265 }
266
267 if(Rn[k-1] > 1)
268 {
269 xj[k] = Rn[k-1];
270 xj[k-1] = 1;
271 mub = k - 1;
272 break;
273 }
274 }
275
276 for(k = mub; k >= 1; k--) xj[k-1] = (2 * k + 1) * xj[k] / X - xj[k+1];
277 double coxj = xj[0] - X * xj[1] * cos(X) + X * xj[0] * sin(X);
278
279 /* --- Computation Dn(alpha) and Dn(alpha*m) (cf MIE's theory)
280 downward recursion - real and imaginary parts */
281
282 double RDnY[10001];
283 double IDnY[10001];
284 double RDnX[10001];
285 RDnY[mu] = 0;
286 IDnY[mu] = 0;
287 RDnX[mu] = 0;
288
289
290 for(k = mu; k >= 1; k--)
291 {
292 RDnX[k-1] = k / X - 1 / (RDnX[k] + k / X);
293 double XnumRDnY = RDnY[k] + Ren * k / X;
294 double XnumIDnY = IDnY[k] + Imn * k / X;
295 double XdenDnY = XnumRDnY * XnumRDnY + XnumIDnY * XnumIDnY;
296 RDnY[k-1] = k * Ren / X - XnumRDnY / XdenDnY;
297 IDnY[k-1] = k * Imn / X + XnumIDnY / XdenDnY;
298 }
299
300 /* --- Initialization of the upward recursions
301 macro to help keep indexing correct, can't be to safe */
302 #define INDEX(X) ((X)+1)
303 double xy[10002];
304 xy[INDEX(-1)] = sin(X) / X;
305 xy[INDEX(0)] = -cos(X) / X;
306
307 double RGnX[10001];
308 double IGnX[10001];
309 RGnX[0] = 0;
310 IGnX[0] = -1;
311 Qsca = 0;
312 Qext = 0;
313
314 double RAn[10001];
315 double IAn[10001];
316 double RBn[10001];
317 double IBn[10001];
318
319 for(k = 1; k <= mu; k++)
320 {
321 if (k <= mub) xj[k] /= coxj;
322 else xj[k] = Rn[k-1] * xj[k-1];
323
324 /* --- Computation of bessel's function y(alpha) */
325 xy[INDEX(k)] = (2 * k - 1) * xy[INDEX(k-1)] / X - xy[INDEX(k-2)];
326 double xJonH = xj[k] / ( xj[k] * xj[k] + xy[INDEX(k)] * xy[INDEX(k)] );
327
328 /* --- Computation of Gn(alpha), Real and Imaginary part */
329 double XdenGNX = (RGnX[k-1] - k/X)*(RGnX[k-1] - k/X) + IGnX[k-1] * IGnX[k-1];
330 RGnX[k] = (k / X - RGnX[k-1])/XdenGNX - k / X;
331 IGnX[k] = IGnX[k-1] / XdenGNX;
332
333 /* --- Computation of An(alpha) and Bn(alpha), Real and Imaginary part */
334 double Xnum1An = RDnY[k] - nr * RDnX[k];
335 double Xnum2An = IDnY[k] + ni * RDnX[k];
336 double Xden1An = RDnY[k] - nr * RGnX[k] - ni * IGnX[k];
337 double Xden2An = IDnY[k] + ni * RGnX[k] - nr * IGnX[k];
338 double XdenAn = Xden1An * Xden1An + Xden2An * Xden2An;
339 double RAnb = (Xnum1An * Xden1An + Xnum2An * Xden2An) / XdenAn;
340 double IAnb = (-Xnum1An * Xden2An + Xnum2An * Xden1An) / XdenAn;
341 RAn[k] = xJonH * (xj[k] * RAnb - xy[INDEX(k)] * IAnb);
342 IAn[k] = xJonH * (xy[INDEX(k)] * RAnb + xj[k] * IAnb);
343
344 double Xnum1Bn = nr * RDnY[k] + ni * IDnY[k] - RDnX[k];
345 double Xnum2Bn = nr * IDnY[k] - ni * RDnY[k];
346 double Xden1Bn = nr * RDnY[k] + ni * IDnY[k] - RGnX[k];
347 double Xden2Bn = nr * IDnY[k] - ni * RDnY[k] - IGnX[k];
348 double XdenBn = Xden1Bn * Xden1Bn + Xden2Bn * Xden2Bn;
349 double RBnb = (Xnum1Bn * Xden1Bn + Xnum2Bn * Xden2Bn) / XdenBn;
350 double IBnb = (-Xnum1Bn * Xden2Bn + Xnum2Bn * Xden1Bn) / XdenBn;
351 RBn[k] = xJonH * (xj[k] * RBnb - xy[INDEX(k)] * IBnb);
352 IBn[k] = xJonH * (xy[INDEX(k)] * RBnb + xj[k] * IBnb);
353
354 /* ---Criterion on the recursion formulas as defined by D. Deirmendjian
355 et al., J. Opt. Soc. Am., 1961, 51, 6, 620-633 */
356 double temp = RAn[k] * RAn[k] + IAn[k] * IAn[k] + RBn[k] * RBn[k] + IBn[k] * IBn[k];
357 if((temp/k) < 1e-14)
358 {
359 mu = k;
360 break;
361 }
362
363 /* --- Computation of the scattering and extinction efficiency factor */
364 double xpond = 2 / X / X * (2 * k + 1);
365 Qsca = Qsca + xpond * temp;
366 Qext = Qext + xpond * (RAn[k] + RBn[k]);
367 }
368
369
370 /* --- Computation of the amplitude functions S1 and S2 (cf MIE's theory)
371 defined by PIn, TAUn, An and Bn with PIn and TAUn related to the
372 Legendre polynomials. */
373 for(int j = 0; j < 83; j++)
374 {
375 double RS1 = 0;
376 double RS2 = 0;
377 double IS1 = 0;
378 double IS2 = 0;
379 double PIn[10001];
380 double TAUn[10001];
381
382 PIn[0] = 0;
383 PIn[1] = 0;
384 TAUn[1] = sixs_sos.cgaus[j];
385
386 for(k = 1; k <= mu; k++)
387 {
388 double co_n = (2.0 * k + 1) / k / (k + 1);
389 RS1 += co_n * (RAn[k] * PIn[k] + RBn[k] * TAUn[k]);
390 RS2 += co_n * (RAn[k] * TAUn[k] + RBn[k] * PIn[k]);
391 IS1 += co_n * (IAn[k] * PIn[k] + IBn[k] * TAUn[k]);
392 IS2 += co_n * (IAn[k] * TAUn[k] + IBn[k] * PIn[k]);
393
394 PIn[k+1] = ((2 * k + 1) * sixs_sos.cgaus[j] * PIn[k] - (k + 1) * PIn[k-1])/k;
395 TAUn[k+1] = (k + 1) * sixs_sos.cgaus[j] * PIn[k + 1] - (k + 2) * PIn[k];
396 }
397
398 /* --- Computation of the scattering intensity efficiency */
399 p11[j] = 2 * (RS1 *RS1 + IS1 * IS1 + RS2 * RS2 + IS2 * IS2)/X/X;
400 }
401 }
402
403
404 /* load parameters from .mie file */
load()405 void AerosolModel::load()
406 {
407 int i;
408 ifstream in(filename.c_str());
409 cin.ignore(numeric_limits<int>::max(),'\n'); /* ignore this line */
410
411 in.ignore(8);
412 for(i = 0; i < 10; i++)
413 {
414 in.ignore(3);
415 in >> sixs_aer.ext[i];
416 in.ignore(6);
417 in >> sca[i];
418 in.ignore(6);
419 in >> sixs_aer.ome[i];
420 in.ignore(6);
421 in >> sixs_aer.gasym[i];
422 in.ignore(3);
423 cin.ignore(numeric_limits<int>::max(),'\n'); /* ignore the rest */
424 }
425
426 /* ignore 3 lines */
427 cin.ignore(numeric_limits<int>::max(),'\n'); /* ignore this line */
428 cin.ignore(numeric_limits<int>::max(),'\n'); /* ignore this line */
429 cin.ignore(numeric_limits<int>::max(),'\n'); /* ignore this line */
430
431 for(i = 0; i < 83; i++)
432 {
433 in.ignore(8);
434 for(int j = 0; j < 10; j++)
435 {
436 in.ignore(1);
437 in >> sixs_sos.phasel[j][i];
438 }
439 cin.ignore(numeric_limits<int>::max(),'\n'); /* ignore the rest */
440 }
441 }
442
443
444 /* do we wish to save this? */
save()445 void AerosolModel::save()
446 {
447 ofstream out(filename.c_str());
448 /* output header */
449 out << " Wlgth Nor_Ext_Co Nor_Sca_Co Sg_Sca_Alb Asymm_Para Extinct_Co Scatter_Co" << endl;
450 int i;
451 for(i = 0; i < 10; i++)
452 {
453 out << setprecision(4); /* set the required precision */
454 out << " " << setw(10) << sixs_disc.wldis[0]
455 << " " << setw(10) << sixs_aer.ext[i]
456 << " " << setw(10) << sca[i]
457 << " " << setw(10) << sixs_aer.ome[i]
458 << " " << setw(10) << sixs_aer.gasym[i]
459 << " " << setw(10) << sixs_aer.ext[i]/nis
460 << " " << setw(10) << sca[i]/nis << endl;
461 }
462
463 out << endl << endl << setw(20) << " " << " Phase Function " << endl;
464 out << " TETA ";
465 for(i = 0; i < 10; i++) out << " " << setw(10) << sixs_disc.wldis[i] << " ";
466 out << endl;
467
468 for(i = 0; i < 83; i++)
469 {
470 out << setprecision(2);
471 out << " " << setw(8) << (180.*acos(sixs_sos.cgaus[i])/M_PI);
472
473 out << setprecision(4);
474 out.setf(ios::scientific, ios::floatfield);
475 for(int j = 0; j < 10; j++) out << " " << setw(14) << sixs_sos.phasel[j][i];
476 out.setf(ios::fixed, ios::floatfield);
477 out << endl;
478 }
479 }
480
481
482 /*
483 To compute the optical scattering parameters (extinction and scattering
484 coefficients, single scattering albedo, phase function, assymetry factor) at the ten discrete
485 wavelengths for the selected model (or created model) from:
486 (1) the characteristics of the basic components of the International Radiation Commission.
487 (1983).
488 dust-like component (D.L., SUBROUTINE DUST)
489 oceanic component (O.C., SUBROUTINE OCEA)
490 water-soluble component (W.S., SUBROUTINE WATE)
491 soot component (S.O., SUBROUTINE SOOT)
492 (2) pre-computed caracteristics,
493 now available are the desertic aerosol model corresponding to background conditions, as
494 described in Shettle(1984), a stratospheric aerosol model as measured Mona Loa (Hawaii)
495 during El Chichon eruption and as described by King et al. (1984), and a biomass burning
496 aerosol model as deduced from measurements taken by sunphotometers in Amazonia.
497 (SUBROUTINES BDM, STM and BBM)
498 (3) computed using the MIE theory with inputs (size distribution, refractive indexes...) given
499 by the user (see SUBROUTINES MIE and EXSCPHASE).
500 These models don't correspond to a mixture of the four basic components.
501 */
aeroso(const double xmud)502 void AerosolModel::aeroso(const double xmud)
503 {
504 /* sra basic components for aerosol model, extinction coefficients are */
505 /* in km-1. */
506 /* dust-like = 1 */
507 /* water-soluble = 2 */
508 /* oceanique = 3 */
509 /* soot = 4 */
510 static const double vi[4] = { 113.983516, 1.13983516e-4, 5.1444150196, 5.977353425e-5 };
511 static const double ni[4] = { 54.734, 1868550., 276.05, 1805820. };
512
513 /* i: 1=dust-like 2=water-soluble 3=oceanic 4=soot */
514 static const double s_ex[4][10] =
515 {
516 {0.1796674e-01,0.1815135e-01,0.1820247e-01,0.1827016e-01,0.1842182e-01,
517 0.1853081e-01,0.1881427e-01,0.1974608e-01,0.1910712e-01,0.1876025e-01},
518 {0.7653460e-06,0.6158538e-06,0.5793444e-06,0.5351736e-06,0.4480091e-06,
519 0.3971033e-06,0.2900993e-06,0.1161433e-06,0.3975192e-07,0.1338443e-07},
520 {0.3499458e-02,0.3574996e-02,0.3596592e-02,0.3622467e-02,0.3676341e-02,
521 0.3708866e-02,0.3770822e-02,0.3692255e-02,0.3267943e-02,0.2801670e-02},
522 {0.8609083e-06,0.6590103e-06,0.6145787e-06,0.5537643e-06,0.4503008e-06,
523 0.3966041e-06,0.2965532e-06,0.1493927e-06,0.1017134e-06,0.6065031e-07}
524 };
525
526 static const double s_sc[4][10] =
527 {
528 {0.1126647e-01,0.1168918e-01,0.1180978e-01,0.1196792e-01,0.1232056e-01,
529 0.1256952e-01,0.1319347e-01,0.1520712e-01,0.1531952e-01,0.1546761e-01},
530 {0.7377123e-06,0.5939413e-06,0.5587120e-06,0.5125148e-06,0.4289210e-06,
531 0.3772760e-06,0.2648252e-06,0.9331806e-07,0.3345499e-07,0.1201109e-07},
532 {0.3499455e-02,0.3574993e-02,0.3596591e-02,0.3622465e-02,0.3676338e-02,
533 0.3708858e-02,0.3770696e-02,0.3677038e-02,0.3233194e-02,0.2728013e-02},
534 {0.2299196e-06,0.1519321e-06,0.1350890e-06,0.1155423e-06,0.8200095e-07,
535 0.6469735e-07,0.3610638e-07,0.6227224e-08,0.1779378e-08,0.3050002e-09}
536 };
537
538 static const double ex2[10] =
539 {
540 43.83631f, 42.12415f, 41.57425f, 40.85399f, 39.1404f,
541 37.89763f, 34.67506f, 24.59f, 17.96726f, 10.57569f
542 };
543
544 static const double sc2[10] =
545 {
546 40.28625f, 39.04473f, 38.6147f, 38.03645f, 36.61054f,
547 35.54456f, 32.69951f, 23.41019f, 17.15375f,10.09731f
548 };
549
550 static const double ex3[10] =
551 {
552 95397.86f, 75303.6f, 70210.64f, 64218.28f, 52430.56f,
553 45577.68f, 31937.77f, 9637.68f, 3610.691f, 810.5614f
554 };
555
556 static const double sc3[10] =
557 {
558 92977.9f, 73397.17f, 68425.49f, 62571.8f, 51049.87f,
559 44348.77f, 31006.21f, 9202.678f, 3344.476f, 664.1915f
560 };
561
562 static const double ex4[10] =
563 {
564 54273040.f, 61981440.f, 63024320.f, 63489470.f, 61467600.f,
565 58179720.f, 46689090.f, 15190620.f, 5133055.f, 899859.4f
566 };
567
568
569 static const double sc4[10] =
570 {
571 54273040.f, 61981440.f, 63024320.f, 63489470.f, 61467600.f,
572 58179720.f, 46689090.f, 15190620.f, 5133055.f, 899859.4f
573 };
574
575 static const double s_asy[4][10] =
576 {
577 {0.896,0.885,0.880,0.877,0.867,0.860,0.845,0.836,0.905,0.871},
578 {0.642,0.633,0.631,0.628,0.621,0.616,0.610,0.572,0.562,0.495},
579 {0.795,0.790,0.788,0.781,0.783,0.782,0.778,0.783,0.797,0.750},
580 {0.397,0.359,0.348,0.337,0.311,0.294,0.253,0.154,0.103,0.055}
581 };
582
583 static const double asy2[10] = { .718f, .712f, .71f, .708f, .704f, .702f, .696f, .68f, .668f, .649f };
584
585 static const double asy3[10] = { .704f, .69f, .686f, .68f, .667f, .659f, .637f, .541f, .437f, .241f };
586 static const double asy4[10] = { .705f, .744f, .751f, .757f, .762f, .759f, .737f, .586f, .372f, .139f };
587
588 /* local */
589 double coef;
590 double sigm;
591 double sumni;
592 double dd[4][10];
593 double pha[5][10][83];
594
595 double ex[4][10];
596 double sc[4][10];
597 double asy[4][10];
598
599 int i; /* crappy VS6 */
600 /* initialize ex, sc & asy */
601 for(i = 0; i < 4; i++)
602 {
603 int j;
604 for(j = 0; j < 10; j++) ex[i][j] = s_ex[i][j];
605 for(j = 0; j < 10; j++) sc[i][j] = s_sc[i][j];
606 for(j = 0; j < 10; j++) asy[i][j] = s_asy[i][j];
607 }
608
609 /* optical properties of aerosol model computed from sra basic comp */
610 for (i = 0; i < 10; ++i)
611 {
612 if(i == 4 && iaer == 0) sixs_aer.ext[i] = 1.f;
613 else sixs_aer.ext[i] = 0.f;
614 sca[i] = 0.f;
615 sixs_aer.ome[i] = 0.f;
616 sixs_aer.gasym[i] = 0.f;
617 sixs_aer.phase[i] = 0.f;
618
619 for (int k = 0; k < 83; ++k) sixs_sos.phasel[i][k] = 0.f;
620 }
621
622 /* return if iear = 0 */
623 if(iaer == 0) return;
624
625 /* look for an interval in cgaus */
626 long int j1 = -1;
627 for (i = 0; i < 82; ++i)
628 if (xmud >= sixs_sos.cgaus[i] && xmud < sixs_sos.cgaus[i+1]) { j1 = i; break; }
629 if(j1 == -1) return; /* unable to find interval */
630
631 coef = -(xmud - sixs_sos.cgaus[j1]) / (sixs_sos.cgaus[j1+1] - sixs_sos.cgaus[j1]);
632
633 switch(iaer)
634 {
635 case 12: /* read from file */
636 {
637 load();
638 for(i = 0; i < 10; i++)
639 sixs_aer.phase[i] = (double)(sixs_sos.phasel[i][j1] +
640 coef*(sixs_sos.phasel[i][j1]-sixs_sos.phasel[i][j1+1]));
641 return;
642 }
643 case 5:
644 {
645 for(i = 0; i < 10; i++)
646 {
647 asy[0][i] = asy2[i];
648 ex[0][i] = ex2[i];
649 sc[0][i] = sc2[i];
650 }
651 break;
652 }
653 case 6:
654 {
655 for(i = 0; i < 10; i++)
656 {
657 asy[0][i] = asy3[i];
658 ex[0][i] = ex3[i];
659 sc[0][i] = sc3[i];
660 }
661 break;
662 }
663 case 7:
664 {
665 for(i = 0; i < 10; i++)
666 {
667 asy[0][i] = asy4[i];
668 ex[0][i] = ex4[i];
669 sc[0][i] = sc4[i];
670 }
671 break;
672 }
673 default:;
674 }
675
676
677 if(iaer >= 5 && iaer <= 11)
678 {
679 /* calling a special aerosol model */
680
681 switch(iaer)
682 {
683 /* (background desert model...) */
684 case 5: bdm(); break;
685 /* (biomass burning model...) */
686 case 6: bbm(); break;
687 /* (stratospherique aerosol model...) */
688 case 7: stm(); break;
689
690 /* (user defined model from size distribution) */
691 case 8:
692 case 9:
693 case 10:
694 case 11: mie (ex, sc, asy); break;
695 }
696
697 for (i = 0; i < 10; i++)
698 {
699 dd[0][i] = (*sixs_aerbas.ph)[i][j1] + coef * ((*sixs_aerbas.ph)[i][j1] - (*sixs_aerbas.ph)[i][j1+1]);
700 for(int k = 0; k < 83; k++) pha[0][i][k] = (*sixs_aerbas.ph)[i][k];
701 }
702
703 mie_in.icp = 1;
704 mie_in.cij[0] = 1.f;
705 /* for normalization of the extinction coefficient */
706 nis = 1. / ex[0][3];
707 }
708 else {
709 /* calling each sra components */
710 mie_in.icp = 4;
711 /* -dust */
712 dust();
713 for(i = 0; i < 10; i++)
714 {
715 dd[0][i] = (*sixs_aerbas.ph)[i][j1] + coef * ((*sixs_aerbas.ph)[i][j1] - (*sixs_aerbas.ph)[i][j1+1]);
716 for(int k = 0; k < 83; k++) pha[0][i][k] = ((*sixs_aerbas.ph))[i][k];
717 }
718
719 /* -water soluble */
720 wate();
721 for(i = 0; i < 10; i++)
722 {
723 dd[1][i] = (*sixs_aerbas.ph)[i][j1]+coef*((*sixs_aerbas.ph)[i][j1]-(*sixs_aerbas.ph)[i][j1+1]);
724 for(int k = 0; k < 83; k++) pha[1][i][k] = (*sixs_aerbas.ph)[i][k];
725 }
726
727
728 /* -oceanic type */
729 ocea();
730 for(i = 0; i < 10; i++)
731 {
732 dd[2][i] = (*sixs_aerbas.ph)[i][j1]+coef*((*sixs_aerbas.ph)[i][j1]-(*sixs_aerbas.ph)[i][j1+1]);
733 for(int k = 0; k < 83; k++) pha[2][i][k] = (*sixs_aerbas.ph)[i][k];
734 }
735
736 /* -soot */
737 soot();
738 for(i = 0; i < 10; i++)
739 {
740 dd[3][i] = (*sixs_aerbas.ph)[i][j1]+coef*((*sixs_aerbas.ph)[i][j1]-(*sixs_aerbas.ph)[i][j1+1]);
741 for(int k = 0; k < 83; k++) pha[3][i][k] = (*sixs_aerbas.ph)[i][k];
742 }
743
744 /* summ of the c/vi calculation */
745 sumni = 0.f;
746 sigm = 0.f;
747
748 for(i = 0; i < 4; i++) sigm+=(double)(c[i]/vi[i]);
749
750 /* cij coefficients calculation */
751 for(i = 0; i < 4; i++) {
752 mie_in.cij[i] = (double)(c[i]/vi[i]/sigm);
753 sumni += mie_in.cij[i]/ni[i];
754 }
755
756 nis = 1. / sumni;
757 }
758
759
760 /* mixing parameters calculation */
761 for(i = 0; i < 10; i++)
762 {
763 for(int j = 0; j < mie_in.icp; j++)
764 {
765 sixs_aer.ext[i] += (double)(ex[j][i] * mie_in.cij[j]);
766 sca[i] += (double)(sc[j][i] * mie_in.cij[j]);
767 sixs_aer.gasym[i] += (double)(sc[j][i] * mie_in.cij[j] * asy[j][i]);
768 sixs_aer.phase[i] += (double)(sc[j][i] * mie_in.cij[j] * dd[j][i]);
769
770 for(int k = 0; k < 83; k++)
771 sixs_sos.phasel[i][k] += (double)(sc[j][i] * mie_in.cij[j] * pha[j][i][k]);
772 }
773
774 sixs_aer.ome[i] = sca[i]/sixs_aer.ext[i];
775 sixs_aer.gasym[i] /= sca[i];
776 sixs_aer.phase[i] /= sca[i];
777
778 for(int k = 0; k < 83; k++) sixs_sos.phasel[i][k] /= sca[i];
779
780 sixs_aer.ext[i] *= (double)nis;
781 sca[i] *= (double)nis;
782 }
783
784 if (filename.size() != 0 && iaer >= 8 && iaer <= 11) save();
785 }
786
parse(const double xmud)787 void AerosolModel::parse(const double xmud)
788 {
789 cin >> iaer;
790 cin.ignore(numeric_limits<int>::max(),'\n');
791
792 /* initialize vars; */
793 mie_in.rmin = 0.f;
794 mie_in.rmax = 0.f;
795 mie_in.icp = 1;
796
797 int i;
798 for(i = 0; i < 4; i++)
799 {
800 mie_in.cij[i] = 0.f;
801
802 mie_in.x1[i] = 0.f;
803 mie_in.x2[i] = 0.f;
804 mie_in.x3[i] = 0.f;
805
806 for(int j = 0; j < 10; j++)
807 {
808 mie_in.rn[j][i] = 0.f;
809 mie_in.ri[j][i] = 0.f;
810 }
811 }
812
813 for(i = 0; i < 50; i++)
814 {
815 mie_in.rsunph[i] = 0.f;
816 mie_in.nrsunph[i] = 0.f;
817 }
818 mie_in.cij[0] = 1.00f;
819
820 switch (iaer)
821 {
822 case 0:
823 case 5:
824 case 6:
825 case 7: break; /* do nothing */
826
827 case 1:
828 {
829 c[0]=0.70f;
830 c[1]=0.29f;
831 c[2]=0.00f;
832 c[3]=0.01f;
833 break;
834 }
835 case 2:
836 {
837 c[0]=0.00f;
838 c[1]=0.05f;
839 c[2]=0.95f;
840 c[3]=0.00f;
841 break;
842 }
843 case 3:
844 {
845 c[0]=0.17f;
846 c[1]=0.61f;
847 c[2]=0.00f;
848 c[3]=0.22f;
849 break;
850 }
851 case 4:
852 {
853 for(i = 0; i < 4; i++) cin >> c[i];
854 cin.ignore(numeric_limits<int>::max(),'\n');
855 break;
856 }
857 case 8:
858 {
859 cin >> mie_in.rmin;
860 cin >> mie_in.rmax;
861 cin >> mie_in.icp;
862 cin.ignore(numeric_limits<int>::max(),'\n');
863
864 if(mie_in.icp >= 4) {
865 G_fatal_error(_("mie_in.icp: %ld > 4, will cause internal buffer overflow"), mie_in.icp);
866 }
867
868 for(i = 0; i < mie_in.icp; i++)
869 {
870 cin >> mie_in.x1[i];
871 cin >> mie_in.x2[i];
872 cin >> mie_in.cij[i];
873 cin.ignore(numeric_limits<int>::max(),'\n');
874
875 int j;
876 for(j = 0; j < 10; j++) cin >> mie_in.rn[j][i];
877 cin.ignore(numeric_limits<int>::max(),'\n');
878
879 for(j = 0; j < 10; j++) cin >> mie_in.ri[j][i];
880 cin.ignore(numeric_limits<int>::max(),'\n');
881 }
882 break;
883 }
884 case 9:
885 {
886 cin >> mie_in.rmin;
887 cin >> mie_in.rmax;
888 cin.ignore(numeric_limits<int>::max(),'\n');
889
890 cin >> mie_in.x1[0];
891 cin >> mie_in.x2[0];
892 cin >> mie_in.x3[0];
893 cin.ignore(numeric_limits<int>::max(),'\n');
894
895 int j;
896 for(j = 0; j < 10; j++) cin >> mie_in.rn[j][0];
897 cin.ignore(numeric_limits<int>::max(),'\n');
898
899
900 for(j = 0; j < 10; j++) cin >> mie_in.ri[j][0];
901 cin.ignore(numeric_limits<int>::max(),'\n');
902
903 break;
904 }
905 case 10:
906 {
907 cin >> mie_in.rmin;
908 cin >> mie_in.rmax;
909 cin.ignore(numeric_limits<int>::max(),'\n');
910
911 cin >> mie_in.x1[0];
912 cin.ignore(numeric_limits<int>::max(),'\n');
913
914 int j;
915 for(j = 0; j < 10; j++) cin >> mie_in.rn[j][0];
916 cin.ignore(numeric_limits<int>::max(),'\n');
917
918 for(j = 0; j < 10; j++) cin >> mie_in.ri[j][0];
919 cin.ignore(numeric_limits<int>::max(),'\n');
920
921 break;
922 }
923 case 11:
924 {
925 cin >> mie_in.irsunph;
926 cin.ignore(numeric_limits<int>::max(),'\n');
927
928 if(mie_in.irsunph >= 50) {
929 G_fatal_error(_("mie_in.irsunph: %ld > 50, will cause internal buffer overflow"), mie_in.irsunph);
930 }
931
932 for(i = 0; i < mie_in.irsunph; i++)
933 {
934 cin >> mie_in.rsunph[i];
935 cin >> mie_in.nrsunph[i];
936 cin.ignore(numeric_limits<int>::max(),'\n');
937
938 double sq = mie_in.rsunph[i]*mie_in.rsunph[i];
939 const double ln10 = 2.3025850929940456840179914546844;
940 mie_in.nrsunph[i] = (double)(mie_in.nrsunph[i]/(sq*sq)/ln10);
941 }
942 mie_in.rmin=mie_in.rsunph[0];
943 mie_in.rmax=mie_in.rsunph[mie_in.irsunph-1]+1e-07f;
944
945 for(i = 0; i < 10; i++) cin >> mie_in.rn[i][0];
946 cin.ignore(numeric_limits<int>::max(),'\n');
947
948 for(i = 0; i < 10; i++) cin >> mie_in.ri[i][0];
949 cin.ignore(numeric_limits<int>::max(),'\n');
950 break;
951 }
952 case 12:
953 { /* read file name */
954 getline(cin,filename);
955 filename = filename.substr(0, filename.find(" "));
956 break;
957 }
958 default: G_warning(_("Unknown aerosol model!"));
959 }
960
961 if(iaer >= 8 && iaer <= 11)
962 {
963 cin >> iaerp;
964 if( iaerp == 1 ) /* read file name */
965 {
966 getline(cin,filename);
967 filename = filename.substr(0, filename.find(" "));
968 filename += ".mie";
969 }
970 }
971
972 aeroso(xmud);
973 }
974
975 /* format 132 */
print132(string s)976 void AerosolModel::print132(string s)
977 {
978 Output::Begin();
979 Output::Repeat(15, ' ');
980 Output::Print(s);
981 Output::Print(" aerosols model");
982 Output::End();
983 }
984
985 /* --- aerosols model ---- */
print()986 void AerosolModel::print()
987 {
988 /* --- aerosols model (type) ---- */
989 Output::Begin();
990 Output::Repeat(10, ' ');
991 Output::Print(" aerosols type identity :");
992 Output::End();
993
994 if(iaer == 4 || (iaer >= 8 && iaer != 11))
995 {
996 Output::Begin();
997 Output::Repeat(15, ' ');
998 Output::Print(" user defined aerosols model ");
999 Output::End();
1000 }
1001
1002 switch(iaer)
1003 {
1004 case 0:
1005 {
1006 Output::Begin();
1007 Output::Repeat(15, ' ');
1008 Output::Print(" no aerosols computed ");
1009 Output::End();
1010 break;
1011 }
1012 case 1: print132(" Continental"); break;
1013 case 2: print132(" Maritime"); break;
1014 case 3: print132(" Urban"); break;
1015 case 4:
1016 {
1017 static const string desc[4] = {
1018 string(" % of dust-like"),
1019 string(" % of water-soluble"),
1020 string(" % of oceanic"),
1021 string(" % of soot")
1022 };
1023
1024 for(int i = 0; i < 4; i++)
1025 {
1026 Output::Begin();
1027 Output::Repeat(26, ' ');
1028 ostringstream s;
1029 s.setf(ios::fixed, ios::floatfield);
1030 s << setprecision(3);
1031 s << c[i] << desc[i] << ends;
1032 Output::Print(s.str());
1033 Output::End();
1034 }
1035 break;
1036 }
1037 case 5: print132(" Desertic"); break;
1038 case 6: print132(" Smoke"); break;
1039 case 7: print132(" Stratospheric"); break;
1040 case 8:
1041 {
1042 Output::Begin();
1043 Output::Repeat(15, ' ');
1044 ostringstream s;
1045 s << "using " << mie_in.icp << " Log-normal size-distribution(s)" << ends;
1046 Output::Print(s.str());
1047 Output::End();
1048
1049 Output::Begin();
1050 Output::Repeat(15, ' ');
1051 Output::Print("Mean radius Stand. Dev. Percent. dencity");
1052 Output::End();
1053
1054 for(int i = 0; i < mie_in.icp; i++)
1055 {
1056 Output::Begin();
1057 Output::Position(41);
1058 ostringstream s1;
1059 s1.setf(ios::fixed, ios::floatfield);
1060 s1 << setprecision(4);
1061 s1 << setw(10) << mie_in.x1[i] << ends;
1062 Output::Print(s1.str());
1063
1064 Output::Position(55);
1065 ostringstream s2;
1066 s2.setf(ios::fixed, ios::floatfield);
1067 s2 << setprecision(3);
1068 s2 << setw(8) << mie_in.x2[i] << ends;
1069 Output::Print(s2.str());
1070
1071 Output::Position(69);
1072 ostringstream s3;
1073 s3.setf(ios::fixed, ios::floatfield);
1074 s3 << setprecision(3);
1075 s3 << setw(11) << mie_in.cij[i] << ends;
1076 Output::Print(s3.str());
1077
1078 Output::End();
1079 }
1080 break;
1081 }
1082 case 9:
1083 {
1084 Output::Begin();
1085 Output::Repeat(15, ' ');
1086 Output::Print("using a Modified Gamma size-distribution");
1087 Output::End();
1088
1089 Output::Begin();
1090 Output::Repeat(19, ' ');
1091 Output::Print("Alpha b Gamma");
1092 Output::End();
1093
1094 Output::Begin();
1095 Output::Position(20);
1096 ostringstream s1;
1097 s1.setf(ios::fixed, ios::floatfield);
1098 s1 << setprecision(3);
1099 s1 << setw(9) << mie_in.x1[0] << ends;
1100 Output::Print(s1.str());
1101
1102 Output::Position(31);
1103 ostringstream s2;
1104 s2.setf(ios::fixed, ios::floatfield);
1105 s2 << setprecision(3);
1106 s2 << setw(9) << mie_in.x2[0] << ends;
1107 Output::Print(s2.str());
1108
1109 Output::Position(47);
1110 ostringstream s3;
1111 s3.setf(ios::fixed, ios::floatfield);
1112 s3 << setprecision(3);
1113 s3 << setw(9) << mie_in.x3[0] << ends;
1114 Output::Print(s3.str());
1115
1116 Output::End();
1117 break;
1118 }
1119 case 10:
1120 {
1121 Output::Begin();
1122 Output::Repeat(15, ' ');
1123 Output::Print("using a Power law size-distribution with alpha=");
1124 ostringstream s;
1125 s.setf(ios::fixed, ios::floatfield);
1126 s << setprecision(1);
1127 s << setw(4) << mie_in.x1[0] << ends;
1128 Output::Print(s.str());
1129 Output::End();
1130
1131 break;
1132 }
1133 case 11: print132(" Sun Photometer"); break;
1134 case 12:
1135 {
1136 Output::Begin();
1137 Output::Repeat(15, ' ');
1138 Output::Print("using data from the file:");
1139 Output::End();
1140
1141 Output::Begin();
1142 Output::Position(25);
1143 Output::Print(filename);
1144 Output::End();
1145 }
1146 }
1147
1148 if(iaer > 7 && iaerp == 1)
1149 {
1150 Output::Begin();
1151 Output::Repeat(15, ' ');
1152 Output::Print(" results saved into the file:");
1153 Output::End();
1154
1155 Output::Begin();
1156 Output::Position(25);
1157 Output::Print(filename);
1158 Output::End();
1159 }
1160 }
1161
Parse(const double xmud)1162 AerosolModel AerosolModel::Parse(const double xmud)
1163 {
1164 AerosolModel aero;
1165 aero.parse(xmud);
1166 return aero;
1167 }
1168
1169