ocssw V2020
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 /* ---------------------------------------------------------------------------------------------- */
187 /* aph_ciotti() - compute aph for input wavelength and size fraction using Cioti */
188 
189 /* ---------------------------------------------------------------------------------------------- */
190 float aph_ciotti(float wave, float sf) {
191  static int firstCall = 1;
192  static float *wtab;
193  static float *ptab, *dptab; // pico
194  static float *mtab, *dmtab; // micro
195  static int ntab = 0;
196 
197  float pico, micro, aph;
198 
199  if (firstCall) {
200 
201  FILE *fp;
202  char *filedir;
203  char filename[FILENAME_MAX];
204  char line [80];
205  int32_t itab;
206  float wv;
207 
208  if ((filedir = getenv("OCDATAROOT")) == NULL) {
209  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
210  exit(1);
211  }
212  strcpy(filename, filedir);
213  strcat(filename, "/common/aph_ciotti_2002_2006.txt");
214 
215  if ((fp = fopen(filename, "r")) == NULL) {
216  fprintf(stderr,
217  "-E- %s line %d: unable to open %s for reading\n", __FILE__, __LINE__, filename);
218  exit(1);
219  }
220 
221  printf("Reading aph* from %s.\n", filename);
222 
223  // number of lines
224  ntab = 0;
225  while (fgets(line, 80, fp))
226  ntab++;
227  rewind(fp);
228 
229  // allocate space
230  if ((wtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
231  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
232  __FILE__, __LINE__);
233  exit(1);
234  }
235  if ((ptab = (float *) calloc(ntab, sizeof (float))) == NULL) {
236  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
237  __FILE__, __LINE__);
238  exit(1);
239  }
240  if ((mtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
241  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
242  __FILE__, __LINE__);
243  exit(1);
244  }
245  if ((dptab = (float *) calloc(ntab, sizeof (float))) == NULL) {
246  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
247  __FILE__, __LINE__);
248  exit(1);
249  }
250  if ((dmtab = (float *) calloc(ntab, sizeof (float))) == NULL) {
251  printf("-E- %s line %d : error allocating memory for Ciotti aph.\n",
252  __FILE__, __LINE__);
253  exit(1);
254  }
255 
256  // read into arrays, skipping header or comment lines
257  itab = 0;
258  while (itab < ntab) {
259  if (fgets(line, 80, fp) == NULL) {
260  fprintf(stderr, "-E- %s line %d: error reading %s at line %d\n", __FILE__, __LINE__, filename, itab);
261  exit(1);
262  }
263  if (line[0] != '/' && line[0] != '!') {
264  sscanf(line, "%f %f %f", &wv, &pico, &micro);
265  wtab [itab] = wv;
266  ptab [itab] = pico * 0.0230 / 0.891;
267  mtab [itab] = micro * 0.0086 / 1.249;
268  itab++;
269  }
270  }
271  ntab = itab;
272 
273  // precompute derivatives for spline interpolation
274  spline(wtab, ptab, ntab, 1e30, 1e30, dptab);
275  spline(wtab, mtab, ntab, 1e30, 1e30, dmtab);
276 
277  firstCall = 0;
278  }
279 
280  if (wave > wtab[ntab - 1]) wave = wtab[ntab - 1];
281 
282  // interpolate coefficients to wavelength
283  splint(wtab, ptab, dptab, ntab, wave, &pico);
284  splint(wtab, mtab, dmtab, ntab, wave, &micro);
285 
286  aph = sf * pico + (1.0 - sf) * micro;
287 
288  return (aph);
289 }
290 
291 
292 /* ---------------------------------------------------------------------------------------------- */
293 /* get_aphstar() - compute aph* for center wavelength, wavelength width, function type, and proxy */
294 
295 /* ---------------------------------------------------------------------------------------------- */
296 float get_aphstar(float wave, int dwave, int ftype, float proxy) {
297  int32_t npts;
298  float wave1, wave2, wv;
299  float aph = 0.0;
300 
301  if (dwave > 1) {
302  npts = dwave + 1;
303  wave1 = wave - dwave / 2;
304  wave2 = wave + dwave / 2;
305  } else {
306  npts = 1;
307  wave1 = wave;
308  wave2 = wave;
309  }
310 
311  for (wv = wave1; wv <= wave2; wv += 1.0) {
312  switch (ftype) {
313  case APHBRICAUD:
314  aph += aph_bricaud(wv, proxy);
315  break;
316  case APHCIOTTI:
317  aph += aph_ciotti(wv, proxy);
318  break;
319  }
320  }
321 
322  aph /= npts;
323 
324  return (aph);
325 }
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:296
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
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
data_t b[NROOTS+1]
Definition: decode_rs.h:77
#define APHCIOTTI
Definition: giop.h:35
Definition: common.h:9
float aph_ciotti(float wave, float sf)
Definition: aph.c:190
#define APHBRICAUD
Definition: giop.h:34
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