OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
qaa.c
Go to the documentation of this file.
1 
91 #include <stdlib.h>
92 #include <math.h>
93 #include <stdarg.h>
94 #include <assert.h>
95 #include "qaa.h"
96 
97 static int idx410 = -1;
98 static int idx440 = -1;
99 static int idx490 = -1;
100 static int idx555 = -1;
101 static int idx670 = -1;
102 static int initialized = 0;
103 static int aph_check = 1;
104 static double S = 0.015;
105 static double acoefs[3];
106 
112 int
114  return initialized;
115 }
116 
131 int
132 qaa_init(int i410, int i440, int i490, int i555, int i670) {
133  idx410 = i410;
134  idx440 = i440;
135  idx490 = i490;
136  idx555 = i555;
137  idx670 = i670;
138 
139  // SeaWiFS coefficients
140 
141  acoefs[0] = -1.146;
142  acoefs[1] = -1.366;
143  acoefs[2] = -0.469;
144 
145  initialized = 1;
146 
147  return 0;
148 }
149 
162 int
163 qaa_set_param(int param, ...) {
164  va_list ap;
165  va_start(ap, param);
166  switch (param) {
167  case QAA_S_PARAM:
168  S = va_arg(ap, double);
169  break;
170  case QAA_COEFS_PARAM:
171  acoefs[0] = va_arg(ap, double);
172  acoefs[1] = va_arg(ap, double);
173  acoefs[2] = va_arg(ap, double);
174  break;
175  case QAA_APH_CHECK:
176  aph_check = va_arg(ap, int);
177  break;
178  }
179  va_end(ap);
180  return 0;
181 }
182 
199 int
200 qaa_v6(int nbands, double *wavel, double *Rrs, double *aw, double *bbw,
201  double *rrs, double *u, double *a, double *bb, unsigned char *flags) {
202 
203  const double g0 = 0.08945;
204  const double g1 = 0.1247;
205 
206  int i, idxref;
207 
208  double rat, aref;
209  double bbpref;
210  double Y;
211  double rho, numer, denom;
212 
213  assert(idx440 >= 0);
214  assert(idx490 >= 0);
215  assert(idx555 >= 0);
216  assert(nbands >= 0);
217 
218  /* pre-test 670 */
219  if ((Rrs[idx670] > 20.0 * pow(Rrs[idx555], 1.5)) ||
220  (Rrs[idx670] < 0.9 * pow(Rrs[idx555], 1.7))) {
221  Rrs[idx670] = 1.27 * pow(Rrs[idx555], 1.47) + 0.00018 * pow(Rrs[idx490] / Rrs[idx555], -3.19);
222  if (flags) *flags |= 0x02;
223  }
224 
225  /* Step 0 */
226  for (i = 0; i < nbands; i++) {
227  rrs[i] = Rrs[i] / (0.52 + 1.7 * Rrs[i]);
228  if (Rrs[i] < 0.0)
229  if (flags) *flags |= 0x01;
230  }
231 
232  /* Step 1 */
233  for (i = 0; i < nbands; i++)
234  u[i] = (sqrt(g0 * g0 + 4.0 * g1 * rrs[i]) - g0) / (2.0 * g1);
235 
236  /* Step 2 */
237  if (Rrs[idx670] >= 0.0015) {
238  aref = aw[idx670] + 0.39*powf(Rrs[idx670]/(Rrs[idx440]+Rrs[idx490]),1.14);
239  idxref = idx670;
240  } else {
241  numer = rrs[idx440] + rrs[idx490];
242  denom = rrs[idx555] + 5.0 * rrs[idx670]*(rrs[idx670] / rrs[idx490]);
243  rho = log10(numer / denom);
244  rho = acoefs[0] + acoefs[1] * rho + acoefs[2] * rho*rho;
245  aref = aw[idx555] + pow(10.0, rho);
246  idxref = idx555;
247  }
248 
249  /* Step 3 */
250  bbpref = ((u[idxref] * aref) / (1.0 - u[idxref])) - bbw[idxref];
251 
252  /* Step 4 */
253  rat = rrs[idx440] / rrs[idx555];
254  Y = 2.0 * (1.0 - 1.2 * exp(-0.9 * rat));
255 
256  /* Step 5 */
257  for (i = 0; i < nbands; i++) {
258  bb[i] = bbpref * pow((wavel[idxref] / wavel[i]), Y) + bbw[i];
259  if (bb[i] < 0.0)
260  if (flags) *flags |= 0x04;
261  }
262 
263  /* Step 6 */
264  for (i = 0; i < nbands; i++) {
265  a[i] = ((1.0 - u[i]) * bb[i]) / u[i];
266  if (a[i] < 0.0)
267  if (flags) *flags |= 0x08;
268  }
269 
270  return 0;
271 }
272 
278 int
279 qaaf_v6(int nbands, float *wavel, float *Rrs, float *aw, float *bbw,
280  float *rrs, float *u, float *a, float *bb,
281  unsigned char *flags) {
282 
283  const float g0 = 0.08945;
284  const float g1 = 0.1247;
285 
286  float rho, numer, denom;
287 
288  float rat, aref;
289  float bbpref;
290  float Y;
291 
292  int i, idxref;
293 
294  assert(idx440 >= 0);
295  assert(idx490 >= 0);
296  assert(idx555 >= 0);
297  assert(idx670 >= 0);
298  assert(nbands >= 0);
299 
300  /* Test for bad Rrs at idx555 */
301  if (Rrs[idx555] <= 0.0)
302  Rrs[idx555] = 0.001;
303 
304  /* pre-test 670 */
305  if ((Rrs[idx670] > 20.0 * powf(Rrs[idx555], 1.5)) ||
306  (Rrs[idx670] < 0.9 * powf(Rrs[idx555], 1.7))) {
307  Rrs[idx670] = 1.27 * powf(Rrs[idx555], 1.47) + 0.00018 * powf(Rrs[idx490] / Rrs[idx555], -3.19);
308  if (flags) *flags |= 0x02;
309  }
310 
311  /* Step 0 */
312  for (i = 0; i < nbands; i++) {
313  rrs[i] = Rrs[i] / (0.52 + 1.7 * Rrs[i]);
314  if (Rrs[i] < 0.0)
315  if (flags) *flags |= 0x01;
316  }
317 
318  /* Step 1 */
319  for (i = 0; i < nbands; i++)
320  u[i] = (sqrt(g0 * g0 + 4.0 * g1 * rrs[i]) - g0) / (2.0 * g1);
321 
322  /* Step 2 */
323  if (Rrs[idx670] >= 0.0015) {
324  aref = aw[idx670] + 0.39*powf(Rrs[idx670]/(Rrs[idx440]+Rrs[idx490]),1.14);
325  idxref = idx670;
326  } else {
327  numer = rrs[idx440] + rrs[idx490];
328  denom = rrs[idx555] + 5.0 * rrs[idx670]*(rrs[idx670] / rrs[idx490]);
329  rho = log10f(numer / denom);
330  rho = acoefs[0] + acoefs[1] * rho + acoefs[2] * rho*rho;
331  aref = aw[idx555] + powf(10.0, rho);
332  idxref = idx555;
333  }
334 
335  /* Step 3 */
336  bbpref = ((u[idxref] * aref) / (1.0 - u[idxref])) - bbw[idxref];
337 
338  /* Step 4 */
339  rat = rrs[idx440] / rrs[idx555];
340  Y = 2.0 * (1.0 - 1.2 * expf(-0.9 * rat));
341 
342  /* Step 5 */
343  for (i = 0; i < nbands; i++) {
344  bb[i] = bbpref * pow((wavel[idxref] / wavel[i]), Y) + bbw[i];
345  if (bb[i] < 0.0)
346  if (flags) *flags |= 0x04;
347  }
348 
349  /* Step 6 */
350  for (i = 0; i < nbands; i++) {
351  a[i] = ((1.0 - u[i]) * bb[i]) / u[i];
352  if (a[i] < 0.0)
353  if (flags) *flags |= 0x08;
354  }
355 
356  return 0;
357 }
358 
379 int
380 qaa_decomp(int nbands, double *wavel, double *rrs, double *a, double *aw,
381  double *adg, double *aph, unsigned char *flags) {
382  int i;
383  double symbol, x1, x2;
384  double zeta, denom, dif1, dif2;
385  double rat, adg440;
386  double Sr;
387 
388  assert(idx410 >= 0);
389  assert(idx440 >= 0);
390  assert(idx555 >= 0);
391  assert(nbands >= 0);
392 
393  /* step 7 */
394  rat = rrs[idx440] / rrs[idx555];
395  symbol = 0.74 + (0.2 / (0.8 + rat));
396 
397  /* step 8 */
398  Sr = S + 0.002 / (0.6 + rat);
399  zeta = exp(Sr * (wavel[idx440] - wavel[idx410]));
400 
401  /* step 9 */
402  denom = zeta - symbol;
403  dif1 = a[idx410] - symbol * a[idx440];
404  dif2 = aw[idx410] - symbol * aw[idx440];
405  adg440 = (dif1 - dif2) / denom;
406 
407  for (i = 0; i < nbands; i++) {
408  adg[i] = adg440 * exp(Sr * (wavel[idx440] - wavel[i]));
409  aph[i] = a[i] - adg[i] - aw[i];
410  }
411 
412  /* check aph443 range */
413 
414  if (aph_check) {
415 
416  x1 = aph[idx440] / a[idx440];
417  if (x1 < 0.15 || x1 > 0.6) {
418 
419  if (flags) *flags |= 0x10;
420 
421  x2 = -0.8 + 1.4 * (a[idx440] - aw[idx440]) / (a[idx410] - aw[idx410]);
422  if (x2 < 0.15) {
423  x2 = 0.15;
424  if (flags) *flags |= 0x20;
425  }
426  if (x2 > 0.6) {
427  x2 = 0.6;
428  if (flags) *flags |= 0x40;
429  }
430 
431  aph[idx440] = a[idx440] * x2;
432  adg440 = a[idx440] - aph[idx440] - aw[idx440];
433 
434  for (i = 0; i < nbands; i++) {
435  adg[i] = adg440 * exp(Sr * (wavel[idx440] - wavel[i]));
436  aph[i] = a[i] - adg[i] - aw[i];
437  }
438 
439  }
440  }
441  return 0;
442 }
443 
449 int
450 qaaf_decomp(int nbands, float *wavel, float *rrs, float *a, float *aw,
451  float *adg, float *aph, unsigned char *flags) {
452  int i;
453  float symbol, x1, x2;
454  float zeta, denom, dif1, dif2;
455  float rat, adg440;
456  float Sr;
457 
458  assert(idx410 >= 0);
459  assert(idx440 >= 0);
460  assert(idx555 >= 0);
461  assert(nbands >= 0);
462 
463  /* step 7 */
464  rat = rrs[idx440] / rrs[idx555];
465  symbol = 0.74 + (0.2 / (0.8 + rat));
466 
467  /* step 8 */
468  Sr = S + 0.002 / (0.6 + rat);
469  zeta = expf(Sr * (wavel[idx440] - wavel[idx410]));
470 
471  /* step 9 */
472  denom = zeta - symbol;
473  dif1 = a[idx410] - symbol * a[idx440];
474  dif2 = aw[idx410] - symbol * aw[idx440];
475  adg440 = (dif1 - dif2) / denom;
476 
477  for (i = 0; i < nbands; i++) {
478  adg[i] = adg440 * expf(Sr * (wavel[idx440] - wavel[i]));
479  aph[i] = a[i] - adg[i] - aw[i];
480  }
481 
482  /* check aph443 range */
483 
484  if (aph_check) {
485 
486  x1 = aph[idx440] / a[idx440];
487  if (x1 < 0.15 || x1 > 0.6) {
488 
489  if (flags) *flags |= 0x10;
490 
491  x2 = -0.8 + 1.4 * (a[idx440] - aw[idx440]) / (a[idx410] - aw[idx410]);
492  if (x2 < 0.15) {
493  x2 = 0.15;
494  if (flags) *flags |= 0x20;
495  }
496  if (x2 > 0.6) {
497  x2 = 0.6;
498  if (flags) *flags |= 0x40;
499  }
500 
501  aph[idx440] = a[idx440] * x2;
502  adg440 = a[idx440] - aph[idx440] - aw[idx440];
503 
504  for (i = 0; i < nbands; i++) {
505  adg[i] = adg440 * expf(Sr * (wavel[idx440] - wavel[i]));
506  aph[i] = a[i] - adg[i] - aw[i];
507  }
508 
509  }
510 
511  }
512 
513  return 0;
514 }
515 
516 #ifdef TEST_QAA
517 
518 #include <stdio.h>
519 
520 /*
521  * to compile:
522  * cc -o qaa -g -DTEST_QAA -I. qaa.c -lm
523  * ./qaa
524  */
525 
526 
527 
528 static void print_out(int n, float *fwl, float *Rrs, float *rrs, float *u,
529  float *a, float *aph, float *adg, float *aw, float *bb, float *bbw) {
530 
531  int i;
532 
533  printf("lamda ");
534  for (i = 0; i < n; i++)
535  printf("%9.0f ", fwl[i]);
536  printf("\n");
537 
538  printf("Rrs : ");
539  for (i = 0; i < n; i++)
540  printf("%9.6f ", Rrs[i]);
541  printf("\n");
542 
543  printf("rrs : ");
544  for (i = 0; i < n; i++)
545  printf("%9.6f ", rrs[i]);
546  printf("\n");
547 
548  printf("u : ");
549  for (i = 0; i < n; i++)
550  printf("%9.6f ", u[i]);
551  printf("\n");
552 
553  printf("a : ");
554  for (i = 0; i < n; i++)
555  printf("%9.6f ", a[i]);
556  printf("\n");
557 
558  printf("aph : ");
559  for (i = 0; i < n; i++)
560  printf("%9.6f ", aph[i]);
561  printf("\n");
562 
563  printf("adg : ");
564  for (i = 0; i < n; i++)
565  printf("%9.6f ", adg[i]);
566  printf("\n");
567 
568  printf("aw : ");
569  for (i = 0; i < n; i++)
570  printf("%9.6f ", aw[i]);
571  printf("\n");
572 
573  printf("bb : ");
574  for (i = 0; i < n; i++)
575  printf("%9.6f ", bb[i]);
576  printf("\n");
577 
578  printf("bbw : ");
579  for (i = 0; i < n; i++)
580  printf("%9.6f ", bbw[i]);
581  printf("\n");
582 
583 }
584 
585 int
586 main(int argc, char *argv[]) {
587 
588  /* the 5th position here will be changed for 640nm later down in code */
589 
590 #define NBANDS 6
591 #define NUM_SPECTRA 24
592  float Rrs_insitu[NUM_SPECTRA][NBANDS] = {
593  /* 412 443 490 510 555 670 */
594  { 0.001919, 0.002297, 0.004420, 0.005547, 0.009138, 0.004110},
595  { 0.001753, 0.002657, 0.004348, 0.005212, 0.007641, 0.003950},
596  { 0.003968, 0.003374, 0.002932, 0.002142, 0.001214, 0.000136},
597  { 0.001332, 0.001414, 0.002116, 0.002464, 0.003891, 0.001669},
598  { 0.001572, 0.001537, 0.002132, 0.002504, 0.004157, 0.001909},
599  { 0.004830, 0.004084, 0.003540, 0.002648, 0.001584, 0.000194},
600  { 0.001145, 0.001457, 0.002501, 0.002945, 0.004180, 0.001524},
601  { 0.000921, 0.001004, 0.001465, 0.001456, 0.001363, 0.000252},
602  { 0.003120, 0.003226, 0.003530, 0.002857, 0.001882, 0.000196},
603  { 0.005170, 0.004682, 0.003794, 0.002535, 0.001343, 0.000131},
604  { 0.004718, 0.004378, 0.003763, 0.002627, 0.001472, 0.000151},
605  { 0.003503, 0.003394, 0.003358, 0.002510, 0.001507, 0.000181},
606  { 0.001005, 0.001131, 0.001879, 0.002219, 0.003316, 0.001100},
607  { 0.007704, 0.006917, 0.005540, 0.003721, 0.002129, 0.000389},
608  { 0.003311, 0.003071, 0.002920, 0.002315, 0.001508, 0.000285},
609  { 0.003476, 0.003285, 0.003073, 0.002347, 0.001436, 0.000200},
610  { 0.004661, 0.005739, 0.008028, 0.007809, 0.006632, 0.001035},
611  { 0.004212, 0.004859, 0.006455, 0.006117, 0.004895, 0.000676},
612  { 0.002090, 0.002280, 0.003091, 0.002780, 0.002177, 0.000444},
613  { 0.003237, 0.003042, 0.002947, 0.002317, 0.001440, 0.000166},
614  { 0.003125, 0.003444, 0.004119, 0.003777, 0.002875, 0.000413},
615  { 0.002637, 0.002896, 0.003477, 0.002987, 0.002085, 0.000251},
616  { 0.003472, 0.003596, 0.003834, 0.003067, 0.001969, 0.000331},
617  { 0.004554, 0.003772, 0.002074, 0.002082, 0.001338, 0.000120}
618  };
619 
620  float fwl[NBANDS] = {412, 443, 490, 510, 555, 670};
621  float aw[NBANDS] = {0.004994, 0.007512, 0.025010, 0.040020, 0.077080, 0.445600};
622  float bbw[NBANDS] = {0.003271, 0.002421, 0.001568, 0.001339, 0.000939, 0.000426};
623 
624  int i, j;
625  float Rrs[NBANDS];
626  float rrs[NBANDS];
627  float u[NBANDS];
628  float a[NBANDS];
629  float bb[NBANDS];
630  float aph[NBANDS];
631  float adg[NBANDS];
632  unsigned char flags;
633 
634  int idx410 = 0;
635  int idx440 = 1;
636  int idx490 = 2;
637  int idx555 = 4;
638  int idx670 = 5;
639  int nbands;
640 
641  // Ping's bbw using 0.0038 * pow((400.0/lambda),4.32);
642  // for ( i = 0; i< NBANDS; i++ )
643  // bbw[i] = 0.0038 * pow(400.0/fwl[i],4.32);
644 
645  qaa_init(idx410, idx440, idx490, idx555, idx670);
647 
648  printf("QAA v6\n");
649  for (j = 0; j < NUM_SPECTRA; j++) {
650 
651  flags = 0;
652  nbands = NBANDS;
653  /* 412 to 670 */
654  for (i = 0; i < nbands; i++)
655  Rrs[i] = Rrs_insitu[j][i];
656 
657  qaaf_v6(nbands, fwl, Rrs, aw, bbw, rrs, u, a, bb, NULL);
658  qaaf_decomp(nbands, fwl, rrs, a, aw, adg, aph, NULL);
659 
660  for (i = 0; i < 6; i++)
661  if (a[i] < aw[i])
662  flags |= 0x08;
663 
664  for (i = 0; i < 6; i++)
665  if (bb[i] < bbw[i])
666  flags |= 0x80;
667 
668  print_out(nbands, fwl, Rrs, rrs, u, a, aph, adg, aw, bb, bbw);
669 
670  if (flags & 0x10)
671  printf("original aph/a ratio was out of range (0.15 to 0.6)\n");
672  if (flags & 0x20)
673  printf(" and was forced to minimum (0.15)\n");
674  if (flags & 0x40)
675  printf(" and was forced to maximum (0.6)\n");
676  printf("\n");
677  }
678 
679  return 0;
680 }
681 #endif
int qaa_init(int i410, int i440, int i490, int i555, int i670)
initalize Quasi-Analytical Algorithm v4
Definition: qaa.c:132
int j
Definition: decode_rs.h:73
@ QAA_S_PARAM
Definition: qaa.h:12
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:99
#define NULL
Definition: decode_rs.h:63
int qaaf_v6(int nbands, float *wavel, float *Rrs, float *aw, float *bbw, float *rrs, float *u, float *a, float *bb, unsigned char *flags)
Quasi-Analytical Algorithm v4.
Definition: qaa.c:279
@ QAA_COEFS_PARAM
Definition: qaa.h:13
int qaa_decomp(int nbands, double *wavel, double *rrs, double *a, double *aw, double *adg, double *aph, unsigned char *flags)
Quasi-Analytical Algorithm - decomposition of total absorption.
Definition: qaa.c:380
int qaa_set_param(int param,...)
set a parameter for Quasi-Analytical Algorithm
Definition: qaa.c:163
int qaa_v6(int nbands, double *wavel, double *Rrs, double *aw, double *bbw, double *rrs, double *u, double *a, double *bb, unsigned char *flags)
Quasi-Analytical Algorithm v6.
Definition: qaa.c:200
flags
Definition: DDAlgorithm.h:22
int16_t * numer[MAXNFILES]
Definition: l2bin.cpp:99
int32_t nbands
data_t u
Definition: decode_rs.h:74
int main(int argc, char *argv[])
Definition: afrt_nc4.cpp:30
@ QAA_APH_CHECK
Definition: qaa.h:14
int qaa_is_initialized(void)
determine if qaa algorithm properly initialized
Definition: qaa.c:113
int i
Definition: decode_rs.h:71
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
int qaaf_decomp(int nbands, float *wavel, float *rrs, float *a, float *aw, float *adg, float *aph, unsigned char *flags)
Quasi-Analytical Algorithm - decomposition of total absorption.
Definition: qaa.c:450
@ NBANDS
Definition: make_L3_v1.1.c:53