ocssw V2020
gsm.c
Go to the documentation of this file.
1 /* =================================================================== */
2 /* module gsm.c - Garver-Seigel-Maritorena 2001 Bio-optical model */
3 /* */
4 /* This module contains the functions to optimize and evaluate the */
5 /* Garver-Seigel-Maritorena 2001 Bio-optical model. The optimization */
6 /* is performed for each set of Rrs values within a scanline, using */
7 /* the Amoeba simplex minimization technique (Numerical Recipes). */
8 /* */
9 /* Reference: */
10 /* Maritorena S., D.A. Siegel & A. Peterson, Optimization of a Semi- */
11 /* Analytical Ocean Color Model for Global Scale Applications, Applied */
12 /* Optics, 41(15): 2705-2714, 2002. */
13 /* */
14 /* Implementation: */
15 /* B. Franz, NASA/OCRPG/SAIC, September 2004 (rewrite of previous) */
16 /* */
17 /* =================================================================== */
18 
19 #include <stdio.h>
20 #include <math.h>
21 #include "l12_proto.h"
22 #include "amoeba.h"
23 #include <gsl/gsl_errno.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_blas.h>
26 #include <gsl/gsl_multifit_nlin.h>
27 
28 #define NPAR 3
29 #define MAXITR 500
30 
31 #define GSMDEFAULT 0
32 #define CHESAPEAKE 1
33 
34 #define AMOEBA 0
35 #define LEVMARQ 1
36 
37 static float badval = BAD_FLT;
38 
39 static int32_t nbandsvis;
40 static float grd1 = 0.0949;
41 static float grd2 = 0.0794;
42 static float *aphstar;
43 static float *adgstar;
44 static float *bbpstar;
45 static float *lambda;
46 static float *aw;
47 static float *bbw;
48 
49 static int32_t coefset = GSMDEFAULT;
50 static int32_t pixday;
51 static float32 pixlon;
52 static float32 pixlat;
53 
54 static double *dFdx;
55 static double *dxda;
56 static double *dxdb;
57 
58 static int32_t GSMRecNum = -1;
59 static double *chl;
60 static double *adg;
61 static double *bbp;
62 static int16 *iter;
63 
64 static float merit;
65 
66 struct datastruct {
67  size_t n;
68  double * y;
69  double * sigma;
70 } data;
71 
72 int gsm_ran(int recnum) {
73  if (recnum == GSMRecNum)
74  return 1;
75  else
76  return 0;
77 }
78 
79 
80 /* ------------------------------------------------------------------- */
81 /* get_cb_model - retrieves model coefficients for Chesapeake region */
82 
83 /* ------------------------------------------------------------------- */
84 void set_cb_model(float wave[], int32_t nwave, float lon, float lat, int32_t day) {
85  static int firstCall = 1;
86  static int *iwtab;
87  static int *interp;
88  static int iseason = -1;
89 
90  static int32_t ntwave = 6;
91  static float twave[6] = {412.0, 443.0, 490.0, 510.0, 555.0, 670.0};
92  static char season[4][7] = {"Spring", "Summer", "Fall", "Winter"};
93 
94  static float tcoefs[5][4][6 + 2] = {
95  { /* Upper Bay */
96  {0.02119, 0.02509, 0.01282, 0.00910, 0.00427, 0.02087, 0.01218, 0.0}, /* Spring */
97  {0.02653, 0.02979, 0.01655, 0.01208, 0.00470, 0.02122, 0.01218, 0.0}, /* Summer */
98  {0.02259, 0.02573, 0.01372, 0.01019, 0.00376, 0.02066, 0.01218, 0.0}, /* Fall */
99  {0.02259, 0.02573, 0.01372, 0.01019, 0.00376, 0.02066, 0.01218, 0.0}
100  }, /* Winter */
101  { /* Mid-Bay */
102  {0.02001, 0.02212, 0.01279, 0.00974, 0.00449, 0.01588, 0.01385, 0.0}, /* Spring */
103  {0.03345, 0.03900, 0.02318, 0.01664, 0.00609, 0.02285, 0.01385, 0.0}, /* Summer */
104  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01385, 0.0}, /* Fall */
105  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01385, 0.0}
106  }, /* Winter */
107  { /* Lower Bay */
108  {0.02001, 0.02212, 0.01279, 0.00974, 0.00449, 0.01588, 0.01330, 0.0}, /* Spring */
109  {0.03345, 0.03900, 0.02318, 0.01664, 0.00609, 0.02285, 0.01330, 0.0}, /* Summer */
110  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01330, 0.0}, /* Fall */
111  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01330, 0.0}
112  }, /* Winter */
113  { /* Inshore MAB */
114  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}, /* Spring */
115  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}, /* Summer */
116  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}, /* Fall */
117  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}
118  }, /* Winter */
119  { /* Offshore MAB */
120  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}, /* Spring */
121  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}, /* Summer */
122  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}, /* Fall */
123  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}
124  } /* Winter */
125 
126  /* aph*412,aph*443,aph*490,aph*510,aph*555,aph*670, Sdg, Sbbp */
127  };
128 
129  float *taphstar, adg_s, bbp_s;
130  int iw, iregion;
131  if (firstCall == 1) {
132  firstCall = 0;
133  if ((iwtab = (int *) calloc(nbandsvis, sizeof (int))) == NULL) {
134  printf("-E- %s line %d : error allocating memory for iwtab in gsm:set_cb_model.\n",
135  __FILE__, __LINE__);
136  exit(1);
137  }
138  if ((interp = (int *) calloc(nbandsvis, sizeof (int))) == NULL) {
139  printf("-E- %s line %d : error allocating memory for interp in gsm:set_cb_model.\n",
140  __FILE__, __LINE__);
141  exit(1);
142  }
143  }
144 
145  if (iseason < 0) {
146 
147  /* first call, determine the season */
148  if (day >= 79 && day < 172) iseason = 0; /* Spring */
149  else if (day >= 172 && day < 265) iseason = 1; /* Summer */
150  else if (day >= 265 && day < 355) iseason = 2; /* Fall */
151  else iseason = 3; /* Winter */
152 
153  printf("\nInitializing Chesapeake GSM model for %s\n", season[iseason]);
154 
155  /* determine if we need to interpolate */
156  for (iw = 0; iw < nwave; iw++) {
157  iwtab[iw] = windex(wave[iw], twave, ntwave);
158  if (fabsf(wave[iw] - twave[iwtab[iw]]) > 0.5) {
159  printf("\nGSM coefficients will be interpolated for %5.1f nm band.\n", wave[iw]);
160  interp[iw] = 1;
161  } else
162  interp[iw] = 0;
163  }
164  }
165 
166  /* determine region */
167  iregion = -1;
168  if (lat < 37.0 && lon > (-(lat + 152) / 2.5))
169  iregion = 4; /* Offshore MAB */
170  else if (lat >= 37.0 && lon > ((lat - 153.31) / 1.5385))
171  iregion = 4; /* Offshore MAB */
172  else if (lon > -75.7)
173  iregion = 3; /* Inshore MAB */
174  else if (lat < 37.6)
175  iregion = 2; /* Lower Bay */
176  else if (lat < 38.6)
177  iregion = 1; /* Mid Bay */
178  else
179  iregion = 0; /* Upper Bay */
180 
181  /* compute coefficients */
182  taphstar = &tcoefs[iregion][iseason][0];
183  adg_s = tcoefs[iregion][iseason][ntwave + 0];
184  bbp_s = tcoefs[iregion][iseason][ntwave + 1];
185  for (iw = 0; iw < nwave; iw++) {
186  if (interp[iw])
187  aphstar[iw] = linterp(twave, taphstar, ntwave, wave[iw]);
188  else
189  aphstar[iw] = taphstar[iwtab[iw]];
190  adgstar[iw] = exp(-adg_s * (wave[iw] - 443.0));
191  bbpstar[iw] = pow((443.0 / wave[iw]), bbp_s);
192  }
193 
194  return;
195 }
196 
197 
198 /* ------------------------------------------------------------------- */
199 /* gsm_amb() GSM model, set-up to work with amoeba optimization */
200 /* */
201 /* Input: */
202 /* auxdata - amoeba interface structure */
203 /* par[] - model parameters */
204 /* */
205 /* Output: */
206 /* function returns chi-squared, which is the quantity minimized. */
207 /* */
208 
209 /* ------------------------------------------------------------------- */
210 double gsm_amb(FITSTRUCT *auxdata, double par[]) {
211 
212  int iw;
213  float x, F, bb, ac;
214  double merit = 0.0;
215  int nwave = auxdata->npnts;
216 
217  if (coefset == CHESAPEAKE) {
218  set_cb_model(lambda, nwave, pixlon, pixlat, pixday);
219  }
220 
221  for (iw = 0; iw < nwave; iw++) {
222 
223  ac = aw [iw] + aphstar[iw] * par[0] + par[1] * adgstar[iw];
224  bb = bbw[iw] + par[2] * bbpstar[iw];
225  x = bb / (ac + bb);
226  F = grd1 * x + grd2 * pow(x, 2);
227 
228  auxdata->yfit[iw] = F;
229 
230  merit += pow((auxdata->y[iw] - F), 2) * auxdata->wgt[iw];
231  }
232 
233  auxdata->merit = merit;
234 
235  return (merit);
236 }
237 
238 
239 /* ------------------------------------------------------------------- */
240 /* fit_gsm_amb() - run amoeba optimization of GSM01 for one spectrum */
241 
242 /* ------------------------------------------------------------------- */
243 int fit_gsm_amb(double Rrs[], double wts[], int32_t npts, double fitparms[],
244  double Rrs_fit[], int16 *itercnt) {
245  static float tol = 1.e-6; /* fractional change in chisqr */
246  static FITSTRUCT auxdata; /* amoeba interface structure */
247  static double init[NPAR * (NPAR + 1)]; /* initial simplex */
248  static double p0[] = {.1, .02, .0001}; /* starting parameters (C,adg,bb)*/
249  static double d0[] = {5., 1.0, 0.01}; /* characteristic length */
250 
251  int i, j;
252  short isml;
253  int status = 1;
254 
255  auxdata.niter = MAXITR; /* max number of iterations allowed */
256  auxdata.nfunc = NPAR; /* number of model parameters */
257  auxdata.npnts = npts; /* number of wavelengths (Rrs values) */
258  auxdata.y = Rrs; /* Input Rrs values */
259  auxdata.wgt = wts; /* Input weights on Rrs values */
260  auxdata.yfit = Rrs_fit; /* Output model predicted Rrs values */
261 
262  /* initialize simplex with first guess model parameters */
263  for (j = 0; j < NPAR + 1; j++) for (i = 0; i < NPAR; i++) init[j * NPAR + i] = p0[i];
264  for (i = 0; i < NPAR; i++) {
265  init[NPAR + i * (NPAR + 1)] += d0[i];
266  fitparms[i] = 0.0;
267  }
268 
269  /* run optimization */
270  isml = amoeba(init, &auxdata, gsm_amb, tol);
271 
272  /* check convergence and record parameter results */
273  if (auxdata.niter > MAXITR)
274  /* failed to converge */
275  status = 0;
276 
277  for (i = 0; i < NPAR; i++)
278  fitparms[i] = init[NPAR * isml + i];
279 
280  *itercnt = auxdata.niter;
281 
282  return (status);
283 }
284 
285 
286 /* ------------------------------------------------------------------- */
287 /* run_gsm_amb() - returns requested GSM product for full scanline */
288 
289 /* ------------------------------------------------------------------- */
290 void run_gsm_amb(l2str *l2rec) {
291 
292  int32_t ip, iw;
293  double model [NPAR];
294  double *Rrs;
295  double *wts;
296  double *Rrs_fit;
297  int16 itercnt;
298  int16 bndcnt;
299  int status;
300 
301  l1str *l1rec = l2rec->l1rec;
302 
303  if ((Rrs = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
304  printf("-E- %s line %d : error allocating memory for Rrs in gsm:run_gsm_amb.\n",
305  __FILE__, __LINE__);
306  exit(1);
307  }
308  if ((wts = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
309  printf("-E- %s line %d : error allocating memory for wts in gsm:run_gsm_amb.\n",
310  __FILE__, __LINE__);
311  exit(1);
312  }
313  if ((Rrs_fit = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
314  printf("-E- %s line %d : error allocating memory for Rrs_fit in gsm:run_gsm_amb.\n",
315  __FILE__, __LINE__);
316  exit(1);
317  }
318  /* if this is a new scanline, fit model for every pixel of the scanline */
319 
320  for (ip = 0; ip < l1rec->npix; ip++) {
321 
322  status = 1;
323  bndcnt = 0;
324  chl [ip] = badval;
325  adg [ip] = badval;
326  bbp [ip] = badval;
327  iter[ip] = 0;
328 
329  if (l1rec->mask[ip] == 0) {
330 
331  /* for CB model selection */
332  pixlon = l1rec->lon[ip];
333  pixlat = l1rec->lat[ip];
334  int16_t year, day;
335  double sec;
336  unix2yds(l1rec->scantime, &year, &day, &sec);
337  pixday = (int32_t) day;
338 
339  for (iw = 0; iw < nbandsvis; iw++) {
340  wts[iw] = 1.0;
341  Rrs[iw] = l2rec->Rrs[ip * l1rec->l1file->nbands + iw];
342  Rrs[iw] = Rrs[iw] / (0.52 + 1.7 * Rrs[iw]);
343  if (Rrs[iw] <= 0.0) status = 0;
344  else bndcnt++;
345  }
346 
347  /* if less than 3 valid radiances, fail pixel */
348  if (bndcnt < 3) {
349  l1rec->flags[ip] |= PRODFAIL;
350  continue;
351  }
352 
353  /* if any negative reflectance, set warning flag */
354  /* but allow processing to continue */
355  if (status == 0) {
356  l1rec->flags[ip] |= PRODWARN;
357  }
358 
359  /* fit model for this pixel */
360  status = fit_gsm_amb(Rrs, wts, nbandsvis, model, Rrs_fit, &itercnt);
361 
362  /* save results and flag failure */
363  if (status != 1 || model[0] <= 0.0) {
364  l1rec->flags[ip] |= PRODFAIL;
365  } else {
366  chl[ip] = model[0];
367  adg[ip] = model[1];
368  bbp[ip] = model[2];
369  }
370  iter[ip] = itercnt;
371  } else {
372  l1rec->flags[ip] |= PRODFAIL;
373  }
374  }
375 
376  GSMRecNum = l1rec->iscan;
377 
378  free(Rrs);
379  free(wts);
380  free(Rrs_fit);
381  return;
382 }
383 
384 
385 /* ------------------------------------------------------------------- */
386 /* gsm_lm_f() GSM01 model, set-up to work with L.M. optimization */
387 /* */
388 /* Input: */
389 /* par[] - model parameters */
390 /* data - data to fit */
391 /* f - gsl_multifit_function */
392 /* */
393 
394 /* ------------------------------------------------------------------- */
395 int gsm_lm_f(const gsl_vector *par, void *data, gsl_vector *f) {
396  static int first = 1;
397 
398  int iw;
399  float x, F, bb, ac;
400  size_t nwave = ((struct datastruct *) data)->n;
401  double *y = ((struct datastruct *) data)->y;
402  double *sigma = ((struct datastruct *) data)->sigma;
403 
404  if (first == 1) {
405 
406  if ((dFdx = (double *) calloc(nwave, sizeof (double))) == NULL) {
407  printf("-E- %s line %d : error allocating memory for dFdx in gsm:gsm_lm_f.\n",
408  __FILE__, __LINE__);
409  exit(1);
410  }
411  if ((dxda = (double *) calloc(nwave, sizeof (double))) == NULL) {
412  printf("-E- %s line %d : error allocating memory for dxda in gsm:gsm_lm_f.\n",
413  __FILE__, __LINE__);
414  exit(1);
415  }
416  if ((dxdb = (double *) calloc(nwave, sizeof (double))) == NULL) {
417  printf("-E- %s line %d : error allocating memory for dxdb in gsm:gsm_lm_f.\n",
418  __FILE__, __LINE__);
419  exit(1);
420  }
421  first = 0;
422  }
423  if (coefset == CHESAPEAKE) {
424  set_cb_model(lambda, nwave, pixlon, pixlat, pixday);
425  }
426 
427  merit = 0;
428  for (iw = 0; iw < nwave; iw++) {
429  ac = aw [iw] +
430  aphstar[iw] * gsl_vector_get(par, 0) + gsl_vector_get(par, 1) * adgstar[iw];
431  bb = bbw[iw] + gsl_vector_get(par, 2) * bbpstar[iw];
432  x = bb / (ac + bb);
433  F = grd1 * x + grd2 * pow(x, 2);
434 
435  dFdx[iw] = grd1 + 2 * grd2*x;
436  dxda[iw] = -x * x / bb;
437  dxdb[iw] = x / bb + dxda[iw];
438 
439  gsl_vector_set(f, iw, (F - y[iw]) / sigma[iw]);
440 
441  merit += pow((F - y[iw]) / sigma[iw], 2);
442  }
443 
444  return GSL_SUCCESS;
445 }
446 
447 
448 /* ------------------------------------------------------------------- */
449 /* gsm_lm_df() GSM01 model, set-up to work with L.M. optimization */
450 /* */
451 /* Input: */
452 /* par[] - model parameters */
453 /* data - data to fit */
454 /* J - Jacobian */
455 /* */
456 
457 /* ------------------------------------------------------------------- */
458 int gsm_lm_df(const gsl_vector *par, void *data, gsl_matrix *J) {
459  int iw;
460  size_t nwave = ((struct datastruct *) data)->n;
461  double *sigma = ((struct datastruct *) data)->sigma;
462 
463  if (coefset == CHESAPEAKE) {
464  set_cb_model(lambda, nwave, pixlon, pixlat, pixday);
465  }
466 
467  for (iw = 0; iw < nwave; iw++) {
468 
469  /* Jacobian matrix J(iw,j) = dfi / dxj, */
470  gsl_matrix_set(J, iw, 0, dFdx[iw] * dxda[iw] * aphstar[iw] / sigma[iw]);
471  gsl_matrix_set(J, iw, 1, dFdx[iw] * dxda[iw] * adgstar[iw] / sigma[iw]);
472  gsl_matrix_set(J, iw, 2, dFdx[iw] * dxdb[iw] * bbpstar[iw] / sigma[iw]);
473  }
474 
475  return GSL_SUCCESS;
476 }
477 
478 int gsm_lm_fdf(const gsl_vector *par, void *data,
479  gsl_vector *f, gsl_matrix *J) {
480  gsm_lm_f(par, data, f);
481  gsm_lm_df(par, data, J);
482 
483  return GSL_SUCCESS;
484 }
485 
486 
487 /* ---------------------------------------------------------------------- */
488 /* run_gsm_lm() - returns requested GSM product for full scanline */
489 
490 /* ---------------------------------------------------------------------- */
491 void run_gsm_lm(l2str *l2rec) {
492  static int firstCall = TRUE;
493 
494  int32_t ip, iw;
495  double *Rrs;
496  double *wts;
497  double *sigma;
498  int16 itercnt;
499  int16 bndcnt;
500  int status;
501  l1str *l1rec = l2rec->l1rec;
502 
503  const gsl_multifit_fdfsolver_type *T;
504  gsl_multifit_fdfsolver *s;
505 
506  static gsl_multifit_function_fdf f;
507  double x_init[3] = {0.0, 0.0, 0.0};
508  gsl_vector_view x = gsl_vector_view_array(x_init, NPAR);
509  if ((Rrs = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
510  printf("-E- %s line %d : error allocating memory for Rrs in gsm:run_gsm_amb.\n",
511  __FILE__, __LINE__);
512  exit(1);
513  }
514  if ((wts = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
515  printf("-E- %s line %d : error allocating memory for wts in gsm:run_gsm_amb.\n",
516  __FILE__, __LINE__);
517  exit(1);
518  }
519  if ((sigma = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
520  printf("-E- %s line %d : error allocating memory for sigma in gsm:run_gsm_amb.\n",
521  __FILE__, __LINE__);
522  exit(1);
523  }
524 
525  if (firstCall == TRUE) {
526 
527  firstCall = FALSE;
528 
529  /* Set up data structure */
530  data.n = nbandsvis;
531  data.y = Rrs;
532  data.sigma = sigma;
533 
534  /* Set up multifit function structure */
535  f.f = &gsm_lm_f;
536  f.df = &gsm_lm_df;
537  f.fdf = &gsm_lm_fdf;
538  f.n = nbandsvis;
539  f.p = NPAR;
540  f.params = &data;
541  }
542 
543  T = gsl_multifit_fdfsolver_lmsder;
544  s = gsl_multifit_fdfsolver_alloc(T, nbandsvis, NPAR);
545 
546  /* if this is a new scanline, fit model for every pixel of the scanline */
547 
548  for (ip = 0; ip < l1rec->npix; ip++) {
549 
550  status = 1;
551  bndcnt = 0;
552  chl [ip] = badval;
553  adg [ip] = badval;
554  bbp [ip] = badval;
555  iter[ip] = 0;
556 
557  if (l1rec->mask[ip] == 0) {
558 
559  /* for CB model selection */
560  pixlon = l1rec->lon[ip];
561  pixlat = l1rec->lat[ip];
562  int16_t year, day;
563  double sec;
564  unix2yds(l1rec->scantime, &year, &day, &sec);
565  pixday = (int32_t) day;
566 
567 
568  for (iw = 0; iw < nbandsvis; iw++) {
569  wts[iw] = 1.0;
570  sigma[iw] = 1 / sqrt(wts[iw]);
571  Rrs[iw] = l2rec->Rrs[ip * l1rec->l1file->nbands + iw];
572  Rrs[iw] = Rrs[iw] / (0.52 + 1.7 * Rrs[iw]);
573  if (Rrs[iw] <= 0.0) status = 0;
574  else bndcnt++;
575  }
576 
577  /* if less than 3 valid radiances, fail pixel */
578  if (bndcnt < 3) {
579  l1rec->flags[ip] |= PRODFAIL;
580  continue;
581  }
582 
583  /* if any negative reflectance, set warning flag */
584  /* but allow processing to continue */
585  if (status == 0) {
586  l1rec->flags[ip] |= PRODWARN;
587  }
588 
589  /* fit model for this pixel */
590  gsl_multifit_fdfsolver_set(s, &f, &x.vector);
591  itercnt = 0;
592  do {
593  itercnt++;
594  status = gsl_multifit_fdfsolver_iterate(s);
595 
596  //printf ("status = %s\n", gsl_strerror (status));
597 
598  if (status)
599  break;
600 
601  status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
602 
603  } while (status == GSL_CONTINUE && itercnt < MAXITR);
604 
605 
606  /* save results and flag failure */
607  if (itercnt == MAXITR || gsl_vector_get(s->x, 0) <= 0.0) {
608  l1rec->flags[ip] |= PRODFAIL;
609  } else {
610  chl[ip] = gsl_vector_get(s->x, 0);
611  adg[ip] = gsl_vector_get(s->x, 1);
612  bbp[ip] = gsl_vector_get(s->x, 2);
613  }
614  iter[ip] = itercnt;
615  } else {
616  l1rec->flags[ip] |= PRODFAIL;
617  }
618  }
619 
620  GSMRecNum = l1rec->iscan;
621  gsl_multifit_fdfsolver_free(s);
622 
623  free(Rrs);
624  free(wts);
625  free(sigma);
626  return;
627 }
628 
629 
630 /* ---------------------------------------------------------------------- */
631 /* run_gsm - runs optimization using requested method */
632 
633 /* ---------------------------------------------------------------------- */
634 void run_gsm(l2str *l2rec) {
635  static int firstCall = TRUE;
636  l1str *l1rec = l2rec->l1rec;
637 
638  if (firstCall == TRUE) {
639 
640  int32_t iw;
641  int32_t taphn = -1;
642  float *taphw = NULL;
643  float *taphs = NULL;
644  float adg_s = -1;
645  float bbp_s = -1;
646 
647  int16_t year, day;
648  double sec;
649  unix2yds(l2rec->l1rec->scantime, &year, &day, &sec);
650  pixday = day;
651 
652  firstCall = FALSE;
653 
654  /* set model coefficient set */
655  coefset = input->gsm_opt;
656 
657  /* build array of sensor wavelengths */
658  for (iw = 0; iw < l1rec->l1file->nbands; iw++) {
659  lambda[iw] = l1rec->l1file->fwave[iw];
660  }
661 
662  /* limit to visible wavelengths */
663  nbandsvis = rdsensorinfo(l1rec->l1file->sensorID, input->evalmask, "NbandsVIS", NULL);
664 
665  /* set static coefficients */
666  if (coefset == GSMDEFAULT) {
667  adg_s = input->gsm_adg_s;
668  bbp_s = input->gsm_bbp_s;
669  taphw = input->gsm_aphw; /* aphstar table wavelengths */
670  taphs = input->gsm_aphs; /* aphstar table coefficients */
671  taphn = 0; /* aphstar table size */
672  for (taphn = 0; taphn < nbandsvis; taphn++)
673  if (taphw[taphn] < 0) break;
674  for (iw = 0; iw < nbandsvis; iw++) {
675  aphstar[iw] = linterp(taphw, taphs, taphn, lambda[iw]);
676  adgstar[iw] = exp(-adg_s * (lambda[iw] - 443.0));
677  bbpstar[iw] = pow((443.0 / lambda[iw]), bbp_s);
678  }
679  }
680 
681  /* save associated irradiances and water absorption/backscatter */
682  /* use band-averaged vales if no nLw out-of-band correction applied */
683 
684  if (input->outband_opt >= 2) {
685  for (iw = 0; iw < nbandsvis; iw++) {
686  aw [iw] = aw_spectra(lambda[iw], BANDW);
687  bbw[iw] = bbw_spectra(lambda[iw], BANDW);
688  }
689  } else {
690  for (iw = 0; iw < nbandsvis; iw++) {
691  aw [iw] = l1rec->l1file->aw[iw];
692  bbw[iw] = l1rec->l1file->bbw[iw];
693  }
694  }
695  }
696 
697  switch (input->gsm_fit) {
698  case AMOEBA:
699  run_gsm_amb(l2rec);
700  break;
701  case LEVMARQ:
702  run_gsm_lm(l2rec);
703  break;
704  default:
705  printf("%s Line %d: Unknown optimization method for GSM %d\n",
706  __FILE__, __LINE__, input->gsm_fit);
707  exit(1);
708  break;
709  }
710 }
711 
712 
713 /* ------------------------------------------------------------------- */
714 /* get_gsm() - returns requested GSM product for full scanline */
715 
716 /* ------------------------------------------------------------------- */
717 void get_gsm(l2str *l2rec, l2prodstr *p, float prod[]) {
718  int prodID = p->cat_ix;
719  int band = p->prod_ix;
720  int32_t ip, iw;
721  static int firstCall = 1;
722  static int last_npix;
723 
724  l1str *l1rec = l2rec->l1rec;
725  filehandle *l1file = l1rec->l1file;
726  int32_t nbands = l1file->nbands;
727 
728  /* return requested product, return -1 for abs/scatter in NIR */
729  if (firstCall) {
730 
731  if ((aphstar = (float *) calloc(nbands, sizeof (float))) == NULL) {
732  printf("-E- %s line %d : error allocating memory for gsm.\n",
733  __FILE__, __LINE__);
734  exit(1);
735  }
736  if ((adgstar = (float *) calloc(nbands, sizeof (double))) == NULL) {
737  printf("-E- %s line %d : error allocating memory for gsm.\n",
738  __FILE__, __LINE__);
739  exit(1);
740  }
741  if ((bbpstar = (float *) calloc(nbands, sizeof (double))) == NULL) {
742  printf("-E- %s line %d : error allocating memory for gsm.\n",
743  __FILE__, __LINE__);
744  exit(1);
745  }
746  if ((lambda = (float *) calloc(nbands, sizeof (float))) == NULL) {
747  printf("-E- %s line %d : error allocating memory for gsm.\n",
748  __FILE__, __LINE__);
749  exit(1);
750  }
751  if ((aw = (float *) calloc(nbands, sizeof (float))) == NULL) {
752  printf("-E- %s line %d : error allocating memory for gsm.\n",
753  __FILE__, __LINE__);
754  exit(1);
755  }
756  if ((bbw = (float *) calloc(nbands, sizeof (float))) == NULL) {
757  printf("-E- %s line %d : error allocating memory for gsm.\n",
758  __FILE__, __LINE__);
759  exit(1);
760  }
761  if ((chl = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
762  printf("-E- %s line %d : error allocating memory for gsm.\n",
763  __FILE__, __LINE__);
764  exit(1);
765  }
766  if ((adg = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
767  printf("-E- %s line %d : error allocating memory for gsm.\n",
768  __FILE__, __LINE__);
769  exit(1);
770  }
771  if ((bbp = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
772  printf("-E- %s line %d : error allocating memory for gsm.\n",
773  __FILE__, __LINE__);
774  exit(1);
775  }
776  if ((iter = (int16 *) calloc(l1rec->npix, sizeof (int16))) == NULL) {
777  printf("-E- %s line %d : error allocating memory for gsm.\n",
778  __FILE__, __LINE__);
779  exit(1);
780  }
781  last_npix = l1rec->npix;
782  firstCall = 0;
783  }
784 
785  if (band >= 0) {
786  if (lambda[band] > 700.0) {
787  for (ip = 0; ip < l1rec->npix; ip++)
788  prod[ip] = -1.0;
789  return;
790  }
791  }
792 
793  iw = p->prod_ix;
794 
795  if (l1rec->npix > last_npix) {
796  free(chl);
797  if ((chl = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
798  printf("-E- %s line %d : error allocating memory for gsm.\n",
799  __FILE__, __LINE__);
800  exit(1);
801  }
802  free(adg);
803  if ((adg = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
804  printf("-E- %s line %d : error allocating memory for gsm.\n",
805  __FILE__, __LINE__);
806  exit(1);
807  }
808  free(bbp);
809  if ((bbp = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
810  printf("-E- %s line %d : error allocating memory for gsm.\n",
811  __FILE__, __LINE__);
812  exit(1);
813  }
814  free(iter);
815  if ((iter = (int16 *) calloc(l1rec->npix, sizeof (int16))) == NULL) {
816  printf("-E- %s line %d : error allocating memory for gsm.\n",
817  __FILE__, __LINE__);
818  exit(1);
819  }
820  last_npix = l1rec->npix;
821  }
822 
823  if (!gsm_ran(l1rec->iscan))
824  run_gsm(l2rec);
825 
826  for (ip = 0; ip < l1rec->npix; ip++) {
827 
828  switch (prodID) {
829 
830  case CAT_chl_gsm:
831  prod[ip] = (float) chl[ip];
832  break;
833 
834  case CAT_adg_gsm:
835  if (band < 0)
836  prod[ip] = (float) adg[ip];
837  else
838  prod[ip] = (float) adg[ip] * adgstar[iw];
839  break;
840 
841  case CAT_bbp_gsm:
842  if (band < 0)
843  prod[ip] = (float) bbp[ip];
844  else
845  prod[ip] = (float) bbp[ip] * bbpstar[iw];
846  break;
847 
848  case CAT_aph_gsm:
849  prod[ip] = (float) aphstar[iw] * chl[ip];
850  break;
851 
852  case CAT_a_gsm:
853  prod[ip] = (float) aw[iw] + aphstar[iw] * chl[ip] + adg[ip] * adgstar[iw];
854  break;
855 
856  case CAT_bb_gsm:
857  prod[ip] = (float) bbw[iw] + bbp[ip] * bbpstar[iw];
858  break;
859 
860  default:
861  printf("-E- %s line %d : erroneous product ID %d passed to gsm.\n",
862  __FILE__, __LINE__, prodID);
863  exit(1);
864  }
865  }
866 
867  return;
868 }
869 
870 
871 /* ------------------------------------------------------------------- */
872 /* get_gsm_iter() - returns iteration count */
873 
874 /* ------------------------------------------------------------------- */
875 int16 *get_iter_gsm(l2str *l2rec) {
876  coefset = input->gsm_opt;
877  if (!gsm_ran(l2rec->l1rec->iscan))
878  run_gsm(l2rec);
879 
880  return iter;
881 }
882 
883 /* Interface to convl12() to return GSM iops */
884 
885 void iops_gsm(l2str *l2rec) {
886  int32_t ib, ip, ipb;
887  int32_t nbands = l2rec->l1rec->l1file->nbands;
888 
889  coefset = input->gsm_opt;
890  if (!gsm_ran(l2rec->l1rec->iscan))
891  run_gsm(l2rec);
892 
893  for (ip = 0; ip < l2rec->l1rec->npix; ip++)
894  for (ib = 0; ib < nbands; ib++) {
895  ipb = ip * nbands + ib;
896  l2rec->a [ipb] = (float) aw[ib] + aphstar[ib] * chl[ip] + adg[ip] * adgstar[ib];
897  l2rec->bb[ipb] = (float) bbw[ib] + bbp[ip] * bbpstar[ib];
898  }
899 
900  return;
901 }
902 
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
#define CAT_chl_gsm
Definition: l2prod.h:66
void interp(double *ephemPtr, int startLoc, double *inTime, int numCoefs, int numCom, int numSets, int velFlag, double *posvel)
integer, parameter int16
Definition: cubeio.f90:3
int j
Definition: decode_rs.h:73
double * y
Definition: gsm.c:68
int32_t day
int status
Definition: l1_czcs_hdf.c:31
#define CAT_aph_gsm
Definition: l2prod.h:98
#define MAXITR
Definition: gsm.c:29
int niter
Definition: amoeba.h:9
#define FALSE
Definition: rice.h:164
#define NULL
Definition: decode_rs.h:63
read l1rec
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
#define PRODWARN
Definition: l2_flags.h:11
#define TRUE
Definition: rice.h:165
float * lat
void run_gsm_lm(l2str *l2rec)
Definition: gsm.c:491
void get_gsm(l2str *l2rec, l2prodstr *p, float prod[])
Definition: gsm.c:717
#define PRODFAIL
Definition: l2_flags.h:39
data_t lambda[NROOTS+1]
Definition: decode_rs.h:75
#define CAT_bb_gsm
Definition: l2prod.h:97
short amoeba(double *pnts, FITSTRUCT *auxdata, double(*func)(FITSTRUCT *, double[]), float tol)
Definition: amoeba.c:14
double merit
Definition: amoeba.h:4
#define BANDW
Definition: l12_parms.h:84
instr * input
const double F
#define BAD_FLT
Definition: jplaeriallib.h:19
struct datastruct data
int16 * get_iter_gsm(l2str *l2rec)
Definition: gsm.c:875
void run_gsm_amb(l2str *l2rec)
Definition: gsm.c:290
void set_cb_model(float wave[], int32_t nwave, float lon, float lat, int32_t day)
Definition: gsm.c:84
int init(int32_t ipr, int32_t jpr, char *efile, char *pfile)
Definition: proj_report.c:51
int J
Definition: Usds.c:62
int fit_gsm_amb(double Rrs[], double wts[], int32_t npts, double fitparms[], double Rrs_fit[], int16 *itercnt)
Definition: gsm.c:243
double precision function f(R1)
Definition: tmd.lp.f:1454
#define NPAR
Definition: gsm.c:28
read recnum
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
int32_t rdsensorinfo(int32_t sensorID, int32_t evalmask, const char *pname, void **pval)
Definition: rdsensorinfo.c:69
double * wgt
Definition: amoeba.h:6
float aw_spectra(int32_t wl, int32_t width)
float bbw_spectra(int32_t wl, int32_t width)
float tol
#define GSMDEFAULT
Definition: gsm.c:31
void unix2yds(double usec, short *year, short *day, double *secs)
float * ac[MAX_BANDS]
double * y
Definition: amoeba.h:5
void run_gsm(l2str *l2rec)
Definition: gsm.c:634
int gsm_ran(int recnum)
Definition: gsm.c:72
double * sigma
Definition: gsm.c:69
#define AMOEBA
Definition: gsm.c:34
short int nfunc
Definition: amoeba.h:2
#define CAT_adg_gsm
Definition: l2prod.h:67
void iops_gsm(l2str *l2rec)
Definition: gsm.c:885
double * yfit
Definition: amoeba.h:8
size_t n
Definition: gsm.c:67
int32_t nbands
Definition: gsm.c:66
float * lon
short int npnts
Definition: amoeba.h:3
data_t s[NROOTS]
Definition: decode_rs.h:75
int gsm_lm_f(const gsl_vector *par, void *data, gsl_vector *f)
Definition: gsm.c:395
int32_t model
Definition: atrem_corl1.h:132
#define CAT_a_gsm
Definition: l2prod.h:96
int gsm_lm_fdf(const gsl_vector *par, void *data, gsl_vector *f, gsl_matrix *J)
Definition: gsm.c:478
double gsm_amb(FITSTRUCT *auxdata, double par[])
Definition: gsm.c:210
int i
Definition: decode_rs.h:71
#define CAT_bbp_gsm
Definition: l2prod.h:68
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:71
int gsm_lm_df(const gsl_vector *par, void *data, gsl_matrix *J)
Definition: gsm.c:458
#define CHESAPEAKE
Definition: gsm.c:32
#define LEVMARQ
Definition: gsm.c:35
float p[MODELMAX]
Definition: atrem_corl1.h:131