NASA Logo
Ocean Color Science Software

ocssw V2022
aph.c
Go to the documentation of this file.
1 /* ============================================================================================== */
2 /* aph.c - functions to compute phytoplankton-specific absorption */
3 /* ============================================================================================== */
4 
5 #include "l12_proto.h"
6 #include "giop.h"
7 
8 /* ---------------------------------------------------------------------------------------------- */
9 /* aph_bricaud() - compute aph for input wavelength and chl using Bricaud */
10 
11 /* ---------------------------------------------------------------------------------------------- */
12 float aph_bricaud_1995(float wave, float chl) {
13  static int firstCall = 1;
14  static float *wtab;
15  static float *atab, *datab;
16  static float *btab, *dbtab;
17  static int ntab = 0;
18 
19  float a, b, aph;
20 
21  if (firstCall) {
22 
23  FILE *fp;
24  char *filedir;
25  char filename[FILENAME_MAX];
26  char line [80];
27  int32_t itab;
28  float a, b, w;
29 
30  if ((filedir = getenv("OCDATAROOT")) == NULL) {
31  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
32  exit(1);
33  }
34  strcpy(filename, filedir);
35  strcat(filename, "/common/aph_bricaud_1995.txt");
36 
37  if ((fp = fopen(filename, "r")) == NULL) {
38  fprintf(stderr,
39  "-E- %s line %d: unable to open %s for reading\n", __FILE__, __LINE__, filename);
40  exit(1);
41  }
42 
43  printf("Reading aph* from %s.\n", filename);
44 
45  // number of lines
46  ntab = 0;
47  while (fgets(line, 80, fp))
48  ntab++;
49  rewind(fp);
50 
51  // allocate space
52  if ((wtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
53  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
54  __FILE__, __LINE__);
55  exit(1);
56  }
57  if ((atab = (float *) calloc(ntab, sizeof (float))) == NULL) {
58  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
59  __FILE__, __LINE__);
60  exit(1);
61  }
62  if ((btab = (float *) calloc(ntab, sizeof (float))) == NULL) {
63  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
64  __FILE__, __LINE__);
65  exit(1);
66  }
67  if ((datab = (float *) calloc(ntab, sizeof (float))) == NULL) {
68  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
69  __FILE__, __LINE__);
70  exit(1);
71  }
72  if ((dbtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
73  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
74  __FILE__, __LINE__);
75  exit(1);
76  }
77 
78  // read into arrays, skipping header or comment lines
79  itab = 0;
80  while (itab < ntab) {
81  if (fgets(line, 80, fp) == NULL) {
82  fprintf(stderr, "-E- %s line %d: error reading %s at line %d\n", __FILE__, __LINE__, filename, itab);
83  exit(1);
84  }
85  if (line[0] != '/' && line[0] != '!') {
86  sscanf(line, "%f,%f,%f", &w, &a, &b);
87  wtab [itab] = w;
88  atab [itab] = a;
89  btab [itab] = b;
90  itab++;
91  }
92  }
93  ntab = itab;
94 
95  // precompute derivatives for spline interpolation
96  spline(wtab, atab, ntab, 1e30, 1e30, datab);
97  spline(wtab, btab, ntab, 1e30, 1e30, dbtab);
98 
99  firstCall = 0;
100  }
101 
102  if (wave > wtab[ntab - 1]) wave = wtab[ntab - 1];
103 
104  // interpolate coefficients to wavelength
105  splint(wtab, atab, datab, ntab, wave, &a);
106  splint(wtab, btab, dbtab, ntab, wave, &b);
107  aph = a * pow(chl, -b);
108 
109  return (aph);
110 }
111 
112 /* ---------------------------------------------------------------------------------------------- */
113 /* aph_bricaud() - compute aph for input wavelength and chl using Bricaud */
114 
115 /* ---------------------------------------------------------------------------------------------- */
116 float aph_bricaud_1998(float wave, float chl) {
117  static int firstCall = 1;
118  static float *table;
119  static float *wtab;
120  static float *aphitab, *daphitab;
121  static float *ephitab, *dephitab;
122  static int ntab = 0;
123 
124  float aph, aphi, ephi;
125 
126  if (firstCall) {
127 
128  char *filedir;
129  char filename[FILENAME_MAX];
130  int ncol, nrow;
131 
132  if ((filedir = getenv("OCDATAROOT")) == NULL) {
133  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
134  exit(1);
135  }
136  strcpy(filename, filedir);
137  strcat(filename, "/common/aph_bricaud_1998.txt");
138  printf("Reading aph* from %s.\n", filename);
139 
140  nrow = table_row_count(filename);
141  if (nrow <= 0) {
142  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, filename);
143  exit(1);
144  }
145 
147  table = table_read_r4(filename, ncol, nrow);
148  ntab = nrow;
149 
150  wtab = &table[nrow * 0];
151  aphitab = &table[nrow * 3];
152  ephitab = &table[nrow * 4];
153 
154  // precompute derivatives for spline interpolation
155  if ((daphitab = (float *) calloc(ntab, sizeof (float))) == NULL) {
156  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
157  __FILE__, __LINE__);
158  exit(1);
159  }
160  if ((dephitab = (float *) calloc(ntab, sizeof (float))) == NULL) {
161  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
162  __FILE__, __LINE__);
163  exit(1);
164  }
165  spline(wtab, aphitab, ntab, 1e30, 1e30, daphitab);
166  spline(wtab, ephitab, ntab, 1e30, 1e30, dephitab);
167 
168  firstCall = 0;
169  }
170 
171  if (wave > wtab[ntab - 1]) wave = wtab[ntab - 1];
172 
173  // interpolate coefficients to wavelength
174  splint(wtab, aphitab, daphitab, ntab, wave, &aphi);
175  splint(wtab, ephitab, dephitab, ntab, wave, &ephi);
176  aph = aphi * pow(chl, ephi - 1);
177 
178  return (aph);
179 }
180 
181 float aph_bricaud(float wave, float chl) {
182  return (aph_bricaud_1998(wave, chl));
183 }
184 
185 /* ---------------------------------------------------------------------------------------------- */
186 /* aph_bricaud() - compute aph for input wavelength and chl using Bricaud */
187 
188 /* ---------------------------------------------------------------------------------------------- */
189 float pderiv_aph_bricaud_1998(float wave, float chl) {
190  static int firstCall = 1;
191  static float *table;
192  static float *wtab;
193  static float *aphitab, *daphitab;
194  static float *ephitab, *dephitab;
195  static int ntab = 0;
196 
197  float aphi, ephi, daphdchl;
198 
199  if (firstCall) {
200 
201  char *filedir;
202  char filename[FILENAME_MAX];
203  int ncol, nrow;
204 
205  if ((filedir = getenv("OCDATAROOT")) == NULL) {
206  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
207  exit(1);
208  }
209  strcpy(filename, filedir);
210  strcat(filename, "/common/aph_bricaud_1998.txt");
211  //printf("Reading aph* from %s.\n", filename);
212 
213  nrow = table_row_count(filename);
214  if (nrow <= 0) {
215  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, filename);
216  exit(1);
217  }
218 
220  table = table_read_r4(filename, ncol, nrow);
221  ntab = nrow;
222 
223  wtab = &table[nrow * 0];
224  aphitab = &table[nrow * 3];
225  ephitab = &table[nrow * 4];
226 
227  // precompute derivatives for spline interpolation
228  if ((daphitab = (float *) calloc(ntab, sizeof (float))) == NULL) {
229  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
230  __FILE__, __LINE__);
231  exit(1);
232  }
233  if ((dephitab = (float *) calloc(ntab, sizeof (float))) == NULL) {
234  printf("-E- %s line %d : error allocating memory for Bricaud aph.\n",
235  __FILE__, __LINE__);
236  exit(1);
237  }
238  spline(wtab, aphitab, ntab, 1e30, 1e30, daphitab);
239  spline(wtab, ephitab, ntab, 1e30, 1e30, dephitab);
240 
241  firstCall = 0;
242  }
243 
244  if (wave > wtab[ntab - 1]) wave = wtab[ntab - 1];
245 
246  // interpolate coefficients to wavelength
247  splint(wtab, aphitab, daphitab, ntab, wave, &aphi);
248  splint(wtab, ephitab, dephitab, ntab, wave, &ephi);
249  //aph = aphi * pow(chl, ephi - 1);
250 
251  daphdchl = aphi*(ephi-1.)* pow(chl,ephi-2.);
252 
253  return (daphdchl);
254 }
255 
256 float aph_bricaud_pderiv(float wave, float chl) {
257  return (pderiv_aph_bricaud_1998(wave, chl));
258 }
259 
260 
261 
262 /* ---------------------------------------------------------------------------------------------- */
263 /* aph_ciotti() - compute aph for input wavelength and size fraction using Cioti */
264 
265 /* ---------------------------------------------------------------------------------------------- */
266 float aph_ciotti(float wave, float sf) {
267  static int firstCall = 1;
268  static float *wtab;
269  static float *ptab, *dptab; // pico
270  static float *mtab, *dmtab; // micro
271  static int ntab = 0;
272 
273  float pico, micro, aph;
274 
275  if (firstCall) {
276 
277  FILE *fp;
278  char *filedir;
279  char filename[FILENAME_MAX];
280  char line [80];
281  int32_t itab;
282  float wv;
283 
284  if ((filedir = getenv("OCDATAROOT")) == NULL) {
285  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
286  exit(1);
287  }
288  strcpy(filename, filedir);
289  strcat(filename, "/common/aph_ciotti_2002_2006.txt");
290 
291  if ((fp = fopen(filename, "r")) == NULL) {
292  fprintf(stderr,
293  "-E- %s line %d: unable to open %s for reading\n", __FILE__, __LINE__, filename);
294  exit(1);
295  }
296 
297  printf("Reading aph* from %s.\n", filename);
298 
299  // number of lines
300  ntab = 0;
301  while (fgets(line, 80, fp))
302  ntab++;
303  rewind(fp);
304 
305  // allocate space
306  if ((wtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
307  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
308  __FILE__, __LINE__);
309  exit(1);
310  }
311  if ((ptab = (float *) calloc(ntab, sizeof (float))) == NULL) {
312  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
313  __FILE__, __LINE__);
314  exit(1);
315  }
316  if ((mtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
317  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
318  __FILE__, __LINE__);
319  exit(1);
320  }
321  if ((dptab = (float *) calloc(ntab, sizeof (float))) == NULL) {
322  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
323  __FILE__, __LINE__);
324  exit(1);
325  }
326  if ((dmtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
327  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
328  __FILE__, __LINE__);
329  exit(1);
330  }
331 
332  // read into arrays, skipping header or comment lines
333  itab = 0;
334  while (itab < ntab) {
335  if (fgets(line, 80, fp) == NULL) {
336  fprintf(stderr, "-E- %s line %d: error reading %s at line %d\n", __FILE__, __LINE__, filename, itab);
337  exit(1);
338  }
339  if (line[0] != '/' && line[0] != '!') {
340  sscanf(line, "%f %f %f", &wv, &pico, &micro);
341  wtab [itab] = wv;
342  ptab [itab] = pico * 0.0230 / 0.891;
343  mtab [itab] = micro * 0.0086 / 1.249;
344  itab++;
345  }
346  }
347  ntab = itab;
348 
349  // precompute derivatives for spline interpolation
350  spline(wtab, ptab, ntab, 1e30, 1e30, dptab);
351  spline(wtab, mtab, ntab, 1e30, 1e30, dmtab);
352 
353  firstCall = 0;
354  }
355 
356  if (wave > wtab[ntab - 1]) wave = wtab[ntab - 1];
357 
358  // interpolate coefficients to wavelength
359  splint(wtab, ptab, dptab, ntab, wave, &pico);
360  splint(wtab, mtab, dmtab, ntab, wave, &micro);
361 
362  aph = sf * pico + (1.0 - sf) * micro;
363 
364  return (aph);
365 }
366 
367 
368 /* ---------------------------------------------------------------------------------------------- */
369 /* get_aphstar() - compute aph* for center wavelength, wavelength width, function type, and proxy */
370 
371 /* ---------------------------------------------------------------------------------------------- */
372 float get_aphstar(float wave, int dwave, int ftype, float proxy) {
373  int32_t npts;
374  float wave1, wave2, wv;
375  float aph = 0.0;
376 
377  if (dwave > 1) {
378  npts = dwave + 1;
379  wave1 = wave - dwave / 2;
380  wave2 = wave + dwave / 2;
381  } else {
382  npts = 1;
383  wave1 = wave;
384  wave2 = wave;
385  }
386 
387  for (wv = wave1; wv <= wave2; wv += 1.0) {
388  switch (ftype) {
389  case APHBRICAUD:
390  aph += aph_bricaud(wv, proxy);
391  break;
392  case APHCIOTTI:
393  aph += aph_ciotti(wv, proxy);
394  break;
395  }
396  }
397 
398  aph /= npts;
399 
400  return (aph);
401 }
402 
403 /* ---------------------------------------------------------------------------------------------- */
404 /* get_uaphstar() - compute aph partial derivative wrt chl for center wavelength, */
405 /* wavelength width type, and proxy */
406 
407 /* ---------------------------------------------------------------------------------------------- */
408 float get_aphstar_pderiv(float wave, int dwave, int ftype, float proxy) {
409  int32_t npts;
410  float wave1, wave2, wv;
411  float daphdchl = 0.0;
412 
413  if (dwave > 1) {
414  npts = dwave + 1;
415  wave1 = wave - dwave / 2;
416  wave2 = wave + dwave / 2;
417  } else {
418  npts = 1;
419  wave1 = wave;
420  wave2 = wave;
421  }
422 
423  for (wv = wave1; wv <= wave2; wv += 1.0) {
424  switch (ftype) {
425  case APHBRICAUD:
426  daphdchl += aph_bricaud_pderiv(wv, proxy);
427  break;
428  case APHCIOTTI:
429  daphdchl += 0; //set to zero at present
430  break;
431  }
432  }
433 
434  daphdchl /= npts;
435 
436  return (daphdchl);
437 }
float * table_read_r4(char *input_filename, int m, int n)
Definition: table_io.cpp:2540
float get_aphstar(float wave, int dwave, int ftype, float proxy)
Definition: aph.c:372
int table_column_count(char *input_filename)
Definition: table_io.cpp:2524
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
#define NULL
Definition: decode_rs.h:63
float aph_bricaud_1998(float wave, float chl)
Definition: aph.c:116
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
int table_row_count(char *input_filename)
Definition: table_io.cpp:2528
float aph_bricaud_pderiv(float wave, float chl)
Definition: aph.c:256
float pderiv_aph_bricaud_1998(float wave, float chl)
Definition: aph.c:189
float get_aphstar_pderiv(float wave, int dwave, int ftype, float proxy)
Definition: aph.c:408
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
ftype
Definition: mapgen.py:234
data_t b[NROOTS+1]
Definition: decode_rs.h:77
#define APHCIOTTI
Definition: giop.h:37
Definition: common.h:10
float aph_ciotti(float wave, float sf)
Definition: aph.c:266
#define APHBRICAUD
Definition: giop.h:36
subroutine splint(xa, ya, y2a, n, x, y)
float aph_bricaud_1995(float wave, float chl)
Definition: aph.c:12
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
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
float aph_bricaud(float wave, float chl)
Definition: aph.c:181