11 int atmocor2(l2str *l2rec, aestr *aerec, int32_t ip) {
12 static int firstCall = 1;
13 static int want_nirRrs = 0;
14 static int want_nirLw = 0;
15 static int want_mumm = 0;
16 static int want_ramp = 1;
17 static int32_t aer_iter_max = 1;
18 static int32_t aer_iter_min = 1;
21 static float radeg =
RADEG;
24 static float p0 =
STDPR;
25 static float df = 0.33;
26 static float cbot = 0.7;
27 static float ctop = 1.3;
28 static float seed_chl = 0.0;
29 static float seed_green = 0.0;
30 static float seed_red = 0.0;
31 static float nir_chg = 0.02;
40 static int32_t swir_s;
41 static int32_t swir_l;
49 l1str *
l1rec = l2rec->l1rec;
53 int32_t brdf_opt =
input->brdf_opt;
54 int32_t aer_opt =
input->aer_opt;
55 int32_t glint_opt =
input->glint_opt;
56 int32_t cirrus_opt =
input->cirrus_opt;
58 float *Fo =
l1rec->Fo;
59 float *Fobar =
l1file->Fobar;
60 float fsol =
l1rec->fsol;
63 float delphi =
l1rec->delphi[ip];
64 float ws =
l1rec->ws [ip];
65 float gc =
l1rec->glint_coef[ip];
67 int32_t *aermodmin = &l2rec->aermodmin[ip];
68 int32_t *aermodmax = &l2rec->aermodmax[ip];
69 float *aermodrat = &l2rec->aerratio [ip];
70 int32_t *aermodmin2 = &l2rec->aermodmin2[ip];
71 int32_t *aermodmax2 = &l2rec->aermodmax2[ip];
72 float *aermodrat2 = &l2rec->aerratio2 [ip];
73 float *eps = &l2rec->eps [ip];
76 float *Lw = &l2rec->Lw [ip *
l1file->nbands];
77 float *nLw = &l2rec->nLw [ip *
l1file->nbands];
78 float *La = &l2rec->La [ip *
l1file->nbands];
79 float *taua = &l2rec->taua [ip *
l1file->nbands];
82 float *brdf = &l2rec->brdf [ip *
l1file->nbands];
83 float *
Rrs = &l2rec->Rrs [ip *
l1file->nbands];
85 float *nLw_unc = &l2rec->nLw_unc[ip *
l1file->nbands];
86 float *Rrs_unc = &l2rec->Rrs_unc[ip *
l1file->nbands];
93 float *tLw_nir, *last_tLw_nir;
96 int32_t iter, iter_max, iter_min, last_iter, iter_reset;
104 float refl_nir = 0, last_refl_nir = 0;
108 if (firstCall == 1) {
111 if ((wave = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
112 printf(
"-E- : Error allocating memory to wave\n");
115 for (ib = 0; ib < nwave; ib++) {
116 wave[ib] =
l1file->fwave[ib];
127 if (nir_s < 0 || nir_l < 0) {
128 printf(
"Aerosol selection bands %d and %d not available for this sensor\n",
133 printf(
"Invalid aerosol selection bands: long (%d) must be greater than short (%d).\n",
137 if (wave[nir_s] < 600) {
138 printf(
"Aerosol selection band(s) must be greater than 600nm");
144 daer =
MAX(nir_l - nir_s, 1);
145 cslp = 1. / (ctop - cbot);
147 printf(
"Aerosol selection bands %d and %d\n",
l1file->iwave[aer_s],
l1file->iwave[aer_l]);
156 aer_iter_max =
input->aer_iter_max;
157 printf(
"NIR correction enabled.\n");
163 aer_iter_max =
input->aer_iter_max;
164 printf(
"MUMM correction enabled.\n");
170 aer_iter_max =
input->aer_iter_max;
173 if (swir_s < 0 || swir_l < 0) {
174 printf(
"Aerosol selection bands %d and %d not available for this sensor\n",
178 printf(
"NIR/SWIR switching correction enabled.\n");
183 aer_iter_max =
input->aer_iter_max;
184 printf(
"NIR correction enabled --> for multi-scattering epsilon.\n");
189 aer_iter_max =
input->aer_iter_max;
190 printf(
"NIR correction enabled --> for spectral matching.\n");
193 if (
input->aer_rrs_short >= 0.0 &&
input->aer_rrs_long >= 0.0) {
196 aer_iter_max =
input->aer_iter_max;
197 printf(
"NIR correction via input Rrs enabled.\n");
201 if (
input->aer_iter_max < 1)
203 if ( want_nirLw || want_nirRrs) {
210 printf(
"%s line %d: can't find red band\n", __FILE__, __LINE__);
218 printf(
"%s line %d: can't find green band\n", __FILE__, __LINE__);
226 if ((taur = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
227 printf(
"-E- : Error allocating memory to taur\n");
230 if ((tLw = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
231 printf(
"-E- : Error allocating memory to tLw\n");
234 if ((
rhown_nir = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
235 printf(
"-E- : Error allocating memory to rhown_nir\n");
238 if ((tLw_nir = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
239 printf(
"-E- : Error allocating memory to tLw_nir\n");
242 if ((last_tLw_nir = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
243 printf(
"-E- : Error allocating memory to last_tLw_nir\n");
246 if ((Ltemp = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
247 printf(
"-E- : Error allocating memory to Ltemp\n");
250 if ((Lunc = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
251 printf(
"-E- : Error allocating memory to Lunc\n");
261 for (ib = 0; ib < nwave; ib++) {
263 ipb = ip * nwave + ib;
272 if (glint_opt != 2) TLg [ib] = 0.0;
276 l2rec->eps[ip] = badval;
278 if (l2rec->l1rec->Lt[ipb] <= 0.0)
279 if (wave[ib] < 1000) nneg++;
298 airmass = 1.0 /
mu0 + 1.0 /
mu;
302 for (ib = 0; ib < nwave; ib++) {
304 ipb = ip * nwave + ib;
310 Ltemp[ib] = l2rec->l1rec->Lt[ipb];
311 Lunc [ib] = l2rec->l1rec->Lt_unc[ipb];
315 Ltemp[ib] = Ltemp[ib] /
l1rec->tg_sol[ipb] /
l1rec->tg_sen[ipb];
316 Lunc [ib] = Lunc [ib] /
l1rec->tg_sol[ipb] /
l1rec->tg_sen[ipb];
321 if (cirrus_opt) Ltemp[ib] -=
l1rec->rho_cirrus[ip] / Ka * Fo[ib] *
mu0 /
pi;
324 Ltemp[ib] /=
l1rec->polcor[ipb];
325 Lunc [ib] /=
l1rec->polcor[ipb];
328 Ltemp[ib] -=
l1rec->tLf[ipb];
331 Ltemp[ib] -=
l1rec->Lr[ipb];
334 if (
input->oxaband_opt == 1) {
335 Ltemp[ib] /=
l1rec->t_o2[ipb];
336 Lunc [ib] /=
l1rec->t_o2[ipb];
342 if (brdf_opt < FOQMOREL && brdf_opt >
NOBRDF)
343 ocbrdf(l2rec, ip, brdf_opt, wave, nwvis,
solz,
senz, delphi, ws, -1.0,
NULL,
NULL, brdf);
349 iter_max = aer_iter_max;
350 iter_min = aer_iter_min;
352 last_refl_nir = 100.;
356 if (glint_opt == 1 &&
l1rec->glint_coef[ip] > glint_min) {
357 iter_max =
MAX(2, iter_max);
358 iter_min =
MAX(2, iter_min);
361 if (glint_opt == 2) {
377 daer =
MAX(aer_l - aer_s, 1);
381 if (want_nirLw || want_nirRrs) {
382 for (ib = 0; ib < nwave; ib++) {
383 last_tLw_nir[ib] = 0.0;
387 Rrs[green] = seed_green;
402 for (ib = 0; ib < nwave; ib++) {
407 if (want_glintcorr) {
410 glint_rad(iter, nwave, aer_s, aer_l,
gc, airmass,
mu0, Fo, taur, taua, tLw, TLg);
413 for (ib = 0; ib < nwave; ib++) {
424 for (ib = aer_s; ib <= aer_l; ib += daer) {
427 tLw_nir[ib] =
rhown_nir[ib] /
pi * Fo[ib] *
mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
430 if (tLw_nir[ib] > tLw[ib] && tLw[ib] > 0.0)
431 tLw_nir[ib] = tLw[ib];
434 tLw[ib] -= tLw_nir[ib];
444 for (ib = aer_s; ib <= aer_l; ib += daer) {
447 tLw_nir[ib] =
rhown_nir[ib] /
pi * Fo[ib] *
mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
450 tLw[ib] -= tLw_nir[ib];
459 get_rhown_eval(
input->fqfile,
Rrs, wave, aer_s, aer_l, nwave, &
l1rec->sw_a_avg[ipb], &
l1rec->sw_bb_avg[ipb], chl,
solz,
senz, delphi,
rhown_nir);
461 for (ib = aer_s; ib <= aer_l; ib += daer) {
464 tLw_nir[ib] =
rhown_nir[ib] /
pi * Fo[ib] *
mu0 * t_sol[ib] * t_sen[ib] / brdf[ib];
467 tLw_nir[ib] = (1.0 - df) * tLw_nir[ib] + df * last_tLw_nir[ib];
471 if (chl > 0.0 && chl <= cbot)
473 else if ((chl > cbot) && (chl < ctop))
474 tLw_nir[ib] *= (cslp * chl + cint);
478 tLw[ib] -= tLw_nir[ib];
485 status =
aerosol(l2rec, aer_opt, aerec, ip, wave, nwave, aer_s, aer_l, Fo, tLw,
486 La, t_sol, t_sen, eps, taua, aermodmin, aermodmax, aermodrat,
487 aermodmin2, aermodmax2, aermodrat2);
489 for (ib = 0; ib < nwave; ib++) {
491 ipb = ip * nwave + ib;
506 for (ib = 0; ib < nwave; ib++) {
509 tLw[ib] = tLw[ib] - La[ib];
510 Lw [ib] = tLw[ib] / t_sen[ib] *
tg_sol[ib];
511 nLw[ib] = Lw [ib] / t_sol[ib] /
tg_sol[ib] /
mu0 / fsol * brdf[ib];
517 for (ib = aer_s; ib <= aer_l; ib += daer)
518 last_tLw_nir[ib] = tLw_nir[ib];
519 for (ib = 0; ib < nwvis; ib++) {
520 Rrs[ib] = nLw[ib] / Fobar[ib];
528 if (chl == badchl && iter_reset == 0 && iter < iter_max) {
539 if (chl == badchl && iter_reset == 1 && iter < iter_max) {
550 for (ib = 0; ib < nwave; ib++) {
561 if (tindx >= 1.05 &&
status == 0) {
562 for (ib = 0; ib < nwvis; ib++) {
563 Rrs[ib] = nLw[ib] / Fobar[ib];
567 if (chl > 0 && (chl <= 1.0 || nLw[nir_l] < 0.08)) {
568 iter_max = aer_iter_max;
571 daer =
MAX(aer_l - aer_s, 1);
584 }
else if (iter < iter_min) {
586 }
else if (want_nirLw && (
fabs(refl_nir - last_refl_nir) <
fabs(nir_chg * refl_nir) || refl_nir < 0.0)) {
588 }
else if (want_mumm || want_nirRrs) {
590 }
else if (iter > iter_max) {
594 last_refl_nir = refl_nir;
598 l2rec->num_iter[ip] = iter;
615 if (want_nirLw || want_nirRrs || want_mumm) {
616 for (ib = aer_s; ib <= aer_l; ib += daer) {
617 tLw[ib] = tLw_nir[ib];
618 Lw [ib] = tLw[ib] / t_sen[ib] *
tg_sol[ib];
619 nLw[ib] = Lw [ib] / t_sol[ib] /
tg_sol[ib] /
mu0 / fsol * brdf[ib];
620 Rrs[ib] = nLw[ib] / Fobar[ib];
633 ocbrdf(l2rec, ip, brdf_opt, wave, nwvis,
solz,
senz, delphi, ws, -1.0, nLw, Fobar, brdf);
634 for (ib = 0; ib < nwvis; ib++) {
635 nLw[ib] = nLw[ib] * brdf[ib];
640 for (ib = 0; ib < nwave; ib++) {
641 if (ib != aer_s && ib != aer_l) {
642 Rrs[ib] = nLw[ib] / Fobar[ib];
643 l2rec->Rrs[ip * nwave + ib] =
Rrs[ib];
645 nLw_unc[ib] = Lunc[ib] / t_sen[ib] / t_sol[ib] /
mu0 / fsol * brdf[ib];
646 Rrs_unc[ib] = nLw_unc[ib] / Fobar[ib];