A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
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 /* (Oct 2023) */
16 /* Output parameter uncertainties computed per McKinna et al (2019) */
17 /* bbp line height model per McKinna et al (2021) */
18 /* bbp chl-based model per Huot et al (2008) */
19 /* */
20 /* Implementation: */
21 /* B. Franz, NASA/OBPG/SAIC, May 2008 */
22 /* */
23 /* Reference: */
24 /* Werdell et al. (2013) Generalized ocean color inversion model for */
25 /* retrieving marine inherent optical properties, Applied Optics, 52, */
26 /* 2019-2037. */
27 /* */
28 /* Notes: */
29 /* This code was written with the intent that it may become a base for */
30 /* multiple algorithms (to be writtenn as wrappers). As such, it can't*/
31 /* be assumed that the starting config is the only config (e.g., */
32 /* number of wavelengths to fit), which leads to some inefficiencies. */
33 /* =================================================================== */
34 
35 #include <stdio.h>
36 #include <math.h>
37 #include "l12_proto.h"
38 #include "giop.h"
39 #include "amoeba.h"
40 #include <gsl/gsl_errno.h>
41 #include <gsl/gsl_vector.h>
42 #include <gsl/gsl_blas.h>
43 #include <gsl/gsl_multifit_nlin.h>
44 #include <gsl/gsl_linalg.h>
45 #include <gsl/gsl_poly.h>
46 #include <gsl/gsl_matrix.h>
47 #include <gsl/gsl_permutation.h>
48 
49 
50 typedef double realtype;
51 
52 typedef struct fit_data_str {
53  double *y;
54  double *w;
55  giopstr *g;
56 } fitstr;
57 
58 static int32_t LastRecNum = -1;
59 static float badval = BAD_FLT;
60 static float bbp_s_max = 5.0;
61 static float bbp_s_min = 0.0;
62 static float adg_s_max = 0.025;
63 static float adg_s_min = 0.01;
64 static float chl_max = 10.0;
65 static float chl_min = 0.03;
66 
67 static float *aw;
68 static float *bbw;
69 static float *foq;
70 
71 static float *aph1;
72 static float *adg1;
73 static float *bbp1;
74 
75 static float *anap1;
76 static float *acdom1;
77 static float *bbph1;
78 static float *bbnap1;
79 
80 
81 // these need to be passed from run_giop if giop is to
82 // support simultaneous products (e.g., a gsm wrapper)
83 
84 static int16 *iter;
85 static int16 *iopf;
86 static float *a;
87 static float *bb;
88 static float *aph;
89 static float *adg;
90 static float *bbp;
91 
92 static float *a_unc;
93 static float *bb_unc;
94 
95 static float *anap;
96 static float *acdom;
97 static float *bbph;
98 static float *bbnap;
99 
100 static float *chl;
101 static float *aph_s;
102 static float *adg_s;
103 static float *uadg_s;
104 static float *bbp_s;
105 static float *ubbp_s;
106 static float *rrsdiff;
107 static float *mRrs;
108 static float **fit_par;
109 static float *chisqr;
110 static int16 max_npar;
111 
112 static int *siop_num;
113 static int allocateRrsRaman = 0;
114 
115 void freeArray(void **a, int32_t m) {
116  int i;
117  for (i = 0; i < m; ++i) {
118  free(a[i]);
119  }
120  free(a);
121 }
122 
123 void freeDArray(double **a, int32_t m) {
124  int i;
125  for (i = 0; i < m; ++i) {
126  free(a[i]);
127  }
128  free(a);
129 }
130 
131 /* ------------------------------------------------------------------- */
132 /* giop_int_tab_file - reads a table of eigenvectors and interpolates */
133 /* to input wavelengths and returns a 2-D array of vectors. */
134 /* Function reads in asscoicated uncertainty vectors from */
135 /* companion file (if avaiable) */
136 
137 /* ------------------------------------------------------------------- */
138 void giop_int_tab_file(char *file, char *ufile, int nx, float *x, float **y, float **uy) {
139 
140  float *table;
141  float *xtab;
142  float *ytab;
143  float *utable;
144  float *uxtab;
145  float *uytab;
146  int ncol, nrow, uncol, unrow;
147  int i, ivec;
148 
149  printf("\nLoading basis vector from %s\n", file);
150 
151 
152  nrow = table_row_count(file);
153  if (nrow <= 0) {
154  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, file);
155  exit(1);
156  }
157 
158  ncol = table_column_count(file);
159  table = table_read_r4(file, ncol, nrow);
160 
161  // Interpolate to input wavelengths (x)
162 
163  xtab = &table[nrow * 0];
164 
165  for (ivec = 0; ivec < ncol - 1; ivec++) {
166  ytab = &table[nrow * (ivec + 1)];
167  for (i = 0; i < nx; i++) {
168  y[ivec][i] = linterp(xtab, ytab, nrow, x[i]);
169  uy[ivec][i] = 0.0;
170  }
171  }
172 
173  table_free_r4(table);
174 
175  //Check if uncertainty file provided
176  if (!ufile) {
177  printf("\nNo uncertainty basis vector file ascociated with %s\n", file);
178  } else {
179 
180  unrow = table_row_count(ufile);
181  if (unrow <= 0) {
182  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, ufile);
183  exit(1);
184  }
185 
186  uncol = table_column_count(ufile);
187  utable = table_read_r4(ufile, uncol, unrow);
188 
189  // Interpolate to input wavelengths (x)
190 
191  uxtab = &utable[unrow * 0];
192 
193  for (ivec = 0; ivec < uncol - 1; ivec++) {
194  uytab = &utable[unrow * (ivec + 1)];
195  for (i = 0; i < nx; i++) {
196  uy[ivec][i] = linterp(uxtab, uytab, unrow, x[i]);
197  }
198  }
199 
200  table_free_r4(utable);
201 
202  }
203 
204 }
205 
206 
207 /* ------------------------------------------------------------------- */
208 /* giop_ctl_start - set starting values for optimization */
209 
210 /* ------------------------------------------------------------------- */
211 void giop_ctl_start(giopstr *g, float chl) {
212 
213  int ipar, ivec;
214 
215  // set starting params based on model elements
216 
217  ipar = 0;
218 
219  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
220  g->par[ipar] = chl / g->aph_nvec;
221  g->len[ipar] = 0.5;
222  ipar++;
223  }
224 
225  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
226  g->par[ipar] = 0.01;
227  g->len[ipar] = 0.1;
228  ipar++;
229  }
230 
231  switch (g->bbp_opt) {
232  case BBPLASFIX:
233  case BBPQAAFIX:
234  break;
235  case BBPLAS:
236  g->par[ipar] = 1.0;
237  g->len[ipar] = 0.01;
238  ipar++;
239  break;
240  default:
241  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
242  g->par[ipar] = 0.001;
243  g->len[ipar] = 0.01;
244  ipar++;
245  }
246  break;
247  }
248 
249  //note: for SVDSIOP ipar > npar
250  if (ipar != g->npar && g->fit_opt != SVDSIOP) {
251  printf("-E- %s Line %d: number of GIOP fit parameters (%d) does not equal number of initialized parameter (%d)\n",
252  __FILE__, __LINE__, g->npar, ipar);
253  exit(1);
254  }
255 }
256 
257 /* ------------------------------------------------------------------- */
258 /* giop_ctl_init - initialize model control structure */
259 /* */
260 /* Provides default values or user-supplied parameters for all static */
261 /* model components. These may later be over-riden by dynamic model */
262 /* parameters (e.g., band-ratio based bbp_s). */
263 /* */
264 
265 /* ------------------------------------------------------------------- */
266 void giop_ctl_init(giopstr *g, int nwave, float wave[],
267  float aw[], float bbw[]) {
268 
269  int iw, iwx;
270 
271  float def_aph_w = 443.0;
272  float def_aph_s = 0.5;
273  float def_adg_w = 443.0;
274  float def_adg_s = 0.02061;
275  float def_bbp_w = 443.0;
276  float def_bbp_s = 1.03373;
277  float def_grd[] = {0.0949, 0.0794};
278 
279  // determine indices of wavelengths to be used (default to all vis)
280 
281  if (input->giop_wave[0] > 0) {
282  for (iw = 0; iw < nwave; iw++) {
283  if (input->giop_wave[iw] <= 0) break;
284  iwx = windex(input->giop_wave[iw], wave, nwave);
285  g->bindx[iw] = iwx;
286  //windex(input->giop_wave[iw],wave,nwave);
287  }
288  g->nwave = iw;
289  } else {
290  g->nwave = MIN(windex(671., wave, nwave) + 1, nwave);
291  for (iw = 0; iw < g->nwave; iw++)
292  g->bindx[iw] = iw;
293  }
294 
295  // set wavelengths and pure-water values
296 
297  for (iw = 0; iw < g->nwave; iw++) {
298  g->wave[iw] = wave[g->bindx[iw]];
299  g->aw [iw] = aw [g->bindx[iw]];
300  g->bbw [iw] = bbw [g->bindx[iw]];
301  }
302 
303  //rrs uncertainties options
304 
305  if (input->giop_rrs_unc_opt >= 0) {
306  g->urrs_opt = input->giop_rrs_unc_opt;
307  } else {
308  g->urrs_opt = URRSNONE;
309  }
310 
311  //Input Rrs
312  /*
313  if (input->giop_rrs_unc[0] > 0) {
314  g->wt_opt = 1;
315  for (iw = 0; iw < nwave; iw++) {
316  if (input->giop_rrs_unc[iw] <= 0) break;
317  g->wts[iw] = 1.0 / pow(input->giop_rrs_unc[iw], 2);
318  }
319  if (iw != g->nwave) {
320  printf("-E- %s line %d: number of giop_rrs_unc (%d) must equal number of giop_wave (%d)",
321  __FILE__, __LINE__, iw, g->nwave);
322  exit(1);
323  }
324  } else {
325  g->wt_opt = 0;
326  for (iw = 0; iw < g->nwave; iw++)
327  g->wts[iw] = 1.0;
328  }
329  */
330 
331  // maximum iterations
332 
333  if (input->giop_maxiter > 0)
334  g->maxiter = input->giop_maxiter;
335  else
336  g->maxiter = 500;
337 
338  // fitting method
339 
340  if (input->giop_fit_opt > 0)
341  g->fit_opt = input->giop_fit_opt;
342  else
343  g->fit_opt = LEVMARQ;
344 
345  // Rrs to bb/(a+bb) method
346 
347  if (input->giop_rrs_opt > 0)
348  g->rrs_opt = input->giop_rrs_opt;
349  else
350  g->rrs_opt = RRSGRD;
351 
352 
353  // set coefficients of Gordon quadratic
354 
355  if (input->giop_grd[0] > -999.0) {
356  g->grd[0] = input->giop_grd[0];
357  g->grd[1] = input->giop_grd[1];
358  } else {
359  g->grd[0] = def_grd[0];
360  g->grd[1] = def_grd[1];
361  }
362 
363  // default basis vectors
364 
365  strcpy(g->aph_tab_file, input->giop_aph_file);
366  strcpy(g->adg_tab_file, input->giop_adg_file);
367  strcpy(g->bbp_tab_file, input->giop_bbp_file);
368  strcpy(g->acdom_tab_file, input->giop_acdom_file);
369  strcpy(g->anap_tab_file, input->giop_anap_file);
370  strcpy(g->bbph_tab_file, input->giop_bbph_file);
371  strcpy(g->bbnap_tab_file, input->giop_bbnap_file);
372 
373  //set defaults
374  g->aph_opt = APHTAB;
375  g->adg_opt = ADGTAB;
376  g->bbp_opt = BBPTAB;
377  g->acdom_opt = ACDOMNONE;
378  g->anap_opt = ANAPNONE;
379  g->bbnap_opt = BBPHNONE;
380  g->bbnap_opt = BBNAPNONE;
381 
382  // aphstar function
383 
384  if (input->giop_aph_opt > 0)
385  g->aph_opt = input->giop_aph_opt;
386 
387  if (input->giop_aph_w > 0.0)
388  g->aph_w = input->giop_aph_w;
389  else
390  g->aph_w = def_aph_w;
391 
392  if (input->giop_aph_s > -999.0)
393  g->aph_s = input->giop_aph_s;
394  else
395  g->aph_s = def_aph_s;
396 
397 
398  // adgstar function
399 
400  //if (input->giop_adg_opt > 0)
401  // g->adg_opt = input->giop_adg_opt;
402  if (input->giop_adg_opt != 1)
403  g->adg_opt = input->giop_adg_opt;
404  else
405  g->adg_opt = ADGS;
406 
407  if (input->giop_adg_w > 0.0)
408  g->adg_w = input->giop_adg_w;
409  else
410  g->adg_w = def_adg_w;
411 
412  if (input->giop_adg_s > -999.0) {
413  g->adg_s = input->giop_adg_s;
414  g->uadg_s = input->giop_uadg_s;
415  } else {
416  g->adg_s = def_adg_s;
417  g->uadg_s = 0.0;
418  }
419 
420  //acdom star function
421  if (input->giop_acdom_opt != 1)
422  g->acdom_opt = input->giop_acdom_opt;
423  else
424  g->acdom_opt = ACDOMNONE;
425 
426  //anap star function
427  if (input->giop_anap_opt != 1)
428  g->anap_opt = input->giop_anap_opt;
429  else
430  g->anap_opt = ANAPNONE;
431 
432  // bbpstar function
433 
434  //if (input->giop_bbp_opt > 0)
435  // g->bbp_opt = input->giop_bbp_opt;
436  if (input->giop_bbp_opt != 1)
437  g->bbp_opt = input->giop_bbp_opt;
438  else
439  g->bbp_opt = BBPS;
440 
441  if (input->giop_bbp_w > 0.0)
442  g->bbp_w = input->giop_bbp_w;
443  else
444  g->bbp_w = def_bbp_w;
445 
446  if (input->giop_bbp_s > -999.0) {
447  g->bbp_s = input->giop_bbp_s;
448  g->ubbp_s = input->giop_ubbp_s;
449  } else {
450  g->bbp_s = def_bbp_s;
451  g->ubbp_s = 0.0;
452  }
453 
454  //bbph star function
455  if (input->giop_bbph_opt != 1)
456  g->bbph_opt = input->giop_bbph_opt;
457  else
458  g->bbph_opt = BBPHNONE;
459 
460  //bbnap star function
461  if (input->giop_bbnap_opt != 1)
462  g->bbnap_opt = input->giop_bbnap_opt;
463  else
464  g->bbnap_opt = BBNAPNONE;
465 
466 
467  // set number of vectors
468 
469  switch (g->aph_opt) {
470  case APHTAB:
471  g->aph_nvec = table_column_count(g->aph_tab_file) - 1;
472  break;
473  default:
474  g->aph_nvec = 1;
475  break;
476  }
477 
478  switch (g->adg_opt) {
479  case ADGTAB:
480  g->adg_nvec = table_column_count(g->adg_tab_file) - 1;
481  break;
482  default:
483  g->adg_nvec = 1;
484  break;
485  }
486 
487  switch (g->acdom_opt) {
488  case ACDOMTAB:
489  g->acdom_nvec = table_column_count(g->acdom_tab_file) - 1;
490  break;
491  default:
492  g->acdom_nvec = 1;
493  break;
494  }
495 
496  switch (g->anap_opt) {
497  case ANAPTAB:
498  g->anap_nvec = table_column_count(g->anap_tab_file) - 1;
499  break;
500  default:
501  g->anap_nvec = 1;
502  break;
503  }
504 
505  switch (g->bbp_opt) {
506  case BBPTAB:
507  g->bbp_nvec = table_column_count(g->bbp_tab_file) - 1;
508  break;
509  case BBPLASFIX:
510  case BBPQAAFIX:
511  case BBPLH:
512  case BBPCHL:
513  g->bbp_nvec = 1;
514  break;
515  default:
516  g->bbp_nvec = 1;
517  break;
518  }
519 
520  switch (g->bbph_opt) {
521  case BBPHTAB:
522  g->bbph_nvec = table_column_count(g->bbph_tab_file) - 1;
523  break;
524  default:
525  g->bbph_nvec = 1;
526  break;
527  }
528 
529  switch (g->bbnap_opt) {
530  case BBNAPTAB:
531  g->bbnap_nvec = table_column_count(g->bbnap_tab_file) - 1;
532  break;
533  default:
534  g->bbnap_nvec = 1;
535  break;
536  }
537 
538  // total number of parameters to be optimized
539  if (g->fit_opt == SVDSIOP) {
540  //only fitting for chl, cdom and nap
541  g->npar = 3;
542  } else {
543  g->npar = g->aph_nvec + g->adg_nvec + g->bbp_nvec;
544  }
545 
546 
547  // allocate space for vectors (one element per sensor wavelength)
548 
549  g->aph_tab_nw = nwave;
550  g->aph_tab_w = (float *) calloc(nwave, sizeof (float));
551  g->aph_tab_s = (float **) allocate2d_float(g->aph_tab_nw, g->aph_nvec);
552  g->uaph_tab_s = (float **) allocate2d_float(g->aph_tab_nw, g->aph_nvec);
553 
554  g->adg_tab_nw = nwave;
555  g->adg_tab_w = (float *) calloc(nwave, sizeof (float));
556  g->adg_tab_s = (float **) allocate2d_float(g->adg_tab_nw, g->adg_nvec);
557  g->uadg_tab_s = (float **) allocate2d_float(g->adg_tab_nw, g->adg_nvec);
558 
559  g->acdom_tab_nw = nwave;
560  g->acdom_tab_w = (float *) calloc(nwave, sizeof (float));
561  g->acdom_tab_s = (float **) allocate2d_float(g->acdom_tab_nw, g->acdom_nvec);
562  g->uacdom_tab_s = (float **) allocate2d_float(g->acdom_tab_nw, g->acdom_nvec);
563 
564  g->anap_tab_nw = nwave;
565  g->anap_tab_w = (float *) calloc(nwave, sizeof (float));
566  g->anap_tab_s = (float **) allocate2d_float(g->anap_tab_nw, g->anap_nvec);
567  g->uanap_tab_s = (float **) allocate2d_float(g->anap_tab_nw, g->anap_nvec);
568 
569  g->bbp_tab_nw = nwave;
570  g->bbp_tab_w = (float *) calloc(nwave, sizeof (float));
571  g->bbp_tab_s = (float **) allocate2d_float(g->bbp_tab_nw, g->bbp_nvec);
572  g->ubbp_tab_s = (float **) allocate2d_float(g->bbp_tab_nw, g->bbp_nvec);
573 
574  g->bbph_tab_nw = nwave;
575  g->bbph_tab_w = (float *) calloc(nwave, sizeof (float));
576  g->bbph_tab_s = (float **) allocate2d_float(g->bbph_tab_nw, g->bbph_nvec);
577  g->ubbph_tab_s = (float **) allocate2d_float(g->bbph_tab_nw, g->bbph_nvec);
578 
579  g->bbnap_tab_nw = nwave;
580  g->bbnap_tab_w = (float *) calloc(nwave, sizeof (float));
581  g->bbnap_tab_s = (float **) allocate2d_float(g->bbnap_tab_nw, g->bbnap_nvec);
582  g->ubbnap_tab_s = (float **) allocate2d_float(g->bbnap_tab_nw, g->bbnap_nvec);
583 
584  // set vector wavelengths (same as ALL sensor for now)
585 
586  for (iw = 0; iw < nwave; iw++) {
587  g->aph_tab_w[iw] = wave[iw];
588  g->adg_tab_w[iw] = wave[iw];
589  g->acdom_tab_w[iw] = wave[iw];
590  g->anap_tab_w[iw] = wave[iw];
591  g->bbp_tab_w[iw] = wave[iw];
592  g->bbph_tab_w[iw] = wave[iw];
593  g->bbnap_tab_w[iw] = wave[iw];
594  }
595 
596  // load tabulated vectors if requested
597 
598  switch (g->aph_opt) {
599  case APHTAB:
600  giop_int_tab_file(g->aph_tab_file, g->uaph_tab_file, nwave, wave, g->aph_tab_s, g->uaph_tab_s);
601  break;
602  }
603 
604  switch (g->adg_opt) {
605  case ADGTAB:
606  giop_int_tab_file(g->adg_tab_file, g->uadg_tab_file,nwave, wave, g->adg_tab_s, g->uadg_tab_s);
607  break;
608  }
609 
610  switch (g->bbp_opt) {
611  case BBPTAB:
612  giop_int_tab_file(g->bbp_tab_file, g->ubbp_tab_file, nwave, wave, g->bbp_tab_s, g->ubbp_tab_s);
613  break;
614  }
615  switch (g->acdom_opt) {
616  case ACDOMTAB:
617  giop_int_tab_file(g->acdom_tab_file,g->uacdom_tab_file, nwave, wave, g->acdom_tab_s, g->uacdom_tab_s);
618  break;
619  }
620  switch (g->anap_opt) {
621  case ANAPTAB:
622  giop_int_tab_file(g->anap_tab_file, g->uanap_tab_file, nwave, wave, g->anap_tab_s, g->uanap_tab_s);
623  break;
624  }
625  switch (g->bbph_opt) {
626  case BBPHTAB:
627  giop_int_tab_file(g->bbph_tab_file, g->ubbph_tab_file, nwave, wave, g->bbph_tab_s, g->ubbph_tab_s);
628  break;
629  }
630  switch (g->bbnap_opt) {
631  case BBNAPTAB:
632  giop_int_tab_file(g->bbnap_tab_file, g->ubbnap_tab_file, nwave, wave, g->bbnap_tab_s, g->ubbnap_tab_s);
633  break;
634  }
635 
636 }
637 
638 
639 /* ------------------------------------------------------------------- */
640 /* */
641 
642 /* ------------------------------------------------------------------- */
643 int giop_ran(int recnum) {
644  if (recnum == LastRecNum)
645  return 1;
646  else
647  return 0;
648 }
649 
650 /* ------------------------------------------------------------------- */
651 /* */
652 
653 /* ------------------------------------------------------------------- */
654 void giop_model(giopstr *g, double par[], int nwave, float wave[], float aw[], float bbw[],
655  float foq[], float aph[], float adg[], float bbp[],
656  double rrs[], double **dfdpar, double **parstar) {
657  // double rrs[],double dfdpar[NBANDS][GIOPMAXPAR],double parstar[NBANDS][GIOPMAXPAR])
658  int iw, iwtab, ivec, ipar;
659  float bb, a;
660  float x;
661  float *aphstar;
662  float *adgstar;
663  float *bbpstar;
664  float *uaphstar;
665  float *uadgstar;
666  float *ubbpstar;
667  float dfdx;
668  float dxda;
669  float dxdb;
670 
671  // note: input wavelength set can vary between optimization and evaluation calls
672  if ((aphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
673  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
674  __FILE__, __LINE__);
675  exit(1);
676  }
677  if ((adgstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
678  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
679  __FILE__, __LINE__);
680  exit(1);
681  }
682  if ((bbpstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
683  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
684  __FILE__, __LINE__);
685  exit(1);
686  }
687  if ((uaphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
688  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
689  __FILE__, __LINE__);
690  exit(1);
691  }
692  if ((uadgstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
693  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
694  __FILE__, __LINE__);
695  exit(1);
696  }
697  if ((ubbpstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
698  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
699  __FILE__, __LINE__);
700  exit(1);
701  }
702 
703  for (iw = 0; iw < nwave; iw++) {
704 
705  // evaluate model for each IOP component at all wavelengths
706  // also compute and store associated uncertainty, based on fitted paraneter
707  // uncertainties as stored in the second half of the par array
708  // note that we're using a simple summation for uncertainties (for now)
709 
710  ipar = 0;
711 
712  switch (g->aph_opt) {
713  case APHGAUSS:
714  aphstar[ipar] = exp(-1. * pow(wave[iw] - g->aph_w, 2) / (2 * pow(g->aph_s, 2)));
715  uaphstar[ipar] = 0;
716  aph[iw] = par[ipar] * aphstar[ipar];
717  aph[iw + nwave] = sqrt(pow(par[ipar + g->npar]*aphstar[ipar] ,2)
718  + pow(par[ipar]*uaphstar[ipar],2));
719  ipar++;
720  break;
721  default:
722  aph[iw] = 0;
723  aph[iw + nwave] = 0;
724  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
725  iwtab = windex(wave[iw], g->aph_tab_w, g->aph_tab_nw);
726  aphstar[ipar] = g->aph_tab_s[ivec][iwtab];
727  uaphstar[ipar] = g->uaph_tab_s[ivec][iwtab];
728  aph[iw] += par[ipar] * aphstar[ipar];
729  aph[iw + nwave] += sqrt(pow(par[ipar + g->npar]*aphstar[ipar],2)
730  + pow(par[ipar]*uaphstar[ipar],2));
731  ipar++;
732  }
733  break;
734  }
735 
736  switch (g->adg_opt) {
737  case ADGS:
738  case ADGSQAA:
739  case ADGSOBPG:
740  adgstar[ipar] = exp(-g->adg_s * (wave[iw] - g->adg_w));
741  uadgstar[ipar] = sqrt(pow(g->uadg_s*(g->adg_w - wave[iw])*exp(-g->adg_s * (wave[iw] - g->adg_w)) ,2));
742  adg[iw] = par[ipar] * adgstar[ipar];
743  adg[iw + nwave] = sqrt(pow(par[ipar + g->npar] * adgstar[ipar],2) + pow(par[ipar]*uadgstar[ipar],2));
744  ipar++;
745  break;
746  default:
747  adg[iw] = 0;
748  adg[iw + nwave] = 0;
749  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
750  iwtab = windex(wave[iw], g->adg_tab_w, g->adg_tab_nw);
751  adgstar[ipar] = g->adg_tab_s[ivec][iwtab];
752  uadgstar[ipar] = g->uadg_tab_s[ivec][iwtab];
753  adg[iw] += par[ipar] * adgstar[ipar];
754  adg[iw + nwave] += sqrt(pow(par[ipar + g->npar] * adgstar[ipar],2)
755  + pow(par[ipar]*uadgstar[ipar],2));
756  ipar++;
757  }
758  break;
759  }
760 
761  switch (g->bbp_opt) {
762  case BBPLASFIX:
763  case BBPQAAFIX:
764  iwtab = windex(wave[iw], g->bbp_tab_w, g->bbp_tab_nw);
765  bbpstar[ipar] = g->bbp_tab_s[0][iwtab];
766  bbp[iw] = bbpstar[ipar];
767  bbp[iw + nwave] = 0.0;
768  break;
769  case BBPLH:
770  case BBPCHL:
771  iwtab = windex(wave[iw], g->bbp_tab_w, g->bbp_tab_nw);
772  bbpstar[ipar] = g->bbp_tab_s[0][iwtab];
773  ubbpstar[ipar] = g->ubbp_tab_s[0][iwtab];
774  bbp[iw] = bbpstar[ipar];
775  bbp[iw + nwave] = ubbpstar[ipar];
776  break;
777  case BBPS:
778  case BBPSLAS:
779  case BBPSHAL:
780  case BBPSQAA:
781  case BBPSCIOTTI:
782  case BBPSMM01:
783  bbpstar[ipar] = pow((g->bbp_w / wave[iw]), g->bbp_s);
784  ubbpstar[ipar] = sqrt(pow(g->ubbp_s*pow(g->bbp_w / wave[iw], g->bbp_s)*log(g->bbp_w / wave[iw]),2));
785  bbp[iw] = par[ipar] * bbpstar[ipar];
786  bbp[iw + nwave] = sqrt(pow(par[ipar + g->npar] * bbpstar[ipar],2) +
787  pow(par[ipar]*ubbpstar[ipar],2));
788  ipar++;
789  break;
790  default:
791  bbp[iw] = 0;
792  bbp[iw + nwave] = 0;
793  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
794  iwtab = windex(wave[iw], g->bbp_tab_w, g->bbp_tab_nw);
795  bbpstar[ipar] = g->bbp_tab_s[ivec][iwtab];
796  ubbpstar[ipar] = g->ubbp_tab_s[ivec][iwtab];
797  bbp[iw] += par[ipar] * bbpstar[ipar];
798  bbp[iw + nwave] += sqrt(pow(par[ipar + g->npar] * bbpstar[ipar],2) +
799  pow(par[ipar]*ubbpstar[ipar],2));
800  ipar++;
801  }
802  break;
803  }
804 
805  a = aw [iw] + aph[iw] + adg[iw];
806  bb = bbw[iw] + bbp[iw];
807  x = bb / (a + bb);
808 
809  switch (g->rrs_opt) {
810  case RRSGRD:
811  rrs[iw] = g->grd[0] * x + g->grd[1] * pow(x, 2);
812  dfdx = g->grd[0] + 2 * g->grd[1] * x;
813  break;
814  case RRSFOQ:
815  rrs[iw] = foq[iw] * x;
816  dfdx = foq[iw];
817  break;
818  }
819 
820  if (dfdpar != NULL) {
821 
822  ipar = 0;
823 
824  dxda = -x * x / bb;
825  dxdb = x / bb + dxda;
826 
827  switch (g->aph_opt) {
828  default:
829  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
830  dfdpar[iw][ipar] = dfdx * dxda * aphstar[ipar];
831  ipar++;
832  }
833  break;
834  }
835 
836  switch (g->adg_opt) {
837  default:
838  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
839  dfdpar[iw][ipar] = dfdx * dxda * adgstar[ipar];
840  ipar++;
841  }
842  break;
843  }
844 
845  switch (g->bbp_opt) {
846  case BBPLASFIX:
847  case BBPQAAFIX:
848  case BBPLH:
849  case BBPCHL:
850  break;
851  default:
852  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
853  dfdpar[iw][ipar] = dfdx * dxdb * bbpstar[ipar];
854  ipar++;
855  }
856  break;
857  }
858  }
859 
860  if (parstar != NULL) {
861 
862  ipar = 0;
863 
864  switch (g->aph_opt) {
865  default:
866  for (ivec = 0; ivec < g->aph_nvec; ivec++) {
867  parstar[iw][ipar] = aphstar[ipar];
868  ipar++;
869  }
870  break;
871  }
872 
873  switch (g->adg_opt) {
874  default:
875  for (ivec = 0; ivec < g->adg_nvec; ivec++) {
876  parstar[iw][ipar] = adgstar[ipar];
877  ipar++;
878  }
879  break;
880  }
881 
882  switch (g->bbp_opt) {
883  case BBPLASFIX:
884  case BBPQAAFIX:
885  case BBPLH:
886  case BBPCHL:
887  break;
888  default:
889  for (ivec = 0; ivec < g->bbp_nvec; ivec++) {
890  parstar[iw][ipar] = bbpstar[ipar];
891  ipar++;
892  }
893  break;
894  }
895  }
896 
897  }
898 
899  free(aphstar);
900  free(adgstar);
901  free(bbpstar);
902  free(uaphstar);
903  free(uadgstar);
904  free(ubbpstar);
905  return;
906 }
907 /* ------------------------------------------------------------------- */
908 /* */
909 
910 /* ------------------------------------------------------------------- */
911 void giop_model_iterate(giopstr *g, double par[], int nwave, float wave[], float aw[], float bbw[],
912  float foq[], float aph[], float adg[], float bbp[], float acdom[], float anap[], float bbph[], float bbnap[],
913  double rrs[], double **dfdpar, double **parstar) {
914  // double rrs[],double dfdpar[NBANDS][GIOPMAXPAR],double parstar[NBANDS][GIOPMAXPAR])
915  int iw, iwtab, idx443;
916 
917  int16 isiop = g->siopIdx;
918 
919  float bb, a;
920  float x;
921  float acdom443;
922 
923  float *aphstar;
924  float *acdomstar;
925  float *anapstar;
926  float *bbphstar;
927  float *bbnapstar;
928 
929  float dfdx;
930  float dxda;
931  float dxdb;
932 
933 
934  /* note: input wavelength set can vary between optimization and evaluation calls*/
935  if ((aphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
936  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
937  __FILE__, __LINE__);
938  exit(1);
939  }
940  if ((acdomstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
941  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
942  __FILE__, __LINE__);
943  exit(1);
944  }
945  if ((anapstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
946  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
947  __FILE__, __LINE__);
948  exit(1);
949  }
950  if ((bbphstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
951  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
952  __FILE__, __LINE__);
953  exit(1);
954  }
955  if ((bbnapstar = (float *) calloc(nwave, sizeof (float))) == NULL) {
956  printf("-E- %s line %d : error allocating memory for GIOP:giop_model.\n",
957  __FILE__, __LINE__);
958  exit(1);
959  }
960 
961 
962  /*Sanity check: do the input tables have the same number of columns?*/
963  switch (g->fit_opt) {
964  case SVDSIOP:
965  if (g->aph_nvec != g->acdom_nvec || g->aph_nvec != g->anap_nvec
966  || g->aph_nvec != g->bbph_nvec || g->aph_nvec != g->bbnap_nvec
967  || g->acdom_nvec != g->anap_nvec || g->acdom_nvec != g->bbph_nvec
968  || g->acdom_nvec != g->bbnap_nvec || g->anap_nvec != g->bbph_nvec
969  || g->anap_nvec != g->bbnap_nvec || g->bbph_nvec != g->bbnap_nvec) {
970 
971  printf("-E- %s line %d: SIOP tables must have same number of columns.\n",
972  __FILE__, __LINE__);
973  exit(1);
974  }
975  break;
976  }
977 
978  idx443 = windex(443., g->acdom_tab_w, g->acdom_tab_nw);
979 
980  for (iw = 0; iw < nwave; iw++) {
981 
982  // evaluate model for each IOP component at all wavelengths
983  // also compute and store associated uncertainty, based on fitted paraneter
984  // uncertainties as stored in the second half of the par array
985  // note that we're using a simple summation for uncertainties (for now)
986 
987  iwtab = windex(wave[iw], g->aph_tab_w, g->aph_tab_nw);
988  aphstar[iw] = g->aph_tab_s[isiop][iwtab];
989  aph[iw] = par[0] * aphstar[iw];
990  aph[iw + nwave] = par[g->npar] * aphstar[iw];
991 
992  /*Esnure that acdom443 is equal to 1.0*/
993  iwtab = windex(wave[iw], g->acdom_tab_w, g->acdom_tab_nw);
994  acdomstar[iw] = g->acdom_tab_s[isiop][iwtab];
995  acdom443 = g->acdom_tab_s[isiop][idx443];
996  acdom[iw] = par[1]*(acdomstar[iw] / acdom443);
997  acdom[iw + nwave] = par[1 + g->npar]*(aphstar[iw] / acdom443);
998 
999  iwtab = windex(wave[iw], g->anap_tab_w, g->anap_tab_nw);
1000  anapstar[iw] = g->anap_tab_s[isiop][iwtab];
1001  anap[iw] = par[2] * anapstar[iw];
1002  anap[iw + nwave] = par[2 + g->npar] * aphstar[iw];
1003 
1004  adg[iw] = par[1] * acdomstar[iw] + par[2] * anapstar[iw];
1005  adg[iw + nwave] = par[1 + g->npar] * acdomstar[iw] + par[2 + g->npar] * anapstar[iw];
1006 
1007  iwtab = windex(wave[iw], g->bbph_tab_w, g->bbph_tab_nw);
1008  bbphstar[iw] = g->bbph_tab_s[isiop][iwtab];
1009  bbph[iw] = par[0] * bbphstar[iw];
1010  bbph[iw + nwave] = par[g->npar] * bbphstar[iw];
1011 
1012  iwtab = windex(wave[iw], g->bbnap_tab_w, g->bbnap_tab_nw);
1013  bbnapstar[iw] = g->bbnap_tab_s[isiop][iwtab];
1014  bbnap[iw] = par[2] * bbnapstar[iw];
1015  bbnap[iw + nwave] = par[2 + g->npar] * bbnapstar[iw];
1016 
1017  bbp[iw] = par[0] * bbphstar[iw] + par[2] * bbnapstar[iw];
1018  bbp[iw + nwave] = par[g->npar] * bbphstar[iw] + par[2 + g->npar] * bbnapstar[iw];
1019 
1020  a = aw [iw] + aph[iw] + adg[iw];
1021  bb = bbw[iw] + bbp[iw];
1022  x = bb / (a + bb);
1023 
1024  switch (g->rrs_opt) {
1025  case RRSGRD:
1026  rrs[iw] = g->grd[0] * x + g->grd[1] * pow(x, 2);
1027  dfdx = g->grd[0] + 2 * g->grd[1] * x;
1028  break;
1029  case RRSFOQ:
1030  rrs[iw] = foq[iw] * x;
1031  dfdx = foq[iw];
1032  break;
1033  }
1034 
1035  if (dfdpar != NULL) {
1036 
1037  dxda = -x * x / bb;
1038  dxdb = x / bb + dxda;
1039 
1040  dfdpar[iw][0] = dfdx * dxda * aphstar[iw];
1041  dfdpar[iw][1] = dfdx * dxda * acdomstar[iw];
1042  dfdpar[iw][2] = dfdx * dxdb * anapstar[iw];
1043  dfdpar[iw][3] = dfdx * dxdb * bbphstar[iw];
1044  dfdpar[iw][4] = dfdx * dxdb * bbnapstar[iw];
1045 
1046 
1047  }
1048 
1049  if (parstar != NULL) {
1050 
1051  parstar[iw][0] = aphstar[iw];
1052  parstar[iw][1] = acdomstar[iw];
1053  parstar[iw][2] = anapstar[iw];
1054  parstar[iw][3] = bbphstar[iw];
1055  parstar[iw][4] = bbnapstar[iw];
1056 
1057  }
1058 
1059  }
1060 
1061  free(aphstar);
1062  free(acdomstar);
1063  free(anapstar);
1064  free(bbphstar);
1065  free(bbnapstar);
1066 
1067  return;
1068 }
1069 
1070 /* ------------------------------------------------------------------- */
1071 /* */
1072 
1073 /* ------------------------------------------------------------------- */
1074 double giop_amb(FITSTRUCT *ambdata, double par[]) {
1075  int iw;
1076 
1077  giopstr *g = (giopstr *) (ambdata->meta);
1078 
1079  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, ambdata->yfit, NULL, NULL);
1080 
1081  ambdata->merit = 0.0;
1082  for (iw = 0; iw < g->nwave; iw++) {
1083  ambdata->merit += pow((ambdata->y[iw] - ambdata->yfit[iw]), 2) * ambdata->wgt[iw];
1084  }
1085 
1086  return (ambdata->merit);
1087 }
1088 
1089 
1090 /* ------------------------------------------------------------------- */
1091 /* */
1092 
1093 /* ------------------------------------------------------------------- */
1094 int fit_giop_amb(giopstr *g, double Rrs[], double wts[], double par[],
1095  double Rrs_fit[], int16 *itercnt) {
1096  int status = 0;
1097 
1098 
1099  static float tol = 1.e-6; /* fractional change in chisqr */
1100  static FITSTRUCT ambdata; /* amoeba interface structure */
1101  static double *init; /* initial simplex */
1102 
1103  static int firstCall = 1;
1104 
1105  int i, j;
1106  short isml;
1107 
1108  ambdata.niter = g->maxiter; /* max number of iterations */
1109  ambdata.nfunc = g->npar; /* number of model parameters */
1110  ambdata.npnts = g->nwave; /* number of wavelengths (Rrs) */
1111  ambdata.y = Rrs; /* Input Rrs values (subsurface) */
1112  ambdata.wgt = wts; /* Input weights on Rrs values */
1113  ambdata.yfit = Rrs_fit; /* Output model predicted Rrs values */
1114  ambdata.meta = g;
1115 
1116  if (firstCall == 1) {
1117  firstCall = 0;
1118  if ((init = (double *) calloc(g->nwave * (g->nwave + 1), sizeof (double))) == NULL) {
1119  printf("-E- %s line %d : error allocating memory for GIOP:gio_model.\n",
1120  __FILE__, __LINE__);
1121  exit(1);
1122  }
1123 
1124  }
1125  /* initialize simplex with first guess model parameters */
1126  for (j = 0; j < g->npar + 1; j++)
1127  for (i = 0; i < g->npar; i++)
1128  init[j * g->npar + i] = g->par[i];
1129 
1130  for (i = 0; i < g->npar; i++) {
1131  init[g->npar + i * (g->npar + 1)] += g->len[i];
1132  par[i] = 0.0;
1133  }
1134 
1135  /* run optimization */
1136  isml = amoeba(init, &ambdata, giop_amb, tol);
1137 
1138  /* check convergence and record parameter results */
1139  if (ambdata.niter >= g->maxiter)
1140  status = 1;
1141 
1142  for (i = 0; i < g->npar; i++) {
1143  par[i] = init[g->npar * isml + i];
1144  }
1145 
1146  *itercnt = ambdata.niter;
1147 
1148  return (status);
1149 }
1150 
1151 
1152 
1153 /* ---------------------------------------------------------------------- */
1154 /* wrapper function for L-M fit to GIOP model */
1155 
1156 /* ---------------------------------------------------------------------- */
1157 int giop_lm_fdf(const gsl_vector *parv, void *data, gsl_vector *f, gsl_matrix *J) {
1158  double *y = ((fitstr *) data)->y;
1159  double *w = ((fitstr *) data)->w;
1160  float *aw = ((fitstr *) data)->g->aw;
1161  float *bbw = ((fitstr *) data)->g->bbw;
1162  float *foq = ((fitstr *) data)->g->foq;
1163  float *wave = ((fitstr *) data)->g->wave;
1164  size_t nwave = ((fitstr *) data)->g->nwave;
1165  size_t npar = ((fitstr *) data)->g->npar;
1166  giopstr *g = ((fitstr *) data)->g;
1167 
1168  double *par;
1169  double *yfit;
1170  double **dydpar;
1171 
1172  double sigma;
1173 
1174  int iw, ipar;
1175 
1176  if ((par = (double *) calloc(nwave, sizeof (double))) == NULL) {
1177  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1178  __FILE__, __LINE__);
1179  exit(1);
1180  }
1181 
1182  /* extract model params from vector */
1183  for (ipar = 0; ipar < npar; ipar++) {
1184  par[ipar] = gsl_vector_get(parv, ipar);
1185  }
1186  if ((yfit = (double *) calloc(nwave, sizeof (double))) == NULL) {
1187  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1188  __FILE__, __LINE__);
1189  exit(1);
1190  }
1191  if ((dydpar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1192  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1193  __FILE__, __LINE__);
1194  exit(1);
1195  }
1196  for (iw = 0; iw < nwave; iw++)
1197  if ((dydpar[iw] = (double *) calloc(nwave, sizeof (double))) == NULL) {
1198  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1199  __FILE__, __LINE__);
1200  exit(1);
1201  }
1202 
1203  /* run model */
1204  giop_model(g, par, nwave, wave, aw, bbw, foq, aph1, adg1, bbp1, yfit, dydpar, NULL);
1205 
1206  /* evaluate and store for lm solver */
1207  for (iw = 0; iw < nwave; iw++) {
1208 
1209  /* silly, but we need sigma and not sigma-squared */
1210  sigma = sqrt(1. / w[iw]);
1211 
1212  /* Function to be minimized */
1213  if (f != NULL)
1214  gsl_vector_set(f, iw, (yfit[iw] - y[iw]) / sigma);
1215 
1216  /* Jacobian matrix */
1217  if (J != NULL)
1218  for (ipar = 0; ipar < npar; ipar++)
1219  gsl_matrix_set(J, iw, ipar, dydpar[iw][ipar] / sigma);
1220  }
1221 
1222  free(yfit);
1223  freeDArray(dydpar, nwave);
1224  free(par);
1225  return GSL_SUCCESS;
1226 }
1227 
1228 int giop_lm_f(const gsl_vector *parv, void *data, gsl_vector *f) {
1229  return (giop_lm_fdf(parv, data, f, NULL));
1230 }
1231 
1232 int giop_lm_df(const gsl_vector *parv, void *data, gsl_matrix *J) {
1233  return (giop_lm_fdf(parv, data, NULL, J));
1234 }
1235 
1236 /* ---------------------------------------------------------------------- */
1237 /*calc_uadg_s() - calculate uncertainty in model ed adg_s slope */
1238 
1239 /* ---------------------------------------------------------------------- */
1240 
1241 float calc_uadg_s(giopstr *g, float Rrs1, float Rrs2, float uRrs1, float uRrs2, float covRrs1Rrs2) {
1242 
1243  float dsdr1,dsdr2;
1244  float uadg_s;
1245 
1246  switch (g->adg_opt) {
1247  case ADGSQAA:
1248  dsdr1 = -0.002 / (0.36*Rrs2 + 1.2*Rrs1 + Rrs1*Rrs1/Rrs2);
1249  dsdr2 = (0.002*Rrs1) / pow(Rrs1 + 0.6*Rrs2,2);
1250  uadg_s = pow(pow(dsdr1*uRrs1,2) + pow(dsdr2*uRrs2,2) + 2*dsdr1*dsdr2*covRrs1Rrs2, 0.5);
1251  break;
1252  case ADGSOBPG:
1253  dsdr1 = 0.0038/(log(10)*Rrs1);
1254  dsdr2 = -0.0038/(log(10)*Rrs2);
1255  uadg_s = pow(pow(dsdr1*uRrs1,2) + pow(dsdr2*uRrs2,2) +2*dsdr1*dsdr2*covRrs1Rrs2, 0.5);
1256  break;
1257  }
1258  return uadg_s;
1259 }
1260 
1261 /* ---------------------------------------------------------------------- */
1262 /*calc_sbbp_unc() - calculate uncertainty in the modelled bbp_s slope */
1263 
1264 /* ---------------------------------------------------------------------- */
1265 
1266 float calc_ubbp_s(giopstr *g, float Rrs1, float Rrs2, float uRrs1, float uRrs2, float covRrs1Rrs2) {
1267 
1268  float v,dvdr1,dvdr2,dsdr1,dsdr2,dsdv;
1269  float dsdchl;
1270  float ubbp_s;
1271 
1272  switch (g->bbp_opt) {
1273  case BBPSHAL:
1274  dsdr1 = 0.8/Rrs1;
1275  dsdr2 = -0.8*Rrs1/pow(Rrs2,2);
1276  ubbp_s = pow( pow(dsdr1*uRrs1,2) + pow(dsdr2*uRrs2,2) + 2*dsdr1*dsdr2*covRrs1Rrs2, 0.5);
1277  break;
1278  case BBPSQAA:
1279  v = -0.9*(Rrs1/Rrs2);
1280  dsdv = -2.4*exp(v);
1281  dvdr1 = -0.9/Rrs1;
1282  dvdr2 = 0.9*(Rrs1/pow(Rrs2,2));
1283  dsdr1 = dsdv*dvdr1;
1284  dsdr2 = dsdv*dvdr2;
1285  ubbp_s = pow(pow(dsdr1*uRrs1,2) + pow(dsdr2*uRrs2,2) + 2*dsdr1*dsdr2*covRrs1Rrs2, 0.5);
1286  break;
1287  case BBPSCIOTTI:
1288  dsdchl = -0.768/(log(10)*g->chl);
1289  ubbp_s = pow(pow( dsdchl*g->uchl,2) ,0.5);
1290  break;
1291  case BBPSMM01:
1292  dsdchl = -0.5/(log(10)*g->chl);
1293  ubbp_s = pow(pow( dsdchl*g->uchl,2) ,0.5);
1294  break;
1295  }
1296 
1297  return ubbp_s;
1298 }
1299 
1300 /* ---------------------------------------------------------------------- */
1301 /*calc_pinv() - calculate the Moore-Penrose pseudo inverse of a matrix */
1302 
1303 /* ---------------------------------------------------------------------- */
1304 
1305 gsl_matrix* calc_pinv(gsl_matrix *A, const realtype rcond) {
1306 
1307  gsl_matrix *V, *Sigma_pinv, *U, *A_pinv;
1308  gsl_matrix *_tmp_mat = NULL;
1309  gsl_vector *_tmp_vec;
1310  gsl_vector *u;
1311  realtype x, cutoff;
1312  size_t i, j;
1313  unsigned int n = A->size1;
1314  unsigned int m = A->size2;
1315  bool was_swapped = false;
1316 
1317 
1318  if (m > n) {
1319  /* libgsl SVD can only handle the case m <= n - transpose matrix */
1320  was_swapped = true;
1321  _tmp_mat = gsl_matrix_alloc(m, n);
1322  gsl_matrix_transpose_memcpy(_tmp_mat, A);
1323  A = _tmp_mat;
1324  i = m;
1325  m = n;
1326  n = i;
1327  }
1328 
1329  /* do SVD */
1330  V = gsl_matrix_alloc(m, m);
1331  u = gsl_vector_alloc(m);
1332  _tmp_vec = gsl_vector_alloc(m);
1333  gsl_linalg_SV_decomp(A, V, u, _tmp_vec);
1334  gsl_vector_free(_tmp_vec);
1335 
1336  /* compute Σ⁻¹ */
1337  Sigma_pinv = gsl_matrix_alloc(m, n);
1338  gsl_matrix_set_zero(Sigma_pinv);
1339  cutoff = rcond * gsl_vector_max(u);
1340 
1341  for (i = 0; i < m; ++i) {
1342  if (gsl_vector_get(u, i) > cutoff) {
1343  x = 1. / gsl_vector_get(u, i);
1344  }
1345  else {
1346  x = 0.;
1347  }
1348  gsl_matrix_set(Sigma_pinv, i, i, x);
1349  }
1350 
1351  /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
1352  U = gsl_matrix_alloc(n, n);
1353  gsl_matrix_set_zero(U);
1354 
1355  for (i = 0; i < n; ++i) {
1356  for (j = 0; j < m; ++j) {
1357  gsl_matrix_set(U, i, j, gsl_matrix_get(A, i, j));
1358  }
1359  }
1360 
1361  if (_tmp_mat != NULL) {
1362  gsl_matrix_free(_tmp_mat);
1363  }
1364 
1365  /* two dot products to obtain pseudoinverse */
1366  _tmp_mat = gsl_matrix_alloc(m, n);
1367  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., _tmp_mat);
1368 
1369  if (was_swapped) {
1370  A_pinv = gsl_matrix_alloc(n, m);
1371  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., U, _tmp_mat, 0., A_pinv);
1372  }
1373  else {
1374  A_pinv = gsl_matrix_alloc(m, n);
1375  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., _tmp_mat, U, 0., A_pinv);
1376  }
1377 
1378  gsl_matrix_free(_tmp_mat);
1379  gsl_matrix_free(U);
1380  gsl_matrix_free(Sigma_pinv);
1381  gsl_vector_free(u);
1382  gsl_matrix_free(V);
1383 
1384  return A_pinv;
1385 }
1386 
1387 
1388 /* ---------------------------------------------------------------------- */
1389 /*par_unc_calc() - compute uncertainties in the GIOP model free parameters*/
1390 /* */
1391 /* Reference: */
1392 /* McKinna et al. (2019) AApproach for Propagating Radiometric Data */
1393 /* Uncertainties Through NASA Ocean Color Algorithms, Front. Earth Sci., */
1394 /* 7, doi: 10.3389/feart.2019.00176. */
1395 /* */
1396 /* Implementation: L. McKinna, GO2Q/NASA GSFC, October 2023 */
1397 
1398 /* ---------------------------------------------------------------------- */
1399 void calc_par_unc(giopstr *g, double par[], double uRrs[], double upar[]) {
1400 
1401  int iw, ix, iy;
1402 
1403  size_t nwave = g->nwave;
1404  int32_t m = g->nwave;
1405  int32_t n = g->npar;
1406 
1407  double *Rrs_f;
1408  double **dydpar;
1409 
1410  gsl_matrix *covRrs = gsl_matrix_alloc(m,m);
1411  gsl_matrix *Jac = gsl_matrix_alloc(m,n);
1412  gsl_matrix *invJac = gsl_matrix_alloc(n,m);
1413  gsl_matrix *tmpMatrix = gsl_matrix_alloc(m,n);
1414  gsl_matrix *covPar = gsl_matrix_alloc(n,n);
1415 
1416  if ((Rrs_f = (double *) calloc(nwave, sizeof (double))) == NULL) {
1417  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1418  __FILE__, __LINE__);
1419  exit(1);
1420  }
1421 
1422  if ((dydpar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1423  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1424  __FILE__, __LINE__);
1425  exit(1);
1426  }
1427 
1428  for (iw = 0; iw < nwave; iw++) {
1429  if ((dydpar[iw] = (double *) calloc(nwave, sizeof (double))) == NULL) {
1430  printf("-E- %s line %d : error allocating memory for GIOP:giop_lm_fdf.\n",
1431  __FILE__, __LINE__);
1432  exit(1);
1433  }
1434  }
1435 
1436 
1437  /* run model */
1438  switch (g->fit_opt) {
1439  case SVDSIOP:
1440  giop_model_iterate(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, acdom1, anap1, bbph1, bbnap1, Rrs_f, dydpar, NULL);
1441  break;
1442  default:
1443  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, Rrs_f, dydpar, NULL);
1444  break;
1445  }
1446 
1447  /* evaluate and store for lm solver */
1448  for (ix = 0; ix < g->nwave; ix++) {
1449  for (iy = 0; iy < g->npar; iy++) {
1450  gsl_matrix_set(Jac, ix, iy, dydpar[ix][iy]);
1451  }
1452  }
1453 
1454  //Set the error-covariance matrix with Rrs variance diagonal terms
1455  for (ix = 0; ix < g->nwave; ix++) {
1456  for (iy = 0; iy < g->nwave; iy++) {
1457  if (ix == iy) {
1458  gsl_matrix_set(covRrs, ix, iy, pow(uRrs[g->bindx[ix]],2));
1459  } else {
1460  gsl_matrix_set(covRrs, ix, iy, 0.0);
1461  }
1462  }
1463  }
1464 
1465  invJac = calc_pinv(Jac,1E-15);
1466 
1467  gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, covRrs, invJac, 0.0, tmpMatrix);
1468  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, invJac, tmpMatrix, 0.0, covPar);
1469 
1470  for (ix = 0; ix < g->npar; ix++){
1471  for (iy = 0; iy < g->npar; iy++){
1472  if (ix == iy) {
1473  upar[ix] = pow(gsl_matrix_get(covPar,ix,iy),0.5);
1474  }
1475  }
1476  }
1477 
1478  free(Rrs_f);
1479  freeDArray(dydpar, nwave);
1480  gsl_matrix_free(Jac);
1481  gsl_matrix_free(covPar);
1482  gsl_matrix_free(invJac);
1483  gsl_matrix_free(covRrs);
1484  gsl_matrix_free(tmpMatrix);
1485 
1486 }
1487 
1488 /* ---------------------------------------------------------------------- */
1489 /* fit_giop_lm() - runs Levenburg-Marquart optimization for one pixel */
1490 
1491 /* ---------------------------------------------------------------------- */
1492 int fit_giop_lm(giopstr *g, double Rrs[], double wts[], double par[], double *chi, int16 *itercnt) {
1493  int status = 0;
1494  int iter;
1495 
1496  static fitstr data;
1497 
1498  size_t npar = g->npar;
1499 
1500  static gsl_multifit_function_fdf func;
1501  const gsl_multifit_fdfsolver_type *t;
1502  gsl_multifit_fdfsolver *s;
1503  gsl_vector_view x;
1504 
1505  gsl_matrix *J = gsl_matrix_alloc(g->nwave, npar);
1506  gsl_matrix *cov = gsl_matrix_alloc(npar, npar);
1507 
1508  double sum, dof, c;
1509  int ipar;
1510 
1511  /* Set up data structure */
1512  data.y = Rrs;
1513  data.w = wts;
1514  data.g = g;
1515 
1516  /* Set up multifit function structure */
1517  func.f = &giop_lm_f;
1518  func.df = &giop_lm_df;
1519  func.fdf = &giop_lm_fdf;
1520  func.n = g->nwave;
1521  func.p = g->npar;
1522  func.params = &data;
1523 
1524  /* Allocate solver space */
1525  t = gsl_multifit_fdfsolver_lmsder;
1526  s = gsl_multifit_fdfsolver_alloc(t, g->nwave, g->npar);
1527 
1528  /* Set start params */
1529  x = gsl_vector_view_array(par, npar);
1530  gsl_multifit_fdfsolver_set(s, &func, &x.vector);
1531 
1532  /* Fit model for this pixel */
1533  status = 1;
1534  for (iter = 0; iter < g->maxiter; iter++) {
1535  gsl_multifit_fdfsolver_iterate(s);
1536  if (gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4) == GSL_SUCCESS) {
1537  status = 0;
1538  break;
1539  }
1540  }
1541  *itercnt = iter;
1542 
1543  /* Compute covariance matrix from Jacobian */
1544  gsl_multifit_fdfsolver_jac(s, J);
1545  gsl_multifit_covar(J, 0.0, cov);
1546 
1547  /* Compute chi-square */
1548  sum = gsl_blas_dnrm2(s->f); // quadrature sum of minimization function
1549  dof = func.n - func.p; // degrees of freedom
1550  *chi = pow(sum, 2.0) / dof; // chi-square per dof
1551 
1552  if (g->wt_opt == 0)
1553  c = sum / sqrt(dof); // use variance in fit to estimate Rrs uncertainty
1554  else
1555  c = 1.0; // Rrs uncertainty included in sum
1556 
1557  /* Retrieve fit params & error estimates */
1558  for (ipar = 0; ipar < npar; ipar++) {
1559  par[ipar] = gsl_vector_get(s->x, ipar);
1560  par[ipar + npar] = sqrt(gsl_matrix_get(cov, ipar, ipar)) * c;
1561  //printf("%d %lf %lf\n",ipar,par[ipar],par[ipar+npar]);
1562  }
1563 
1564  gsl_multifit_fdfsolver_free(s);
1565  gsl_matrix_free(cov);
1566  gsl_matrix_free(J);
1567 
1568  return (status);
1569 }
1570 
1571 
1572 /* ---------------------------------------------------------------------- */
1573 /* fit_giop_svd() - constrained matrix solution */
1574 
1575 /* ---------------------------------------------------------------------- */
1576 int fit_giop_svd(giopstr *g, double rrs[], double wts[], double par[]) {
1577  int status = 0; /* init to success */
1578 
1579  double *rrs_fit;
1580  double **parstar;
1581 
1582  size_t nwave = g->nwave;
1583  size_t npar = g->npar;
1584 
1585  int iw, ipar;
1586 
1587  double a0, a1, a2;
1588  double u0, u1, u;
1589 
1590  gsl_matrix *A = gsl_matrix_alloc(nwave, npar);
1591  gsl_matrix *V = gsl_matrix_alloc(npar, npar);
1592  gsl_vector *S = gsl_vector_alloc(npar);
1593  gsl_vector *W = gsl_vector_alloc(npar);
1594  gsl_vector *x = gsl_vector_alloc(npar);
1595  gsl_vector *b = gsl_vector_alloc(nwave);
1596 
1597  if ((rrs_fit = (double *) calloc(nwave, sizeof (double))) == NULL) {
1598  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1599  __FILE__, __LINE__);
1600  exit(1);
1601  }
1602  if ((parstar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1603  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1604  __FILE__, __LINE__);
1605  exit(1);
1606  }
1607  for (iw = 0; iw < nwave; iw++)
1608  if ((parstar[iw] = (double *) calloc(nwave, sizeof (double))) == NULL) {
1609  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1610  __FILE__, __LINE__);
1611  exit(1);
1612  }
1613 
1614  /* run model to get parstar wavelength-dependence terms */
1615  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, rrs_fit, NULL, parstar);
1616 
1617  /* clear fit parameters */
1618  for (ipar = 0; ipar < g->npar; ipar++)
1619  par[ipar] = BAD_FLT;
1620 
1621  /* build A matrix and b vector for A x = b */
1622 
1623  for (iw = 0; iw < g->nwave; iw++) {
1624 
1625  /* get u = bb/(a+bb) */
1626  switch (g->rrs_opt) {
1627  case RRSGRD:
1628  a2 = g->grd[1];
1629  a1 = g->grd[0];
1630  a0 = -rrs[iw];
1631  if ((gsl_poly_solve_quadratic(a2, a1, a0, &u0, &u1) == 2) && u1 > 0.0)
1632  u = u1;
1633  else {
1634  status = 1;
1635  goto cleanup;
1636  }
1637  break;
1638  case RRSFOQ:
1639  u = rrs[iw] / g->foq[iw];
1640  break;
1641  }
1642 
1643  gsl_vector_set(b, iw, -(u * g->aw[iw] + (u - 1.0) * g->bbw[iw]));
1644 
1645  for (ipar = 0; ipar < g->npar; ipar++) {
1646  if (ipar < (g->aph_nvec + g->adg_nvec))
1647  gsl_matrix_set(A, iw, ipar, parstar[iw][ipar] * u); // a
1648  else
1649  gsl_matrix_set(A, iw, ipar, parstar[iw][ipar]*(u - 1.0)); // bb
1650  }
1651  }
1652 
1653  /* solve A x = b for x */
1654  status = gsl_linalg_SV_decomp(A, V, S, W);
1655  status = gsl_linalg_SV_solve(A, V, S, b, x);
1656 
1657  /* extract fitted parameters */
1658  if (status == 0) {
1659  for (ipar = 0; ipar < g->npar; ipar++) {
1660  par[ipar] = gsl_vector_get(x, ipar);
1661  }
1662  }
1663 
1664 cleanup:
1665 
1666  gsl_matrix_free(A);
1667  gsl_matrix_free(V);
1668  gsl_vector_free(S);
1669  gsl_vector_free(x);
1670  gsl_vector_free(b);
1671 
1672  free(rrs_fit);
1673  freeDArray(parstar, nwave);
1674  return (status);
1675 }
1676 
1677 /* ---------------------------------------------------------------------- */
1678 /* fit_giop_svd_siop() - adaptive linear matrix inversion solution using */
1679 /* specific inherent optical properties (SIOPS) and */
1680 /* SVD to solve system of equations */
1681 /* */
1682 /* Reference: */
1683 /* Brando et al. (2012) Adaptive semianalytical inversion of ocean color */
1684 /* radiometry in optically complex waters, Applied Optics, 51(15), */
1685 /* 2808-2833. */
1686 /* */
1687 /* Implementation: L. McKinna, SAIC/NASA GSFC, November 2015 */
1688 /* */
1689 
1690 /* ---------------------------------------------------------------------- */
1691 int fit_giop_svd_siop(giopstr *g, double rrs[], double wts[], double par[], double *chi) {
1692 
1693  int16 status = 0; /* init to success */
1694 
1695  double *rrs_fit;
1696  double **parstar;
1697  double *parArr;
1698  double *parFit;
1699  double *rmse;
1700  int *badSolution;
1701 
1702 
1703  size_t nwave = g->nwave;
1704  size_t npar = g->npar; //Number of parameters fixed at 3
1705  size_t nvec = g->aph_nvec;
1706 
1707  int iw, ipar, iv, smlIdx;
1708 
1709  double a0, a1, a2;
1710  double u0, u1, u;
1711 
1712  double diffSq, diffSqSum, smlRmse, sumRrs;
1713 
1714  gsl_matrix *A = gsl_matrix_alloc(nwave, npar);
1715  gsl_matrix *V = gsl_matrix_alloc(npar, npar);
1716  gsl_vector *S = gsl_vector_alloc(npar);
1717  gsl_vector *W = gsl_vector_alloc(npar);
1718  gsl_vector *x = gsl_vector_alloc(npar);
1719  gsl_vector *b = gsl_vector_alloc(nwave);
1720 
1721 
1722  if ((rrs_fit = (double *) calloc(nwave, sizeof (double))) == NULL) {
1723  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1724  __FILE__, __LINE__);
1725  exit(1);
1726  }
1727  if ((parstar = (double **) calloc(nwave, sizeof (double *))) == NULL) {
1728  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1729  __FILE__, __LINE__);
1730  exit(1);
1731  }
1732  for (iw = 0; iw < nwave; iw++) {
1733  if ((parstar[iw] = (double *) calloc(5, sizeof (double))) == NULL) {
1734  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1735  __FILE__, __LINE__);
1736  exit(1);
1737  }
1738  }
1739  if ((parArr = (double *) calloc(npar * nvec, sizeof (double *))) == NULL) {
1740  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1741  __FILE__, __LINE__);
1742  exit(1);
1743  }
1744  if ((rmse = (double *) calloc(nvec, sizeof (double *))) == NULL) {
1745  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1746  __FILE__, __LINE__);
1747  exit(1);
1748  }
1749  if ((badSolution = (int *) calloc(nvec, sizeof (int *))) == NULL) {
1750  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1751  __FILE__, __LINE__);
1752  exit(1);
1753  }
1754  if ((parFit = (double *) calloc(npar, sizeof (double *))) == NULL) {
1755  printf("-E- %s line %d : error allocating memory for GIOP:fit_giop_svd.\n",
1756  __FILE__, __LINE__);
1757  exit(1);
1758  }
1759 
1760  /* intialize fit parameters */
1761  for (ipar = 0; ipar < g->npar; ipar++)
1762  par[ipar] = BAD_FLT;
1763 
1764  for (ipar = 0; ipar < 3; ipar++)
1765  parFit[ipar] = BAD_FLT;
1766 
1767 
1768  /* clear fit parameters */
1769  for (iv = 0; iv < nvec; iv++)
1770  badSolution[iv] = BAD_INT;
1771 
1772 
1773  /*Iterate over siop combinations*/
1774  for (iv = 0; iv < nvec; iv++) {
1775 
1776  g->siopIdx = iv;
1777 
1778  /* run model to get parstar wavelength-dependence terms */
1779  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);
1780 
1781  /* build A matrix and b vector for A x = b */
1782  for (iw = 0; iw < g->nwave; iw++) {
1783 
1784  /* get u = bb/(a+bb) */
1785  switch (g->rrs_opt) {
1786  case RRSGRD:
1787  a2 = g->grd[1];
1788  a1 = g->grd[0];
1789  a0 = -rrs[iw];
1790  if ((gsl_poly_solve_quadratic(a2, a1, a0, &u0, &u1) == 2) && u1 > 0.0)
1791  u = u1;
1792  else {
1793  status = 1;
1794  goto cleanup;
1795  }
1796  break;
1797  case RRSFOQ:
1798  u = rrs[iw] / g->foq[iw];
1799  break;
1800  }
1801 
1802  gsl_vector_set(b, iw, -(g->aw[iw] * u + g->bbw[iw]*(u - 1.0)));
1803 
1804  // SIOP implementation Aij in eq 12 of Brando et al 2012
1805  // A is 3 components + water
1806  // i=0 (phy): ipar=0 for aphy* and ipar 3 for bbphy*
1807  // i=1 (CDOM): ipar=1 for aCDOM*
1808  // i=2 (NAP): ipar=2 for anap* and ipar 4 for bbnap*
1809  gsl_matrix_set(A, iw, 0, parstar[iw][0] * u + parstar[iw][3]*(u - 1.0)); // a+ bb
1810  gsl_matrix_set(A, iw, 1, parstar[iw][1] * u); // ONLY a
1811  gsl_matrix_set(A, iw, 2, parstar[iw][2] * u + parstar[iw][4]*(u - 1.0)); // a+ bb
1812 
1813  }
1814 
1815  /* solve A x = b for x */
1816  /*SVD decomposition*/
1817  status = gsl_linalg_SV_decomp(A, V, S, W);
1818  status = gsl_linalg_SV_solve(A, V, S, b, x);
1819 
1820 
1821  for (ipar = 0; ipar < npar; ipar++) {
1822 
1823  /*Test for non-physical (negative) concentrations or matrix inversion*/
1824  /*failure. For such circumstances set to BAD_FLT*/
1825  if (gsl_vector_get(x, ipar) < 0.0 || status != 0) {
1826  badSolution[iv] = 1;
1827  for (ipar = 0; ipar < npar; ipar++) {
1828  parArr[iv * npar + ipar] = BAD_FLT;
1829  }
1830  break;
1831 
1832  } else if (gsl_vector_get(x, ipar) >= 0.0) {
1833  badSolution[iv] = 0;
1834  parArr[iv * npar + ipar] = gsl_vector_get(x, ipar);
1835  }
1836  }
1837  }
1838 
1839  /*Find optimal set of parameters*/
1840  g->siopIdx = 0;
1841  diffSq = 0;
1842  diffSqSum = 0;
1843 
1844  /*Compute the root mean square error (rmse) metric for all SIOP combinations. */
1845  /*The smallest rmse valid values is returned as the optimal solution.*/
1846  for (iv = 0; iv < nvec; iv++) {
1847 
1848  g->siopIdx = iv;
1849 
1850  for (ipar = 0; ipar < npar; ipar++) {
1851  //if (badSolution[iv]) {
1852  // parFit[ipar] = BAD_FLT;
1853  //} else {
1854  parFit[ipar] = parArr[iv * npar + ipar];
1855  //}
1856  }
1857 
1858  /*Calculate fitted Rrs*/
1859  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);
1860 
1861  /*reset the sum of squares values*/
1862  diffSq = 0;
1863  diffSqSum = 0;
1864  sumRrs = 0;
1865 
1866  for (iw = 0; iw < g->nwave; iw++) {
1867  diffSq = pow((rrs_fit[iw] - rrs[iw]), 2);
1868  diffSqSum += diffSq;
1869  sumRrs += rrs[iw];
1870  }
1871 
1872 
1873  if (badSolution[iv]) {
1874  rmse[iv] = 999;
1875  } else {
1876  /*Save value to relative RMSE array*/
1877  rmse[iv] = pow((diffSqSum / (g->nwave - 1)), 0.5); // sumRrs/(g->nwave - 1) );
1878  }
1879  }
1880 
1881  /*Find the smallest relative RMSE value*/
1882  /*Initialize with zeroth values*/
1883  smlRmse = rmse[0];
1884  smlIdx = 0;
1885 
1886  for (iv = 0; iv < nvec; iv++) {
1887  if (rmse[iv] < smlRmse) {
1888  smlRmse = rmse[iv];
1889  smlIdx = iv;
1890  }
1891  }
1892 
1893 
1894  //Save index of best SIOP combination to giop record
1895  g->siopIdx = smlIdx;
1896 
1897  /*Get optimal concentrations values from array of all concentrations*/
1898  /*If there is still a negative solution, set it to zero*/
1899  if (rmse[smlIdx] != 999 && badSolution[smlIdx] != 1) {
1900 
1901  for (ipar = 0; ipar < npar; ipar++) {
1902  par[ipar] = parArr[smlIdx * npar + ipar];
1903  }
1904  *chi = smlRmse;
1905  status = 0;
1906 
1907 
1908  } else {
1909  //No solution status
1910  for (ipar = 0; ipar < npar; ipar++) {
1911  par[ipar] = BAD_FLT;
1912  }
1913  *chi = BAD_FLT;
1914  status = -99;
1915 
1916  }
1917 
1918 cleanup:
1919  gsl_matrix_free(A);
1920  gsl_matrix_free(V);
1921  gsl_vector_free(S);
1922  gsl_vector_free(x);
1923  gsl_vector_free(b);
1924 
1925 
1926  free(rrs_fit);
1927  free(parArr);
1928  free(rmse);
1929  free(badSolution);
1930  freeDArray(parstar, nwave);
1931  return (status);
1932 }
1933 
1934 /* ------------------------------------------------------------------- */
1935 /* giop_chl() - returns magnitude of aph* as chl */
1936 
1937 /* ------------------------------------------------------------------- */
1938 float32 giop_chl(giopstr *g, int16 iopf, double par[], float *uchl) {
1939  float32 chl = badval;
1940  int ipar;
1941 
1942  float uchl2_temp = badval;
1943  *uchl = badval;
1944 
1945  switch (g->aph_opt) {
1946  case APHGAUSS:
1947  break;
1948  default:
1949  if ((iopf & IOPF_FAILED) == 0) {
1950  chl = 0.0;
1951  uchl2_temp = 0.0;
1952  for (ipar = 0; ipar < g->aph_nvec; ipar++) {
1953  chl += par[ipar];
1954  uchl2_temp += pow(par[ipar + g->npar],2);
1955  }
1956  *uchl = pow(uchl2_temp,0.5);
1957  }
1958  break;
1959  }
1960 
1961 
1962  return (chl);
1963 }
1964 
1965 /* ---------------------------------------------------------------------- */
1966 /* Convert Rrs[0+] to Rrs[0-] */
1967 
1968 /* ---------------------------------------------------------------------- */
1969 float rrs_above_to_below(float Rrs) {
1970  return (Rrs / (0.52 + 1.7 * Rrs));
1971 
1972 }
1973 
1974 /* ---------------------------------------------------------------------- */
1975 /* Uncertainty estimate convert Rrs[0+] to Rrs[0-] */
1976 
1977 /* ---------------------------------------------------------------------- */
1978 float rrs_above_to_below_unc(float Rrs, float uRrs) {
1979  float dydr;
1980  dydr = 0.52 / pow(0.52 + 1.7*Rrs, 2.);
1981  return pow(pow(dydr*uRrs,2.),0.5);
1982 
1983 }
1984 
1985 /* ---------------------------------------------------------------------- */
1986 /* Convert Rrs[0-] to Rrs[0+] */
1987 
1988 /* ---------------------------------------------------------------------- */
1989 float rrs_below_to_above(float rrs_s) {
1990  return ( (0.52 * rrs_s) / (1.0 - 1.7 * rrs_s));
1991 
1992 }
1993 
1994 /* ---------------------------------------------------------------------- */
1995 /* Uncertainty estimate convert Rrs[0-] to Rrs[0+] */
1996 
1997 /* ---------------------------------------------------------------------- */
1998 float rrs_below_to_above_unc(float rrs_s, float urrs_s) {
1999  float dydr;
2000  dydr = 0.52/pow(1. - 1.7*rrs_s,2.);
2001  return pow(pow(dydr*urrs_s,2.), 0.5);
2002 
2003 }
2004 
2005 /* ---------------------------------------------------------------------- */
2006 /* get_bbp_lh computes spectral backscattering coefficient using the */
2007 /* reflectance line height model of McKinna et al (2021) */
2008 /* -----------------------------------------------------------------------*/
2009 int get_bbp_lh(l2str *l2rec, giopstr *g, int ipb2, float tab_wave[],
2010  float tab_bbp[], float tab_ubbp[], int tab_nwave) {
2011 
2012  int iw;
2013  int i1, i2, i3;
2014  float Rrs1, Rrs2, Rrs3;
2015  float uRrs1, uRrs2, uRrs3, covRrs1Rrs2, covRrs1Rrs3, covRrs2Rrs3;
2016  float LH, uLH;
2017  float *uRrs2_new;
2018  uRrs2_new = (float*) calloc(1, sizeof(float));
2019 
2020  float bbp555_lh = badval;
2021  float ubbp555_lh = 0;
2022 
2023  float *wave = l2rec->l1rec->l1file->fwave;
2024  int32 nwave = l2rec->l1rec->l1file->nbands;
2025 
2026  static int status = -1;
2027 
2028  //intialize arrays
2029  for (iw = 0; iw < tab_nwave; iw++) {
2030  g->bbp_tab_w[iw] = wave[iw];
2031  tab_ubbp[iw] = badval;
2032  }
2033 
2034  //LH backscattering coefficients
2035  float alh[2] = {-2.5770 , 281.27};
2036  //LH backscattering uncertainties and covariance
2037  float ualh[3] = {2.4819E-2, 20.777, 0.24852};
2038 
2039  i1 = windex(490.0, wave, nwave);
2040  i2 = windex(555.0, wave, nwave);
2041  i3 = windex(670.0, wave, nwave);
2042 
2043  Rrs1 = l2rec->Rrs[ipb2 + i1];
2044  Rrs2 = l2rec->Rrs[ipb2 + i2];
2045  Rrs3 = l2rec->Rrs[ipb2 + i3];
2046 
2047  uRrs1 = g->uRrs_a[i1];
2048  uRrs2 = g->uRrs_a[i2];
2049  uRrs3 = g->uRrs_a[i3];
2050 
2051  covRrs1Rrs2 = 0.0;
2052  covRrs1Rrs3 = 0.0;
2053  covRrs2Rrs3 = 0.0;
2054 
2055  Rrs2 = conv_rrs_to_555(Rrs2, l2rec->l1rec->l1file->fwave[i2], uRrs2, uRrs2_new);
2056  //conv_rrs_to_555(float Rrs, float wave, float uRrs_in, float *uRrs_out )
2057  //uRrs2 = uconv_rrs_to_555(Rrs2, uRrs2, l2rec->l1rec->l1file->fwave[i2]);
2058 
2059  LH = Rrs2 - (0.64*Rrs1 + 0.36*Rrs3 );
2060  uLH = sqrt( pow(0.64*uRrs1,2) + pow(*uRrs2_new,2) +
2061  pow(0.36*uRrs3,2) - 2*0.64*covRrs1Rrs2 -
2062  2*0.64*0.36*covRrs1Rrs3 - 0.36*covRrs2Rrs3);
2063 
2064  //Compute bbp(555) and uncertainty using an empirical model
2065  bbp555_lh = pow(10.,alh[0] + alh[1]*LH);
2066 
2067  ubbp555_lh = sqrt(pow(log(10)*ualh[0]*bbp555_lh,2) +
2068  pow(LH*log(10)*ualh[1]*bbp555_lh,2) +
2069  pow(uLH*alh[1]*log(10)*bbp555_lh,2) +
2070  2*ualh[2]*pow(log(10)*bbp555_lh,2));
2071 
2072  if (bbp555_lh < 0) {
2073  return(0);
2074  } else {
2075  status = 1;
2076  }
2077 
2078  //Compute bbp_s according to Lee et al (2002)
2079  // update power-law exponent based on QAA band-ratio relationship
2080  i1 = windex(443.0, wave, nwave);
2081  i2 = windex(550.0, wave, nwave);
2082  //QAA bbp_s relationship derived from NOMAD data.
2083  //No Raman scattering correction applied to Rrs.
2084  Rrs1 = l2rec->Rrs[ipb2 + i1];
2085  Rrs2 = l2rec->Rrs[ipb2 + i2];
2086  uRrs1 = g->uRrs_a[i1];
2087  uRrs2 = g->uRrs_a[i2];
2088  covRrs1Rrs2 = 0.0;
2089 
2090  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2091  Rrs1 = rrs_above_to_below(Rrs1);
2092  Rrs2 = rrs_above_to_below(Rrs2);
2093  uRrs1 = rrs_above_to_below_unc(Rrs1, uRrs1);
2094  uRrs2 = rrs_above_to_below_unc(Rrs2, uRrs2);
2095  g->bbp_s = MAX(MIN(2.0 * (1.0 - 1.2 * exp(-0.9 * Rrs1 / Rrs2)), bbp_s_max), bbp_s_min);
2096  if (g->bbp_s != bbp_s_min || g->bbp_s != bbp_s_max) {
2097  g->ubbp_s = calc_ubbp_s(g, Rrs1, Rrs2, uRrs1, *uRrs2_new, covRrs1Rrs2);
2098  } else {
2099  g->ubbp_s = 0.0;
2100  }
2101  } else if (Rrs2 > 0.0) {
2102  g->bbp_s = bbp_s_min;
2103  } else {
2104  g->bbp_s = BAD_FLT;
2105  status = 0;
2106  }
2107 
2108  if (status == 0) {
2109  for (iw = 0; iw < tab_nwave; iw++) {
2110  g->bbp_tab_w[iw] = wave[iw];
2111  tab_bbp[iw] = 0;
2112  tab_ubbp[iw] = 0;
2113  }
2114  } else {
2115  for (iw = 0; iw < tab_nwave; iw++) {
2116  g->bbp_tab_w[iw] = wave[iw];
2117  tab_bbp[iw] = bbp555_lh*pow(555./tab_wave[iw], g->bbp_s);
2118  tab_ubbp[iw] = sqrt(
2119  pow(bbp555_lh*g->ubbp_s*pow(555./tab_wave[iw], g->bbp_s)*log(555./tab_wave[iw]),2)
2120  + pow(ubbp555_lh*pow(555./tab_wave[iw], g->bbp_s),2));
2121  }
2122  }
2123 
2124  free(uRrs2_new);
2125  return status;
2126 }
2127 
2128 /* ---------------------------------------------------------------------- */
2129 /* get_bbp_chl computes spectral backscattering coefficient using the */
2130 /* chlorophyll-dependent model of Huot et al (2008) */
2131 /* ----------------------------------------------------------------------- */
2132 int get_bbp_chl(l2str *l2rec, giopstr *g, int ipb2, float tab_wave[],
2133  float tab_bbp[], float tab_ubbp[], int tab_nwave) {
2134 
2135  int iw;
2136  int i1, i2;
2137  float Rrs1, Rrs2;
2138  float uRrs1, uRrs2,covRrs1Rrs2;
2139 
2140  float bbp555_chl = badval;
2141  float ubbp555_chl = 0;
2142 
2143  float *wave = l2rec->l1rec->l1file->fwave;
2144  int32 nwave = l2rec->l1rec->l1file->nbands;
2145 
2146  static int status = -1;
2147 
2148  //intialize arrays
2149  for (iw = 0; iw < tab_nwave; iw++) {
2150  g->bbp_tab_w[iw] = wave[iw];
2151  tab_ubbp[iw] = badval;
2152  }
2153 
2154  //Huot backscattering coefficients
2155  float achl[2] = {0, 0};
2156  //Huot backscattering uncertainties and covariance
2157  float uachl[2] = {1E-4, 0.02};
2158 
2159  i1 = windex(550.0, wave, nwave);
2160 
2161  achl[0] = 2.267E-3 - 5.058E-6 * (l2rec->l1rec->l1file->fwave[i1] - 550.);
2162  achl[1] = 0.565 - 4.86E-4 * (l2rec->l1rec->l1file->fwave[i1] - 550.);
2163 
2164  //Compute bbp(555) and uncertainty estimate using an empirical model
2165  bbp555_chl = achl[0]*pow(g->chl,achl[1]);
2166  ubbp555_chl = sqrt( g->uchl*pow(achl[0]*achl[1]*pow(g->chl,(achl[1]-1)),2) +
2167  pow(uachl[0]*pow(g->chl,achl[1]),2) +
2168  pow(uachl[1]* log(g->chl)*achl[0]*pow(g->chl,achl[1]),2));
2169 
2170 
2171  if (bbp555_chl < 0) {
2172  return(0);
2173  } else {
2174  status = 1;
2175  }
2176 
2177  //Compute bbp_s according to Lee et al (2002)
2178  // update power-law exponent based on QAA band-ratio relationship
2179  i1 = windex(443.0, wave, nwave);
2180  i2 = windex(550.0, wave, nwave);
2181  //QAA bbp_s relationship derived from NOMAD data.
2182  //No Raman scattering correction applied to Rrs.
2183  Rrs1 = l2rec->Rrs[ipb2 + i1];
2184  Rrs2 = l2rec->Rrs[ipb2 + i2];
2185  uRrs1 = g->uRrs_a[i1];
2186  uRrs2 = g->uRrs_a[i2];
2187  covRrs1Rrs2 = 0.0;
2188 
2189  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2190  Rrs1 = rrs_above_to_below(Rrs1);
2191  Rrs2 = rrs_above_to_below(Rrs2);
2192  uRrs1 = rrs_above_to_below_unc(Rrs1, uRrs1);
2193  uRrs2 = rrs_above_to_below_unc(Rrs2, uRrs2);
2194  g->bbp_s = MAX(MIN(2.0 * (1.0 - 1.2 * exp(-0.9 * Rrs1 / Rrs2)), bbp_s_max), bbp_s_min);
2195  if (g->bbp_s != bbp_s_min || g->bbp_s != bbp_s_max) {
2196  g->ubbp_s = calc_ubbp_s(g, Rrs1, Rrs2, uRrs1, uRrs2, covRrs1Rrs2);
2197  } else {
2198  g->ubbp_s = 0.0;
2199  }
2200  } else if (Rrs2 > 0.0) {
2201  g->bbp_s = bbp_s_min;
2202  } else {
2203  g->bbp_s = BAD_FLT;
2204  status = 0;
2205  }
2206 
2207  if (status == 0) {
2208  for (iw = 0; iw < tab_nwave; iw++) {
2209  g->bbp_tab_w[iw] = wave[iw];
2210  tab_bbp[iw] = 0;
2211  tab_ubbp[iw] = 0;
2212  }
2213  } else {
2214  for (iw = 0; iw < tab_nwave; iw++) {
2215  g->bbp_tab_w[iw] = wave[iw];
2216  tab_bbp[iw] = bbp555_chl*pow(555./tab_wave[iw], g->bbp_s);
2217  tab_ubbp[iw] = sqrt(
2218  pow(bbp555_chl*g->ubbp_s*pow(555./tab_wave[iw], g->bbp_s)*log(555./tab_wave[iw]),2)
2219  + pow(ubbp555_chl*pow(555./tab_wave[iw], g->bbp_s),2));
2220  }
2221  }
2222 
2223  return status;
2224 }
2225 
2226 /* ---------------------------------------------------------------------- */
2227 /* run_giop - runs optimization using requested method */
2228 
2229 /* ---------------------------------------------------------------------- */
2230 void run_giop(l2str *l2rec) {
2231  static int firstCall = 1;
2232  static giopstr giopctl;
2233  static giopstr *g = &giopctl;
2234 
2235  double *par;
2236  double *upar;
2237  double *Rrs_a; /* above water, per fit band */
2238  double *uRrs_a; /* uncertainty above water, per fit band */
2239  double *Rrs_b; /* below water, per fit band */
2240  double *uRrs_b; /* uncertainty below water, per fit band */
2241  double *rrs_diff; /*observed-modeled rrs diff , per fit band */
2242  double *Rrs_f; /* modeled Rrs, per fit band */
2243  double *wts; /* weights, per fit band */
2244  double *uRrs;
2245 
2246  float Rrs1, Rrs2;
2247  /*uncertainties in Rrs1, Rrs2, and covariance*/
2248  float uRrs1, uRrs2, covRrs1Rrs2;
2249  int32 i1, i2;
2250 
2251  int16 itercnt;
2252  int16 bndcnt;
2253  int16 status;
2254 
2255  int32 ipar, ip, iw, ib, ipb, ipb2, ierr;
2256  double chi = BAD_FLT;
2257 
2258  l1str *l1rec = l2rec->l1rec;
2259  uncertainty_t *uncertainty=l1rec->uncertainty;
2260 
2261  float *wave = l1rec->l1file->fwave;
2262  int32 nwave = l1rec->l1file->nbands;
2263  static int32 npix;
2264 
2265 
2266  float aph_norm = BAD_FLT;
2267  float dudtau, dtaudchl,dudchl,dvdchl,daphdchl;
2268 
2269  if (firstCall) {
2270  npix = l1rec->npix;
2271 
2272  firstCall = 0;
2273 
2274  // initialize control structure (to get npar)
2275  if ((bbw = calloc(nwave, sizeof (float))) == NULL) {
2276  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2277  __FILE__, __LINE__);
2278  exit(1);
2279  }
2280  if ((aw = calloc(nwave, sizeof (float))) == NULL) {
2281  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2282  __FILE__, __LINE__);
2283  exit(1);
2284  }
2285  if ((foq = calloc(nwave, sizeof (float))) == NULL) {
2286  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2287  __FILE__, __LINE__);
2288  exit(1);
2289  }
2290  if ((aph1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2291  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2292  __FILE__, __LINE__);
2293  exit(1);
2294  }
2295  if ((acdom1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2296  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2297  __FILE__, __LINE__);
2298  exit(1);
2299  }
2300  if ((anap1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2301  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2302  __FILE__, __LINE__);
2303  exit(1);
2304  }
2305  if ((adg1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2306  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2307  __FILE__, __LINE__);
2308  exit(1);
2309  }
2310  if ((bbph1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2311  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2312  __FILE__, __LINE__);
2313  exit(1);
2314  }
2315  if ((bbnap1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2316  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2317  __FILE__, __LINE__);
2318  exit(1);
2319  }
2320  if ((bbp1 = calloc(2 * nwave, sizeof (float))) == NULL) {
2321  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2322  __FILE__, __LINE__);
2323  exit(1);
2324  }
2325  if ((g->wave = (float *) calloc(nwave, sizeof (float))) == NULL) {
2326  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2327  __FILE__, __LINE__);
2328  exit(1);
2329  }
2330  if ((g->bindx = (int *) calloc(nwave, sizeof (int))) == NULL) {
2331  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2332  __FILE__, __LINE__);
2333  exit(1);
2334  }
2335  if ((g->aw = (float *) calloc(nwave, sizeof (float))) == NULL) {
2336  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2337  __FILE__, __LINE__);
2338  exit(1);
2339  }
2340  if ((g->bbw = (float *) calloc(nwave, sizeof (float))) == NULL) {
2341  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2342  __FILE__, __LINE__);
2343  exit(1);
2344  }
2345  if ((g->wts = (float *) calloc(nwave, sizeof (float))) == NULL) {
2346  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2347  __FILE__, __LINE__);
2348  exit(1);
2349  }
2350  if ((g->foq = (float *) calloc(nwave, sizeof (float))) == NULL) {
2351  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2352  __FILE__, __LINE__);
2353  exit(1);
2354  }
2355  if ((g->par = (double *) calloc(nwave, sizeof (double))) == NULL) {
2356  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2357  __FILE__, __LINE__);
2358  exit(1);
2359  }
2360  if ((g->len = (double *) calloc(nwave, sizeof (double))) == NULL) {
2361  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2362  __FILE__, __LINE__);
2363  exit(1);
2364  }
2365  if ((g->Rrs_a = (float *) calloc(nwave, sizeof (float))) == NULL) {
2366  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2367  __FILE__, __LINE__);
2368  exit(1);
2369  }
2370  if ((g->uRrs_a = (float *) calloc(nwave, sizeof (float))) == NULL) {
2371  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
2372  __FILE__, __LINE__);
2373  exit(1);
2374  }
2375 
2376 
2377  giop_ctl_init(g, nwave, wave, aw, bbw);
2378 
2379  // allocate static storage for one scanline
2380 
2381  if ((iter = calloc(npix, sizeof (int16))) == NULL) {
2382  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2383  __FILE__, __LINE__);
2384  exit(1);
2385  }
2386  if ((iopf = calloc(npix, sizeof (int16))) == NULL) {
2387  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2388  __FILE__, __LINE__);
2389  exit(1);
2390  }
2391  if ((mRrs = calloc(npix * nwave, sizeof (float))) == NULL) {
2392  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2393  __FILE__, __LINE__);
2394  exit(1);
2395  }
2396  if ((a = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2397  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2398  __FILE__, __LINE__);
2399  exit(1);
2400  }
2401  if ((a_unc = calloc(npix * nwave, sizeof (float))) == NULL) {
2402  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2403  __FILE__, __LINE__);
2404  exit(1);
2405  }
2406  if ((aph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2407  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2408  __FILE__, __LINE__);
2409  exit(1);
2410  }
2411  if ((acdom = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2412  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2413  __FILE__, __LINE__);
2414  exit(1);
2415  }
2416  if ((anap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2417  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2418  __FILE__, __LINE__);
2419  exit(1);
2420  }
2421  if ((adg = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2422  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2423  __FILE__, __LINE__);
2424  exit(1);
2425  }
2426  if ((bbph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2427  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2428  __FILE__, __LINE__);
2429  exit(1);
2430  }
2431  if ((bbnap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2432  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2433  __FILE__, __LINE__);
2434  exit(1);
2435  }
2436  if ((bbp = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2437  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2438  __FILE__, __LINE__);
2439  exit(1);
2440  }
2441  if ((bb = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2442  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2443  __FILE__, __LINE__);
2444  exit(1);
2445  }
2446  if ((bb_unc = calloc(npix * nwave, sizeof (float))) == NULL) {
2447  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2448  __FILE__, __LINE__);
2449  exit(1);
2450  }
2451  if ((chl = calloc(npix * 2, sizeof (float))) == NULL) {
2452  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2453  __FILE__, __LINE__);
2454  exit(1);
2455  }
2456  if ((rrsdiff = calloc(npix, sizeof (float))) == NULL) {
2457  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2458  __FILE__, __LINE__);
2459  exit(1);
2460  }
2461  if ((aph_s = calloc(npix, sizeof (float))) == NULL) {
2462  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2463  __FILE__, __LINE__);
2464  exit(1);
2465  }
2466  if ((adg_s = calloc(npix, sizeof (float))) == NULL) {
2467  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2468  __FILE__, __LINE__);
2469  exit(1);
2470  }
2471  if ((uadg_s = calloc(npix, sizeof (float))) == NULL) {
2472  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2473  __FILE__, __LINE__);
2474  exit(1);
2475  }
2476  if ((bbp_s = calloc(npix, sizeof (float))) == NULL) {
2477  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2478  __FILE__, __LINE__);
2479  exit(1);
2480  }
2481  if ((ubbp_s = calloc(npix, sizeof (float))) == NULL) {
2482  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2483  __FILE__, __LINE__);
2484  exit(1);
2485  }
2486  if ((fit_par = allocate2d_float(npix, g->npar)) == NULL) {
2487  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2488  __FILE__, __LINE__);
2489  exit(1);
2490  }
2491  if ((siop_num = calloc(npix, sizeof (int))) == NULL) {
2492  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2493  __FILE__, __LINE__);
2494  exit(1);
2495  }
2496  max_npar = g->npar;
2497  if ((chisqr = calloc(npix, sizeof (float))) == NULL) {
2498  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2499  __FILE__, __LINE__);
2500  exit(1);
2501  }
2502 
2503  /*Note: Calculation of Raman scattering contribution to Rrs is at present*/
2504  /*only supported in l2gen. Where Raman Rrs is not calculated (i.e. l3gen smi)*/
2505  /*set Rrs_raman to zeros */
2506  if (l2rec->Rrs_raman == NULL) {
2507 
2508  printf("\n");
2509  printf("No Raman scattering correction applied to Rrs. \n");
2510  printf("\n");
2511 
2512  allocateRrsRaman = 1;
2513  l2rec->Rrs_raman = (float*) allocateMemory(npix *
2514  l1rec->l1file->nbands * sizeof (float), "Rrs_ram");
2515  }
2516 
2517  }
2518 
2519  if ((par = calloc(2 * nwave, sizeof (double))) == NULL) {
2520  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2521  __FILE__, __LINE__);
2522  exit(1);
2523  }
2524 
2525  if ((upar = calloc(g->npar, sizeof (double))) == NULL) {
2526  printf("-E- %s line %d : error allocating memory for GIOP:run_giop.\n",
2527  __FILE__, __LINE__);
2528  exit(1);
2529  }
2530 
2531 
2532 
2533  // reallocate if npix changes - l3gen-ism...
2534  if (l1rec->npix > npix) {
2535  npix = l1rec->npix;
2536  free(iter);
2537  free(iopf);
2538  free(mRrs);
2539  free(a);
2540  free(a_unc);
2541  free(aph);
2542  free(acdom);
2543  free(anap);
2544  free(adg);
2545  free(bbph);
2546  free(bbnap);
2547  free(bbp);
2548  free(bb);
2549  free(bb_unc);
2550  free(chl);
2551  free(rrsdiff);
2552  free(aph_s);
2553  free(adg_s);
2554  free(uadg_s);
2555  free(bbp_s);
2556  free(ubbp_s);
2557  free(fit_par);
2558  free(siop_num);
2559  free(chisqr);
2560  if (allocateRrsRaman) {
2561  free(l2rec->Rrs_raman);
2562  l2rec->Rrs_raman = (float*) allocateMemory(npix *
2563  l1rec->l1file->nbands * sizeof (float), "Rrs_ram");
2564  }
2565 
2566  if ((iter = calloc(npix, sizeof (int16))) == NULL) {
2567  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2568  __FILE__, __LINE__);
2569  exit(1);
2570  }
2571  if ((iopf = calloc(npix, sizeof (int16))) == NULL) {
2572  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2573  __FILE__, __LINE__);
2574  exit(1);
2575  }
2576  if ((mRrs = calloc(npix * nwave, sizeof (float))) == NULL) {
2577  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2578  __FILE__, __LINE__);
2579  exit(1);
2580  }
2581  if ((a = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2582  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2583  __FILE__, __LINE__);
2584  exit(1);
2585  }
2586  if ((a_unc = calloc(npix * nwave, sizeof (float))) == NULL) {
2587  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2588  __FILE__, __LINE__);
2589  exit(1);
2590  }
2591  if ((aph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2592  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2593  __FILE__, __LINE__);
2594  exit(1);
2595  }
2596  if ((acdom = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2597  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2598  __FILE__, __LINE__);
2599  exit(1);
2600  }
2601  if ((anap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2602  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2603  __FILE__, __LINE__);
2604  exit(1);
2605  }
2606  if ((adg = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2607  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2608  __FILE__, __LINE__);
2609  exit(1);
2610  }
2611  if ((bbph = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2612  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2613  __FILE__, __LINE__);
2614  exit(1);
2615  }
2616  if ((bbnap = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2617  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2618  __FILE__, __LINE__);
2619  exit(1);
2620  }
2621  if ((bbp = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2622  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2623  __FILE__, __LINE__);
2624  exit(1);
2625  }
2626  if ((bb = calloc(npix * nwave * 2, sizeof (float))) == NULL) {
2627  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2628  __FILE__, __LINE__);
2629  exit(1);
2630  }
2631  if ((bb_unc = calloc(npix * nwave, sizeof (float))) == NULL) {
2632  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2633  __FILE__, __LINE__);
2634  exit(1);
2635  }
2636  if ((chl = calloc(npix * 2, sizeof (float))) == NULL) {
2637  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2638  __FILE__, __LINE__);
2639  exit(1);
2640  }
2641  if ((rrsdiff = calloc(npix, sizeof (float))) == NULL) {
2642  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2643  __FILE__, __LINE__);
2644  exit(1);
2645  }
2646  if ((aph_s = calloc(npix, sizeof (float))) == NULL) {
2647  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2648  __FILE__, __LINE__);
2649  exit(1);
2650  }
2651  if ((adg_s = calloc(npix, sizeof (float))) == NULL) {
2652  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2653  __FILE__, __LINE__);
2654  exit(1);
2655  }
2656  if ((uadg_s = calloc(npix, sizeof (float))) == NULL) {
2657  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2658  __FILE__, __LINE__);
2659  exit(1);
2660  }
2661  if ((bbp_s = calloc(npix, sizeof (float))) == NULL) {
2662  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2663  __FILE__, __LINE__);
2664  exit(1);
2665  }
2666  if ((ubbp_s = calloc(npix, sizeof (float))) == NULL) {
2667  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2668  __FILE__, __LINE__);
2669  exit(1);
2670  }
2671  if ((fit_par = allocate2d_float(npix, g->npar)) == NULL) {
2672  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2673  __FILE__, __LINE__);
2674  exit(1);
2675  }
2676  if ((siop_num = calloc(npix, sizeof (int))) == NULL) {
2677  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2678  __FILE__, __LINE__);
2679  exit(1);
2680  }
2681  max_npar = g->npar;
2682  if ((chisqr = calloc(npix, sizeof (float))) == NULL) {
2683  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2684  __FILE__, __LINE__);
2685  exit(1);
2686  }
2687 
2688  }
2689 
2690  if ((Rrs_a = calloc(nwave, sizeof (double))) == NULL) {
2691  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2692  __FILE__, __LINE__);
2693  exit(1);
2694  }
2695  if ((uRrs_a = calloc(nwave, sizeof (double))) == NULL) {
2696  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2697  __FILE__, __LINE__);
2698  exit(1);
2699  }
2700  if ((Rrs_b = calloc(nwave, sizeof (double))) == NULL) {
2701  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2702  __FILE__, __LINE__);
2703  exit(1);
2704  }
2705  if ((uRrs_b = calloc(nwave, sizeof (double))) == NULL) {
2706  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2707  __FILE__, __LINE__);
2708  exit(1);
2709  }
2710  if ((rrs_diff = calloc(nwave, sizeof (double))) == NULL) {
2711  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2712  __FILE__, __LINE__);
2713  exit(1);
2714  }
2715  if ((Rrs_f = calloc(nwave, sizeof (double))) == NULL) {
2716  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2717  __FILE__, __LINE__);
2718  exit(1);
2719  }
2720  if ((wts = calloc(nwave, sizeof (double))) == NULL) {
2721  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2722  __FILE__, __LINE__);
2723  exit(1);
2724  }
2725  if ((uRrs = calloc(nwave, sizeof (double))) == NULL) {
2726  printf("-E- %s line %d : error allocating memory for GIOP.\n",
2727  __FILE__, __LINE__);
2728  exit(1);
2729  }
2730 
2731  // re-initialize control structure when giop_iterate=1 (under evaluation)
2732  // PJW 9 Jan 2013
2733  if (input->giop_iterate > 0) {
2734  g->adg_s = input->giop_adg_s;
2735  g->bbp_s = input->giop_bbp_s;
2736  }
2737 
2738 
2739  for (ip = 0; ip < l1rec->npix; ip++) {
2740 
2741  ipb = ip*nwave;
2742  ipb2 = ip * l1rec->l1file->nbands;
2743  ierr = l1rec->npix * nwave + ipb;
2744 
2745  // initialize output arrays and counters
2746 
2747  for (ib = 0; ib < nwave; ib++) {
2748  a [ipb + ib] = badval;
2749  a_unc [ipb + ib] = 0.0;
2750  bb [ipb + ib] = badval;
2751  bb_unc [ipb + ib] = 0.0;
2752  aph [ipb + ib] = badval;
2753  adg [ipb + ib] = badval;
2754  acdom [ipb + ib] = badval;
2755  anap[ipb + ib] = badval;
2756  bbp [ipb + ib] = badval;
2757  bbph[ipb + ib] = badval;
2758  bbnap [ipb + ib] = badval;
2759  mRrs[ipb + ib] = badval;
2760  a [ierr + ib] = 0.0;
2761  bb [ierr + ib] = 0.0;
2762  aph [ierr + ib] = 0.0;
2763  adg [ierr + ib] = 0.0;
2764  acdom [ierr + ib] = 0.0;
2765  anap [ierr + ib] = 0.0;
2766  bbp [ierr + ib] = 0.0;
2767  bbph [ierr + ib] = 0.0;
2768  bbnap [ierr + ib] = 0.0;
2769  }
2770  chl [ip] = badval;
2771  iter[ip] = 0;
2772  iopf[ip] = 0;
2773  status = 0;
2774  bndcnt = 0;
2775  itercnt = 0;
2776  rrsdiff[ip] = badval;
2777  chisqr[ip] = badval;
2778 
2779  aph_s[ip] = badval;
2780  adg_s[ip] = badval;
2781  uadg_s[ip] = 0.0;
2782  bbp_s[ip] = badval;
2783  ubbp_s[ip] = 0.0;
2784  siop_num[ip] = BAD_INT;
2785 
2786  for (ipar = 0; ipar < g->npar; ipar++)
2787  fit_par[ip][ipar] = badval;
2788 
2789 
2790  // flag and skip if pixel already masked
2791 
2792  if (l1rec->mask[ip]) {
2793  iopf[ip] |= IOPF_ISMASKED;
2794  iopf[ip] |= IOPF_FAILED;
2795  continue;
2796  }
2797 
2798  // set values that depend on sea state
2799 
2800  for (iw = 0; iw < nwave; iw++) {
2801  aw [iw] = l1rec->sw_a [ipb2 + iw];
2802  bbw[iw] = l1rec->sw_bb[ipb2 + iw];
2803  }
2804  for (iw = 0; iw < g->nwave; iw++) {
2805  g->aw [iw] = aw [g->bindx[iw]];
2806  g->bbw[iw] = bbw [g->bindx[iw]];
2807  }
2808 
2809  //define GIOP Rrs and accsociated uncertainties
2810  for (iw = 0; iw < nwave; iw++) {
2811  g->Rrs_a[iw] = l2rec->Rrs[ipb2 + iw];
2812 
2813  switch (g->urrs_opt) {
2814  case URRSNONE:
2815  g->uRrs_a[iw] = 0.0;
2816  break;
2817  case URRSCALC:
2818  if (uncertainty) {
2819  g->uRrs_a[iw] = l2rec->Rrs_unc[ipb2 + iw];
2820  } else {
2821  g->uRrs_a[iw] = 0.0;
2822  }
2823  break;
2824  case URRSREL:
2825  g->uRrs_a[iw] = g->Rrs_a[iw]*input->giop_rrs_unc[iw];
2826  break;
2827  case URRSABS:
2828  g->uRrs_a[iw] = input->giop_rrs_unc[iw];
2829  break;
2830  default:
2831  g->uRrs_a[iw] = 0.0;
2832  break;
2833  }
2834  }
2835 
2836 
2837  //define giop model fit weights for select bands
2838  for (iw = 0; iw < g->nwave; iw++) {
2839  ib = g->bindx[iw];
2840  switch (g->urrs_opt) {
2841  case URRSCALC:
2842  if (uncertainty) {
2843  g->wts[iw] = 1./pow(g->uRrs_a[g->bindx[iw]], 2);
2844  } else {
2845  g->wts[iw] = 1.0;
2846  }
2847  break;
2848  case URRSREL:
2849  case URRSABS:
2850  g->wt_opt = 1;
2851  g->wts[iw] = 1./pow(g->uRrs_a[g->bindx[iw]], 2);
2852  break;
2853  default:
2854  g->wt_opt = 0;
2855  g->wts[iw] = 1.0;
2856  break;
2857  }
2858  }
2859 
2860  // set dynamic, non-optimized model parameters
2861  //
2862 
2863  aph_s[ip] = g->aph_s;
2864  adg_s[ip] = g->adg_s;
2865  bbp_s[ip] = g->bbp_s;
2866  uadg_s[ip] = g->uadg_s;
2867  ubbp_s[ip] = g->ubbp_s;
2868 
2869 
2870  // get starting chlorophyll from default algorithm
2871 
2872  g->chl = MAX(MIN(l2rec->chl[ip], chl_max), chl_min);
2873 
2874  //Get starting chlorophyll uncertainties for default algorithm
2875  //If not computed, set to zero
2876  if (uncertainty) {
2877  g->uchl = l2rec->chl_unc[ip];
2878  } else {
2879  g->uchl = 0.0;
2880  }
2881 
2882  switch (g->rrs_opt) {
2883  case RRSFOQ:
2884  foqint_morel(input->fqfile, wave, nwave, 0.0, 0.0, 0.0, g->chl, foq);
2885  for (iw = 0; iw < g->nwave; iw++) {
2886  g->foq[iw] = foq[g->bindx[iw]];
2887  }
2888  break;
2889  }
2890 
2891  // aph function
2892 
2893  switch (g->aph_opt) {
2894  case APHBRICAUD:
2895  // replaces default tabbed values with Bricaud chl-based values
2896  for (iw = 0; iw < nwave; iw++) {
2897  g->aph_tab_w[iw] = wave[iw];
2898  // get the aph* normalization factor - the 0.055 is a "typical"
2899  // aph* value for 443nm
2900  aph_norm = 0.055 / get_aphstar(443., BANDW, APHBRICAUD, g->chl);
2901  g->aph_tab_s[0][iw] = aph_norm * get_aphstar(wave[iw], BANDW, APHBRICAUD, g->chl);
2902 
2903  //Estimate aph* Bricaud standard uncertainty
2904  dudtau = -0.055/pow(get_aphstar(443., BANDW, APHBRICAUD, g->chl),2);
2905  dtaudchl = get_aphstar_pderiv(443., BANDW, APHBRICAUD, g->chl);
2906  dudchl = dudtau*dtaudchl;
2907  dvdchl = get_aphstar_pderiv(wave[iw], BANDW, APHBRICAUD, g->chl);
2908  daphdchl = dudchl*g->aph_tab_s[0][iw] + dvdchl*aph_norm;
2909  g->uaph_tab_s[0][iw] = sqrt(pow(daphdchl*g->uchl,2));
2910  }
2911  g->aph_tab_nw = nwave;
2912  g->aph_nvec = 1;
2913  break;
2914  case APHCIOTTI:
2915  // replaces default tabbed values with Ciotti size-fraction-based values
2916  for (iw = 0; iw < nwave; iw++) {
2917  g->aph_tab_w[iw] = wave[iw];
2918  g->aph_tab_s[0][iw] = get_aphstar(wave[iw], BANDW, APHCIOTTI, g->aph_s);
2919  g->uaph_tab_s[0][iw] = 0;
2920  }
2921  g->aph_tab_nw = nwave;
2922  g->aph_nvec = 1;
2923  break;
2924  }
2925 
2926  // adg function
2927 
2928  switch (g->adg_opt) {
2929  case ADGSQAA:
2930  // update exponential based on QAA band-ratio relationship
2931  i1 = windex(443.0, wave, nwave);
2932  i2 = windex(550.0, wave, nwave);
2933  //QAA Sdg relationship derived from NOMAD data.
2934  //No Raman scattering correction applied to Rrs.
2935  Rrs1 = l2rec->Rrs[ipb2 + i1];
2936  Rrs2 = l2rec->Rrs[ipb2 + i2];
2937  uRrs1 = g->uRrs_a[i1];
2938  uRrs2 = g->uRrs_a[i2];
2939  covRrs1Rrs2 = 0.0;
2940  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2941  Rrs1 = rrs_above_to_below(Rrs1);
2942  Rrs2 = rrs_above_to_below(Rrs2);
2943  uRrs1 = rrs_above_to_below_unc(Rrs1,uRrs2);
2944  uRrs2 = rrs_above_to_below_unc(Rrs2,uRrs2);
2945  covRrs1Rrs2 = 0.0;
2946  g->adg_s = MAX(MIN(0.015 + 0.002 / (0.6 + Rrs1 / Rrs2), adg_s_max), adg_s_min);
2947  if (g->adg_s != adg_s_min || g->adg_s != adg_s_max) {
2948  g->uadg_s = calc_uadg_s(g, Rrs1, Rrs2, uRrs1, uRrs2, covRrs1Rrs2);
2949  } else {
2950  g->uadg_s = 0.0;
2951  }
2952  } else if (Rrs2 > 0.0) {
2953  g->adg_s = adg_s_min;
2954  g->uadg_s = 0.0;
2955  } else {
2956  g->adg_s = BAD_FLT;
2957  iopf[ip] |= IOPF_BADRRS;
2958  status = 1;
2959  }
2960  break;
2961  case ADGSOBPG:
2962  // update exponential based on OBPG band-ratio relationship
2963  i1 = windex(412.0, wave, nwave);
2964  i2 = windex(550.0, wave, nwave);
2965  //QAA Sdg relationship derived from NOMAD data.
2966  //No Raman scattering correction applied to Rrs.
2967  Rrs1 = l2rec->Rrs[ipb2 + i1];
2968  Rrs2 = l2rec->Rrs[ipb2 + i2];
2969  uRrs1 = g->uRrs_a[i1];
2970  uRrs2 = g->uRrs_a[i2];
2971  covRrs1Rrs2 = 0.0;
2972  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
2973  g->adg_s = MAX(MIN(0.015 + 0.0038 * log10(Rrs1 / Rrs2), adg_s_max), adg_s_min);
2974  if (g->adg_s != adg_s_min || g->adg_s != adg_s_max) {
2975  g->uadg_s = calc_uadg_s(g, Rrs1, Rrs2, uRrs1, uRrs2, covRrs1Rrs2);
2976  } else {
2977  g->uadg_s = 0.0;
2978  }
2979  } else if (Rrs2 > 0.0) {
2980  g->adg_s = adg_s_min;
2981  g->uadg_s = 0.0;
2982  } else {
2983  g->adg_s = BAD_FLT;
2984  iopf[ip] |= IOPF_BADRRS;
2985  status = 1;
2986  }
2987  break;
2988  }
2989 
2990  // bbp function
2991 
2992  switch (g->bbp_opt) {
2993  case BBPSHAL:
2994  // update power-law exponent based on HAL band-ratio relationship
2995  i1 = windex(490.0, wave, nwave);
2996  i2 = windex(550.0, wave, nwave);
2997  Rrs1 = l2rec->Rrs[ipb2 + i1] - l2rec->Rrs_raman[ipb2 + i1]; //CHECK IF RAMAN SHOULD BE HERE!!
2998  Rrs2 = l2rec->Rrs[ipb2 + i2] - l2rec->Rrs_raman[ipb2 + i2];
2999  uRrs1 = g->uRrs_a[i1];
3000  uRrs2 = g->uRrs_a[i2];
3001  covRrs1Rrs2 = 0.0;
3002  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
3003  g->bbp_s = MAX(MIN(0.8 * (Rrs1 / Rrs2) + 0.2, bbp_s_max), bbp_s_min);
3004  if (g->bbp_s != bbp_s_min || g->bbp_s != bbp_s_max) {
3005  g->ubbp_s = calc_ubbp_s(g, Rrs1, Rrs2, uRrs1, uRrs2, covRrs1Rrs2);
3006  } else {
3007  g->ubbp_s = 0.0;
3008  }
3009  } else if (Rrs2 > 0.0) {
3010  g->bbp_s = bbp_s_min;
3011  } else {
3012  g->bbp_s = BAD_FLT;
3013  iopf[ip] |= IOPF_BADRRS;
3014  status = 1;
3015  }
3016  break;
3017  case BBPSQAA:
3018  // update power-law exponent based on QAA band-ratio relationship
3019  i1 = windex(443.0, wave, nwave);
3020  i2 = windex(550.0, wave, nwave);
3021  //QAA bbp_s relationship derived from NOMAD data.
3022  //No Raman scattering correction applied to Rrs.
3023  Rrs1 = l2rec->Rrs[ipb2 + i1];
3024  Rrs2 = l2rec->Rrs[ipb2 + i2];
3025  uRrs1 = g->uRrs_a[i1];
3026  uRrs2 = g->uRrs_a[i2];
3027  covRrs1Rrs2 = 0.0;
3028 
3029  if (Rrs1 > 0.0 && Rrs2 > 0.0) {
3030  Rrs1 = rrs_above_to_below(Rrs1);
3031  Rrs2 = rrs_above_to_below(Rrs2);
3032  uRrs1 = rrs_above_to_below_unc(Rrs1, uRrs1);
3033  uRrs2 = rrs_above_to_below_unc(Rrs2, uRrs2);
3034  g->bbp_s = MAX(MIN(2.0 * (1.0 - 1.2 * exp(-0.9 * Rrs1 / Rrs2)), bbp_s_max), bbp_s_min);
3035  if (g->bbp_s != bbp_s_min || g->bbp_s != bbp_s_max) {
3036  g->ubbp_s = calc_ubbp_s(g, Rrs1, Rrs2, uRrs1, uRrs2, covRrs1Rrs2);
3037  } else {
3038  g->ubbp_s = 0.0;
3039  }
3040  } else if (Rrs2 > 0.0) {
3041  g->bbp_s = bbp_s_min;
3042  } else {
3043  g->bbp_s = BAD_FLT;
3044  iopf[ip] |= IOPF_BADRRS;
3045  status = 1;
3046  }
3047  break;
3048  case BBPSCIOTTI:
3049  // update power-law exponent based on Ciotti chl relationship
3050  if (l2rec->chl[ip] > 0.0) {
3051  g->bbp_s = MAX(MIN(1.0 - 0.768 * log10(g->chl), bbp_s_max), bbp_s_min);
3052  if (g->bbp_s != bbp_s_min || g->bbp_s != bbp_s_max) {
3053  g->ubbp_s = calc_ubbp_s(g, 0, 0, 0, 0, 0);
3054  } else {
3055  g->ubbp_s = 0.0;
3056  }
3057  } else {
3058  g->bbp_s = BAD_FLT;
3059  iopf[ip] |= IOPF_BADRRS;
3060  status = 1;
3061  }
3062  break;
3063  case BBPSMM01:
3064  // update power-law exponent based on Morel chl relationship
3065  if (l2rec->chl[ip] > 0.0) {
3066  g->bbp_s = MAX(MIN(0.5 * (0.3 - log10(g->chl)), bbp_s_max), bbp_s_min);
3067  if (g->bbp_s != bbp_s_min || g->bbp_s != bbp_s_max) {
3068  g->ubbp_s = calc_ubbp_s(g, 0, 0, 0, 0, 0);
3069  } else {
3070  g->ubbp_s = 0.0;
3071  }
3072  } else {
3073  g->bbp_s = BAD_FLT;
3074  iopf[ip] |= IOPF_BADRRS;
3075  status = 1;
3076  }
3077  break;
3078  case BBPSLAS:
3079  // update power-law exponent based on Loisel & Stramski model
3080  if ((g->bbp_s = get_bbp_las_eta(l2rec, ip)) == BAD_FLT) {
3081  iopf[ip] |= IOPF_BADRRS;
3082  status = 1;
3083  }
3084  g->ubbp_s = 0.0;
3085  printf("GIOP message: bbp slope uncertainties not computed during interface to LAS model.\n");
3086  printf(" ubbp_s uncertainty values set to: %f.\n", g->ubbp_s);
3087  g->bbp_nvec = 1;
3088  break;
3089  case BBPLAS:
3090  case BBPLASFIX:
3091  // replaces default tabbed values with results from Loisel & Stramski model
3092  if (get_bbp_las(l2rec, ip, g->bbp_tab_w, g->bbp_tab_s[0], g->bbp_tab_nw) == 0) {
3093  iopf[ip] |= IOPF_BADRRS;
3094  status = 1;
3095  }
3096  g->ubbp_s = 0.0;
3097  printf("GIOP message: bbp slope uncertainties not computed during interface to LAS model.\n");
3098  printf(" ubbp_s uncertainty values set to: %f.\n", g->ubbp_s);
3099  g->bbp_nvec = 1;
3100  break;
3101  case BBPQAAFIX:
3102  // replaces default tabbed values with results from QAA model
3103  if (get_bbp_qaa(l2rec, ip, g->bbp_tab_w, g->bbp_tab_s[0], g->bbp_tab_nw) == 0) {
3104  iopf[ip] |= IOPF_BADRRS;
3105  status = 1;
3106  }
3107  g->ubbp_s = 0.0;
3108  printf("GIOP message: bbp slope uncertainties not computed during interface to QAA model.\n");
3109  printf(" ubbp_s uncertainty values set to: %f.\n", g->ubbp_s);
3110  g->bbp_nvec = 1;
3111  break;
3112  case BBPLH:
3113  if (get_bbp_lh(l2rec, g, ipb2, g->bbp_tab_w, g->bbp_tab_s[0], g->ubbp_tab_s[0], g->bbp_tab_nw) == 0) {
3114  iopf[ip] |= IOPF_BADRRS;
3115  status = 1;
3116  }
3117  g->bbp_nvec = 1;
3118  break;
3119  case BBPCHL:
3120  if (get_bbp_chl(l2rec, g, ipb2, g->bbp_tab_w, g->bbp_tab_s[0], g->ubbp_tab_s[0], g->bbp_tab_nw) == 0) {
3121  iopf[ip] |= IOPF_BADRRS;
3122  status = 1;
3123  }
3124  g->bbp_nvec = 1;
3125  break;
3126  }
3127 
3128  // save updated model config, flag and skip if unable to compute
3129 
3130  aph_s[ip] = g->aph_s;
3131  adg_s[ip] = g->adg_s;
3132  bbp_s[ip] = g->bbp_s;
3133  uadg_s[ip] = g->uadg_s;
3134  ubbp_s[ip] = g->ubbp_s;
3135 
3136  if (status != 0) {
3137  iopf[ip] |= IOPF_BADRRS;
3138  continue;
3139  }
3140 
3141  // convert to subsurface reflectance
3142 
3143  for (iw = 0; iw < g->nwave; iw++) {
3144  ib = g->bindx[iw];
3145  Rrs_a[iw] = l2rec->Rrs[ipb2 + ib] - l2rec->Rrs_raman[ipb2 + ib];
3146  Rrs_b[iw] = rrs_above_to_below(Rrs_a[iw]);
3147  uRrs_b[iw] = rrs_above_to_below_unc(Rrs_a[iw],g->uRrs_a[iw]);
3148 
3149  wts [iw] = g->wts[iw]; // /rrs_above_to_below(Rrs_a[iw])*sqrt(l2rec->nobs[ip]);
3150  //wts[iw] = 1.0 / pow(uRrs_b[iw], 2); // testing with fitting weights
3151  if (Rrs_b[iw] > 0.0) {
3152  bndcnt++;
3153  }
3154  }
3155 
3156  // if less than npar valid reflectances, flag and skip pixel */
3157 
3158  if (bndcnt < g->npar) {
3159  iopf[ip] |= IOPF_BADRRS;
3160  continue;
3161  }
3162 
3163  // initialize model parameters
3164 
3165  giop_ctl_start(g, g->chl);
3166 
3167  for (ipar = 0; ipar < g->npar; ipar++) {
3168  par[ipar] = g->par[ipar]; // model params
3169  par[ipar + g->npar] = 0.0; // model param errors
3170  upar[ipar] = 0.0; // uncertainty model param errors
3171  }
3172 
3173  // run model optimization for this pixel
3174 
3175  switch (g->fit_opt) {
3176  case AMOEBA:
3177  status = fit_giop_amb(g, Rrs_b, wts, par, Rrs_f, &itercnt);
3178  break;
3179  case LEVMARQ:
3180  status = fit_giop_lm(g, Rrs_b, wts, par, &chi, &itercnt);
3181  break;
3182  case SVDFIT:
3183  status = fit_giop_svd(g, Rrs_b, wts, par);
3184  break;
3185  case SVDSIOP:
3186  status = fit_giop_svd_siop(g, Rrs_b, wts, par, &chi);
3187  /*If an optimal svd_siop solution not reached*/
3188  if (status == -99) {
3189  siop_num[ip] = -1;
3190  }
3191  break;
3192  default:
3193  printf("%s Line %d: Unknown optimization method for GIOP %d\n",
3194  __FILE__, __LINE__, g->fit_opt);
3195  exit(1);
3196  break;
3197  }
3198 
3199  //Compute fit parameter uncertainties
3200 
3201  if (status == 0) {
3202  calc_par_unc(g, par, uRrs_b, upar);
3203  }
3204 
3205  //if no solution, flag as IOP prod failure
3206 
3207  if (status != 0) {
3208  iopf[ip] |= IOPF_FAILED;
3209  if (itercnt >= g->maxiter)
3210  iopf[ip] |= IOPF_MAXITER;
3211  continue;
3212  }
3213 
3214  //save parameter uncertainties to par array
3215 
3216  for (ipar = 0; ipar < g->npar; ipar++) {
3217  if (isfinite(par[ipar + g->npar]))
3218  //estimate combined model misfit uncertainty + radiometric uncertainty
3219  par[ipar+g->npar] = sqrt(pow(upar[ipar],2) +pow(par[ipar+g->npar],2));
3220  }
3221 
3222  // save final params
3223 
3224  for (ipar = 0; ipar < g->npar; ipar++) {
3225  if (isfinite(par[ipar]))
3226  fit_par[ip][ipar] = par[ipar];
3227  }
3228 
3229  chisqr[ip] = chi;
3230 
3231  // evaluate model at fitted bands
3232 
3233  switch (g->fit_opt) {
3234  case SVDSIOP:
3235  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);
3236  break;
3237  default:
3238  giop_model(g, par, g->nwave, g->wave, g->aw, g->bbw, g->foq, aph1, adg1, bbp1, Rrs_f, NULL, NULL);
3239  break;
3240  }
3241 
3242 
3243  // bogus evaluation, flag and skip
3244 
3245  for (iw = 0; iw < g->nwave; iw++) {
3246  if (!isfinite(Rrs_f[iw])) {
3247  iopf[ip] |= IOPF_NAN;
3248  break;
3249  }
3250  }
3251  if (iopf[ip] & IOPF_NAN) {
3252  iopf[ip] |= IOPF_FAILED;
3253  continue;
3254  }
3255 
3256  // check goodness of fit
3257 
3258  rrsdiff[ip] = 0.0;
3259  for (iw = 0; iw < g->nwave; iw++) {
3260  if (g->wave[iw] >= 400 && g->wave[iw] <= 600) {
3261  if (fabs(Rrs_b[iw]) > 1e-7 && fabs(Rrs_f[iw] - Rrs_b[iw]) > 1e-5)
3262  rrsdiff[ip] += fabs(Rrs_f[iw] - Rrs_b[iw]) / fabs(Rrs_b[iw]);
3263  }
3264  }
3265  rrsdiff[ip] /= g->nwave;
3266  if (rrsdiff[ip] > input->giop_rrs_diff) {
3267  iopf[ip] |= IOPF_RRSDIFF;
3268  }
3269 
3270  //if (status == 0) {
3271  // calc_par_unc(g, par, rrs_diff, upar);
3272  //}
3273 
3274  // store in static globals
3275 
3276  switch (g->fit_opt) {
3277  case SVDSIOP:
3278  giop_model_iterate(g, par, nwave, wave, aw, bbw, foq, aph1, adg1, bbp1, acdom1, anap1, bbph1, bbnap1, Rrs_f, NULL, NULL);
3279 
3280  mRrs[ipb + ib] = rrs_below_to_above(Rrs_f[ib]) + l2rec->Rrs_raman[ipb2 + ib];
3281 
3282  for (ib = 0; ib < nwave; ib++) {
3283 
3284  if (isfinite(aph1[ib])) {
3285  aph[ipb + ib] = aph1[ib];
3286  aph[ierr + ib] = aph1[ib + nwave];
3287  } else {
3288  iopf[ip] |= IOPF_NAN;
3289  }
3290  if (isfinite(acdom1[ib])) {
3291  acdom[ipb + ib] = acdom1[ib];
3292  acdom[ierr + ib] = acdom1[ib + nwave];
3293  } else {
3294  iopf[ip] |= IOPF_NAN;
3295  }
3296  if (isfinite(anap1[ib])) {
3297  anap[ipb + ib] = anap1[ib];
3298  anap[ierr + ib] = anap1[ib + nwave];
3299  } else {
3300  iopf[ip] |= IOPF_NAN;
3301  }
3302  if (isfinite(aph1[ib]) && isfinite(acdom1[ib]) && isfinite(anap1[ib])) {
3303  adg[ipb + ib] = acdom1[ib] + anap1[ib];
3304  adg[ierr + ib] = pow(pow(adg1[ib + nwave],2) + pow(acdom1[ib + nwave],2),0.5);
3305  a[ipb + ib] = aw[ib] + aph1[ib] + acdom1[ib] + anap1[ib];
3306  a[ierr + ib] = pow( pow(aph1[ib + nwave],2) + pow(adg1[ib + nwave],2) + pow(acdom1[ib + nwave],2) ,0.5 );
3307  } else {
3308  iopf[ip] |= IOPF_NAN;
3309  }
3310  if (isfinite(bbnap1[ib])) {
3311  bbnap [ipb + ib] = bbnap1[ib];
3312  bbnap [ierr + ib] = bbnap1[ib + nwave];
3313  } else {
3314  iopf[ip] |= IOPF_NAN;
3315  }
3316  if (isfinite(bbph1[ib])) {
3317  bbph [ipb + ib] = bbph1[ib];
3318  bbph [ierr + ib] = bbph1[ib + nwave];
3319  } else {
3320  iopf[ip] |= IOPF_NAN;
3321  }
3322  if (isfinite(bbph1[ib]) && isfinite(bbnap1[ib])) {
3323  bbp[ipb + ib] = bbph1[ib] + bbnap1[ib];
3324  bbp[ierr + ib] = pow(pow(bbph1[ib + nwave],2) + pow(bbnap1[ib + nwave],2),0.5);
3325  bb [ipb + ib] = bbw[ib] + bbph1[ib] + bbnap1[ib];
3326  bb [ierr + ib] = pow(pow(bbph1[ib + nwave],2) + pow(bbnap[ib + nwave],2),0.5);
3327  } else {
3328  iopf[ip] |= IOPF_NAN;
3329  }
3330  }
3331 
3332  /* check IOP ranges and set flags */
3333  set_iop_flag(wave, nwave, &a[ipb], &aph[ipb], &adg[ipb],
3334  &bb[ipb], &bbp[ipb], &iopf[ip]);
3335 
3336  /*Optimal SIOP combination (for iterative aLMI solution)*/
3337  siop_num[ip] = 1 + g->siopIdx;
3338 
3339  /* aLMI compute chlorophyll */
3340  chl [ip] = par[0];
3341 
3342  break;
3343 
3344  default:
3345  //Fit model to full visible bandset
3346  giop_model(g, par, nwave, wave, aw, bbw, foq, aph1, adg1, bbp1, Rrs_f, NULL, NULL);
3347 
3348  for (ib = 0; ib < nwave; ib++) {
3349 
3350  mRrs[ipb + ib] = rrs_below_to_above(Rrs_f[ib]) + l2rec->Rrs_raman[ipb2 + ib];
3351 
3352  if (isfinite(aph1[ib])) {
3353  aph[ipb + ib] = aph1[ib];
3354  aph[ierr + ib] = aph1[ib + nwave];
3355  } else {
3356  iopf[ip] |= IOPF_NAN;
3357  }
3358  if (isfinite(adg1[ib])) {
3359  adg[ipb + ib] = adg1[ib];
3360  adg[ierr + ib] = adg1[ib + nwave];
3361  } else {
3362  iopf[ip] |= IOPF_NAN;
3363  }
3364  if (isfinite(aph1[ib]) && isfinite(adg1[ib])) {
3365  a[ipb + ib] = aw[ib] + aph1[ib] + adg1[ib];
3366  a[ierr + ib] = pow( pow(aph1[ib + nwave],2) + pow(adg1[ib + nwave],2), 0.5 );
3367  a_unc[ipb + ib] = pow( pow(aph1[ib + nwave],2) + pow(adg1[ib + nwave],2), 0.5 );
3368 
3369  } else {
3370  iopf[ip] |= IOPF_NAN;
3371  }
3372  if (isfinite(bbp1[ib])) {
3373  bbp[ipb + ib] = bbp1[ib];
3374  bbp[ierr + ib] = bbp1[ib + nwave];
3375  bb [ipb + ib] = bbw[ib] + bbp1[ib];
3376  bb [ierr + ib] = bbp1[ib + nwave];
3377  bb_unc[ipb + ib] = bbp1[ib + nwave];
3378  } else {
3379  iopf[ip] |= IOPF_NAN;
3380  }
3381  }
3382  // check IOP ranges and set flags
3383 
3384  set_iop_flag(wave, nwave, &a[ipb], &aph[ipb], &adg[ipb],
3385  &bb[ipb], &bbp[ipb], &iopf[ip]);
3386 
3387  // compute chlorophyll
3388 
3389  chl [ip] = giop_chl(g, iopf[ip], par, &chl[ip + l1rec->npix]);
3390  break;
3391  }
3392 
3393 
3394  iter[ip] = itercnt;
3395  }
3396 
3397 
3398 
3399  // fail pixels where any flags were set
3400 
3401  for (ip = 0; ip < l1rec->npix; ip++)
3402  if (iopf[ip] != 0) l1rec->flags[ip] |= PRODFAIL;
3403 
3404 
3405  LastRecNum = l1rec->iscan;
3406  free(Rrs_a);
3407  free(Rrs_b);
3408  free(Rrs_f);
3409  free(wts);
3410  free(par);
3411  free(upar);
3412  free(uRrs_a);
3413  free(uRrs_b);
3414 }
3415 
3416 
3417 /*----------------------------------------------------------------------*/
3418 /* extract_band_3d() - utility function to extract selected band subset */
3419 
3420 /*----------------------------------------------------------------------*/
3421 static void extract_band_3d(float *in_buf, float *out_buf, int numPixels, int numBands) {
3422  float * out_ptr = out_buf;
3423  for (int pix = 0; pix < numPixels; pix++) {
3424  float *in_ptr = in_buf + pix * numBands;
3425  for (int band_3d = 0; band_3d < input->nwavelengths_3d; band_3d++) {
3426  int band = input->wavelength_3d_index[band_3d];
3427  *out_ptr = in_ptr[band];
3428  out_ptr++;
3429  }
3430  }
3431 }
3432 
3433 /*----------------------------------------------------------------------*/
3434 /* extract_band_3d_unc() - utility function to extract selected band */
3435 /* subset of spectral product uncertainties */
3436 
3437 /*----------------------------------------------------------------------*/
3438 static void extract_band_3d_unc(float *in_buf, float *out_buf, int numPixels, int numBands) {
3439  float * out_ptr = out_buf;
3440  for (int pix = 0; pix < numPixels; pix++) {
3441  float *in_ptr = in_buf + (pix * numBands + numPixels*numBands) ;
3442  for (int band_3d = 0; band_3d < input->nwavelengths_3d; band_3d++) {
3443  int band = input->wavelength_3d_index[band_3d];
3444  *out_ptr = in_ptr[band];
3445  out_ptr++;
3446  }
3447  }
3448 }
3449 
3450 /* ------------------------------------------------------------------- */
3451 /* get_giop() - returns requested GIOP product for full scanline */
3452 
3453 /* ------------------------------------------------------------------- */
3454 void get_giop(l2str *l2rec, l2prodstr *p, float prod[]) {
3455  int prodID = p->cat_ix;
3456  int ib = p->prod_ix;
3457 
3458  int32_t ip, ipb, ierr;
3459  l1str *l1rec = l2rec->l1rec;
3460 
3461  /*Before running GIOP, check if valid output products selected*/
3462  /*Note: aCDOM, aNAP, bbPh, bbNAP are only valid choices for giop_fit_opt SVDSIOP*/
3463  switch (input->giop_fit_opt) {
3464  case SVDSIOP:
3465  break;
3466  default:
3467  switch (prodID) {
3468  case CAT_acdom_giop:
3469  case CAT_anap_giop:
3470  case CAT_bbph_giop:
3471  case CAT_bbnap_giop:
3472  case CAT_acdom_unc_giop:
3473  case CAT_anap_unc_giop:
3474  case CAT_bbph_unc_giop:
3475  case CAT_bbnap_unc_giop:
3476  printf("-E- %s line %d : products acdom, anap, bbph and bbnap are only applicable with giop_fit_opt=SVDSIOP.\n",
3477  __FILE__, __LINE__);
3478  exit(1);
3479  break;
3480  }
3481  }
3482 
3483  if (!giop_ran(l1rec->iscan))
3484  run_giop(l2rec);
3485 
3486  if(p->rank == 3) {
3487  switch (prodID) {
3488 
3489  case CAT_aph_giop:
3490  extract_band_3d(aph, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3491  break;
3492 
3493  case CAT_adg_giop:
3494  extract_band_3d(adg, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3495  break;
3496 
3497  case CAT_acdom_giop:
3498  extract_band_3d(acdom, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3499  break;
3500 
3501  case CAT_anap_giop:
3502  extract_band_3d(anap, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3503  break;
3504 
3505  case CAT_bbp_giop:
3506  extract_band_3d(bbp, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3507  break;
3508 
3509  case CAT_a_giop:
3510  extract_band_3d(a, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3511  break;
3512 
3513  case CAT_bb_giop:
3514  extract_band_3d(bb, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3515  break;
3516 
3517  case CAT_bbph_giop:
3518  extract_band_3d(bbph, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3519  break;
3520 
3521  case CAT_bbnap_giop:
3522  extract_band_3d(bbph, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3523  break;
3524 
3525  case CAT_a_unc_giop:
3526  extract_band_3d_unc(a, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3527  break;
3528 
3529  case CAT_aph_unc_giop:
3530  extract_band_3d_unc(aph, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3531  break;
3532 
3533  case CAT_adg_unc_giop:
3534  extract_band_3d_unc(adg, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3535  break;
3536 
3537  case CAT_acdom_unc_giop:
3538  extract_band_3d_unc(acdom, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3539  break;
3540 
3541  case CAT_anap_unc_giop:
3542  extract_band_3d_unc(anap, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3543  break;
3544 
3545  case CAT_bb_unc_giop:
3546  extract_band_3d_unc(bb, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3547  break;
3548 
3549  case CAT_bbp_unc_giop:
3550  extract_band_3d_unc(bbp, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3551  break;
3552 
3553  case CAT_bbph_unc_giop:
3554  extract_band_3d_unc(bbph, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3555  break;
3556 
3557  case CAT_bbnap_unc_giop:
3558  extract_band_3d_unc(bbnap, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
3559  break;
3560 
3561  default:
3562  printf("-E- %s line %d : product ID %d passed to GIOP does not support 3D data sets.\n",
3563  __FILE__, __LINE__, prodID);
3564  exit(1);
3565  } // switch prodID
3566 
3567  } else {
3568 
3569  for (ip = 0; ip < l1rec->npix; ip++) {
3570 
3571  ipb = ip * l1rec->l1file->nbands + ib;
3572  ierr = l1rec->npix * l1rec->l1file->nbands + ipb;
3573 
3574  switch (prodID) {
3575 
3576  case CAT_mRrs_giop:
3577  prod[ip] = (float) mRrs[ipb];
3578  break;
3579 
3580  case CAT_aph_giop:
3581  prod[ip] = (float) aph[ipb];
3582  break;
3583 
3584  case CAT_adg_giop:
3585  prod[ip] = (float) adg[ipb];
3586  break;
3587 
3588  case CAT_bbp_giop:
3589  prod[ip] = (float) bbp[ipb];
3590  break;
3591 
3592  case CAT_a_giop:
3593  prod[ip] = (float) a[ipb];
3594  break;
3595 
3596  case CAT_bb_giop:
3597  prod[ip] = (float) bb[ipb];
3598  break;
3599 
3600  case CAT_acdom_giop:
3601  prod[ip] = (float) acdom[ipb];
3602  break;
3603 
3604  case CAT_anap_giop:
3605  prod[ip] = (float) anap[ipb];
3606  break;
3607 
3608  case CAT_bbph_giop:
3609  prod[ip] = (float) bbph[ipb];
3610  break;
3611 
3612  case CAT_bbnap_giop:
3613  prod[ip] = (float) bbnap[ipb];
3614  break;
3615 
3616  case CAT_chl_giop:
3617  prod[ip] = (float) chl[ip];
3618  break;
3619 
3620  case CAT_opt_siop_giop:
3621  prod[ip] = (int) siop_num[ip];
3622  break;
3623 
3624  case CAT_aph_unc_giop:
3625  prod[ip] = (float) aph[ierr];
3626  break;
3627 
3628  case CAT_adg_unc_giop:
3629  prod[ip] = (float) adg[ierr];
3630  break;
3631 
3632  case CAT_acdom_unc_giop:
3633  prod[ip] = (float) adg[ierr];
3634  break;
3635 
3636  case CAT_anap_unc_giop:
3637  prod[ip] = (float) adg[ierr];
3638  break;
3639 
3640  case CAT_bbp_unc_giop:
3641  prod[ip] = (float) bbp[ierr];
3642  break;
3643 
3644  case CAT_bbph_unc_giop:
3645  prod[ip] = (float) bbph[ierr];
3646  break;
3647 
3648  case CAT_bbnap_unc_giop:
3649  prod[ip] = (float) bbnap[ierr];
3650  break;
3651 
3652  case CAT_a_unc_giop:
3653  prod[ip] = (float) a[ierr];
3654  break;
3655 
3656  case CAT_bb_unc_giop:
3657  prod[ip] = (float) bb[ierr];
3658  break;
3659 
3660  case CAT_chl_unc_giop:
3661  prod[ip] = (float) chl[ip + l1rec->npix];
3662  break;
3663 
3664  case CAT_aphs_giop:
3665  prod[ip] = (float) aph_s[ip];
3666  break;
3667 
3668  case CAT_adgs_giop:
3669  prod[ip] = (float) adg_s[ip];
3670  break;
3671 
3672  case CAT_adgs_unc_giop:
3673  prod[ip] = (float) uadg_s[ip];
3674  break;
3675 
3676  case CAT_bbps_giop:
3677  prod[ip] = (float) bbp_s[ip];
3678  break;
3679 
3680  case CAT_bbps_unc_giop:
3681  prod[ip] = (float) ubbp_s[ip];
3682  break;
3683 
3684  case CAT_rrsdiff_giop:
3685  prod[ip] = (float) rrsdiff[ip];
3686  break;
3687 
3688  case CAT_chisqr_giop:
3689  prod[ip] = (float) chisqr[ip];
3690  break;
3691 
3692  case CAT_fitpar_giop:
3693  if (ib >= max_npar) {
3694  printf("-E- %s line %d : output request for GIOP fit parameter %d exceeds number of fit parameters %d.\n",
3695  __FILE__, __LINE__, ib, max_npar);
3696  exit(1);
3697  }
3698  prod[ip] = (float) fit_par[ip][ib];
3699  break;
3700 
3701  default:
3702  printf("-E- %s line %d : erroneous product ID %d passed to GIOP.\n",
3703  __FILE__, __LINE__, prodID);
3704  exit(1);
3705  }
3706  } // for ip
3707  } // if rank != 3
3708 
3709  return;
3710 }
3711 
3712 
3713 /* ------------------------------------------------------------------- */
3714 /* get_iter_giop() - returns iteration count */
3715 
3716 /* ------------------------------------------------------------------- */
3717 int16 *get_iter_giop(l2str *l2rec) {
3718  if (!giop_ran(l2rec->l1rec->iscan))
3719  run_giop(l2rec);
3720 
3721  return iter;
3722 }
3723 
3724 
3725 /* ------------------------------------------------------------------- */
3726 /* get_flags_giop() - returns iteration count */
3727 
3728 /* ------------------------------------------------------------------- */
3729 int16 *get_flags_giop(l2str *l2rec) {
3730  if (!giop_ran(l2rec->l1rec->iscan))
3731  run_giop(l2rec);
3732 
3733  return iopf;
3734 }
3735 
3736 
3737 /* ------------------------------------------------------------------- */
3738 /* Interface to convl12() to return GIOP iops */
3739 
3740 /* ------------------------------------------------------------------- */
3741 void iops_giop(l2str *l2rec) {
3742  int32_t ib, ip, ipb, ipb2;
3743 
3744  int32_t nbands = l2rec->l1rec->l1file->nbands;
3745  int32_t npix = l2rec->l1rec->npix;
3746 
3747  if (!giop_ran(l2rec->l1rec->iscan))
3748  run_giop(l2rec);
3749 
3750  for (ip = 0; ip < npix; ip++) {
3751  for (ib = 0; ib < nbands; ib++) {
3752  ipb2 = ip * nbands + ib;
3753  ipb = ip * nbands + ib;
3754  l2rec->a [ipb2] = (float) a[ipb];
3755  l2rec->bb[ipb2] = (float) bb[ipb];
3756  }
3757  }
3758 
3759  return;
3760 }
3761 
3762 
3763 // PJW 9 Jan 2013
3764 /* ------------------------------------------------------------------- */
3765 /* give external access to local *chl */
3766 
3767 /* ------------------------------------------------------------------- */
3768 float* giop_get_chl_pointer() {
3769  return chl;
3770 }
3771 /* ------------------------------------------------------------------- */
3772 /* give external access to local *adg */
3773 
3774 /* ------------------------------------------------------------------- */
3775 float* giop_get_adg_pointer() {
3776  return adg;
3777 }
3778 /* ------------------------------------------------------------------- */
3779 /* give external access to local *bbp */
3780 
3781 /* ------------------------------------------------------------------- */
3782 float* giop_get_bbp_pointer() {
3783  return bbp;
3784 }
3785 /* ------------------------------------------------------------------- */
3786 /* give external access to local *aph */
3787 
3788 /* ------------------------------------------------------------------- */
3789 float* giop_get_aph_pointer() {
3790  return aph;
3791 }
3792 /* ------------------------------------------------------------------- */
3793 /* give external access to local **fit_par */
3794 
3795 /* ------------------------------------------------------------------- */
3796 float** giop_get_fitpar_pointer() {
3797  return fit_par;
3798 }
3799 
3800 /* ------------------------------------------------------------------- */
3801 /* give external access to local *bbp_s */
3802 
3803 /* ------------------------------------------------------------------- */
3804 float* giop_get_bbp_s_pointer() {
3805  return bbp_s;
3806 }
3807 
3808 /* ------------------------------------------------------------------- */
3809 float* giop_get_ubbp_s_pointer() {
3810  return ubbp_s;
3811 }
3812 
3813 /* ------------------------------------------------------------------- */
3814 float* giop_get_a_unc_pointer() {
3815  return a_unc;
3816 }
3817 
3818 /* ------------------------------------------------------------------- */
3819 float* giop_get_bb_unc_pointer() {
3820  return bb_unc;
3821 }
freeArray
void freeArray(void **a, int32_t m)
Definition: giop.c:115
ANAPTAB
#define ANAPTAB
Definition: giop.h:45
MAX
#define MAX(A, B)
Definition: swl0_utils.h:25
giop_ctl_init
void giop_ctl_init(giopstr *g, int nwave, float wave[], float aw[], float bbw[])
Definition: giop.c:266
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
BBPCHL
#define BBPCHL
Definition: giop.h:31
g
e g
Definition: HOWTO_Add_a_product.txt:5
giop_get_ubbp_s_pointer
float * giop_get_ubbp_s_pointer()
Definition: giop.c:3809
j
int j
Definition: decode_rs.h:73
get_aphstar
float get_aphstar(float wave, int dwave, int ftype, float proxy)
Definition: aph.c:372
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:1576
extract_band_3d
float * extract_band_3d(float *in_buf, float *out_buf)
Definition: prodgen.c:67
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:58
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
calc_uadg_s
float calc_uadg_s(giopstr *g, float Rrs1, float Rrs2, float uRrs1, float uRrs2, float covRrs1Rrs2)
Definition: giop.c:1241
giop_get_aph_pointer
float * giop_get_aph_pointer()
Definition: giop.c:3789
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
BBPSMM01
#define BBPSMM01
Definition: giop.h:25
BBPQAAFIX
#define BBPQAAFIX
Definition: giop.h:29
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
calc_par_unc
void calc_par_unc(giopstr *g, double par[], double uRrs[], double upar[])
Definition: giop.c:1399
color_dtdb.v
v
Definition: color_dtdb.py:430
fit_giop_amb
int fit_giop_amb(giopstr *g, double Rrs[], double wts[], double par[], double Rrs_fit[], int16 *itercnt)
Definition: giop.c:1094
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:53
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:1157
SVDSIOP
#define SVDSIOP
Definition: giop.h:9
CAT_anap_giop
#define CAT_anap_giop
Definition: l2prod.h:357
ACDOMNONE
#define ACDOMNONE
Definition: giop.h:42
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
giop_get_adg_pointer
float * giop_get_adg_pointer()
Definition: giop.c:3775
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
giop_int_tab_file
void giop_int_tab_file(char *file, char *ufile, int nx, float *x, float **y, float **uy)
Definition: giop.c:138
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:3768
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:1969
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:54
CAT_aph_unc_giop
#define CAT_aph_unc_giop
Definition: l2prod.h:201
APHGAUSS
#define APHGAUSS
Definition: giop.h:35
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:55
giop_amb
double giop_amb(FITSTRUCT *ambdata, double par[])
Definition: giop.c:1074
A
const float A
Definition: calc_par.c:100
rrs_above_to_below_unc
float rrs_above_to_below_unc(float Rrs, float uRrs)
Definition: giop.c:1978
ADGSQAA
#define ADGSQAA
Definition: giop.h:14
fit_data_str::w
double * w
Definition: giop.c:54
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:643
giop_get_a_unc_pointer
float * giop_get_a_unc_pointer()
Definition: giop.c:3814
rrs_below_to_above
float rrs_below_to_above(float rrs_s)
Definition: giop.c:1989
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:34
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
get_aphstar_pderiv
float get_aphstar_pderiv(float wave, int dwave, int ftype, float proxy)
Definition: aph.c:408
hico.value_locate.x
list x
Definition: value_locate.py:64
rrs_below_to_above_unc
float rrs_below_to_above_unc(float rrs_s, float urrs_s)
Definition: giop.c:1998
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:911
calc_ubbp_s
float calc_ubbp_s(giopstr *g, float Rrs1, float Rrs2, float uRrs1, float uRrs2, float covRrs1Rrs2)
Definition: giop.c:1266
fit_giop_lm
int fit_giop_lm(giopstr *g, double Rrs[], double wts[], double par[], double *chi, int16 *itercnt)
Definition: giop.c:1492
tol
float tol
Definition: pml_iop_tables.c:45
conv_rrs_to_555
float conv_rrs_to_555(float Rrs, float wave, float uRrs_in, float *uRrs_out)
Definition: convert_band.c:17
create_images.Rrs
def Rrs
Definition: create_images.py:92
URRSREL
#define URRSREL
Definition: giop.h:65
CAT_opt_siop_giop
#define CAT_opt_siop_giop
Definition: l2prod.h:364
giop_chl
float32 giop_chl(giopstr *g, int16 iopf, double par[], float *uchl)
Definition: giop.c:1938
URRSCALC
#define URRSCALC
Definition: giop.h:63
get_iter_giop
int16 * get_iter_giop(l2str *l2rec)
Definition: giop.c:3717
BBPHNONE
#define BBPHNONE
Definition: giop.h:50
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:3782
giop_ctl_start
void giop_ctl_start(giopstr *g, float chl)
Definition: giop.c:211
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:1232
APHCIOTTI
#define APHCIOTTI
Definition: giop.h:37
RRSFOQ
#define RRSFOQ
Definition: giop.h:59
realtype
double realtype
Definition: giop.c:50
freeDArray
void freeDArray(double **a, int32_t m)
Definition: giop.c:123
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:654
get_bbp_qaa
int get_bbp_qaa(l2str *l2rec, int ip, float tab_wave[], float tab_bbp[], int tab_nwave)
Definition: get_qaa.c:353
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
calc_pinv
gsl_matrix * calc_pinv(gsl_matrix *A, const realtype rcond)
Definition: giop.c:1305
fabs
#define fabs(a)
Definition: misc.h:93
APHBRICAUD
#define APHBRICAUD
Definition: giop.h:36
FITSTRUCT::meta
void * meta
Definition: amoeba.h:10
get_flags_giop
int16 * get_flags_giop(l2str *l2rec)
Definition: giop.c:3729
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
BBPLASFIX
#define BBPLASFIX
Definition: giop.h:28
giop_get_bb_unc_pointer
float * giop_get_bb_unc_pointer()
Definition: giop.c:3819
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
get_bbp_lh
int get_bbp_lh(l2str *l2rec, giopstr *g, int ipb2, float tab_wave[], float tab_bbp[], float tab_ubbp[], int tab_nwave)
Definition: giop.c:2009
giop_get_fitpar_pointer
float ** giop_get_fitpar_pointer()
Definition: giop.c:3796
FITSTRUCT::npnts
short int npnts
Definition: amoeba.h:3
run_giop
void run_giop(l2str *l2rec)
Definition: giop.c:2230
s
data_t s[NROOTS]
Definition: decode_rs.h:75
URRSABS
#define URRSABS
Definition: giop.h:64
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:123
giop_lm_f
int giop_lm_f(const gsl_vector *parv, void *data, gsl_vector *f)
Definition: giop.c:1228
BANDW
#define BANDW
Definition: l1.h:53
fit_giop_svd_siop
int fit_giop_svd_siop(giopstr *g, double rrs[], double wts[], double par[], double *chi)
Definition: giop.c:1691
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:3804
BBPLH
#define BBPLH
Definition: giop.h:30
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_adgs_unc_giop
#define CAT_adgs_unc_giop
Definition: l2prod.h:545
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:53
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
URRSNONE
#define URRSNONE
Definition: giop.h:62
fit_data_str
Definition: giop.c:52
get_bbp_chl
int get_bbp_chl(l2str *l2rec, giopstr *g, int ipb2, float tab_wave[], float tab_bbp[], float tab_ubbp[], int tab_nwave)
Definition: giop.c:2132
modis.modis_utils.cleanup
def cleanup(lut)
Definition: modis_utils.py:219
ANAPNONE
#define ANAPNONE
Definition: giop.h:46
ADGS
#define ADGS
Definition: giop.h:13
iops_giop
void iops_giop(l2str *l2rec)
Definition: giop.c:3741
CAT_bbps_unc_giop
#define CAT_bbps_unc_giop
Definition: l2prod.h:546
npix
int npix
Definition: get_cmp.c:28
ACDOMTAB
#define ACDOMTAB
Definition: giop.h:41
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
BBPHTAB
#define BBPHTAB
Definition: giop.h:49
get_giop
void get_giop(l2str *l2rec, l2prodstr *p, float prod[])
Definition: giop.c:3454
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

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster