OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
las_iop_ksm.c
Go to the documentation of this file.
1 /* --------------------------------------------------------------------------------------- */
2 /* Loisel & Stramski IOP model */
3 /* */
4 /* References: */
5 /* */
6 /* H. Loisel, J.M. Nicolas, A. Sciandra, D. Stramski, and A. Poteau. 2006. The */
7 /* spectral dependency of optical backscattering by marine particles from satellite */
8 /* remote sensing of the global ocean. Journal of Geophysical Research, 111, C09024, */
9 /* doi: 10.1029/2005JC003367. */
10 /* */
11 /* H. Loisel, D. Stramski, B. G. Mitchell, F. Fell, V. Fournier-Sicre, B. Lemasle */
12 /* and M. Babin. 2001. Comparison of the ocean inherent optical properties obtained */
13 /* from measurements and inverse modeling. Applied Optics.40: 2384-2397. */
14 /* */
15 /* H. Loisel, E. Bosc D. Stramski, K. Oubelker and P-Y. Deschamps. 2001. Seasonal */
16 /* variability of the backscattering coefficients in the Mediterranean Sea based */
17 /* on Satellite SeaWIFS imagery. Geophysical Research letters. 28: 4203-4206. */
18 /* */
19 /* H. Loisel, J.M. Nicolas, P.Y.Deschamps, and R. Frouin. 2002. Seasonal and */
20 /* inter-annual variability of the particulate matter in the global ocean. */
21 /* Geophysical Research letters 29: 2996. */
22 /* */
23 /* H. Loisel and D. Stramski. 2000. Estimation of the inherent optical properties */
24 /* of natural waters from irradiance attenuation coefficient and reflectance in */
25 /* the presence of Raman scattering. Applied Optics. 39: 3001-3011. */
26 /* */
27 /* Fortran code provided by H. Loisel, May 2008. */
28 /* */
29 /* Implementation: */
30 /* */
31 /* B. Franz, NASA/OBPG, June 2008. */
32 /* --------------------------------------------------------------------------------------- */
33 
34 #include <stdlib.h>
35 #include <math.h>
36 #include "l12_proto.h"
37 #include "l2prod.h"
38 #include "amoeba.h"
39 
40 #define LASNSOL 3
41 #define LASNKC 3
42 #define LASNRC 2
43 
44 typedef float r_array[LASNSOL][LASNRC];
45 typedef float k_array[LASNSOL][LASNKC];
46 
47 typedef struct las_table_struc {
48  int nwav;
49  int nsol;
50  int nrc;
51  int nkc;
52  float *wav;
53  float sol[LASNSOL];
54  r_array *rc;
55  k_array *kc;
56 } lastabstr;
57 
58 static lastabstr tab; // R(0-) and <Kd>1 coefficient table
59 
60 static int LastRecNum = -1;
61 static float badData = BAD_FLT;
62 
63 static int ib490 = -1; // Index of 490 wavelength //RJH
64 static int nwave = -1; // number of wavelengths to fit
65 static int nwave_ksm = -1; // number of wavelengths to fit //RJH
66 static float *wave; // sensor wavelengths to fit
67 static float *wave_ksm; // sensor wavelengths to fit //RJH
68 static float *a; // total absorption coefficient per band and pixel
69 static float *b; // total scattering coefficient per band and pixel
70 static float *c; // beam attenuation coefficient per band and pixel
71 static float *bb; // backscatter coefficient per band and pixel
72 static float *bbp; // particulate backscatter coefficient per band and pixel
73 static float *aw; // pure-water total absorption per band
74 static float *bbw; // pure-water backscattering per band
75 static float *bbw_ksm; //RJH
76 
77 
78 /* ------------------------------------------------------------------------------- */
79 /* read_las_tables() - load Loisel & Stramski R0 and Kd coefficient tables */
80 
81 /* ------------------------------------------------------------------------------- */
82 void read_las_tables(int sensorID) {
83  FILE *fp;
84  char *tmp_str;
85  char fname[FILENAME_MAX];
86  char sensor_name[FILENAME_MAX];
87  int iwav, isol;
88  static int firstCall = 1;
89 
90 
91  tab.nsol = LASNSOL;
92  tab.nkc = LASNKC;
93  tab.nrc = LASNRC;
94 
95  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
96  printf("OCDATAROOT environment variable is not defined.\n");
97  exit(1);
98  }
99 
100  /* sensor specific LUTs used - only available SeaWiFS */
101  switch (sensorID) {
102  case SEAWIFS:
103  sprintf(sensor_name, "seawifs");
104  tab.nwav = 5;
105  break;
106  default: /* interpolate from seawifs wavelengths */
107  sprintf(sensor_name, "seawifs");
108  tab.nwav = 5;
109  break;
110  }
111 
112  if (tab.rc == NULL) {
113  if ((tab.rc = (r_array *) calloc(tab.nwav, sizeof (r_array))) == NULL) {
114  printf("-E- : Error allocating memory to read_las_tables\n");
115  exit(FATAL_ERROR);
116  }
117  if ((tab.kc = (k_array *) calloc(tab.nwav, sizeof (k_array))) == NULL) {
118  printf("-E- : Error allocating memory to read_las_tables\n");
119  exit(FATAL_ERROR);
120  }
121  if ((tab.wav = (float *) calloc(tab.nwav, sizeof (float))) == NULL) {
122  printf("-E- : Error allocating memory to read_las_tables\n");
123  exit(FATAL_ERROR);
124  }
125  }
126  printf("Loading Loisel & Stramski IOP model tables:\n");
127 
128  sprintf(fname, "%s/%s/iop/las/%s", tmp_str, sensor_name, "LUT_RRS");
129  printf(" %s\n", fname);
130  if ((fp = fopen(fname, "r")) == NULL) {
131  printf("Error opening %s\n", fname);
132  exit(1);
133  }
134  for (isol = 0; isol < tab.nsol; isol++) for (iwav = 0; iwav < tab.nwav; iwav++)
135  fscanf(fp, "%f %f %f %f", &tab.wav[iwav], &tab.sol[isol], &tab.rc[iwav][isol][0], &tab.rc[iwav][isol][1]);
136  fclose(fp);
137 
138  sprintf(fname, "%s/%s/iop/las/%s", tmp_str, sensor_name, "LUT_KD");
139  printf(" %s\n", fname);
140  if ((fp = fopen(fname, "r")) == NULL) {
141  printf("Error opening %s\n", fname);
142  exit(1);
143  }
144  for (isol = 0; isol < tab.nsol; isol++) for (iwav = 0; iwav < tab.nwav; iwav++)
145  fscanf(fp, "%f %f %f %f %f", &tab.wav[iwav],
146  &tab.sol[isol], &tab.kc[iwav][isol][0], &tab.kc[iwav][isol][1], &tab.kc[iwav][isol][2]);
147  fclose(fp);
148 }
149 
150 
151 /* ------------------------------------------------------------------------------- */
152 /* computes R(0-) and <Kd>1 at one wavelength from Rrs and Rrs(443)/Rrs(555) */
153 
154 /* ------------------------------------------------------------------------------- */
155 void get_R0_Kd(float rat, float Rrs, float wave, float solz, float *R0, float *Kd) {
156  float log_q = log10(rat);
157  int iwav = windex(wave, tab.wav, tab.nwav);
158  float dwav = wave - tab.wav[iwav];
159 
160  float R01, R02;
161  float Kd1, Kd2;
162 
163  if (fabs(dwav) > 1.0) {
164 
165  // if the input wavelength is not one of the table wavelengths, call the
166  // function recursively for bounding table wavelengths and interpolate
167 
168  int iwav1, iwav2;
169  if (dwav > 0.0) {
170  iwav2 = MIN(iwav + 1, tab.nwav - 1);
171  iwav1 = iwav2 - 1;
172  } else {
173  iwav1 = MAX(iwav - 1, 0);
174  iwav2 = iwav1 + 1;
175  }
176  get_R0_Kd(rat, Rrs, tab.wav[iwav1], solz, &R01, &Kd1);
177  get_R0_Kd(rat, Rrs, tab.wav[iwav2], solz, &R02, &Kd2);
178  *R0 = R01 + (wave - tab.wav[iwav1]) / (tab.wav[iwav2] - tab.wav[iwav1])*(R02 - R01);
179  *Kd = Kd1 + (wave - tab.wav[iwav1]) / (tab.wav[iwav2] - tab.wav[iwav1])*(Kd2 - Kd1);
180 
181  } else {
182 
183  // compute R(0-) and <Kd>1 for bounding solar zenith angles of the
184  // coefficient table and interpolate to input solar zenith angle
185 
186  int isol1, isol2;
187  float solz1, solz2;
188 
189  isol2 = 0;
190  while (tab.sol[isol2] < solz) isol2++;
191  isol2 = MAX(MIN(isol2, tab.nsol - 1), 1);
192  isol1 = isol2 - 1;
193 
194  solz1 = tab.sol[isol1];
195  solz2 = tab.sol[isol2];
196 
197  R01 = tab.rc[iwav][isol1][0] * pow(Rrs, tab.rc[iwav][isol1][1]);
198  R02 = tab.rc[iwav][isol2][0] * pow(Rrs, tab.rc[iwav][isol2][1]);
199 
200  Kd1 = pow(10.0, (tab.kc[iwav][isol1][0] * log_q + tab.kc[iwav][isol1][1]) / (tab.kc[iwav][isol1][2] + log_q));
201  Kd2 = pow(10.0, (tab.kc[iwav][isol2][0] * log_q + tab.kc[iwav][isol2][1]) / (tab.kc[iwav][isol2][2] + log_q));
202 
203  *R0 = R01 + (solz - solz1)*(R02 - R01) / (solz2 - solz1);
204  *Kd = Kd1 + (solz - solz1)*(Kd2 - Kd1) / (solz2 - solz1);
205  }
206 
207 }
208 
209 
210 /* ------------------------------------------------------------------------------- */
211 /* have we run for this scan line? */
212 
213 /* ------------------------------------------------------------------------------- */
214 int las_ran(int recnum) {
215  if (recnum == LastRecNum)
216  return 1;
217  else
218  return 0;
219 }
220 
221 /* ------------------------------------------------------------------------------- */
222 /* allocates private arrays */
223 
224 /* ------------------------------------------------------------------------------- */
225 void alloc_las(int nbands) {
226  a = (float*) calloc(nbands, sizeof (float));
227  b = (float*) calloc(nbands, sizeof (float));
228  c = (float*) calloc(nbands, sizeof (float));
229  bb = (float*) calloc(nbands, sizeof (float));
230  bbp = (float*) calloc(nbands, sizeof (float));
231  aw = (float*) calloc(nbands, sizeof (float));
232  bbw = (float*) calloc(nbands, sizeof (float));
233  bbw_ksm = (float*) calloc(nbands, sizeof (float)); //RJH
234  wave_ksm = (float*) calloc(nbands, sizeof (float)); //RJH
235 }
236 
237 
238 /* ------------------------------------------------------------------------------- */
239 /* computes eta = bw/b */
240 
241 /* ------------------------------------------------------------------------------- */
242 float eta_func(float bbw, float bb) {
243  float bw = bbw * 2.0;
244  float bbp, b;
245 
246  if (bb <= bbw)
247  bbp = bbw / 10.0;
248  else
249  bbp = bb - bbw;
250 
251  b = bw + bbp / 0.0183;
252 
253  return (bw / b);
254 }
255 
256 
257 /* ------------------------------------------------------------------------------- */
258 /* computes alpha exponent of bb function */
259 
260 /* ------------------------------------------------------------------------------- */
261 float alpha_func(float solz, float eta) {
262  float a = -0.0007 * solz + 0.1024;
263  float b = -0.0042 * solz + 0.3594;
264  return (a * log10(eta) + b);
265 }
266 
267 
268 /* ------------------------------------------------------------------------------- */
269 /* computes delta exponent of bb function */
270 
271 /* ------------------------------------------------------------------------------- */
272 float delta_func(float eta) {
273  return (0.871 + 0.40 * eta - 1.83 * pow(eta, 2.0));
274 }
275 
276 /* ------------------------------------------------------------------------------- */
277 /* run the model, collect the output */
278 
279 /* ------------------------------------------------------------------------------- */
280 void run_las(l2str *l2rec, int ip) {
281  static int firstCall = 1;
282  static int ib443 = -1;
283  static int ib555 = -1;
284  static int maxiter = 4;
285  static float nw = 1.334;
286  float *Rrs, ratio;
287  float R0, Kd;
288  float solz, mu;
289  float alpha, delta, eta;
290  float h, m, minh, maxh, ec50;
291  int status;
292  int ib, it;
293 
294  l1str *l1rec = l2rec->l1rec;
295  int32_t nbands = l1rec->l1file->nbands;
296 
297  if (firstCall) {
298 
299  firstCall = 0;
300 
301  alloc_las(nbands);
302  wave = l1rec->l1file->fwave;
303  ib443 = windex(443, wave, nbands);
304  ib490 = windex(490, wave, nbands);
305  ib555 = windex(555, wave, nbands);
306  nwave = ib555 + 1;
307  nwave_ksm = ib555 - ib490 + 1; //RJH
308 
309 
310  read_las_tables(l1rec->l1file->sensorID);
311 
312  get_aw_bbw(l2rec, wave, nwave, aw, bbw); //RJH committed
313 
314  for (ib = 0; ib < nwave_ksm; ib++) { //RJH
315  bbw_ksm[ib] = bbw[ib + ib490];
316  wave_ksm[ib] = wave[ib + ib490];
317  }
318  }
319 
320  if ((Rrs = (float *) calloc(nwave, sizeof (float))) == NULL) {
321  printf("-E- : Error allocating memory to run_las\n");
322  exit(FATAL_ERROR);
323  }
324 
325  status = 0;
326 
327  // check for negative radiances
328 
329  for (ib = 0; ib < nwave; ib++) {
330  Rrs[ib] = l2rec->Rrs[ip * nbands + ib];
331  if (Rrs[ib] <= 0.0) status = 1;
332  }
333 
334  // if no negatives and pixel not already masked, run model
335 
336  if (status == 0 && !l1rec->mask[ip]) {
337 
338  if ((input->brdf_opt & FOQMOREL) != 0)
339  solz = 0.0;
340  else
341  solz = MIN(l1rec->solz[ip], 60.0);
342 
343  ratio = Rrs[ib443] / Rrs[ib555];
344  mu = cos(asin(sin(solz / RADEG) / nw));
345 
346  for (ib = 0; ib < nwave; ib++) {
347 
348  get_R0_Kd(ratio, Rrs[ib], wave[ib], solz, &R0, &Kd);
349 
350  eta = 0.05;
351  for (it = 0; it < maxiter; it++) {
352  alpha = alpha_func(solz, eta);
353  delta = delta_func(eta);
354  bb[ib] = Kd * pow(10, alpha) * pow(R0, delta);
355  eta = eta_func(bbw[ib], bb[ib]);
356  }
357 
358  minh = -0.0034 * solz + 1.248;
359  maxh = -0.0084 * solz + 1.7223;
360  ec50 = 0.0033 * solz - 1.9837;
361 
362  m = 0.0034 * pow(solz, 2) - 0.0674 * solz - 9.34;
363  h = pow(10.0, minh + (maxh - minh) / (1 + pow(10, (ec50 - log10(eta)) * m)));
364 
365  a[ib] = mu * Kd / sqrt(1 + h * R0 / (1 - R0));
366  b[ib] = bbw[ib]*2.0 / eta;
367  c[ib] = a[ib] + b[ib];
368  bbp[ib] = bb[ib] - bbw[ib];
369 
370  if (!isfinite(a[ib]) || !isfinite(bb[ib])) {
371  a [ib] = badData;
372  b [ib] = badData;
373  c [ib] = badData;
374  bb [ib] = badData;
375  bbp[ib] = badData;
376  }
377  }
378 
379  } else {
380 
381  for (ib = 0; ib < nwave; ib++) {
382  a [ib] = badData;
383  b [ib] = badData;
384  c [ib] = badData;
385  bb [ib] = badData;
386  bbp[ib] = badData;
387  }
388  LastRecNum = l1rec->iscan;
389  free(Rrs);
390  return;
391  }
392 
393  LastRecNum = l1rec->iscan;
394  free(Rrs);
395  return;
396 }
397 
398 
399 /* ------------------------------------------------------------------- */
400 /* */
401 
402 /* ------------------------------------------------------------------- */
403 double las_eta_amb(FITSTRUCT *ambdata, double par[]) {
404  int iw;
405 
406  ambdata->merit = 0.0;
407  for (iw = 0; iw < nwave; iw++) {
408  ambdata->yfit[iw] = par[0] * pow(wave[0] / wave[iw], par[1]);
409  ambdata->merit += pow((ambdata->y[iw] - ambdata->yfit[iw]), 2) * ambdata->wgt[iw];
410  }
411 
412  return (ambdata->merit);
413 }
414 
415 double las_eta_amb_ksm(FITSTRUCT *ambdata, double par[]) {
416  int iw;
417 
418  ambdata->merit = 0.0;
419  for (iw = 0; iw < nwave_ksm; iw++) {
420  ambdata->yfit[iw] = par[0] * pow(wave_ksm[0] / wave_ksm[iw], par[1]); //RJH
421  ambdata->merit += pow((ambdata->y[iw] - ambdata->yfit[iw]), 2) * ambdata->wgt[iw];
422  }
423 
424  return (ambdata->merit);
425 }
426 
427 
428 
429 /* ------------------------------------------------------------------- */
430 /* */
431 /* ------------------------------------------------------------------- */
432 #define NPARETALAS 2
433 
434 float fit_las_eta_amb(float *bbp) {
435  static int firstCall = 1;
436  static float tol = 1.e-6; /* fractional change in chisqr */
437  static FITSTRUCT ambdata; /* amoeba interface structure */
438  static double init[NPARETALAS * (NPARETALAS + 1)]; /* initial simplex */
439  static double par [NPARETALAS];
440  static double len [NPARETALAS];
441  static double *wts;
442  static double *y;
443  static double *yfit;
444  static int npar = NPARETALAS;
445  static int maxiter = 100;
446 
447  int i, j;
448  short isml;
449  int status = 0;
450 
451  if (firstCall) {
452  firstCall = 0;
453  if ((wts = (double *) calloc(nwave, sizeof (double))) == NULL) {
454  printf("-E- : Error allocating memory to fit_las_eta_amb\n");
455  exit(FATAL_ERROR);
456  }
457  if ((y = (double *) calloc(nwave, sizeof (double))) == NULL) {
458  printf("-E- : Error allocating memory to fit_las_eta_amb\n");
459  exit(FATAL_ERROR);
460  }
461  if ((yfit = (double *) calloc(nwave, sizeof (double))) == NULL) {
462  printf("-E- : Error allocating memory to fit_las_eta_amb\n");
463  exit(FATAL_ERROR);
464  }
465  for (i = 0; i < nwave; i++) {
466  wts[i] = 1.0;
467  }
468  ambdata.nfunc = npar; /* number of model parameters */
469  ambdata.npnts = nwave; /* number of wavelengths */
470  ambdata.wgt = wts; /* Input weights on values */
471  ambdata.meta = NULL;
472 
473  }
474  ambdata.niter = maxiter;
475 
476  for (i = 0; i < nwave; i++) {
477  y[i] = (double) bbp[i];
478  }
479 
480  ambdata.y = y; /* Input values */
481  ambdata.yfit = yfit; /* Output model predicted values */
482 
483  par[0] = y[0];
484  len[0] = 0.01;
485  par[1] = 0.0;
486  len[1] = 2.0;
487 
488  /* initialize simplex with first guess model parameters */
489  for (j = 0; j < npar + 1; j++)
490  for (i = 0; i < npar; i++)
491  init[j * npar + i] = par[i];
492 
493  for (i = 0; i < npar; i++) {
494  init[npar + i * (npar + 1)] += len[i];
495  par[i] = 0.0;
496  }
497 
498  /* run optimization */
499  isml = amoeba(init, &ambdata, las_eta_amb, tol);
500 
501  for (i = 0; i < npar; i++) {
502  par[i] = init[npar * isml + i];
503  }
504 
505  /* check convergence and record parameter results */
506  if (ambdata.niter >= maxiter || par[0] == BAD_FLT || !isfinite(par[1]))
507  return (BAD_FLT);
508  else
509  return (par[1]);
510 }
511 
512 float fit_las_eta_amb_ksm(float *bbp) {
513  static int firstCall = 1;
514  static float tol = 1.e-6; /* fractional change in chisqr */
515  static FITSTRUCT ambdata; /* amoeba interface structure */
516  static double init[NPARETALAS * (NPARETALAS + 1)]; /* initial simplex */
517  static double par [NPARETALAS];
518  static double len [NPARETALAS];
519  static double *wts;
520  static double *y;
521  static double *yfit;
522  static int npar = NPARETALAS;
523  static int maxiter = 100;
524 
525  int i, j;
526  short isml;
527  int status = 0;
528 
529  if (firstCall) {
530  firstCall = 0;
531  if ((wts = (double *) calloc(nwave_ksm, sizeof (double))) == NULL) {
532  printf("-E- : Error allocating memory to fit_las_eta_amb\n");
533  exit(FATAL_ERROR);
534  }
535  if ((y = (double *) calloc(nwave_ksm, sizeof (double))) == NULL) {
536  printf("-E- : Error allocating memory to fit_las_eta_amb\n");
537  exit(FATAL_ERROR);
538  }
539  if ((yfit = (double *) calloc(nwave_ksm, sizeof (double))) == NULL) {
540  printf("-E- : Error allocating memory to fit_las_eta_amb\n");
541  exit(FATAL_ERROR);
542  }
543  for (i = 0; i < nwave_ksm; i++) {
544  wts[i] = 1.0;
545  }
546  ambdata.nfunc = npar; /* number of model parameters */
547  ambdata.npnts = nwave_ksm; /* number of wavelengths */
548  ambdata.wgt = wts; /* Input weights on values */
549  ambdata.meta = NULL;
550 
551  }
552  ambdata.niter = maxiter;
553 
554  for (i = 0; i < nwave_ksm; i++) {
555  y[i] = (double) bbp[i + ib490];
556  }
557 
558  ambdata.y = y; /* Input values */
559  ambdata.yfit = yfit; /* Output model predicted values */
560 
561  par[0] = y[0];
562  len[0] = 0.01;
563  par[1] = 0.0;
564  len[1] = 2.0;
565 
566  /* initialize simplex with first guess model parameters */
567  for (j = 0; j < npar + 1; j++)
568  for (i = 0; i < npar; i++)
569  init[j * npar + i] = par[i];
570 
571  for (i = 0; i < npar; i++) {
572  init[npar + i * (npar + 1)] += len[i];
573  par[i] = 0.0;
574  }
575 
576  /* run optimization */
577  isml = amoeba(init, &ambdata, las_eta_amb_ksm, tol);
578 
579  for (i = 0; i < npar; i++) {
580  par[i] = init[npar * isml + i];
581  }
582 
583  /* check convergence and record parameter results */
584  if (ambdata.niter >= maxiter || par[0] == BAD_FLT || !isfinite(par[1]))
585  return (BAD_FLT);
586  else
587  return (par[1]);
588 }
589 
590 
591 
592 /* ------------------------------------------------------------------------------- */
593 /* interface to l2_hdf_generic() */
594 
595 /* ------------------------------------------------------------------------------- */
596 void get_las(l2str *l2rec, l2prodstr *p, float prod[]) {
597  int prodID = p->cat_ix;
598  int ib = p->prod_ix;
599  int ip;
600 
601  l1str *l1rec = l2rec->l1rec;
602 
603 
604  for (ip = 0; ip < l1rec->npix; ip++) {
605 
606  run_las(l2rec, ip);
607 
608  if (ib >= nwave) {
609  prod[ip] = p->badData;
610  l1rec->flags[ip] |= PRODFAIL;
611  continue;
612  }
613 
614  switch (prodID) {
615 
616  case CAT_a_las:
617  if (a[ib] > badData)
618  prod[ip] = a[ib];
619  else {
620  prod[ip] = p->badData;
621  l1rec->flags[ip] |= PRODFAIL;
622  }
623  break;
624 
625  case CAT_b_las:
626  if (b[ib] > badData)
627  prod[ip] = b[ib];
628  else {
629  prod[ip] = p->badData;
630  l1rec->flags[ip] |= PRODFAIL;
631  }
632  break;
633 
634  case CAT_c_las:
635  if (c[ib] > badData)
636  prod[ip] = c[ib];
637  else {
638  prod[ip] = p->badData;
639  l1rec->flags[ip] |= PRODFAIL;
640  }
641  break;
642 
643  case CAT_bb_las:
644  if (bb[ib] > badData)
645  prod[ip] = bb[ib];
646  else {
647  prod[ip] = p->badData;
648  l1rec->flags[ip] |= PRODFAIL;
649  }
650  break;
651 
652  case CAT_bbp_las:
653  if (bbp[ib] > badData)
654  prod[ip] = bbp[ib];
655  else {
656  prod[ip] = p->badData;
657  l1rec->flags[ip] |= PRODFAIL;
658  }
659  break;
660 
661  case CAT_bbps_las:
662  prod[ip] = get_bbp_las_eta(l2rec, ip);
663  if (prod[ip] == BAD_FLT) {
664  prod[ip] = p->badData;
665  l1rec->flags[ip] |= PRODFAIL;
666  }
667  break;
668 
669  default:
670  printf("-E- %s line %d : erroneous product ID %d passed to Loisel & Stramski model().\n",
671  __FILE__, __LINE__, prodID);
672  exit(1);
673  }
674  }
675 
676  return;
677 }
678 
679 
680 
681 /* ------------------------------------------------------------------------------- */
682 /* interface to convl12() to return iops */
683 
684 /* ------------------------------------------------------------------------------- */
685 void iops_las(l2str *l2rec) {
686  int ib, ip;
687  int32_t nbands = l2rec->l1rec->l1file->nbands; //RJH committed
688 
689  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
690  run_las(l2rec, ip);
691  for (ib = 0; ib < nwave; ib++) {
692  l2rec->a [ip * nbands + ib] = a [ib];
693  l2rec->bb[ip * nbands + ib] = bb[ib];
694  }
695  for (ib = nwave; ib < nbands; ib++) {
696  l2rec->a [ip * nbands + ib] = badData;
697  l2rec->bb[ip * nbands + ib] = badData;
698  }
699  }
700 
701  return;
702 }
703 
704 
705 /* ------------------------------------------------------------------------------- */
706 /* interface to giop() */
707 
708 /* ------------------------------------------------------------------------------- */
709 int get_bbp_las(l2str *l2rec, int ip, float tab_wave[], float tab_bbp[], int tab_nwave) {
710  int iw;
711 
712  run_las(l2rec, ip);
713 
714  for (iw = 0; iw < nwave; iw++) {
715  if (bbp[iw] < 0)
716  return (0);
717  }
718 
719  for (iw = 0; iw < tab_nwave; iw++) {
720  tab_bbp[iw] = linterp(wave, &bbp[0], nwave, tab_wave[iw]);
721  }
722 
723  return (1);
724 }
725 
726 /* ------------------------------------------------------------------------------- */
727 /* interface to giop() */
728 
729 /* ------------------------------------------------------------------------------- */
730 float get_bbp_las_eta(l2str *l2rec, int ip) {
731  int iw;
732  float las;
733 
734  run_las(l2rec, ip);
735 
736  for (iw = 0; iw < nwave; iw++) {
737  if (bbp[iw] < 0)
738  return (BAD_FLT);
739  }
740 
741  las = fit_las_eta_amb(&bbp[0]);
742  return (las);
743 }
744 
745 float get_bbp_las_eta_ksm(l2str *l2rec, int ip) {
746  int iw;
747  float las;
748 
749  run_las(l2rec, ip);
750 
751  for (iw = 0; iw < nwave_ksm; iw++) {
752  if (bbp[iw + ib490] < 0)
753  return (BAD_FLT);
754  }
755 
756  las = fit_las_eta_amb_ksm(&bbp[0]);
757  return (las);
758 }
759 
get_las
void get_las(l2str *l2rec, l2prodstr *p, float prod[])
Definition: las_iop_ksm.c:596
las_table_struc::wav
float * wav
Definition: las_iop.c:52
MAX
#define MAX(A, B)
Definition: swl0_utils.h:26
MIN
#define MIN(x, y)
Definition: rice.h:169
NPARETALAS
#define NPARETALAS
Definition: las_iop_ksm.c:432
j
int j
Definition: decode_rs.h:73
status
int status
Definition: l1_czcs_hdf.c:32
get_bbp_las_eta
float get_bbp_las_eta(l2str *l2rec, int ip)
Definition: las_iop_ksm.c:730
alloc_las
void alloc_las(int nbands)
Definition: las_iop_ksm.c:225
k_array
float k_array[LASNSOL][LASNKC]
Definition: las_iop_ksm.c:45
FITSTRUCT::niter
int niter
Definition: amoeba.h:9
FOQMOREL
#define FOQMOREL
Definition: l12_parms.h:53
LASNSOL
#define LASNSOL
Definition: las_iop_ksm.c:40
NULL
#define NULL
Definition: decode_rs.h:63
las_eta_amb_ksm
double las_eta_amb_ksm(FITSTRUCT *ambdata, double par[])
Definition: las_iop_ksm.c:415
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
las_table_struc
Definition: las_iop.c:47
CAT_bb_las
#define CAT_bb_las
Definition: l2prod.h:189
mu
float mu
Definition: atrem_corl1v2.h:324
h
float h[MODELMAX]
Definition: atrem_corl1.h:131
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
get_bbp_las
int get_bbp_las(l2str *l2rec, int ip, float tab_wave[], float tab_bbp[], int tab_nwave)
Definition: las_iop_ksm.c:709
amoeba
short amoeba(double *pnts, FITSTRUCT *auxdata, double(*func)(FITSTRUCT *, double[]), float tol)
Definition: amoeba.c:14
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
las_ran
int las_ran(int recnum)
Definition: las_iop_ksm.c:214
FITSTRUCT::merit
double merit
Definition: amoeba.h:4
input
instr * input
Definition: main_l2binmatch.cpp:27
l12_proto.h
get_R0_Kd
void get_R0_Kd(float rat, float Rrs, float wave, float solz, float *R0, float *Kd)
Definition: las_iop_ksm.c:155
LASNKC
#define LASNKC
Definition: las_iop_ksm.c:41
CAT_a_las
#define CAT_a_las
Definition: l2prod.h:186
get_bbp_las_eta_ksm
float get_bbp_las_eta_ksm(l2str *l2rec, int ip)
Definition: las_iop_ksm.c:745
init
int init(int32_t ipr, int32_t jpr, char *efile, char *pfile)
Definition: proj_report.c:51
las_table_struc::kc
k_array * kc
Definition: las_iop.c:55
fit_las_eta_amb_ksm
float fit_las_eta_amb_ksm(float *bbp)
Definition: las_iop_ksm.c:512
recnum
read recnum
Definition: HOWTO_Add_a_sensor.txt:72
linterp
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
run_las
void run_las(l2str *l2rec, int ip)
Definition: las_iop_ksm.c:280
alpha_func
float alpha_func(float solz, float eta)
Definition: las_iop_ksm.c:261
fit_las_eta_amb
float fit_las_eta_amb(float *bbp)
Definition: las_iop_ksm.c:434
FITSTRUCT::wgt
double * wgt
Definition: amoeba.h:6
las_table_struc::rc
r_array * rc
Definition: las_iop.c:54
color_dtdb.y
y
Definition: color_dtdb.py:430
LASNRC
#define LASNRC
Definition: las_iop_ksm.c:42
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
tol
float tol
Definition: pml_iop_tables.c:45
delta
const double delta
Definition: VcstCmnConsts.h:523
create_images.Rrs
def Rrs
Definition: create_images.py:92
read_las_tables
void read_las_tables(int sensorID)
Definition: las_iop_ksm.c:82
las_table_struc::sol
float sol[LASNSOL]
Definition: las_iop.c:53
get_aw_bbw
void get_aw_bbw(l2str *l2rec, float wave[], int nwave, float *aw, float *bbw)
Definition: water_spectra.c:179
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
r_array
float r_array[LASNSOL][LASNRC]
Definition: las_iop_ksm.c:44
las_table_struc::nkc
int nkc
Definition: las_iop.c:51
FITSTRUCT::y
double * y
Definition: amoeba.h:5
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
eta_func
float eta_func(float bbw, float bb)
Definition: las_iop_ksm.c:242
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
FITSTRUCT::nfunc
short int nfunc
Definition: amoeba.h:2
las_table_struc::nrc
int nrc
Definition: las_iop.c:50
FITSTRUCT::yfit
double * yfit
Definition: amoeba.h:8
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
las_table_struc::nwav
int nwav
Definition: las_iop.c:48
las_table_struc::nsol
int nsol
Definition: las_iop.c:49
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
fabs
#define fabs(a)
Definition: misc.h:93
FITSTRUCT::meta
void * meta
Definition: amoeba.h:10
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
CAT_bbp_las
#define CAT_bbp_las
Definition: l2prod.h:190
delta_func
float delta_func(float eta)
Definition: las_iop_ksm.c:272
FITSTRUCT::npnts
short int npnts
Definition: amoeba.h:3
color_dtdb.alpha
alpha
Definition: color_dtdb.py:352
CAT_bbps_las
#define CAT_bbps_las
Definition: l2prod.h:216
c
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
CAT_b_las
#define CAT_b_las
Definition: l2prod.h:187
CAT_c_las
#define CAT_c_las
Definition: l2prod.h:188
SEAWIFS
#define SEAWIFS
Definition: sensorDefs.h:12
i
int i
Definition: decode_rs.h:71
iops_las
void iops_las(l2str *l2rec)
Definition: las_iop_ksm.c:685
a
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
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:97
badData
Every product should define the units l2prod badData
Definition: HOWTO_Add_a_product.txt:98
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
las_eta_amb
double las_eta_amb(FITSTRUCT *ambdata, double par[])
Definition: las_iop_ksm.c:403


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster