OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
giop.c
Go to the documentation of this file.
1 /* =================================================================== */
2 /* module giop.c: generic iop model */
3 /* */
4 /* (May 2003) */
5 /* The module contains the functions to optimize and evaluate a user */
6 /* controlled IOP model, which currently defaults to the GSM 2002. */
7 /* */
8 /* (May 2015) */
9 /* Raman scattering correction now applied to Rrs */
10 /* */
11 /* (Dec 2015) */
12 /* Includes iterative adaptive SIOP matrix inversion that iterates */
13 /* through a number of IOP spectral shapes. */
14 /* */
15 /* Implementation: */
16 /* B. Franz, NASA/OBPG/SAIC, May 2008 */
17 /* */
18 /* Reference: */
19 /* Werdell et al. (2013) Generalized ocean color inversion model for */
20 /* retrieving marine inherent optical properties, Applied Optics, 52, */
21 /* 2019-2037. */
22 /* */
23 /* Notes: */
24 /* This code was written with the intent that it may become a base for */
25 /* multiple algorithms (to be writtenn as wrappers). As such, it can't*/
26 /* be assumed that the starting config is the only config (e.g., */
27 /* number of wavelengths to fit), which leads to some inefficiencies. */
28 /* =================================================================== */
29 
30 #include <stdio.h>
31 #include <math.h>
32 #include "l12_proto.h"
33 #include "giop.h"
34 #include "amoeba.h"
35 #include <gsl/gsl_errno.h>
36 #include <gsl/gsl_vector.h>
37 #include <gsl/gsl_blas.h>
38 #include <gsl/gsl_multifit_nlin.h>
39 #include <gsl/gsl_linalg.h>
40 #include <gsl/gsl_poly.h>
41 
42 typedef struct fit_data_str {
43  double *y;
44  double *w;
45  giopstr *g;
46 } fitstr;
47 
48 static int32_t LastRecNum = -1;
49 static float badval = BAD_FLT;
50 static float bbp_s_max = 5.0;
51 static float bbp_s_min = 0.0;
52 static float adg_s_max = 0.025;
53 static float adg_s_min = 0.01;
54 static float chl_max = 10.0;
55 static float chl_min = 0.03;
56 
57 static float *aw;
58 static float *bbw;
59 static float *foq;
60 
61 static float *aph1;
62 static float *adg1;
63 static float *bbp1;
64 
65 static float *anap1;
66 static float *acdom1;
67 static float *bbph1;
68 static float *bbnap1;
69 
70 
71 // these need to be passed from run_giop if giop is to
72 // support simultaneous products (e.g., a gsm wrapper)
73 
74 static int16 *iter;
75 static int16 *iopf;
76 static float *a;
77 static float *bb;
78 static float *aph;
79 static float *adg;
80 static float *bbp;
81 
82 static float *anap;
83 static float *acdom;
84 static float *bbph;
85 static float *bbnap;
86 
87 static float *chl;
88 static float *aph_s;
89 static float *adg_s;
90 static float *bbp_s;
91 static float *rrsdiff;
92 static float *mRrs;
93 static float **fit_par;
94 static float *chisqr;
95 static int16 max_npar;
96 
97 static int *siop_num;
98 static int allocateRrsRaman = 0;
99 
100 void freeArray(void **a, int32_t m) {
101  int i;
102  for (i = 0; i < m; ++i) {
103  free(a[i]);
104  }
105  free(a);
106 }
107 
108 void freeDArray(double **a, int32_t m) {
109  int i;
110  for (i = 0; i < m; ++i) {
111  free(a[i]);
112  }
113  free(a);
114 }
115 
116 /* ------------------------------------------------------------------- */
117 /* giop_int_tab_file - reads a table of eigenvectors and interpolates */
118 /* to input wavelengths and returns a 2-D array of vectors, */
119 
120 /* ------------------------------------------------------------------- */
121 void giop_int_tab_file(char *file, int nx, float *x, float **y) {
122 
123  float *table;
124  float *xtab;
125  float *ytab;
126  int ncol, nrow;
127  int i, ivec;
128 
129  printf("\nLoading default basis vector from %s\n", file);
130 
131 
132  nrow = table_row_count(file);
133  if (nrow <= 0) {
134  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, file);
135  exit(1);
136  }
137 
138  ncol = table_column_count(file);
139  table = table_read_r4(file, ncol, nrow);
140 
141  // Interpolate to input wavelengths (x)
142 
143  xtab = &table[nrow * 0];
144 
145  for (ivec = 0; ivec < ncol - 1; ivec++) {
146  ytab = &table[nrow * (ivec + 1)];
147  for (i = 0; i < nx; i++) {
148  y[ivec][i] = linterp(xtab, ytab, nrow, x[i]);
149  }
150  }
151 
152  table_free_r4(table);
153 }
154 
155 
156 /* ------------------------------------------------------------------- */
157 /* giop_ctl_start - set starting values for optimization */
158 
159 /* ------------------------------------------------------------------- */
160 void giop_ctl_start(giopstr *g, float chl) {
161 
162  int ipar, ivec;
163 
164  // set starting params based on model elements
165 
166  ipar = 0;
167 
168  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
169  g->par[ipar] = chl / g->aph_nvec;
170  g->len[ipar] = 0.5;
171  ipar++;
172  }
173 
174  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
175  g->par[ipar] = 0.01;
176  g->len[ipar] = 0.1;
177  ipar++;
178  }
179 
180  switch (g->bbp_opt) {
181  case BBPLASFIX:
182  case BBPQAAFIX:
183  break;
184  case BBPLAS:
185  g->par[ipar] = 1.0;
186  g->len[ipar] = 0.01;
187  ipar++;
188  break;
189  default:
190  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
191  g->par[ipar] = 0.001;
192  g->len[ipar] = 0.01;
193  ipar++;
194  }
195  break;
196  }
197 
198  //note: for SVDSIOP ipar > npar
199  if (ipar != g->npar && g->fit_opt != SVDSIOP) {
200  printf("-E- %s Line %d: number of GIOP fit parameters (%d) does not equal number of initialized parameter (%d)\n",
201  __FILE__, __LINE__, g->npar, ipar);
202  exit(1);
203  }
204 }
205 
206 /* ------------------------------------------------------------------- */
207 /* giop_ctl_init - initialize model control structure */
208 /* */
209 /* Provides default values or user-supplied parameters for all static */
210 /* model components. These may later be over-riden by dynamic model */
211 /* parameters (e.g., band-ratio based bbp_s). */
212 /* */
213 
214 /* ------------------------------------------------------------------- */
215 void giop_ctl_init(giopstr *g, int nwave, float wave[],
216  float aw[], float bbw[]) {
217 
218  int iw, iwx;
219 
220  float def_aph_w = 443.0;
221  float def_aph_s = 0.5;
222  float def_adg_w = 443.0;
223  float def_adg_s = 0.02061;
224  float def_bbp_w = 443.0;
225  float def_bbp_s = 1.03373;
226  float def_grd[] = {0.0949, 0.0794};
227 
228  // determine indices of wavelengths to be used (default to all vis)
229 
230  if (input->giop_wave[0] > 0) {
231  for (iw = 0; iw < nwave; iw++) {
232  if (input->giop_wave[iw] <= 0) break;
233  iwx = windex(input->giop_wave[iw], wave, nwave);
234  g->bindx[iw] = iwx;
235  //windex(input->giop_wave[iw],wave,nwave);
236  }
237  g->nwave = iw;
238  } else {
239  g->nwave = MIN(windex(671., wave, nwave) + 1, nwave);
240  for (iw = 0; iw < g->nwave; iw++)
241  g->bindx[iw] = iw;
242  }
243 
244  // set wavelengths and pure-water values
245 
246  for (iw = 0; iw < g->nwave; iw++) {
247  g->wave[iw] = wave[g->bindx[iw]];
248  g->aw [iw] = aw [g->bindx[iw]];
249  g->bbw [iw] = bbw [g->bindx[iw]];
250  }
251 
252  // set input Rrs wts based on specified uncertainties
253 
254  if (input->giop_rrs_unc[0] > 0) {
255  g->wt_opt = 1;
256  for (iw = 0; iw < nwave; iw++) {
257  if (input->giop_rrs_unc[iw] <= 0) break;
258  g->wts[iw] = 1.0 / pow(input->giop_rrs_unc[iw], 2);
259  }
260  if (iw != g->nwave) {
261  printf("-E- %s line %d: number of giop_rrs_unc (%d) must equal number of giop_wave (%d)",
262  __FILE__, __LINE__, iw, g->nwave);
263  exit(1);
264  }
265  } else {
266  g->wt_opt = 0;
267  for (iw = 0; iw < g->nwave; iw++)
268  g->wts[iw] = 1.0;
269  }
270 
271 
272  // maximum iterations
273 
274  if (input->giop_maxiter > 0)
275  g->maxiter = input->giop_maxiter;
276  else
277  g->maxiter = 500;
278 
279  // fitting method
280 
281  if (input->giop_fit_opt > 0)
282  g->fit_opt = input->giop_fit_opt;
283  else
284  g->fit_opt = LEVMARQ;
285 
286  // Rrs to bb/(a+bb) method
287 
288  if (input->giop_rrs_opt > 0)
289  g->rrs_opt = input->giop_rrs_opt;
290  else
291  g->rrs_opt = RRSGRD;
292 
293 
294  // set coefficients of Gordon quadratic
295 
296  if (input->giop_grd[0] > -999.0) {
297  g->grd[0] = input->giop_grd[0];
298  g->grd[1] = input->giop_grd[1];
299  } else {
300  g->grd[0] = def_grd[0];
301  g->grd[1] = def_grd[1];
302  }
303 
304  // default basis vectors
305 
306  strcpy(g->aph_tab_file, input->giop_aph_file);
307  strcpy(g->adg_tab_file, input->giop_adg_file);
308  strcpy(g->bbp_tab_file, input->giop_bbp_file);
309  strcpy(g->acdom_tab_file, input->giop_acdom_file);
310  strcpy(g->anap_tab_file, input->giop_anap_file);
311  strcpy(g->bbph_tab_file, input->giop_bbph_file);
312  strcpy(g->bbnap_tab_file, input->giop_bbnap_file);
313 
314  //set defaults
315  g->aph_opt = APHTAB;
316  g->adg_opt = ADGTAB;
317  g->bbp_opt = BBPTAB;
318  g->acdom_opt = ACDOMNONE;
319  g->anap_opt = ANAPNONE;
320  g->bbnap_opt = BBPHNONE;
321  g->bbnap_opt = BBNAPNONE;
322 
323  // aphstar function
324 
325  if (input->giop_aph_opt > 0)
326  g->aph_opt = input->giop_aph_opt;
327 
328  if (input->giop_aph_w > 0.0)
329  g->aph_w = input->giop_aph_w;
330  else
331  g->aph_w = def_aph_w;
332 
333  if (input->giop_aph_s > -999.0)
334  g->aph_s = input->giop_aph_s;
335  else
336  g->aph_s = def_aph_s;
337 
338 
339  // adgstar function
340 
341  //if (input->giop_adg_opt > 0)
342  // g->adg_opt = input->giop_adg_opt;
343  if (input->giop_adg_opt != 1)
344  g->adg_opt = input->giop_adg_opt;
345  else
346  g->adg_opt = ADGS;
347 
348  if (input->giop_adg_w > 0.0)
349  g->adg_w = input->giop_adg_w;
350  else
351  g->adg_w = def_adg_w;
352 
353  if (input->giop_adg_s > -999.0)
354  g->adg_s = input->giop_adg_s;
355  else
356  g->adg_s = def_adg_s;
357 
358  //acdom star function
359  if (input->giop_acdom_opt != 1)
360  g->acdom_opt = input->giop_acdom_opt;
361  else
362  g->adg_opt = ACDOMNONE;
363 
364  //anap star function
365  if (input->giop_anap_opt != 1)
366  g->anap_opt = input->giop_anap_opt;
367  else
368  g->anap_opt = ANAPNONE;
369 
370  // bbpstar function
371 
372  //if (input->giop_bbp_opt > 0)
373  // g->bbp_opt = input->giop_bbp_opt;
374  if (input->giop_bbp_opt != 1)
375  g->bbp_opt = input->giop_bbp_opt;
376  else
377  g->bbp_opt = BBPS;
378 
379  if (input->giop_bbp_w > 0.0)
380  g->bbp_w = input->giop_bbp_w;
381  else
382  g->bbp_w = def_bbp_w;
383 
384  if (input->giop_bbp_s > -999.0)
385  g->bbp_s = input->giop_bbp_s;
386  else
387  g->bbp_s = def_bbp_s;
388 
389  //bbph star function
390  if (input->giop_bbph_opt != 1)
391  g->bbph_opt = input->giop_bbph_opt;
392  else
393  g->bbph_opt = BBPHNONE;
394 
395  //bbnap star function
396  if (input->giop_bbnap_opt != 1)
397  g->bbnap_opt = input->giop_bbnap_opt;
398  else
399  g->bbnap_opt = BBNAPNONE;
400 
401 
402  // set number of vectors
403 
404  switch (g->aph_opt) {
405  case APHTAB:
406  g->aph_nvec = table_column_count(g->aph_tab_file) - 1;
407  break;
408  default:
409  g->aph_nvec = 1;
410  break;
411  }
412 
413  switch (g->adg_opt) {
414  case ADGTAB:
415  g->adg_nvec = table_column_count(g->adg_tab_file) - 1;
416  break;
417  default:
418  g->adg_nvec = 1;
419  break;
420  }
421 
422  switch (g->acdom_opt) {
423  case ACDOMTAB:
424  g->acdom_nvec = table_column_count(g->acdom_tab_file) - 1;
425  break;
426  default:
427  g->acdom_nvec = 1;
428  break;
429  }
430 
431  switch (g->anap_opt) {
432  case ANAPTAB:
433  g->anap_nvec = table_column_count(g->anap_tab_file) - 1;
434  break;
435  default:
436  g->anap_nvec = 1;
437  break;
438  }
439 
440  switch (g->bbp_opt) {
441  case BBPTAB:
442  g->bbp_nvec = table_column_count(g->bbp_tab_file) - 1;
443  break;
444  case BBPLASFIX:
445  case BBPQAAFIX:
446  g->bbp_nvec = 1;
447  break;
448  default:
449  g->bbp_nvec = 1;
450  break;
451  }
452 
453  switch (g->bbph_opt) {
454  case BBPHTAB:
455  g->bbph_nvec = table_column_count(g->bbph_tab_file) - 1;
456  break;
457  default:
458  g->bbph_nvec = 1;
459  break;
460  }
461 
462  switch (g->bbnap_opt) {
463  case BBNAPTAB:
464  g->bbnap_nvec = table_column_count(g->bbnap_tab_file) - 1;
465  break;
466  default:
467  g->bbnap_nvec = 1;
468  break;
469  }
470 
471  // total number of parameters to be optimized
472  if (g->fit_opt == SVDSIOP) {
473  //only fitting for chl, cdom and nap
474  g->npar = 3;
475  } else {
476  g->npar = g->aph_nvec + g->adg_nvec + g->bbp_nvec;
477  }
478 
479 
480  // allocate space for vectors (one element per sensor wavelength)
481 
482  g->aph_tab_nw = nwave;
483  g->aph_tab_w = (float *) calloc(nwave, sizeof (float));
484  g->aph_tab_s = (float **) allocate2d_float(g->aph_tab_nw, g->aph_nvec);
485 
486  g->adg_tab_nw = nwave;
487  g->adg_tab_w = (float *) calloc(nwave, sizeof (float));
488  g->adg_tab_s = (float **) allocate2d_float(g->adg_tab_nw, g->adg_nvec);
489 
490  g->acdom_tab_nw = nwave;
491  g->acdom_tab_w = (float *) calloc(nwave, sizeof (float));
492  g->acdom_tab_s = (float **) allocate2d_float(g->acdom_tab_nw, g->acdom_nvec);
493 
494  g->anap_tab_nw = nwave;
495  g->anap_tab_w = (float *) calloc(nwave, sizeof (float));
496  g->anap_tab_s = (float **) allocate2d_float(g->anap_tab_nw, g->anap_nvec);
497 
498  g->bbp_tab_nw = nwave;
499  g->bbp_tab_w = (float *) calloc(nwave, sizeof (float));
500  g->bbp_tab_s = (float **) allocate2d_float(g->bbp_tab_nw, g->bbp_nvec);
501 
502  g->bbph_tab_nw = nwave;
503  g->bbph_tab_w = (float *) calloc(nwave, sizeof (float));
504  g->bbph_tab_s = (float **) allocate2d_float(g->bbph_tab_nw, g->bbph_nvec);
505 
506  g->bbnap_tab_nw = nwave;
507  g->bbnap_tab_w = (float *) calloc(nwave, sizeof (float));
508  g->bbnap_tab_s = (float **) allocate2d_float(g->bbnap_tab_nw, g->bbnap_nvec);
509 
510  // set vector wavelengths (same as ALL sensor for now)
511 
512  for (iw = 0; iw < nwave; iw++) {
513  g->aph_tab_w[iw] = wave[iw];
514  g->adg_tab_w[iw] = wave[iw];
515  g->acdom_tab_w[iw] = wave[iw];
516  g->anap_tab_w[iw] = wave[iw];
517  g->bbp_tab_w[iw] = wave[iw];
518  g->bbph_tab_w[iw] = wave[iw];
519  g->bbnap_tab_w[iw] = wave[iw];
520  }
521 
522  // load tabulated vectors if requested
523 
524  switch (g->aph_opt) {
525  case APHTAB:
526  giop_int_tab_file(g->aph_tab_file, nwave, wave, g->aph_tab_s);
527  break;
528  }
529 
530  switch (g->adg_opt) {
531  case ADGTAB:
532  giop_int_tab_file(g->adg_tab_file, nwave, wave, g->adg_tab_s);
533  break;
534  }
535 
536  switch (g->bbp_opt) {
537  case BBPTAB:
538  giop_int_tab_file(g->bbp_tab_file, nwave, wave, g->bbp_tab_s);
539  break;
540  }
541  switch (g->acdom_opt) {
542  case ACDOMTAB:
543  giop_int_tab_file(g->acdom_tab_file, nwave, wave, g->acdom_tab_s);
544  break;
545  }
546  switch (g->anap_opt) {
547  case ANAPTAB:
548  giop_int_tab_file(g->anap_tab_file, nwave, wave, g->anap_tab_s);
549  break;
550  }
551  switch (g->bbph_opt) {
552  case BBPHTAB:
553  giop_int_tab_file(g->bbph_tab_file, nwave, wave, g->bbph_tab_s);
554  break;
555  }
556  switch (g->bbnap_opt) {
557  case BBNAPTAB:
558  giop_int_tab_file(g->bbnap_tab_file, nwave, wave, g->bbnap_tab_s);
559  break;
560  }
561 
562 }
563 
564 
565 /* ------------------------------------------------------------------- */
566 /* */
567 
568 /* ------------------------------------------------------------------- */
569 int giop_ran(int recnum) {
570  if (recnum == LastRecNum)
571  return 1;
572  else
573  return 0;
574 }
575 
576 /* ------------------------------------------------------------------- */
577 /* */
578 
579 /* ------------------------------------------------------------------- */
580 void giop_model(giopstr *g, double par[], int nwave, float wave[], float aw[], float bbw[],
581  float foq[], float aph[], float adg[], float bbp[],
582  double rrs[], double **dfdpar, double **parstar) {
583  // double rrs[],double dfdpar[NBANDS][GIOPMAXPAR],double parstar[NBANDS][GIOPMAXPAR])
584  int iw, iwtab, ivec, ipar;
585  float bb, a;
586  float x;
587  float *aphstar;
588  float *adgstar;
589  float *bbpstar;
590  float dfdx;
591  float dxda;
592  float dxdb;
593 
594  // note: input wavelength set can vary between optimization and evaluation calls
595  if ((aphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
596  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
597  __FILE__, __LINE__);
598  exit(1);
599  }
600  if ((adgstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
601  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
602  __FILE__, __LINE__);
603  exit(1);
604  }
605  if ((bbpstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
606  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
607  __FILE__, __LINE__);
608  exit(1);
609  }
610 
611  for (iw = 0; iw < nwave; iw++) {
612 
613  // evaluate model for each IOP component at all wavelengths
614  // also compute and store associated uncertainty, based on fitted paraneter
615  // uncertainties as stored in the second half of the par array
616  // note that we're using a simple summation for uncertainties (for now)
617 
618  ipar = 0;
619 
620  switch (g->aph_opt) {
621  case APHGAUSS:
622  aphstar[ipar] = exp(-1. * pow(wave[iw] - g->aph_w, 2) / (2 * pow(g->aph_s, 2)));
623  aph[iw] = par[ipar] * aphstar[ipar];
624  aph[iw + nwave] = par[ipar + g->npar] * aphstar[ipar];
625  ipar++;
626  break;
627  default:
628  aph[iw] = 0;
629  aph[iw + nwave] = 0;
630  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
631  iwtab = windex(wave[iw], g->aph_tab_w, g->aph_tab_nw);
632  aphstar[ipar] = g->aph_tab_s[ivec][iwtab];
633  aph[iw] += par[ipar] * aphstar[ipar];
634  aph[iw + nwave] += par[ipar + g->npar] * aphstar[ipar];
635  ipar++;
636  }
637  break;
638  }
639 
640  switch (g->adg_opt) {
641  case ADGS:
642  case ADGSQAA:
643  case ADGSOBPG:
644  adgstar[ipar] = exp(-g->adg_s * (wave[iw] - g->adg_w));
645  adg[iw] = par[ipar] * adgstar[ipar];
646  adg[iw + nwave] = par[ipar + g->npar] * adgstar[ipar];
647  ipar++;
648  break;
649  default:
650  adg[iw] = 0;
651  adg[iw + nwave] = 0;
652  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
653  iwtab = windex(wave[iw], g->adg_tab_w, g->adg_tab_nw);
654  adgstar[ipar] = g->adg_tab_s[ivec][iwtab];
655  adg[iw] += par[ipar] * adgstar[ipar];
656  adg[iw + nwave] += par[ipar + g->npar] * adgstar[ipar]; // uncertainty (wag)
657  ipar++;
658  }
659  break;
660  }
661 
662  switch (g->bbp_opt) {
663  case BBPLASFIX:
664  case BBPQAAFIX:
665  iwtab = windex(wave[iw], g->bbp_tab_w, g->bbp_tab_nw);
666  bbpstar[ipar] = g->bbp_tab_s[0][iwtab];
667  bbp[iw] = bbpstar[ipar];
668  break;
669  case BBPS:
670  case BBPSLAS:
671  case BBPSHAL:
672  case BBPSQAA:
673  case BBPSCIOTTI:
674  case BBPSMM01:
675  bbpstar[ipar] = pow((g->bbp_w / wave[iw]), g->bbp_s);
676  bbp[iw] = par[ipar] * bbpstar[ipar];
677  bbp[iw + nwave] = par[ipar + g->npar] * bbpstar[ipar];
678  ipar++;
679  break;
680  default:
681  bbp[iw] = 0;
682  bbp[iw + nwave] = 0;
683  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
684  iwtab = windex(wave[iw], g->bbp_tab_w, g->bbp_tab_nw);
685  bbpstar[ipar] = g->bbp_tab_s[ivec][iwtab];
686  bbp[iw] += par[ipar] * bbpstar[ipar];
687  bbp[iw + nwave] += par[ipar + g->npar] * bbpstar[ipar];
688  ipar++;
689  }
690  break;
691  }
692 
693  a = aw [iw] + aph[iw] + adg[iw];
694  bb = bbw[iw] + bbp[iw];
695  x = bb / (a + bb);
696 
697  switch (g->rrs_opt) {
698  case RRSGRD:
699  rrs[iw] = g->grd[0] * x + g->grd[1] * pow(x, 2);
700  dfdx = g->grd[0] + 2 * g->grd[1] * x;
701  break;
702  case RRSFOQ:
703  rrs[iw] = foq[iw] * x;
704  dfdx = foq[iw];
705  break;
706  }
707 
708  if (dfdpar != NULL) {
709 
710  ipar = 0;
711 
712  dxda = -x * x / bb;
713  dxdb = x / bb + dxda;
714 
715  switch (g->aph_opt) {
716  default:
717  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
718  dfdpar[iw][ipar] = dfdx * dxda * aphstar[ipar];
719  ipar++;
720  }
721  break;
722  }
723 
724  switch (g->adg_opt) {
725  default:
726  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
727  dfdpar[iw][ipar] = dfdx * dxda * adgstar[ipar];
728  ipar++;
729  }
730  break;
731  }
732 
733  switch (g->bbp_opt) {
734  case BBPLASFIX:
735  case BBPQAAFIX:
736  break;
737  default:
738  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
739  dfdpar[iw][ipar] = dfdx * dxdb * bbpstar[ipar];
740  ipar++;
741  }
742  break;
743  }
744  }
745 
746  if (parstar != NULL) {
747 
748  ipar = 0;
749 
750  switch (g->aph_opt) {
751  default:
752  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
753  parstar[iw][ipar] = aphstar[ipar];
754  ipar++;
755  }
756  break;
757  }
758 
759  switch (g->adg_opt) {
760  default:
761  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
762  parstar[iw][ipar] = adgstar[ipar];
763  ipar++;
764  }
765  break;
766  }
767 
768  switch (g->bbp_opt) {
769  case BBPLASFIX:
770  case BBPQAAFIX:
771  break;
772  default:
773  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
774  parstar[iw][ipar] = bbpstar[ipar];
775  ipar++;
776  }
777  break;
778  }
779  }
780 
781  }
782 
783  free(aphstar);
784  free(adgstar);
785  free(bbpstar);
786  return;
787 }
788 /* ------------------------------------------------------------------- */
789 /* */
790 
791 /* ------------------------------------------------------------------- */
792 void giop_model_iterate(giopstr *g, double par[], int nwave, float wave[], float aw[], float bbw[],
793  float foq[], float aph[], float adg[], float bbp[], float acdom[], float anap[], float bbph[], float bbnap[],
794  double rrs[], double **dfdpar, double **parstar) {
795  // double rrs[],double dfdpar[NBANDS][GIOPMAXPAR],double parstar[NBANDS][GIOPMAXPAR])
796  int iw, iwtab, idx443;
797 
798  int16 isiop = g->siopIdx;
799 
800  float bb, a;
801  float x;
802  float acdom443;
803 
804  float *aphstar;
805  float *acdomstar;
806  float *anapstar;
807  float *bbphstar;
808  float *bbnapstar;
809 
810  float dfdx;
811  float dxda;
812  float dxdb;
813 
814 
815  /* note: input wavelength set can vary between optimization and evaluation calls*/
816  if ((aphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
817  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
818  __FILE__, __LINE__);
819  exit(1);
820  }
821  if ((acdomstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
822  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
823  __FILE__, __LINE__);
824  exit(1);
825  }
826  if ((anapstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
827  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
828  __FILE__, __LINE__);
829  exit(1);
830  }
831  if ((bbphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
832  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
833  __FILE__, __LINE__);
834  exit(1);
835  }
836  if ((bbnapstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
837  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
838  __FILE__, __LINE__);
839  exit(1);
840  }
841 
842 
843  /*Sanity check: do the input tables have the same number of columns?*/
844  switch (g->fit_opt) {
845  case SVDSIOP:
846  if (g->aph_nvec != g->acdom_nvec || g->aph_nvec != g->anap_nvec
847  || g->aph_nvec != g->bbph_nvec || g->aph_nvec != g->bbnap_nvec
848  || g->acdom_nvec != g->anap_nvec || g->acdom_nvec != g->bbph_nvec
849  || g->acdom_nvec != g->bbnap_nvec || g->anap_nvec != g->bbph_nvec
850  || g->anap_nvec != g->bbnap_nvec || g->bbph_nvec != g->bbnap_nvec) {
851 
852  printf("-E- %s line %d: SIOP tables must have same number of columns.\n",
853  __FILE__, __LINE__);
854  exit(1);
855  }
856  break;
857  }
858 
859  idx443 = windex(443., g->acdom_tab_w, g->acdom_tab_nw);
860 
861  for (iw = 0; iw < nwave; iw++) {
862 
863  // evaluate model for each IOP component at all wavelengths
864  // also compute and store associated uncertainty, based on fitted paraneter
865  // uncertainties as stored in the second half of the par array
866  // note that we're using a simple summation for uncertainties (for now)
867 
868  iwtab = windex(wave[iw], g->aph_tab_w, g->aph_tab_nw);
869  aphstar[iw] = g->aph_tab_s[isiop][iwtab];
870  aph[iw] = par[0] * aphstar[iw];
871  aph[iw + nwave] = par[g->npar] * aphstar[iw];
872 
873  /*Esnure that acdom443 is equal to 1.0*/
874  iwtab = windex(wave[iw], g->acdom_tab_w, g->acdom_tab_nw);
875  acdomstar[iw] = g->acdom_tab_s[isiop][iwtab];
876  acdom443 = g->acdom_tab_s[isiop][idx443];
877  acdom[iw] = par[1]*(acdomstar[iw] / acdom443);
878  acdom[iw + nwave] = par[1 + g->npar]*(aphstar[iw] / acdom443);
879 
880  iwtab = windex(wave[iw], g->anap_tab_w, g->anap_tab_nw);
881  anapstar[iw] = g->anap_tab_s[isiop][iwtab];
882  anap[iw] = par[2] * anapstar[iw];
883  anap[iw + nwave] = par[2 + g->npar] * aphstar[iw];
884 
885  adg[iw] = par[1] * acdomstar[iw] + par[2] * anapstar[iw];
886  adg[iw + nwave] = par[1 + g->npar] * acdomstar[iw] + par[2 + g->npar] * anapstar[iw];
887 
888  iwtab = windex(wave[iw], g->bbph_tab_w, g->bbph_tab_nw);
889  bbphstar[iw] = g->bbph_tab_s[isiop][iwtab];
890  bbph[iw] = par[0] * bbphstar[iw];
891  bbph[iw + nwave] = par[g->npar] * bbphstar[iw];
892 
893  iwtab = windex(wave[iw], g->bbnap_tab_w, g->bbnap_tab_nw);
894  bbnapstar[iw] = g->bbnap_tab_s[isiop][iwtab];
895  bbnap[iw] = par[2] * bbnapstar[iw];
896  bbnap[iw + nwave] = par[2 + g->npar] * bbnapstar[iw];
897 
898  bbp[iw] = par[0] * bbphstar[iw] + par[2] * bbnapstar[iw];
899  bbp[iw + nwave] = par[g->npar] * bbphstar[iw] + par[2 + g->npar] * bbnapstar[iw];
900 
901  a = aw [iw] + aph[iw] + adg[iw];
902  bb = bbw[iw] + bbp[iw];
903  x = bb / (a + bb);
904 
905  switch (g->rrs_opt) {
906  case RRSGRD:
907  rrs[iw] = g->grd[0] * x + g->grd[1] * pow(x, 2);
908  dfdx = g->grd[0] + 2 * g->grd[1] * x;
909  break;
910  case RRSFOQ:
911  rrs[iw] = foq[iw] * x;
912  dfdx = foq[iw];
913  break;
914  }
915 
916  if (dfdpar != NULL) {
917 
918  dxda = -x * x / bb;
919  dxdb = x / bb + dxda;
920 
921  dfdpar[iw][0] = dfdx * dxda * aphstar[iw];
922  dfdpar[iw][1] = dfdx * dxda * acdomstar[iw];
923  dfdpar[iw][2] = dfdx * dxdb * anapstar[iw];
924  dfdpar[iw][3] = dfdx * dxdb * bbphstar[iw];
925  dfdpar[iw][4] = dfdx * dxdb * bbnapstar[iw];
926 
927 
928  }
929 
930  if (parstar != NULL) {
931 
932  parstar[iw][0] = aphstar[iw];
933  parstar[iw][1] = acdomstar[iw];
934  parstar[iw][2] = anapstar[iw];
935  parstar[iw][3] = bbphstar[iw];
936  parstar[iw][4] = bbnapstar[iw];
937 
938  }
939 
940  }
941 
942  free(aphstar);
943  free(acdomstar);
944  free(anapstar);
945  free(bbphstar);
946  free(bbnapstar);
947 
948  return;
949 }
950 
951 /* ------------------------------------------------------------------- */
952 /* */
953 
954 /* ------------------------------------------------------------------- */
955 double giop_amb(FITSTRUCT *ambdata, double par[]) {
956  int iw;
957 
958  giopstr *g = (giopstr *) (ambdata->meta);
959 
960  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, ambdata->yfit, NULL, NULL);
961 
962  ambdata->merit = 0.0;
963  for (iw = 0; iw < g->nwave; iw++) {
964  ambdata->merit += pow((ambdata->y[iw] - ambdata->yfit[iw]), 2) * ambdata->wgt[iw];
965  }
966 
967  return (ambdata->merit);
968 }
969 
970 
971 /* ------------------------------------------------------------------- */
972 /* */
973 
974 /* ------------------------------------------------------------------- */
975 int fit_giop_amb(giopstr *g, double Rrs[], double wts[], double par[],
976  double Rrs_fit[], int16 *itercnt) {
977  int status = 0;
978 
979 
980  static float tol = 1.e-6; /* fractional change in chisqr */
981  static FITSTRUCT ambdata; /* amoeba interface structure */
982  static double *init; /* initial simplex */
983 
984  static int firstCall = 1;
985 
986  int i, j;
987  short isml;
988 
989  ambdata.niter = g->maxiter; /* max number of iterations */
990  ambdata.nfunc = g->npar; /* number of model parameters */
991  ambdata.npnts = g->nwave; /* number of wavelengths (Rrs) */
992  ambdata.y = Rrs; /* Input Rrs values (subsurface) */
993  ambdata.wgt = wts; /* Input weights on Rrs values */
994  ambdata.yfit = Rrs_fit; /* Output model predicted Rrs values */
995  ambdata.meta = g;
996 
997  if (firstCall == 1) {
998  firstCall = 0;
999  if ((init = (double *) calloc(g->nwave * (g->nwave + 1), sizeof (double))) == NULL) {
1000  printf("-E- %s line %d : error allocating memory for GIOP:gio_model.\n",
1001  __FILE__, __LINE__);
1002  exit(1);
1003  }
1004 
1005  }
1006  /* initialize simplex with first guess model parameters */
1007  for (j = 0; j < g->npar + 1; j++)
1008  for (i = 0; i < g->npar; i++)
1009  init[j * g->npar + i] = g->par[i];
1010 
1011  for (i = 0; i < g->npar; i++) {
1012  init[g->npar + i * (g->npar + 1)] += g->len[i];
1013  par[i] = 0.0;
1014  }
1015 
1016  /* run optimization */
1017  isml = amoeba(init, &ambdata, giop_amb, tol);
1018 
1019  /* check convergence and record parameter results */
1020  if (ambdata.niter >= g->maxiter)
1021  status = 1;
1022 
1023  for (i = 0; i < g->npar; i++) {
1024  par[i] = init[g->npar * isml + i];
1025  }
1026 
1027  *itercnt = ambdata.niter;
1028 
1029  return (status);
1030 }
1031 
1032 
1033 
1034 /* ---------------------------------------------------------------------- */
1035 /* wrapper function for L-M fit to GIOP model */
1036 
1037 /* ---------------------------------------------------------------------- */
1038 int giop_lm_fdf(const gsl_vector *parv, void *data, gsl_vector *f, gsl_matrix *J) {
1039  double *y = ((fitstr *) data)->y;
1040  double *w = ((fitstr *) data)->w;
1041  float *aw = ((fitstr *) data)->g->aw;
1042  float *bbw = ((fitstr *) data)->g->bbw;
1043  float *foq = ((fitstr *) data)->g->foq;
1044  float *wave = ((fitstr *) data)->g->wave;
1045  size_t nwave = ((fitstr *) data)->g->nwave;
1046  size_t npar = ((fitstr *) data)->g->npar;
1047  giopstr *g = ((fitstr *) data)->g;
1048 
1049  double *par;
1050  double *yfit;
1051  double **dydpar;
1052 
1053  double sigma;
1054 
1055  int iw, ipar;
1056 
1057  if ((par = (double *) calloc(nwave, sizeof (double))) == NULL) {
1058  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1059  __FILE__, __LINE__);
1060  exit(1);
1061  }
1062 
1063  /* extract model params from vector */
1064  for (ipar = 0; ipar < npar; ipar++) {
1065  par[ipar] = gsl_vector_get(parv, ipar);
1066  }
1067  if ((yfit = (double *) calloc(nwave, sizeof (double))) == NULL) {
1068  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1069  __FILE__, __LINE__);
1070  exit(1);
1071  }
1072  if ((dydpar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1073  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1074  __FILE__, __LINE__);
1075  exit(1);
1076  }
1077  for (iw = 0; iw < nwave; iw++)
1078  if ((dydpar[iw] = (double *) calloc(nwave, sizeof (double))) == NULL) {
1079  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1080  __FILE__, __LINE__);
1081  exit(1);
1082  }
1083 
1084  /* run model */
1085  giop_model(g, par, nwave, wave, aw, bbw, foq, aph1, adg1, bbp1, yfit, dydpar, NULL);
1086 
1087  /* evaluate and store for lm solver */
1088  for (iw = 0; iw < nwave; iw++) {
1089 
1090  /* silly, but we need sigma and not sigma-squared */
1091  sigma = sqrt(1. / w[iw]);
1092 
1093  /* Function to be minimized */
1094  if (f != NULL)
1095  gsl_vector_set(f, iw, (yfit[iw] - y[iw]) / sigma);
1096 
1097  /* Jacobian matrix */
1098  if (J != NULL)
1099  for (ipar = 0; ipar < npar; ipar++)
1100  gsl_matrix_set(J, iw, ipar, dydpar[iw][ipar] / sigma);
1101  }
1102 
1103  free(yfit);
1104  freeDArray(dydpar, nwave);
1105  free(par);
1106  return GSL_SUCCESS;
1107 }
1108 
1109 int giop_lm_f(const gsl_vector *parv, void *data, gsl_vector *f) {
1110  return (giop_lm_fdf(parv, data, f, NULL));
1111 }
1112 
1113 int giop_lm_df(const gsl_vector *parv, void *data, gsl_matrix *J) {
1114  return (giop_lm_fdf(parv, data, NULL, J));
1115 }
1116 
1117 
1118 
1119 /* ---------------------------------------------------------------------- */
1120 /* fit_giop_lm() - runs Levenburg-Marquart optimization for one pixel */
1121 
1122 /* ---------------------------------------------------------------------- */
1123 int fit_giop_lm(giopstr *g, double Rrs[], double wts[], double par[], double *chi, int16 *itercnt) {
1124  int status = 0;
1125  int iter;
1126 
1127  static fitstr data;
1128 
1129  size_t npar = g->npar;
1130 
1131  static gsl_multifit_function_fdf func;
1132  const gsl_multifit_fdfsolver_type *t;
1133  gsl_multifit_fdfsolver *s;
1134  gsl_vector_view x;
1135 
1136  gsl_matrix *J = gsl_matrix_alloc(g->nwave, npar);
1137  gsl_matrix *cov = gsl_matrix_alloc(npar, npar);
1138 
1139  double sum, dof, c;
1140  int ipar;
1141 
1142  /* Set up data structure */
1143  data.y = Rrs;
1144  data.w = wts;
1145  data.g = g;
1146 
1147  /* Set up multifit function structure */
1148  func.f = &giop_lm_f;
1149  func.df = &giop_lm_df;
1150  func.fdf = &giop_lm_fdf;
1151  func.n = g->nwave;
1152  func.p = g->npar;
1153  func.params = &data;
1154 
1155  /* Allocate solver space */
1156  t = gsl_multifit_fdfsolver_lmsder;
1157  s = gsl_multifit_fdfsolver_alloc(t, g->nwave, g->npar);
1158 
1159  /* Set start params */
1160  x = gsl_vector_view_array(par, npar);
1161  gsl_multifit_fdfsolver_set(s, &func, &x.vector);
1162 
1163  /* Fit model for this pixel */
1164  status = 1;
1165  for (iter = 0; iter < g->maxiter; iter++) {
1166  gsl_multifit_fdfsolver_iterate(s);
1167  if (gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4) == GSL_SUCCESS) {
1168  status = 0;
1169  break;
1170  }
1171  }
1172  *itercnt = iter;
1173 
1174  /* Compute covariance matrix from Jacobian */
1175  gsl_multifit_fdfsolver_jac(s, J);
1176  gsl_multifit_covar(J, 0.0, cov);
1177 
1178  /* Compute chi-square */
1179  sum = gsl_blas_dnrm2(s->f); // quadrature sum of minimization function
1180  dof = func.n - func.p; // degrees of freedom
1181  *chi = pow(sum, 2.0) / dof; // chi-square per dof
1182 
1183  if (g->wt_opt == 0)
1184  c = sum / sqrt(dof); // use variance in fit to estimate Rrs uncertainty
1185  else
1186  c = 1.0; // Rrs uncertainty included in sum
1187 
1188  /* Retrieve fit params & error estimates */
1189  for (ipar = 0; ipar < npar; ipar++) {
1190  par[ipar] = gsl_vector_get(s->x, ipar);
1191  par[ipar + npar] = sqrt(gsl_matrix_get(cov, ipar, ipar)) * c;
1192  //printf("%d %lf %lf\n",ipar,par[ipar],par[ipar+npar]);
1193  }
1194 
1195  gsl_multifit_fdfsolver_free(s);
1196  gsl_matrix_free(cov);
1197  gsl_matrix_free(J);
1198 
1199  return (status);
1200 }
1201 
1202 
1203 /* ---------------------------------------------------------------------- */
1204 /* fit_giop_svd() - constrained matrix solution */
1205 
1206 /* ---------------------------------------------------------------------- */
1207 int fit_giop_svd(giopstr *g, double rrs[], double wts[], double par[]) {
1208  int status = 0; /* init to success */
1209 
1210  double *rrs_fit;
1211  double **parstar;
1212 
1213  size_t nwave = g->nwave;
1214  size_t npar = g->npar;
1215 
1216  int iw, ipar;
1217 
1218  double a0, a1, a2;
1219  double u0, u1, u;
1220 
1221  gsl_matrix *A = gsl_matrix_alloc(nwave, npar);
1222  gsl_matrix *V = gsl_matrix_alloc(npar, npar);
1223  gsl_vector *S = gsl_vector_alloc(npar);
1224  gsl_vector *W = gsl_vector_alloc(npar);
1225  gsl_vector *x = gsl_vector_alloc(npar);
1226  gsl_vector *b = gsl_vector_alloc(nwave);
1227 
1228  if ((rrs_fit = (double *) calloc(nwave, sizeof (double))) == NULL) {
1229  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1230  __FILE__, __LINE__);
1231  exit(1);
1232  }
1233  if ((parstar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1234  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1235  __FILE__, __LINE__);
1236  exit(1);
1237  }
1238  for (iw = 0; iw < nwave; iw++)
1239  if ((parstar[iw] = (double *) calloc(nwave, sizeof (double))) == NULL) {
1240  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1241  __FILE__, __LINE__);
1242  exit(1);
1243  }
1244 
1245  /* run model to get parstar wavelength-dependence terms */
1246  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, rrs_fit, NULL, parstar);
1247 
1248  /* clear fit parameters */
1249  for (ipar = 0; ipar < g->npar; ipar++)
1250  par[ipar] = BAD_FLT;
1251 
1252  /* build A matrix and b vector for A x = b */
1253 
1254  for (iw = 0; iw < g->nwave; iw++) {
1255 
1256  /* get u = bb/(a+bb) */
1257  switch (g->rrs_opt) {
1258  case RRSGRD:
1259  a2 = g->grd[1];
1260  a1 = g->grd[0];
1261  a0 = -rrs[iw];
1262  if ((gsl_poly_solve_quadratic(a2, a1, a0, &u0, &u1) == 2) && u1 > 0.0)
1263  u = u1;
1264  else {
1265  status = 1;
1266  goto cleanup;
1267  }
1268  break;
1269  case RRSFOQ:
1270  u = rrs[iw] / g->foq[iw];
1271  break;
1272  }
1273 
1274  gsl_vector_set(b, iw, -(u * g->aw[iw] + (u - 1.0) * g->bbw[iw]));
1275 
1276  for (ipar = 0; ipar < g->npar; ipar++) {
1277  if (ipar < (g->aph_nvec + g->adg_nvec))
1278  gsl_matrix_set(A, iw, ipar, parstar[iw][ipar] * u); // a
1279  else
1280  gsl_matrix_set(A, iw, ipar, parstar[iw][ipar]*(u - 1.0)); // bb
1281  }
1282  }
1283 
1284  /* solve A x = b for x */
1285  status = gsl_linalg_SV_decomp(A, V, S, W);
1286  status = gsl_linalg_SV_solve(A, V, S, b, x);
1287 
1288  /* extract fitted parameters */
1289  if (status == 0) {
1290  for (ipar = 0; ipar < g->npar; ipar++) {
1291  par[ipar] = gsl_vector_get(x, ipar);
1292  }
1293  }
1294 
1295 cleanup:
1296 
1297  gsl_matrix_free(A);
1298  gsl_matrix_free(V);
1299  gsl_vector_free(S);
1300  gsl_vector_free(x);
1301  gsl_vector_free(b);
1302 
1303  free(rrs_fit);
1304  freeDArray(parstar, nwave);
1305  return (status);
1306 }
1307 
1308 /* ---------------------------------------------------------------------- */
1309 /* fit_giop_svd_siop() - adaptive linear matrix inversion solution using */
1310 /* specific inherent optical properties (SIOPS) and */
1311 /* SVD to solve system of equations */
1312 /* */
1313 /* Reference: */
1314 /* Brando et al. (2012) Adaptive semianalytical inversion of ocean color */
1315 /* radiometry in optically complex waters, Applied Optics, 51(15), */
1316 /* 2808-2833. */
1317 /* */
1318 /* Implementation: L. McKinna, SAIC/NASA GSFC, November 2015 */
1319 /* */
1320 
1321 /* ---------------------------------------------------------------------- */
1322 int fit_giop_svd_siop(giopstr *g, double rrs[], double wts[], double par[], double *chi) {
1323 
1324  int16 status = 0; /* init to success */
1325 
1326  double *rrs_fit;
1327  double **parstar;
1328  double *parArr;
1329  double *parFit;
1330  double *rmse;
1331  int *badSolution;
1332 
1333 
1334  size_t nwave = g->nwave;
1335  size_t npar = g->npar; //Number of parameters fixed at 3
1336  size_t nvec = g->aph_nvec;
1337 
1338  int iw, ipar, iv, smlIdx;
1339 
1340  double a0, a1, a2;
1341  double u0, u1, u;
1342 
1343  double diffSq, diffSqSum, smlRmse, sumRrs;
1344 
1345  gsl_matrix *A = gsl_matrix_alloc(nwave, npar);
1346  gsl_matrix *V = gsl_matrix_alloc(npar, npar);
1347  gsl_vector *S = gsl_vector_alloc(npar);
1348  gsl_vector *W = gsl_vector_alloc(npar);
1349  gsl_vector *x = gsl_vector_alloc(npar);
1350  gsl_vector *b = gsl_vector_alloc(nwave);
1351 
1352 
1353  if ((rrs_fit = (double *) calloc(nwave, sizeof (double))) == NULL) {
1354  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1355  __FILE__, __LINE__);
1356  exit(1);
1357  }
1358  if ((parstar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1359  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1360  __FILE__, __LINE__);
1361  exit(1);
1362  }
1363  for (iw = 0; iw < nwave; iw++) {
1364  if ((parstar[iw] = (double *) calloc(5, sizeof (double))) == NULL) {
1365  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1366  __FILE__, __LINE__);
1367  exit(1);
1368  }
1369  }
1370  if ((parArr = (double *) calloc(npar * nvec, sizeof (double *))) == NULL) {
1371  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1372  __FILE__, __LINE__);
1373  exit(1);
1374  }
1375  if ((rmse = (double *) calloc(nvec, sizeof (double *))) == NULL) {
1376  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1377  __FILE__, __LINE__);
1378  exit(1);
1379  }
1380  if ((badSolution = (int *) calloc(nvec, sizeof (int *))) == NULL) {
1381  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1382  __FILE__, __LINE__);
1383  exit(1);
1384  }
1385  if ((parFit = (double *) calloc(npar, sizeof (double *))) == NULL) {
1386  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1387  __FILE__, __LINE__);
1388  exit(1);
1389  }
1390 
1391  /* intialize fit parameters */
1392  for (ipar = 0; ipar < g->npar; ipar++)
1393  par[ipar] = BAD_FLT;
1394 
1395  for (ipar = 0; ipar < 3; ipar++)
1396  parFit[ipar] = BAD_FLT;
1397 
1398 
1399  /* clear fit parameters */
1400  for (iv = 0; iv < nvec; iv++)
1401  badSolution[iv] = BAD_INT;
1402 
1403 
1404  /*Iterate over siop combinations*/
1405  for (iv = 0; iv < nvec; iv++) {
1406 
1407  g->siopIdx = iv;
1408 
1409  /* run model to get parstar wavelength-dependence terms */
1410  giop_model_iterate(g, parFit, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, acdom1, anap1, bbph1, bbnap1, rrs_fit, NULL, parstar);
1411 
1412  /* build A matrix and b vector for A x = b */
1413  for (iw = 0; iw < g->nwave; iw++) {
1414 
1415  /* get u = bb/(a+bb) */
1416  switch (g->rrs_opt) {
1417  case RRSGRD:
1418  a2 = g->grd[1];
1419  a1 = g->grd[0];
1420  a0 = -rrs[iw];
1421  if ((gsl_poly_solve_quadratic(a2, a1, a0, &u0, &u1) == 2) && u1 > 0.0)
1422  u = u1;
1423  else {
1424  status = 1;
1425  goto cleanup;
1426  }
1427  break;
1428  case RRSFOQ:
1429  u = rrs[iw] / g->foq[iw];
1430  break;
1431  }
1432 
1433  gsl_vector_set(b, iw, -(g->aw[iw] * u + g->bbw[iw]*(u - 1.0)));
1434 
1435  // SIOP implementation Aij in eq 12 of Brando et al 2012
1436  // A is 3 components + water
1437  // i=0 (phy): ipar=0 for aphy* and ipar 3 for bbphy*
1438  // i=1 (CDOM): ipar=1 for aCDOM*
1439  // i=2 (NAP): ipar=2 for anap* and ipar 4 for bbnap*
1440  gsl_matrix_set(A, iw, 0, parstar[iw][0] * u + parstar[iw][3]*(u - 1.0)); // a+ bb
1441  gsl_matrix_set(A, iw, 1, parstar[iw][1] * u); // ONLY a
1442  gsl_matrix_set(A, iw, 2, parstar[iw][2] * u + parstar[iw][4]*(u - 1.0)); // a+ bb
1443 
1444  }
1445 
1446  /* solve A x = b for x */
1447  /*SVD decomposition*/
1448  status = gsl_linalg_SV_decomp(A, V, S, W);
1449  status = gsl_linalg_SV_solve(A, V, S, b, x);
1450 
1451 
1452  for (ipar = 0; ipar < npar; ipar++) {
1453 
1454  /*Test for non-physical (negative) concentrations or matrix inversion*/
1455  /*failure. For such circumstances set to BAD_FLT*/
1456  if (gsl_vector_get(x, ipar) < 0.0 || status != 0) {
1457  badSolution[iv] = 1;
1458  for (ipar = 0; ipar < npar; ipar++) {
1459  parArr[iv * npar + ipar] = BAD_FLT;
1460  }
1461  break;
1462 
1463  } else if (gsl_vector_get(x, ipar) >= 0.0) {
1464  badSolution[iv] = 0;
1465  parArr[iv * npar + ipar] = gsl_vector_get(x, ipar);
1466  }
1467  }
1468  }
1469 
1470  /*Find optimal set of parameters*/
1471  g->siopIdx = 0;
1472  diffSq = 0;
1473  diffSqSum = 0;
1474 
1475  /*Compute the root mean square error (rmse) metric for all SIOP combinations. */
1476  /*The smallest rmse valid values is returned as the optimal solution.*/
1477  for (iv = 0; iv < nvec; iv++) {
1478 
1479  g->siopIdx = iv;
1480 
1481  for (ipar = 0; ipar < npar; ipar++) {
1482  //if (badSolution[iv]) {
1483  // parFit[ipar] = BAD_FLT;
1484  //} else {
1485  parFit[ipar] = parArr[iv * npar + ipar];
1486  //}
1487  }
1488 
1489  /*Calculate fitted Rrs*/
1490  giop_model_iterate(g, parFit, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, acdom1, anap1, bbph1, bbnap1, rrs_fit, NULL, parstar);
1491 
1492  /*reset the sum of squares values*/
1493  diffSq = 0;
1494  diffSqSum = 0;
1495  sumRrs = 0;
1496 
1497  for (iw = 0; iw < g->nwave; iw++) {
1498  diffSq = pow((rrs_fit[iw] - rrs[iw]), 2);
1499  diffSqSum += diffSq;
1500  sumRrs += rrs[iw];
1501  }
1502 
1503 
1504  if (badSolution[iv]) {
1505  rmse[iv] = 999;
1506  } else {
1507  /*Save value to relative RMSE array*/
1508  rmse[iv] = pow((diffSqSum / (g->nwave - 1)), 0.5); // sumRrs/(g->nwave - 1) );
1509  }
1510  }
1511 
1512  /*Find the smallest relative RMSE value*/
1513  /*Initialize with zeroth values*/
1514  smlRmse = rmse[0];
1515  smlIdx = 0;
1516 
1517  for (iv = 0; iv < nvec; iv++) {
1518  if (rmse[iv] < smlRmse) {
1519  smlRmse = rmse[iv];
1520  smlIdx = iv;
1521  }
1522  }
1523 
1524 
1525  //Save index of best SIOP combination to giop record
1526  g->siopIdx = smlIdx;
1527 
1528  /*Get optimal concentrations values from array of all concentrations*/
1529  /*If there is still a negative solution, set it to zero*/
1530  if (rmse[smlIdx] != 999 && badSolution[smlIdx] != 1) {
1531 
1532  for (ipar = 0; ipar < npar; ipar++) {
1533  par[ipar] = parArr[smlIdx * npar + ipar];
1534  }
1535  *chi = smlRmse;
1536  status = 0;
1537 
1538 
1539  } else {
1540  //No solution status
1541  for (ipar = 0; ipar < npar; ipar++) {
1542  par[ipar] = BAD_FLT;
1543  }
1544  *chi = BAD_FLT;
1545  status = -99;
1546 
1547  }
1548 
1549 cleanup:
1550  gsl_matrix_free(A);
1551  gsl_matrix_free(V);
1552  gsl_vector_free(S);
1553  gsl_vector_free(x);
1554  gsl_vector_free(b);
1555 
1556 
1557  free(rrs_fit);
1558  free(parArr);
1559  free(rmse);
1560  free(badSolution);
1561  freeDArray(parstar, nwave);
1562  return (status);
1563 }
1564 
1565 /* ------------------------------------------------------------------- */
1566 /* giop_chl() - returns magnitude of aph* as chl */
1567 
1568 /* ------------------------------------------------------------------- */
1569 float32 giop_chl(giopstr *g, int16 iopf, double par[], float *chl_err) {
1570  float32 chl = badval;
1571  int ipar;
1572 
1573  *chl_err = badval;
1574 
1575  switch (g->aph_opt) {
1576  case APHGAUSS:
1577  break;
1578  default:
1579  if ((iopf & IOPF_FAILED) == 0) {
1580  chl = 0.0;
1581  *chl_err = 0.0;
1582  for (ipar = 0; ipar < g->aph_nvec; ipar++) {
1583  chl += par[ipar];
1584  *chl_err += par[ipar + g->npar];
1585  }
1586  }
1587  break;
1588  }
1589 
1590 
1591  return (chl);
1592 }
1593 
1594 /* ---------------------------------------------------------------------- */
1595 /* Convert Rrs[0+] to Rrs[0-] */
1596 
1597 /* ---------------------------------------------------------------------- */
1598 float rrs_above_to_below(float Rrs) {
1599  return (Rrs / (0.52 + 1.7 * Rrs));
1600 
1601 }
1602 
1603 /* ---------------------------------------------------------------------- */
1604 /* Convert Rrs[0-] to Rrs[0+] */
1605 
1606 /* ---------------------------------------------------------------------- */
1607 float rrs_below_to_above(float rrs_s) {
1608  return ( (0.52 * rrs_s) / (1.0 - 1.7 * rrs_s));
1609 
1610 }
1611 
1612 /* ---------------------------------------------------------------------- */
1613 /* run_giop - runs optimization using requested method */
1614 
1615 /* ---------------------------------------------------------------------- */
1616 void run_giop(l2str *l2rec) {
1617  static int firstCall = 1;
1618  static giopstr giopctl;
1619  static giopstr *g = &giopctl;
1620 
1621  double *par;
1622  double *Rrs_a; /* above water, per fit band */
1623  double *Rrs_b; /* below water, per fit band */
1624  double *Rrs_f; /* modeled Rrs, per fit band */
1625  double *wts; /* weights, per fit band */
1626 
1627  float Rrs1, Rrs2;
1628  int32 i1, i2;
1629 
1630  int16 itercnt;
1631  int16 bndcnt;
1632  int16 status;
1633 
1634  int32 ipar, ip, iw, ib, ipb, ipb2, ierr;
1635  double chi = BAD_FLT;
1636 
1637  l1str *l1rec = l2rec->l1rec;
1638 
1639  float *wave = l1rec->l1file->fwave;
1640  int32 nwave = l1rec->l1file->nbands;
1641  static int32 npix;
1642 
1643 
1644  float aph_norm = BAD_FLT;
1645 
1646  if (firstCall) {
1647  npix = l1rec->npix;
1648 
1649  firstCall = 0;
1650 
1651  // initialize control structure (to get npar)
1652  if ((bbw = calloc(nwave, sizeof (float))) == NULL) {
1653  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1654  __FILE__, __LINE__);
1655  exit(1);
1656  }
1657  if ((aw = calloc(nwave, sizeof (float))) == NULL) {
1658  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1659  __FILE__, __LINE__);
1660  exit(1);
1661  }
1662  if ((foq = calloc(nwave, sizeof (float))) == NULL) {
1663  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1664  __FILE__, __LINE__);
1665  exit(1);
1666  }
1667  if ((aph1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1668  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1669  __FILE__, __LINE__);
1670  exit(1);
1671  }
1672  if ((acdom1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1673  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1674  __FILE__, __LINE__);
1675  exit(1);
1676  }
1677  if ((anap1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1678  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1679  __FILE__, __LINE__);
1680  exit(1);
1681  }
1682  if ((adg1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1683  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1684  __FILE__, __LINE__);
1685  exit(1);
1686  }
1687  if ((bbph1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1688  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1689  __FILE__, __LINE__);
1690  exit(1);
1691  }
1692  if ((bbnap1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1693  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1694  __FILE__, __LINE__);
1695  exit(1);
1696  }
1697  if ((bbp1 = calloc(2 * nwave, sizeof (float))) == NULL) {
1698  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1699  __FILE__, __LINE__);
1700  exit(1);
1701  }
1702  if ((g->wave = (float *) calloc(nwave, sizeof (float))) == NULL) {
1703  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1704  __FILE__, __LINE__);
1705  exit(1);
1706  }
1707  if ((g->bindx = (int *) calloc(nwave, sizeof (int))) == NULL) {
1708  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1709  __FILE__, __LINE__);
1710  exit(1);
1711  }
1712  if ((g->aw = (float *) calloc(nwave, sizeof (float))) == NULL) {
1713  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1714  __FILE__, __LINE__);
1715  exit(1);
1716  }
1717  if ((g->bbw = (float *) calloc(nwave, sizeof (float))) == NULL) {
1718  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1719  __FILE__, __LINE__);
1720  exit(1);
1721  }
1722  if ((g->wts = (float *) calloc(nwave, sizeof (float))) == NULL) {
1723  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1724  __FILE__, __LINE__);
1725  exit(1);
1726  }
1727  if ((g->foq = (float *) calloc(nwave, sizeof (float))) == NULL) {
1728  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1729  __FILE__, __LINE__);
1730  exit(1);
1731  }
1732  if ((g->par = (double *) calloc(nwave, sizeof (double))) == NULL) {
1733  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1734  __FILE__, __LINE__);
1735  exit(1);
1736  }
1737  if ((g->len = (double *) calloc(nwave, sizeof (double))) == NULL) {
1738  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
1739  __FILE__, __LINE__);
1740  exit(1);
1741  }
1742 
1743 
1744  giop_ctl_init(g, nwave, wave, aw, bbw);
1745 
1746  // allocate static storage for one scanline
1747 
1748  if ((iter = calloc(npix, sizeof (int16))) == NULL) {
1749  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1750  __FILE__, __LINE__);
1751  exit(1);
1752  }
1753  if ((iopf = calloc(npix, sizeof (int16))) == NULL) {
1754  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1755  __FILE__, __LINE__);
1756  exit(1);
1757  }
1758  if ((mRrs = calloc(npix * nwave, sizeof (float))) == NULL) {
1759  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1760  __FILE__, __LINE__);
1761  exit(1);
1762  }
1763  if ((a = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1764  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1765  __FILE__, __LINE__);
1766  exit(1);
1767  }
1768  if ((aph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1769  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1770  __FILE__, __LINE__);
1771  exit(1);
1772  }
1773  if ((acdom = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1774  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1775  __FILE__, __LINE__);
1776  exit(1);
1777  }
1778  if ((anap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1779  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1780  __FILE__, __LINE__);
1781  exit(1);
1782  }
1783  if ((adg = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1784  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1785  __FILE__, __LINE__);
1786  exit(1);
1787  }
1788  if ((bbph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1789  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1790  __FILE__, __LINE__);
1791  exit(1);
1792  }
1793  if ((bbnap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1794  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1795  __FILE__, __LINE__);
1796  exit(1);
1797  }
1798  if ((bbp = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1799  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1800  __FILE__, __LINE__);
1801  exit(1);
1802  }
1803  if ((bb = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1804  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1805  __FILE__, __LINE__);
1806  exit(1);
1807  }
1808  if ((chl = calloc(npix * 2, sizeof (float))) == NULL) {
1809  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1810  __FILE__, __LINE__);
1811  exit(1);
1812  }
1813  if ((rrsdiff = calloc(npix, sizeof (float))) == NULL) {
1814  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1815  __FILE__, __LINE__);
1816  exit(1);
1817  }
1818  if ((aph_s = calloc(npix, sizeof (float))) == NULL) {
1819  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1820  __FILE__, __LINE__);
1821  exit(1);
1822  }
1823  if ((adg_s = calloc(npix, sizeof (float))) == NULL) {
1824  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1825  __FILE__, __LINE__);
1826  exit(1);
1827  }
1828  if ((bbp_s = calloc(npix, sizeof (float))) == NULL) {
1829  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1830  __FILE__, __LINE__);
1831  exit(1);
1832  }
1833  if ((fit_par = allocate2d_float(npix, g->npar)) == NULL) {
1834  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1835  __FILE__, __LINE__);
1836  exit(1);
1837  }
1838  if ((siop_num = calloc(npix, sizeof (int))) == NULL) {
1839  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1840  __FILE__, __LINE__);
1841  exit(1);
1842  }
1843  max_npar = g->npar;
1844  if ((chisqr = calloc(npix, sizeof (float))) == NULL) {
1845  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1846  __FILE__, __LINE__);
1847  exit(1);
1848  }
1849 
1850  /*Note: Calculation of Raman scattering contribution to Rrs is at present*/
1851  /*only supported in l2gen. Where Raman Rrs is not calculated (i.e. l3gen smi)*/
1852  /*set Rrs_raman to zeros */
1853  if (l2rec->Rrs_raman == NULL) {
1854 
1855  printf("\n");
1856  printf("No Raman scattering correction applied to Rrs. \n");
1857  printf("\n");
1858 
1859  allocateRrsRaman = 1;
1860  l2rec->Rrs_raman = (float*) allocateMemory(npix *
1861  l1rec->l1file->nbands * sizeof (float), "Rrs_ram");
1862  }
1863 
1864  }
1865 
1866  if ((par = calloc(2 * nwave, sizeof (double))) == NULL) {
1867  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
1868  __FILE__, __LINE__);
1869  exit(1);
1870  }
1871 
1872  // reallocate if npix changes - l3gen-ism...
1873  if (l1rec->npix > npix) {
1874  npix = l1rec->npix;
1875  free(iter);
1876  free(iopf);
1877  free(mRrs);
1878  free(a);
1879  free(aph);
1880  free(acdom);
1881  free(anap);
1882  free(adg);
1883  free(bbph);
1884  free(bbnap);
1885  free(bbp);
1886  free(bb);
1887  free(chl);
1888  free(rrsdiff);
1889  free(aph_s);
1890  free(adg_s);
1891  free(bbp_s);
1892  free(fit_par);
1893  free(siop_num);
1894  free(chisqr);
1895  if (allocateRrsRaman) {
1896  free(l2rec->Rrs_raman);
1897  l2rec->Rrs_raman = (float*) allocateMemory(npix *
1898  l1rec->l1file->nbands * sizeof (float), "Rrs_ram");
1899  }
1900 
1901  if ((iter = calloc(npix, sizeof (int16))) == NULL) {
1902  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1903  __FILE__, __LINE__);
1904  exit(1);
1905  }
1906  if ((iopf = calloc(npix, sizeof (int16))) == NULL) {
1907  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1908  __FILE__, __LINE__);
1909  exit(1);
1910  }
1911  if ((mRrs = calloc(npix * nwave, sizeof (float))) == NULL) {
1912  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1913  __FILE__, __LINE__);
1914  exit(1);
1915  }
1916  if ((a = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1917  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1918  __FILE__, __LINE__);
1919  exit(1);
1920  }
1921  if ((aph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1922  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1923  __FILE__, __LINE__);
1924  exit(1);
1925  }
1926  if ((acdom = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1927  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1928  __FILE__, __LINE__);
1929  exit(1);
1930  }
1931  if ((anap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1932  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1933  __FILE__, __LINE__);
1934  exit(1);
1935  }
1936  if ((adg = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1937  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1938  __FILE__, __LINE__);
1939  exit(1);
1940  }
1941  if ((bbph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1942  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1943  __FILE__, __LINE__);
1944  exit(1);
1945  }
1946  if ((bbnap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1947  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1948  __FILE__, __LINE__);
1949  exit(1);
1950  }
1951  if ((bbp = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1952  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1953  __FILE__, __LINE__);
1954  exit(1);
1955  }
1956  if ((bb = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
1957  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1958  __FILE__, __LINE__);
1959  exit(1);
1960  }
1961  if ((chl = calloc(npix * 2, sizeof (float))) == NULL) {
1962  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1963  __FILE__, __LINE__);
1964  exit(1);
1965  }
1966  if ((rrsdiff = calloc(npix, sizeof (float))) == NULL) {
1967  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1968  __FILE__, __LINE__);
1969  exit(1);
1970  }
1971  if ((aph_s = calloc(npix, sizeof (float))) == NULL) {
1972  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1973  __FILE__, __LINE__);
1974  exit(1);
1975  }
1976  if ((adg_s = calloc(npix, sizeof (float))) == NULL) {
1977  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1978  __FILE__, __LINE__);
1979  exit(1);
1980  }
1981  if ((bbp_s = calloc(npix, sizeof (float))) == NULL) {
1982  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1983  __FILE__, __LINE__);
1984  exit(1);
1985  }
1986  if ((fit_par = allocate2d_float(npix, g->npar)) == NULL) {
1987  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1988  __FILE__, __LINE__);
1989  exit(1);
1990  }
1991  if ((siop_num = calloc(npix, sizeof (int))) == NULL) {
1992  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1993  __FILE__, __LINE__);
1994  exit(1);
1995  }
1996  max_npar = g->npar;
1997  if ((chisqr = calloc(npix, sizeof (float))) == NULL) {
1998  printf("-E- %s line %d : error allocating memory for GIOP.\n",
1999  __FILE__, __LINE__);
2000  exit(1);
2001  }
2002 
2003  }
2004 
2005  if ((Rrs_a = calloc(nwave, sizeof (double))) == NULL) {
2006  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2007  __FILE__, __LINE__);
2008  exit(1);
2009  }
2010  if ((Rrs_b = calloc(nwave, sizeof (double))) == NULL) {
2011  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2012  __FILE__, __LINE__);
2013  exit(1);
2014  }
2015  if ((Rrs_f = calloc(nwave, sizeof (double))) == NULL) {
2016  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2017  __FILE__, __LINE__);
2018  exit(1);
2019  }
2020  if ((wts = calloc(nwave, sizeof (double))) == NULL) {
2021  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2022  __FILE__, __LINE__);
2023  exit(1);
2024  }
2025 
2026  // re-initialize control structure when giop_iterate=1 (under evaluation)
2027  // PJW 9 Jan 2013
2028  if (input->giop_iterate > 0) {
2029  g->adg_s = input->giop_adg_s;
2030  g->bbp_s = input->giop_bbp_s;
2031  }
2032 
2033 
2034  for (ip = 0; ip < l1rec->npix; ip++) {
2035 
2036  ipb = ip*nwave;
2037  ipb2 = ip * l1rec->l1file->nbands;
2038  ierr = l1rec->npix * nwave + ipb;
2039 
2040  // initialize output arrays and counters
2041 
2042  for (ib = 0; ib < nwave; ib++) {
2043  a [ipb + ib] = badval;
2044  bb [ipb + ib] = badval;
2045  aph [ipb + ib] = badval;
2046  adg [ipb + ib] = badval;
2047  acdom [ipb + ib] = badval;
2048  anap[ipb + ib] = badval;
2049  bbp [ipb + ib] = badval;
2050  bbph[ipb + ib] = badval;
2051  bbnap [ipb + ib] = badval;
2052  mRrs[ipb + ib] = badval;
2053  a [ierr + ib] = 0.0;
2054  bb [ierr + ib] = 0.0;
2055  aph [ierr + ib] = 0.0;
2056  adg [ierr + ib] = 0.0;
2057  acdom [ierr + ib] = 0.0;
2058  anap [ierr + ib] = 0.0;
2059  bbp [ierr + ib] = 0.0;
2060  bbph [ierr + ib] = 0.0;
2061  bbnap [ierr + ib] = 0.0;
2062  }
2063  chl [ip] = badval;
2064  iter[ip] = 0;
2065  iopf[ip] = 0;
2066  status = 0;
2067  bndcnt = 0;
2068  itercnt = 0;
2069  rrsdiff[ip] = badval;
2070  chisqr[ip] = badval;
2071 
2072  aph_s[ip] = badval;
2073  adg_s[ip] = badval;
2074  bbp_s[ip] = badval;
2075  siop_num[ip] = BAD_INT;
2076 
2077  for (ipar = 0; ipar < g->npar; ipar++)
2078  fit_par[ip][ipar] = badval;
2079 
2080 
2081  // flag and skip if pixel already masked
2082 
2083  if (l1rec->mask[ip]) {
2084  iopf[ip] |= IOPF_ISMASKED;
2085  iopf[ip] |= IOPF_FAILED;
2086  continue;
2087  }
2088 
2089  // set values that depend on sea state
2090 
2091  for (iw = 0; iw < nwave; iw++) {
2092  aw [iw] = l1rec->sw_a [ipb2 + iw];
2093  bbw[iw] = l1rec->sw_bb[ipb2 + iw];
2094  }
2095  for (iw = 0; iw < g->nwave; iw++) {
2096  g->aw [iw] = aw [g->bindx[iw]];
2097  g->bbw[iw] = bbw [g->bindx[iw]];
2098  }
2099 
2100  //
2101  // set dynamic, non-optimized model parameters
2102  //
2103 
2104  aph_s[ip] = g->aph_s;
2105  adg_s[ip] = g->adg_s;
2106  bbp_s[ip] = g->bbp_s;
2107 
2108  // get starting chlorophyll from default algorithm
2109 
2110  g->chl = MAX(MIN(l2rec->chl[ip], chl_max), chl_min);
2111 
2112  // Rrs function
2113 
2114  switch (g->rrs_opt) {
2115  case RRSFOQ:
2116  foqint_morel(input->fqfile, wave, nwave, 0.0, 0.0, 0.0, g->chl, foq);
2117  for (iw = 0; iw < g->nwave; iw++) {
2118  g->foq[iw] = foq[g->bindx[iw]];
2119  }
2120  break;
2121  }
2122 
2123  // aph function
2124 
2125  switch (g->aph_opt) {
2126  case APHBRICAUD:
2127  // replaces default tabbed values with Bricaud chl-based values
2128  for (iw = 0; iw < nwave; iw++) {
2129  g->aph_tab_w[iw] = wave[iw];
2130  // get the aph* normalization factor - the 0.055 is a "typical"
2131  // aph* value for 443nm
2132  aph_norm = 0.055 / get_aphstar(443., BANDW, APHBRICAUD, g->chl);
2133  g->aph_tab_s[0][iw] = aph_norm * get_aphstar(wave[iw], BANDW, APHBRICAUD, g->chl);
2134  }
2135  g->aph_tab_nw = nwave;
2136  g->aph_nvec = 1;
2137  break;
2138  case APHCIOTTI:
2139  // replaces default tabbed values with Ciotti size-fraction-based values
2140  for (iw = 0; iw < nwave; iw++) {
2141  g->aph_tab_w[iw] = wave[iw];
2142  g->aph_tab_s[0][iw] = get_aphstar(wave[iw], BANDW, APHCIOTTI, g->aph_s);
2143  }
2144  g->aph_tab_nw = nwave;
2145  g->aph_nvec = 1;
2146  break;
2147  }
2148 
2149  // adg function
2150 
2151  switch (g->adg_opt) {
2152  case ADGSQAA:
2153  // update exponential based on QAA band-ratio relationship
2154  i1 = windex(443.0, wave, nwave);
2155  i2 = windex(550.0, wave, nwave);
2156  Rrs1 = l2rec->Rrs[ipb2 + i1] - l2rec->Rrs_raman[ipb2 + i1];
2157  Rrs2 = l2rec->Rrs[ipb2 + i2] - l2rec->Rrs_raman[ipb2 + i2];
2158  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2159  Rrs1 = rrs_above_to_below(Rrs1);
2160  Rrs2 = rrs_above_to_below(Rrs2);
2161  g->adg_s = MAX(MIN(0.015 + 0.002 / (0.6 + Rrs1 / Rrs2), adg_s_max), adg_s_min);
2162  } else if (Rrs2 > 0.0) {
2163  g->adg_s = adg_s_min;
2164  } else {
2165  g->adg_s = BAD_FLT;
2166  iopf[ip] |= IOPF_BADRRS;
2167  status = 1;
2168  }
2169  break;
2170  case ADGSOBPG:
2171  // update exponential based on OBPG band-ratio relationship
2172  i1 = windex(412.0, wave, nwave);
2173  i2 = windex(550.0, wave, nwave);
2174  Rrs1 = l2rec->Rrs[ipb2 + i1] - l2rec->Rrs_raman[ipb2 + i1];
2175  Rrs2 = l2rec->Rrs[ipb2 + i2] - l2rec->Rrs_raman[ipb2 + i2];
2176  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2177  g->adg_s = MAX(MIN(0.015 + 0.0038 * log10(Rrs1 / Rrs2), adg_s_max), adg_s_min);
2178  } else if (Rrs2 > 0.0) {
2179  g->adg_s = adg_s_min;
2180  } else {
2181  g->adg_s = BAD_FLT;
2182  iopf[ip] |= IOPF_BADRRS;
2183  status = 1;
2184  }
2185  break;
2186  }
2187 
2188  // bbp function
2189 
2190  switch (g->bbp_opt) {
2191  case BBPSHAL:
2192  // update power-law exponent based on HAL band-ratio relationship
2193  i1 = windex(490.0, wave, nwave);
2194  i2 = windex(550.0, wave, nwave);
2195  Rrs1 = l2rec->Rrs[ipb2 + i1] - l2rec->Rrs_raman[ipb2 + i1];
2196  Rrs2 = l2rec->Rrs[ipb2 + i2] - l2rec->Rrs_raman[ipb2 + i2];
2197  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2198  g->bbp_s = MAX(MIN(0.8 * (Rrs1 / Rrs2) + 0.2, bbp_s_max), bbp_s_min);
2199  } else if (Rrs2 > 0.0) {
2200  g->bbp_s = bbp_s_min;
2201  } else {
2202  g->bbp_s = BAD_FLT;
2203  iopf[ip] |= IOPF_BADRRS;
2204  status = 1;
2205  }
2206  break;
2207  case BBPSQAA:
2208  // update power-law exponent based on QAA band-ratio relationship
2209  i1 = windex(443.0, wave, nwave);
2210  i2 = windex(550.0, wave, nwave);
2211  Rrs1 = l2rec->Rrs[ipb2 + i1] - l2rec->Rrs_raman[ipb2 + i1];
2212  Rrs2 = l2rec->Rrs[ipb2 + i2] - l2rec->Rrs_raman[ipb2 + i2];
2213  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2214  Rrs1 = rrs_above_to_below(Rrs1);
2215  Rrs2 = rrs_above_to_below(Rrs2);
2216  g->bbp_s = MAX(MIN(2.0 * (1.0 - 1.2 * exp(-0.9 * Rrs1 / Rrs2)), bbp_s_max), bbp_s_min);
2217  } else if (Rrs2 > 0.0) {
2218  g->bbp_s = bbp_s_min;
2219  } else {
2220  g->bbp_s = BAD_FLT;
2221  iopf[ip] |= IOPF_BADRRS;
2222  status = 1;
2223  }
2224  break;
2225  case BBPSCIOTTI:
2226  // update power-law exponent based on Ciotti chl relationship
2227  if (l2rec->chl[ip] > 0.0) {
2228  g->bbp_s = MAX(MIN(1.0 - 0.768 * log10(g->chl), bbp_s_max), bbp_s_min);
2229  } else {
2230  g->bbp_s = BAD_FLT;
2231  iopf[ip] |= IOPF_BADRRS;
2232  status = 1;
2233  }
2234  break;
2235  case BBPSMM01:
2236  // update power-law exponent based on Morel chl relationship
2237  if (l2rec->chl[ip] > 0.0) {
2238  g->bbp_s = MAX(MIN(0.5 * (0.3 - log10(g->chl)), bbp_s_max), bbp_s_min);
2239  } else {
2240  g->bbp_s = BAD_FLT;
2241  iopf[ip] |= IOPF_BADRRS;
2242  status = 1;
2243  }
2244  break;
2245  case BBPSLAS:
2246  // update power-law exponent based on Loisel & Stramski model
2247  if ((g->bbp_s = get_bbp_las_eta(l2rec, ip)) == BAD_FLT) {
2248  iopf[ip] |= IOPF_BADRRS;
2249  status = 1;
2250  }
2251  g->bbp_nvec = 1;
2252  break;
2253  case BBPLAS:
2254  case BBPLASFIX:
2255  // replaces default tabbed values with results from Loisel & Stramski model
2256  if (get_bbp_las(l2rec, ip, g->bbp_tab_w, g->bbp_tab_s[0], g->bbp_tab_nw) == 0) {
2257  iopf[ip] |= IOPF_BADRRS;
2258  status = 1;
2259  }
2260  g->bbp_nvec = 1;
2261  break;
2262  case BBPQAAFIX:
2263  // replaces default tabbed values with results from QAA model
2264  if (get_bbp_qaa(l2rec, ip, g->bbp_tab_w, g->bbp_tab_s[0], g->bbp_tab_nw) == 0) {
2265  iopf[ip] |= IOPF_BADRRS;
2266  status = 1;
2267  }
2268  g->bbp_nvec = 1;
2269  break;
2270  }
2271 
2272  // save updated model config, flag and skip if unable to compute
2273 
2274  aph_s[ip] = g->aph_s;
2275  adg_s[ip] = g->adg_s;
2276  bbp_s[ip] = g->bbp_s;
2277 
2278  if (status != 0) {
2279  iopf[ip] |= IOPF_BADRRS;
2280  continue;
2281  }
2282 
2283  // convert to subsurface reflectance
2284 
2285  for (iw = 0; iw < g->nwave; iw++) {
2286  ib = g->bindx[iw];
2287  Rrs_a[iw] = l2rec->Rrs[ipb2 + ib] - l2rec->Rrs_raman[ipb2 + ib];
2288  Rrs_b[iw] = rrs_above_to_below(Rrs_a[iw]);
2289  wts [iw] = g->wts[iw]; // /rrs_above_to_below(Rrs_a[iw])*sqrt(l2rec->nobs[ip]);
2290  if (Rrs_b[iw] > 0.0) {
2291  bndcnt++;
2292  }
2293  }
2294 
2295  // if less than npar valid reflectances, flag and skip pixel */
2296 
2297  if (bndcnt < g->npar) {
2298  iopf[ip] |= IOPF_BADRRS;
2299  continue;
2300  }
2301 
2302  // initialize model parameters
2303 
2304  giop_ctl_start(g, g->chl);
2305 
2306  for (ipar = 0; ipar < g->npar; ipar++) {
2307  par[ipar] = g->par[ipar]; // model params
2308  par[ipar + g->npar] = 0.0; // model param errors
2309  }
2310 
2311  // run model optimization for this pixel
2312 
2313  switch (g->fit_opt) {
2314  case AMOEBA:
2315  status = fit_giop_amb(g, Rrs_b, wts, par, Rrs_f, &itercnt);
2316  break;
2317  case LEVMARQ:
2318  status = fit_giop_lm(g, Rrs_b, wts, par, &chi, &itercnt);
2319  break;
2320  case SVDFIT:
2321  status = fit_giop_svd(g, Rrs_b, wts, par);
2322  break;
2323  case SVDSIOP:
2324  status = fit_giop_svd_siop(g, Rrs_b, wts, par, &chi);
2325  /*If an optimal svd_siop solution not reached*/
2326  if (status == -99) {
2327  siop_num[ip] = -1;
2328  }
2329  break;
2330  default:
2331  printf("%s Line %d: Unknown optimization method for GIOP %d\n",
2332  __FILE__, __LINE__, g->fit_opt);
2333  exit(1);
2334  break;
2335  }
2336 
2337 
2338  // if inversion failed, flag and skip
2339 
2340  if (status != 0) {
2341  iopf[ip] |= IOPF_FAILED;
2342  if (itercnt >= g->maxiter)
2343  iopf[ip] |= IOPF_MAXITER;
2344  continue;
2345  }
2346 
2347 
2348  /* save final params */
2349  for (ipar = 0; ipar < g->npar; ipar++) {
2350  if (isfinite(par[ipar]))
2351  fit_par[ip][ipar] = par[ipar];
2352  }
2353 
2354  chisqr[ip] = chi;
2355 
2356  // evaluate model at fitted bands
2357 
2358  switch (g->fit_opt) {
2359  case SVDSIOP:
2360  giop_model_iterate(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, acdom1, anap1, bbph1, bbnap1, Rrs_f, NULL, NULL);
2361  break;
2362  default:
2363  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, Rrs_f, NULL, NULL);
2364  break;
2365  }
2366 
2367 
2368  // bogus evaluation, flag and skip
2369 
2370  for (iw = 0; iw < g->nwave; iw++) {
2371  if (!isfinite(Rrs_f[iw])) {
2372  iopf[ip] |= IOPF_NAN;
2373  break;
2374  }
2375  }
2376  if (iopf[ip] & IOPF_NAN) {
2377  iopf[ip] |= IOPF_FAILED;
2378  continue;
2379  }
2380 
2381  // check goodness of fit
2382 
2383  rrsdiff[ip] = 0.0;
2384  for (iw = 0; iw < g->nwave; iw++) {
2385  if (g->wave[iw] >= 400 && g->wave[iw] <= 600) {
2386  if (fabs(Rrs_b[iw]) > 1e-7 && fabs(Rrs_f[iw] - Rrs_b[iw]) > 1e-5)
2387  rrsdiff[ip] += fabs(Rrs_f[iw] - Rrs_b[iw]) / fabs(Rrs_b[iw]);
2388  }
2389  }
2390  rrsdiff[ip] /= g->nwave;
2391  if (rrsdiff[ip] > input->giop_rrs_diff) {
2392  iopf[ip] |= IOPF_RRSDIFF;
2393  }
2394 
2395 
2396  // store in static globals
2397 
2398  switch (g->fit_opt) {
2399  case SVDSIOP:
2400  giop_model_iterate(g, par, nwave, wave, aw, bbw, foq, aph1, adg1, bbp1, acdom1, anap1, bbph1, bbnap1, Rrs_f, NULL, NULL);
2401 
2402  mRrs[ipb + ib] = rrs_below_to_above((Rrs_f[ib] + l2rec->Rrs_raman[ipb2 + ib]));
2403 
2404  for (ib = 0; ib < nwave; ib++) {
2405 
2406  if (isfinite(aph1[ib])) {
2407  aph[ipb + ib] = aph1[ib];
2408  aph[ierr + ib] = aph1[ib + nwave];
2409  } else {
2410  iopf[ip] |= IOPF_NAN;
2411  }
2412  if (isfinite(acdom1[ib])) {
2413  acdom[ipb + ib] = acdom1[ib];
2414  acdom[ierr + ib] = acdom1[ib + nwave];
2415  } else {
2416  iopf[ip] |= IOPF_NAN;
2417  }
2418  if (isfinite(anap1[ib])) {
2419  anap[ipb + ib] = anap1[ib];
2420  anap[ierr + ib] = anap1[ib + nwave];
2421  } else {
2422  iopf[ip] |= IOPF_NAN;
2423  }
2424  if (isfinite(aph1[ib]) && isfinite(acdom1[ib]) && isfinite(anap1[ib])) {
2425  adg[ipb + ib] = acdom1[ib] + anap1[ib];
2426  adg[ierr + ib] = adg1[ib + nwave] + acdom1[ib + nwave];
2427  a[ipb + ib] = aw[ib] + aph1[ib] + acdom1[ib] + anap1[ib];
2428  a[ierr + ib] = aph1[ib + nwave] + adg1[ib + nwave] + acdom1[ib + nwave];
2429  } else {
2430  iopf[ip] |= IOPF_NAN;
2431  }
2432  if (isfinite(bbnap1[ib])) {
2433  bbnap [ipb + ib] = bbnap1[ib];
2434  bbnap [ierr + ib] = bbnap1[ib + nwave];
2435  } else {
2436  iopf[ip] |= IOPF_NAN;
2437  }
2438  if (isfinite(bbph1[ib])) {
2439  bbph [ipb + ib] = bbph1[ib];
2440  bbph [ierr + ib] = bbph1[ib + nwave];
2441  } else {
2442  iopf[ip] |= IOPF_NAN;
2443  }
2444  if (isfinite(bbph1[ib]) && isfinite(bbnap1[ib])) {
2445  bbp[ipb + ib] = bbph1[ib] + bbnap1[ib];
2446  bbp[ierr + ib] = bbph1[ib + nwave] + bbnap1[ib + nwave];
2447  bb [ipb + ib] = bbw[ib] + bbph1[ib] + bbnap1[ib];
2448  bb [ierr + ib] = bbw[ib + nwave] + bbph1[ib + nwave] + bbnap[ib + nwave];
2449  } else {
2450  iopf[ip] |= IOPF_NAN;
2451  }
2452  }
2453 
2454  /* check IOP ranges and set flags */
2455  set_iop_flag(wave, nwave, &a[ipb], &aph[ipb], &adg[ipb],
2456  &bb[ipb], &bbp[ipb], &iopf[ip]);
2457 
2458  /*Optimal SIOP combination (for iterative aLMI solution)*/
2459  siop_num[ip] = 1 + g->siopIdx;
2460 
2461  /* aLMI compute chlorophyll */
2462  chl [ip] = par[0];
2463 
2464  break;
2465 
2466  default:
2467  //Fit model to full visible bandset
2468  giop_model(g, par, nwave, wave, aw, bbw, foq, aph1, adg1, bbp1, Rrs_f, NULL, NULL);
2469 
2470  for (ib = 0; ib < nwave; ib++) {
2471 
2472  mRrs[ipb + ib] = rrs_below_to_above((Rrs_f[ib] + l2rec->Rrs_raman[ipb2 + ib]));
2473 
2474  if (isfinite(aph1[ib])) {
2475  aph[ipb + ib] = aph1[ib];
2476  aph[ierr + ib] = aph1[ib + nwave];
2477  } else {
2478  iopf[ip] |= IOPF_NAN;
2479  }
2480  if (isfinite(adg1[ib])) {
2481  adg[ipb + ib] = adg1[ib];
2482  adg[ierr + ib] = adg1[ib + nwave];
2483  } else {
2484  iopf[ip] |= IOPF_NAN;
2485  }
2486  if (isfinite(aph1[ib]) && isfinite(adg1[ib])) {
2487  a[ipb + ib] = aw[ib] + aph1[ib] + adg1[ib];
2488  a[ierr + ib] = aph1[ib + nwave] + adg1[ib + nwave];
2489 
2490  } else {
2491  iopf[ip] |= IOPF_NAN;
2492  }
2493  if (isfinite(bbp1[ib])) {
2494  bbp[ipb + ib] = bbp1[ib];
2495  bbp[ierr + ib] = bbp1[ib + nwave];
2496  bb [ipb + ib] = bbw[ib] + bbp1[ib];
2497  bb [ierr + ib] = bbp1[ib + nwave];
2498  } else {
2499  iopf[ip] |= IOPF_NAN;
2500  }
2501  }
2502  // check IOP ranges and set flags
2503 
2504  set_iop_flag(wave, nwave, &a[ipb], &aph[ipb], &adg[ipb],
2505  &bb[ipb], &bbp[ipb], &iopf[ip]);
2506 
2507  // compute chlorophyll
2508 
2509  chl [ip] = giop_chl(g, iopf[ip], par, &chl[ip + l1rec->npix]);
2510  break;
2511  }
2512 
2513 
2514  iter[ip] = itercnt;
2515  }
2516 
2517 
2518 
2519  // fail pixels where any flags were set
2520 
2521  for (ip = 0; ip < l1rec->npix; ip++)
2522  if (iopf[ip] != 0) l1rec->flags[ip] |= PRODFAIL;
2523 
2524 
2525  LastRecNum = l1rec->iscan;
2526  free(Rrs_a);
2527  free(Rrs_b);
2528  free(Rrs_f);
2529  free(wts);
2530  free(par);
2531 }
2532 
2533 
2534 /* ------------------------------------------------------------------- */
2535 /* get_giop() - returns requested GIOP product for full scanline */
2536 
2537 /* ------------------------------------------------------------------- */
2538 void get_giop(l2str *l2rec, l2prodstr *p, float prod[]) {
2539  int prodID = p->cat_ix;
2540  int ib = p->prod_ix;
2541 
2542  int32_t ip, ipb, ierr;
2543  l1str *l1rec = l2rec->l1rec;
2544 
2545  /*Before running GIOP, check if valid output products selected*/
2546  /*Note: aCDOM, aNAP, bbPh, bbNAP are only valid choices for giop_fit_opt SVDSIOP*/
2547  switch (input->giop_fit_opt) {
2548  case SVDSIOP:
2549  break;
2550  default:
2551  switch (prodID) {
2552  case CAT_acdom_giop:
2553  case CAT_anap_giop:
2554  case CAT_bbph_giop:
2555  case CAT_bbnap_giop:
2556  case CAT_acdom_unc_giop:
2557  case CAT_anap_unc_giop:
2558  case CAT_bbph_unc_giop:
2559  case CAT_bbnap_unc_giop:
2560  printf("-E- %s line %d : products acdom, anap, bbph and bbnap are only applicable with giop_fit_opt=SVDSIOP.\n",
2561  __FILE__, __LINE__);
2562  exit(1);
2563  break;
2564  }
2565  }
2566 
2567  if (!giop_ran(l1rec->iscan))
2568  run_giop(l2rec);
2569 
2570 
2571  for (ip = 0; ip < l1rec->npix; ip++) {
2572 
2573  ipb = ip * l1rec->l1file->nbands + ib;
2574  ierr = l1rec->npix * l1rec->l1file->nbands + ipb;
2575 
2576  switch (prodID) {
2577 
2578  case CAT_mRrs_giop:
2579  prod[ip] = (float) mRrs[ipb];
2580  break;
2581 
2582  case CAT_aph_giop:
2583  prod[ip] = (float) aph[ipb];
2584  break;
2585 
2586  case CAT_adg_giop:
2587  prod[ip] = (float) adg[ipb];
2588  break;
2589 
2590  case CAT_bbp_giop:
2591  prod[ip] = (float) bbp[ipb];
2592  break;
2593 
2594  case CAT_a_giop:
2595  prod[ip] = (float) a[ipb];
2596  break;
2597 
2598  case CAT_bb_giop:
2599  prod[ip] = (float) bb[ipb];
2600  break;
2601 
2602  case CAT_acdom_giop:
2603  prod[ip] = (float) acdom[ipb];
2604  break;
2605 
2606  case CAT_anap_giop:
2607  prod[ip] = (float) anap[ipb];
2608  break;
2609 
2610  case CAT_bbph_giop:
2611  prod[ip] = (float) bbph[ipb];
2612  break;
2613 
2614  case CAT_bbnap_giop:
2615  prod[ip] = (float) bbnap[ipb];
2616  break;
2617 
2618  case CAT_chl_giop:
2619  prod[ip] = (float) chl[ip];
2620  break;
2621 
2622  case CAT_opt_siop_giop:
2623  prod[ip] = (int) siop_num[ip];
2624  break;
2625 
2626  case CAT_aph_unc_giop:
2627  prod[ip] = (float) aph[ierr];
2628  break;
2629 
2630  case CAT_adg_unc_giop:
2631  prod[ip] = (float) adg[ierr];
2632  break;
2633 
2634  case CAT_acdom_unc_giop:
2635  prod[ip] = (float) adg[ierr];
2636  break;
2637 
2638  case CAT_anap_unc_giop:
2639  prod[ip] = (float) adg[ierr];
2640  break;
2641 
2642  case CAT_bbp_unc_giop:
2643  prod[ip] = (float) bbp[ierr];
2644  break;
2645 
2646  case CAT_bbph_unc_giop:
2647  prod[ip] = (float) bbph[ierr];
2648  break;
2649 
2650  case CAT_bbnap_unc_giop:
2651  prod[ip] = (float) bbnap[ierr];
2652  break;
2653 
2654  case CAT_a_unc_giop:
2655  prod[ip] = (float) a[ierr];
2656  break;
2657 
2658  case CAT_bb_unc_giop:
2659  prod[ip] = (float) bb[ierr];
2660  break;
2661 
2662  case CAT_chl_unc_giop:
2663  prod[ip] = (float) chl[ip + l1rec->npix];
2664  break;
2665 
2666  case CAT_aphs_giop:
2667  prod[ip] = (float) aph_s[ip];
2668  break;
2669 
2670  case CAT_adgs_giop:
2671  prod[ip] = (float) adg_s[ip];
2672  break;
2673 
2674  case CAT_bbps_giop:
2675  prod[ip] = (float) bbp_s[ip];
2676  break;
2677 
2678  case CAT_rrsdiff_giop:
2679  prod[ip] = (float) rrsdiff[ip];
2680  break;
2681 
2682  case CAT_chisqr_giop:
2683  prod[ip] = (float) chisqr[ip];
2684  break;
2685 
2686  case CAT_fitpar_giop:
2687  if (ib >= max_npar) {
2688  printf("-E- %s line %d : output request for GIOP fit parameter %d exceeds number of fit parameters %d.\n",
2689  __FILE__, __LINE__, ib, max_npar);
2690  exit(1);
2691  }
2692  prod[ip] = (float) fit_par[ip][ib];
2693  break;
2694 
2695  default:
2696  printf("-E- %s line %d : erroneous product ID %d passed to GIOP.\n",
2697  __FILE__, __LINE__, prodID);
2698  exit(1);
2699  }
2700  }
2701 
2702  return;
2703 }
2704 
2705 
2706 /* ------------------------------------------------------------------- */
2707 /* get_iter_giop() - returns iteration count */
2708 
2709 /* ------------------------------------------------------------------- */
2710 int16 *get_iter_giop(l2str *l2rec) {
2711  if (!giop_ran(l2rec->l1rec->iscan))
2712  run_giop(l2rec);
2713 
2714  return iter;
2715 }
2716 
2717 
2718 /* ------------------------------------------------------------------- */
2719 /* get_flags_giop() - returns iteration count */
2720 
2721 /* ------------------------------------------------------------------- */
2722 int16 *get_flags_giop(l2str *l2rec) {
2723  if (!giop_ran(l2rec->l1rec->iscan))
2724  run_giop(l2rec);
2725 
2726  return iopf;
2727 }
2728 
2729 
2730 /* ------------------------------------------------------------------- */
2731 /* Interface to convl12() to return GIOP iops */
2732 
2733 /* ------------------------------------------------------------------- */
2734 void iops_giop(l2str *l2rec) {
2735  int32_t ib, ip, ipb, ipb2;
2736 
2737  int32_t nbands = l2rec->l1rec->l1file->nbands;
2738  int32_t npix = l2rec->l1rec->npix;
2739 
2740  if (!giop_ran(l2rec->l1rec->iscan))
2741  run_giop(l2rec);
2742 
2743  for (ip = 0; ip < npix; ip++) {
2744  for (ib = 0; ib < nbands; ib++) {
2745  ipb2 = ip * nbands + ib;
2746  ipb = ip * nbands + ib;
2747  l2rec->a [ipb2] = (float) a[ipb];
2748  l2rec->bb[ipb2] = (float) bb[ipb];
2749  }
2750  }
2751 
2752  return;
2753 }
2754 
2755 
2756 // PJW 9 Jan 2013
2757 /* ------------------------------------------------------------------- */
2758 /* give external access to local *chl */
2759 
2760 /* ------------------------------------------------------------------- */
2761 float* giop_get_chl_pointer() {
2762  return chl;
2763 }
2764 /* ------------------------------------------------------------------- */
2765 /* give external access to local *adg */
2766 
2767 /* ------------------------------------------------------------------- */
2768 float* giop_get_adg_pointer() {
2769  return adg;
2770 }
2771 /* ------------------------------------------------------------------- */
2772 /* give external access to local *bbp */
2773 
2774 /* ------------------------------------------------------------------- */
2775 float* giop_get_bbp_pointer() {
2776  return bbp;
2777 }
2778 /* ------------------------------------------------------------------- */
2779 /* give external access to local *aph */
2780 
2781 /* ------------------------------------------------------------------- */
2782 float* giop_get_aph_pointer() {
2783  return aph;
2784 }
2785 /* ------------------------------------------------------------------- */
2786 /* give external access to local **fit_par */
2787 
2788 /* ------------------------------------------------------------------- */
2789 float** giop_get_fitpar_pointer() {
2790  return fit_par;
2791 }
2792 
2793 /* ------------------------------------------------------------------- */
2794 /* give external access to local *bbp_s */
2795 
2796 /* ------------------------------------------------------------------- */
2797 float* giop_get_bbp_s_pointer() {
2798  return bbp_s;
2799 }
2800 
freeArray
void freeArray(void **a, int32_t m)
Definition: giop.c:100
ANAPTAB
#define ANAPTAB
Definition: giop.h:43
MAX
#define MAX(A, B)
Definition: swl0_utils.h:26
giop_ctl_init
void giop_ctl_init(giopstr *g, int nwave, float wave[], float aw[], float bbw[])
Definition: giop.c:215
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
CAT_adgs_giop
#define CAT_adgs_giop
Definition: l2prod.h:205
MIN
#define MIN(x, y)
Definition: rice.h:169
t
data_t t[NROOTS+1]
Definition: decode_rs.h:77
table_read_r4
float * table_read_r4(char *input_filename, int m, int n)
Definition: table_io.cpp:2540
g
e g
Definition: HOWTO_Add_a_product.txt:5
j
int j
Definition: decode_rs.h:73
get_aphstar
float get_aphstar(float wave, int dwave, int ftype, float proxy)
Definition: aph.c:296
BBPSCIOTTI
#define BBPSCIOTTI
Definition: giop.h:24
status
int status
Definition: l1_czcs_hdf.c:32
CAT_mRrs_giop
#define CAT_mRrs_giop
Definition: l2prod.h:212
ADGSOBPG
#define ADGSOBPG
Definition: giop.h:15
table_column_count
int table_column_count(char *input_filename)
Definition: table_io.cpp:2524
CAT_chl_giop
#define CAT_chl_giop
Definition: l2prod.h:197
CAT_bbp_unc_giop
#define CAT_bbp_unc_giop
Definition: l2prod.h:200
FITSTRUCT::niter
int niter
Definition: amoeba.h:9
fit_giop_svd
int fit_giop_svd(giopstr *g, double rrs[], double wts[], double par[])
Definition: giop.c:1207
IOPF_ISMASKED
#define IOPF_ISMASKED
Definition: flags_iop.h:8
CAT_aphs_giop
#define CAT_aphs_giop
Definition: l2prod.h:204
table
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT table(e.g. m1) to be used for interpolation. The table values will be linearly interpolated using the tables corresponding to the node times bracketing the granule time. If the granule time falls before the time of the first node or after the time of the last node
SVDFIT
#define SVDFIT
Definition: giop.h:8
RRSGRD
#define RRSGRD
Definition: giop.h:56
get_bbp_las_eta
float get_bbp_las_eta(l2str *l2rec, int ip)
Definition: las_iop.c:627
allocateMemory
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
CAT_adg_giop
#define CAT_adg_giop
Definition: l2prod.h:196
CAT_aph_giop
#define CAT_aph_giop
Definition: l2prod.h:195
BBPSLAS
#define BBPSLAS
Definition: giop.h:26
NULL
#define NULL
Definition: decode_rs.h:63
BBPTAB
#define BBPTAB
Definition: giop.h:19
giop_get_aph_pointer
float * giop_get_aph_pointer()
Definition: giop.c:2782
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
BBPSMM01
#define BBPSMM01
Definition: giop.h:25
BBPQAAFIX
#define BBPQAAFIX
Definition: giop.h:29
giop_int_tab_file
void giop_int_tab_file(char *file, int nx, float *x, float **y)
Definition: giop.c:121
fit_giop_amb
int fit_giop_amb(giopstr *g, double Rrs[], double wts[], double par[], double Rrs_fit[], int16 *itercnt)
Definition: giop.c:975
IOPF_MAXITER
#define IOPF_MAXITER
Definition: flags_iop.h:10
CAT_bbp_giop
#define CAT_bbp_giop
Definition: l2prod.h:194
CAT_rrsdiff_giop
#define CAT_rrsdiff_giop
Definition: l2prod.h:208
CAT_acdom_giop
#define CAT_acdom_giop
Definition: l2prod.h:356
BBNAPTAB
#define BBNAPTAB
Definition: giop.h:51
IOPF_RRSDIFF
#define IOPF_RRSDIFF
Definition: flags_iop.h:13
giop_lm_fdf
int giop_lm_fdf(const gsl_vector *parv, void *data, gsl_vector *f, gsl_matrix *J)
Definition: giop.c:1038
SVDSIOP
#define SVDSIOP
Definition: giop.h:9
CAT_anap_giop
#define CAT_anap_giop
Definition: l2prod.h:357
ACDOMNONE
#define ACDOMNONE
Definition: giop.h:40
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
giop_get_adg_pointer
float * giop_get_adg_pointer()
Definition: giop.c:2768
CAT_bbps_giop
#define CAT_bbps_giop
Definition: l2prod.h:206
MDN.metrics.rmse
def rmse(y, y_hat)
Definition: metrics.py:75
CAT_anap_unc_giop
#define CAT_anap_unc_giop
Definition: l2prod.h:361
amoeba
short amoeba(double *pnts, FITSTRUCT *auxdata, double(*func)(FITSTRUCT *, double[]), float tol)
Definition: amoeba.c:14
table_row_count
int table_row_count(char *input_filename)
Definition: table_io.cpp:2528
giop_get_chl_pointer
float * giop_get_chl_pointer()
Definition: giop.c:2761
CAT_bb_giop
#define CAT_bb_giop
Definition: l2prod.h:193
rrs_above_to_below
float rrs_above_to_below(float Rrs)
Definition: giop.c:1598
FITSTRUCT::merit
double merit
Definition: amoeba.h:4
file
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 file
Definition: HISTORY.txt:413
input
instr * input
Definition: main_l2binmatch.cpp:27
IOPF_NAN
#define IOPF_NAN
Definition: flags_iop.h:12
BBNAPNONE
#define BBNAPNONE
Definition: giop.h:52
CAT_aph_unc_giop
#define CAT_aph_unc_giop
Definition: l2prod.h:201
APHGAUSS
#define APHGAUSS
Definition: giop.h:33
table_free_r4
void table_free_r4(float *table)
Definition: table_io.cpp:2552
l12_proto.h
fit_data_str::g
giopstr * g
Definition: giop.c:45
giop_amb
double giop_amb(FITSTRUCT *ambdata, double par[])
Definition: giop.c:955
ADGSQAA
#define ADGSQAA
Definition: giop.h:14
fit_data_str::w
double * w
Definition: giop.c:44
init
int init(int32_t ipr, int32_t jpr, char *efile, char *pfile)
Definition: proj_report.c:51
J
int J
Definition: Usds.c:62
giop_ran
int giop_ran(int recnum)
Definition: giop.c:569
rrs_below_to_above
float rrs_below_to_above(float rrs_s)
Definition: giop.c:1607
ADGTAB
#define ADGTAB
Definition: giop.h:12
f
double precision function f(R1)
Definition: tmd.lp.f:1454
recnum
read recnum
Definition: HOWTO_Add_a_sensor.txt:72
CAT_chisqr_giop
#define CAT_chisqr_giop
Definition: l2prod.h:209
linterp
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
APHTAB
#define APHTAB
Definition: giop.h:32
CAT_a_unc_giop
#define CAT_a_unc_giop
Definition: l2prod.h:198
AMOEBA
#define AMOEBA
Definition: giop.h:6
CAT_fitpar_giop
#define CAT_fitpar_giop
Definition: l2prod.h:210
hico.value_locate.x
list x
Definition: value_locate.py:64
FITSTRUCT::wgt
double * wgt
Definition: amoeba.h:6
color_dtdb.y
y
Definition: color_dtdb.py:430
CAT_chl_unc_giop
#define CAT_chl_unc_giop
Definition: l2prod.h:203
func
subroutine func(x, conec, n, bconecno, bn, units, u, inno, i, outno, o, Input, Targ, p, sqerr)
Definition: ffnet.f:287
giop_model_iterate
void giop_model_iterate(giopstr *g, double par[], int nwave, float wave[], float aw[], float bbw[], float foq[], float aph[], float adg[], float bbp[], float acdom[], float anap[], float bbph[], float bbnap[], double rrs[], double **dfdpar, double **parstar)
Definition: giop.c:792
fit_giop_lm
int fit_giop_lm(giopstr *g, double Rrs[], double wts[], double par[], double *chi, int16 *itercnt)
Definition: giop.c:1123
tol
float tol
Definition: pml_iop_tables.c:45
create_images.Rrs
def Rrs
Definition: create_images.py:92
CAT_opt_siop_giop
#define CAT_opt_siop_giop
Definition: l2prod.h:364
get_iter_giop
int16 * get_iter_giop(l2str *l2rec)
Definition: giop.c:2710
BBPHNONE
#define BBPHNONE
Definition: giop.h:48
FITSTRUCT::y
double * y
Definition: amoeba.h:5
CAT_bb_unc_giop
#define CAT_bb_unc_giop
Definition: l2prod.h:199
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
BBPLAS
#define BBPLAS
Definition: giop.h:27
foqint_morel
void foqint_morel(char *file, float wave[], int32_t nwave, float solz, float senzp, float phi, float chl, float brdf[])
Definition: brdf.c:314
giop_get_bbp_pointer
float * giop_get_bbp_pointer()
Definition: giop.c:2775
giop_ctl_start
void giop_ctl_start(giopstr *g, float chl)
Definition: giop.c:160
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
FITSTRUCT::nfunc
short int nfunc
Definition: amoeba.h:2
set_iop_flag
void set_iop_flag(float wave[], int32_t nwave, float a[], float aph[], float adg[], float bb[], float32 bbp[], int16_t *flag)
Definition: flags_iop.c:152
CAT_bbph_unc_giop
#define CAT_bbph_unc_giop
Definition: l2prod.h:362
FITSTRUCT::yfit
double * yfit
Definition: amoeba.h:8
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
CAT_adg_unc_giop
#define CAT_adg_unc_giop
Definition: l2prod.h:202
giop_lm_df
int giop_lm_df(const gsl_vector *parv, void *data, gsl_matrix *J)
Definition: giop.c:1113
APHCIOTTI
#define APHCIOTTI
Definition: giop.h:35
giop_chl
float32 giop_chl(giopstr *g, int16 iopf, double par[], float *chl_err)
Definition: giop.c:1569
RRSFOQ
#define RRSFOQ
Definition: giop.h:57
freeDArray
void freeDArray(double **a, int32_t m)
Definition: giop.c:108
CAT_bbnap_unc_giop
#define CAT_bbnap_unc_giop
Definition: l2prod.h:363
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
u
data_t u
Definition: decode_rs.h:74
giop_model
void giop_model(giopstr *g, double par[], int nwave, float wave[], float aw[], float bbw[], float foq[], float aph[], float adg[], float bbp[], double rrs[], double **dfdpar, double **parstar)
Definition: giop.c:580
get_bbp_qaa
int get_bbp_qaa(l2str *l2rec, int ip, float tab_wave[], float tab_bbp[], int tab_nwave)
Definition: get_qaa.c:367
CAT_a_giop
#define CAT_a_giop
Definition: l2prod.h:192
LEVMARQ
#define LEVMARQ
Definition: giop.h:7
CAT_bbph_giop
#define CAT_bbph_giop
Definition: l2prod.h:358
fabs
#define fabs(a)
Definition: misc.h:93
APHBRICAUD
#define APHBRICAUD
Definition: giop.h:34
FITSTRUCT::meta
void * meta
Definition: amoeba.h:10
get_flags_giop
int16 * get_flags_giop(l2str *l2rec)
Definition: giop.c:2722
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
BBPLASFIX
#define BBPLASFIX
Definition: giop.h:28
BAD_INT
#define BAD_INT
Definition: genutils.h:23
get_bbp_las
int get_bbp_las(l2str *l2rec, int ip, float tab_wave[], float tab_bbp[], int tab_nwave)
Definition: las_iop.c:605
giop_get_fitpar_pointer
float ** giop_get_fitpar_pointer()
Definition: giop.c:2789
FITSTRUCT::npnts
short int npnts
Definition: amoeba.h:3
run_giop
void run_giop(l2str *l2rec)
Definition: giop.c:1616
s
data_t s[NROOTS]
Definition: decode_rs.h:75
IOPF_BADRRS
#define IOPF_BADRRS
Definition: flags_iop.h:11
allocate2d_float
float ** allocate2d_float(size_t h, size_t w)
Allocate a two-dimensional array of type float of a given size.
Definition: allocate2d.c:125
giop_lm_f
int giop_lm_f(const gsl_vector *parv, void *data, gsl_vector *f)
Definition: giop.c:1109
BANDW
#define BANDW
Definition: l1.h:52
fit_giop_svd_siop
int fit_giop_svd_siop(giopstr *g, double rrs[], double wts[], double par[], double *chi)
Definition: giop.c:1322
for
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
giop_get_bbp_s_pointer
float * giop_get_bbp_s_pointer()
Definition: giop.c:2797
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
IOPF_FAILED
#define IOPF_FAILED
Definition: flags_iop.h:9
CAT_acdom_unc_giop
#define CAT_acdom_unc_giop
Definition: l2prod.h:360
CAT_bbnap_giop
#define CAT_bbnap_giop
Definition: l2prod.h:359
BBPSQAA
#define BBPSQAA
Definition: giop.h:22
BBPSHAL
#define BBPSHAL
Definition: giop.h:21
i
int i
Definition: decode_rs.h:71
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
fit_data_str::y
double * y
Definition: giop.c:43
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
fit_data_str
Definition: giop.c:42
modis.modis_utils.cleanup
def cleanup(lut)
Definition: modis_utils.py:219
ANAPNONE
#define ANAPNONE
Definition: giop.h:44
ADGS
#define ADGS
Definition: giop.h:13
iops_giop
void iops_giop(l2str *l2rec)
Definition: giop.c:2734
npix
int npix
Definition: get_cmp.c:27
ACDOMTAB
#define ACDOMTAB
Definition: giop.h:39
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
BBPHTAB
#define BBPHTAB
Definition: giop.h:47
get_giop
void get_giop(l2str *l2rec, l2prodstr *p, float prod[])
Definition: giop.c:2538
BBPS
#define BBPS
Definition: giop.h:20


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