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
get_cal_swf.c
Go to the documentation of this file.
1 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 
3  Function: get_cal_swf
4 
5  Reads and applies calibration table
6 
7  Return:
8  calibrated L1B radiances
9 
10  Sean Bailey, Futuretech Corporation, 2 Jun 2008
11 
12  Modification history:
13  Gene Eplee, SAIC, 17 June 2010 Changed fitting function type
14  descriptor from sds attribute to
15  sds fitting parameter variable
16 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 
18 #include <stdlib.h>
19 #include "hdf4utils.h"
20 #include <timeutils.h>
21 #include "cal_l1a.h"
22 #include "call1a_proto.h"
23 #include "genutils.h"
24 #include "getcal_proto.h"
25 
26 #include <hdf.h>
27 #include <mfhdf.h>
28 
29 
30 float ***alloc3d_float(int m, int n, int o) {
31  int i, j;
32  float ***p;
33  p = malloc(m * sizeof (float *));
34 
35  for (i = 0; i < m; i++) {
36  p[ i ] = malloc(n * sizeof (float **));
37  for (j = 0; j < n; j++) {
38  p[ i ][ j ] = malloc(o * sizeof (float ***));
39  }
40  }
41 
42  return (p);
43 }
44 
45 float ****alloc4d_float(int l, int m, int n, int o) {
46  int i, j, k;
47  float ****p;
48  p = malloc(l * sizeof (float *));
49 
50  for (i = 0; i < l; i++) {
51  p[ i ] = malloc(m * sizeof (float **));
52  for (j = 0; j < m; j++) {
53  p[ i ][ j ] = malloc(n * sizeof (float ***));
54  for (k = 0; k < n; k++) {
55  p[ i ][ j ][ k ] = malloc(o * sizeof (float ****));
56 
57  }
58 
59  }
60  }
61 
62  return (p);
63 }
64 
65 int32_t ***alloc3d_long(int m, int n, int o) {
66  int i, j;
67  int32_t ***p;
68  p = malloc(m * sizeof (int32_t *));
69 
70  for (i = 0; i < m; i++) {
71  p[ i ] = malloc(n * sizeof (int32_t **));
72  for (j = 0; j < n; j++) {
73  p[ i ][ j ] = malloc(o * sizeof (int32_t ***));
74  }
75  }
76 
77  return (p);
78 }
79 
80 
81 #define NBAND 8
82 #define NTEMP 256
83 // #define NCOEF 5
84 #define NCOEF 6
85 #define NDET 4
86 #define NGAIN 4
87 #define NTDI 256
88 #define NKNEE 5
89 #define K1 0
90 #define K3 1
91 #define K4 2
92 #define MSIDE 3
93 #define DARK 5
94 #define GR 4
95 #define NK 4
96 
97 static short ftype[NK];
98 
99 static short gain3corr;
100 static float ***radcor;
101 static int32_t *k1_epoch;
102 // static int32_t *k1_ftype;
103 static int num_k1_epoch;
104 
105 static float ***g3corr;
106 static int32_t *gr_epoch;
107 // static int32_t *gr_ftype;
108 static int num_gr_epoch;
109 
110 static float ***cnts2rad;
111 static int32_t ***det_off;
112 
113 static float ***fptempcor;
114 static int32_t *k3_epoch;
115 // static int32_t *k3_ftype;
116 static int num_k3_epoch;
117 
118 static float ***scanmod;
119 static int32_t *k4_epoch;
120 // static int32_t *k4_ftype;
121 static int num_k4_epoch;
122 
123 static float ****msidecor;
124 static int32_t *ms_epoch;
125 // static int32_t *ms_ftype;
126 static int num_ms_epoch;
127 
128 static float ***dark_restore;
129 static int32_t *dkrest_epoch;
130 static int num_dkrest_epoch;
131 static int haveDarkRestore = 0;
132 
133 static float fp_temp[NBAND][NTEMP];
134 static short tdi_list[NDET][NTDI];
135 static float g_f[NBAND][1024];
136 
137 static int32_t ref_year;
138 static int32_t ref_day;
139 static int32_t ref_min;
140 static short prev_tdi[NBAND] = {-1, -1, -1, -1, -1, -1, -1, -1};
141 static short prev_gain[NBAND] = {-1, -1, -1, -1, -1, -1, -1, -1};
142 static int32_t prev_syear;
143 static int32_t prev_sday;
144 static int32_t prev_smsec;
145 static double ref_jsec;
146 
147 /* -------------------------------------------------------------------------- */
148 /* read_caltable() - called once to load cal table static arrays */
149 
150 /* -------------------------------------------------------------------------- */
151 
152 void read_caltable(char *cal_path) {
153 
154  static int firstCall = 1;
155  char name [H4_MAX_NC_NAME] = "";
156  char sdsname[H4_MAX_NC_NAME] = "";
157 
158  int32 sd_id;
159  int32 sds_id;
160  int32 rank;
161  int32 nt;
162  int32 dims[H4_MAX_VAR_DIMS];
163  int32 nattrs, attr_index, count, num_type;
164  char attr_name[64];
165  int32 start4[4] = {0, 0, 0, 0};
166  int32 start[3] = {0, 0, 0};
167  int32 status;
168  float *r;
169  int32_t *dr;
170 
171  int32_t i, j, k, l;
172  int m = 1;
173 
174  if (!firstCall) return;
175 
176  /* Open the cal cal_path... */
177  sd_id = SDstart(cal_path, DFACC_RDONLY);
178  if (sd_id == FAIL) {
179  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
180  __FILE__, __LINE__, cal_path, DFACC_RDONLY);
181  exit(1);
182  }
183 
184  /* Functional type index */
185  strcpy(sdsname, "ftype");
186  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
187  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
188  status = SDreaddata(sds_id, start, NULL, dims, ftype);
189  if (status != 0) {
190  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
191  __FILE__, __LINE__, sdsname, cal_path);
192  exit(1);
193  } else {
194  status = SDendaccess(sds_id);
195  }
196 
197  /* Radiometric coefficients */
198  strcpy(sdsname, "radiometric_coef");
199  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
200  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
201  radcor = alloc3d_float(dims[0], dims[1], dims[2]);
202  m = 1;
203  for (i = 0; i < rank; i++)
204  m *= dims[i];
205  r = (float *) malloc(m * sizeof (float *));
206  /* Get epoch start times */
207  attr_index = SDfindattr(sds_id, "k1_epochs");
208  status = SDattrinfo(sds_id, attr_index, attr_name, &num_type, &count);
209  num_k1_epoch = count;
210  k1_epoch = (int32_t *) malloc(count * sizeof (int32_t *));
211  status = SDreadattr(sds_id, attr_index, k1_epoch);
212  /* Get functional type index*/
213  /* attr_index = SDfindattr(sds_id,"ftype index");
214  status = SDattrinfo(sds_id,attr_index,attr_name,&num_type, &count);
215  k1_ftype = (int32_t *) malloc(count * sizeof(int32_t *));
216  status = SDreadattr(sds_id,attr_index,k1_ftype); */
217 
218  /* Get the data...*/
219  status = SDreaddata(sds_id, start, NULL, dims, r);
220  if (status != 0) {
221  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
222  __FILE__, __LINE__, sdsname, cal_path);
223  exit(1);
224  } else {
225  status = SDendaccess(sds_id);
226  for (i = 0; i < dims[0]; i++) {
227  for (j = 0; j < dims[1]; j++) {
228  for (k = 0; k < dims[2]; k++) {
229  m = i * dims[1] * dims[2] + j * dims[2] + k;
230  radcor[i][j][k] = r[m];
231  }
232  }
233  }
234  }
235 
236  free(r);
237 
238  /* Get the dark counts */
239  strcpy(sdsname, "dark_counts");
240  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
241  if (sds_id >= 0) {
242  haveDarkRestore = 1;
243  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
244  dark_restore = alloc3d_float(dims[0], dims[1], dims[2]);
245  m = 1;
246  for (i = 0; i < rank; i++)
247  m *= dims[i];
248  r = (float *) malloc(m * sizeof (float *));
249  /* Get epoch start times */
250  attr_index = SDfindattr(sds_id, "dn0_epochs");
251  status = SDattrinfo(sds_id, attr_index, attr_name, &num_type, &count);
252  num_dkrest_epoch = count;
253  dkrest_epoch = (int32_t *) malloc(count * sizeof (int32_t *));
254  status = SDreadattr(sds_id, attr_index, dkrest_epoch);
255  status = SDreaddata(sds_id, start, NULL, dims, r);
256  if (status != 0) {
257  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
258  __FILE__, __LINE__, sdsname, cal_path);
259  exit(1);
260  } else {
261  status = SDendaccess(sds_id);
262  for (i = 0; i < dims[0]; i++) {
263  for (j = 0; j < dims[1]; j++) {
264  for (k = 0; k < dims[2]; k++) {
265  m = i * dims[1] * dims[2] + j * dims[2] + k;
266  dark_restore[i][j][k] = r[m];
267  }
268  }
269  }
270  }
271  free(r);
272 
273  }
274 
275 
276  /* Gain 3:1 correction coefficients */
277  strcpy(sdsname, "gainratio_coef");
278  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
279  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
280  if (status == 0) {
281  if (want_verbose)
282  printf("\nGain 3 correction applied.\n");
283  gain3corr = 1;
284  g3corr = alloc3d_float(dims[0], dims[1], dims[2]);
285  m = 1;
286  for (i = 0; i < rank; i++)
287  m *= dims[i];
288  r = (float *) malloc(m * sizeof (float *));
289  /* Get epoch start times */
290  attr_index = SDfindattr(sds_id, "gr_epochs");
291  status = SDattrinfo(sds_id, attr_index, attr_name, &num_type, &count);
292  num_gr_epoch = count;
293  gr_epoch = (int32_t *) malloc(count * sizeof (int32_t *));
294  status = SDreadattr(sds_id, attr_index, gr_epoch);
295  /* Get functional type index*/
296  /* attr_index = SDfindattr(sds_id,"ftype index");
297  status = SDattrinfo(sds_id,attr_index,attr_name,&num_type, &count);
298  gr_ftype = (int32_t *) malloc(count * sizeof(int32_t *));
299  status = SDreadattr(sds_id,attr_index,gr_ftype); */
300  /* Get the data...*/
301  status = SDreaddata(sds_id, start, NULL, dims, r);
302  if (status != 0) {
303  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
304  __FILE__, __LINE__, sdsname, cal_path);
305  exit(1);
306  } else {
307  status = SDendaccess(sds_id);
308  for (i = 0; i < dims[0]; i++) {
309  for (j = 0; j < dims[1]; j++) {
310  for (k = 0; k < dims[2]; k++) {
311  m = i * dims[1] * dims[2] + j * dims[2] + k;
312  g3corr[i][j][k] = r[m];
313  }
314  }
315  }
316  }
317 
318  free(r);
319  } else {
320  gain3corr = 0;
321  }
322 
323  /* Counts to Radiance coefficients */
324  strcpy(sdsname, "counts_to_radiance");
325  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
326  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
327  cnts2rad = alloc3d_float(dims[2], dims[1], dims[0]);
328  m = 1;
329  for (i = 0; i < rank; i++)
330  m *= dims[i];
331  r = (float *) malloc(m * sizeof (float *));
332  /* Get the data...*/
333  status = SDreaddata(sds_id, start, NULL, dims, r);
334  if (status != 0) {
335  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
336  __FILE__, __LINE__, sdsname, cal_path);
337  exit(1);
338  } else {
339  status = SDendaccess(sds_id);
340  for (i = 0; i < dims[0]; i++) {
341  for (j = 0; j < dims[1]; j++) {
342  for (k = 0; k < dims[2]; k++) {
343  m = i * dims[1] * dims[2] + j * dims[2] + k;
344  cnts2rad[k][j][i] = r[m];
345  }
346  }
347  }
348  }
349 
350  free(r);
351 
352 
353  /* Detector offsets */
354  strcpy(sdsname, "detector_offsets");
355  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
356  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
357  det_off = alloc3d_long(dims[2], dims[1], dims[0]);
358  m = 1;
359  for (i = 0; i < rank; i++)
360  m *= dims[i];
361  dr = (int32_t *) malloc(m * sizeof (int32_t *));
362  /* Get the data...*/
363  status = SDreaddata(sds_id, start, NULL, dims, dr);
364  if (status != 0) {
365  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
366  __FILE__, __LINE__, sdsname, cal_path);
367  exit(1);
368  } else {
369  status = SDendaccess(sds_id);
370  for (i = 0; i < dims[0]; i++) {
371  for (j = 0; j < dims[1]; j++) {
372  for (k = 0; k < dims[2]; k++) {
373  m = i * dims[1] * dims[2] + j * dims[2] + k;
374  det_off[k][j][i] = dr[m];
375  }
376  }
377  }
378  }
379 
380  free(dr);
381 
382  /* Temperature coefficients */
383  strcpy(sdsname, "temperature_coef");
384  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
385  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
386  fptempcor = alloc3d_float(dims[0], dims[1], dims[2]);
387  m = 1;
388  for (i = 0; i < rank; i++)
389  m *= dims[i];
390  r = (float *) malloc(m * sizeof (float *));
391  /* Get epoch start times */
392  attr_index = SDfindattr(sds_id, "k3_epochs");
393  status = SDattrinfo(sds_id, attr_index, attr_name, &num_type, &count);
394  num_k3_epoch = count;
395  k3_epoch = (int32_t *) malloc(count * sizeof (int32_t *));
396  status = SDreadattr(sds_id, attr_index, k3_epoch);
397  /* Get functional type index*/
398  /* attr_index = SDfindattr(sds_id,"ftype index");
399  status = SDattrinfo(sds_id,attr_index,attr_name,&num_type, &count);
400  k3_ftype = (int32_t *) malloc(count * sizeof(int32_t *));
401  status = SDreadattr(sds_id,attr_index,k3_ftype); */
402  /* Get the data...*/
403  status = SDreaddata(sds_id, start, NULL, dims, r);
404  if (status != 0) {
405  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
406  __FILE__, __LINE__, sdsname, cal_path);
407  exit(1);
408  } else {
409  status = SDendaccess(sds_id);
410  for (i = 0; i < dims[0]; i++) {
411  for (j = 0; j < dims[1]; j++) {
412  for (k = 0; k < dims[2]; k++) {
413  m = i * dims[1] * dims[2] + j * dims[2] + k;
414  fptempcor[i][j][k] = r[m];
415  }
416  }
417  }
418  }
419 
420  free(r);
421 
422  /* Scan Modulation correction*/
423  strcpy(sdsname, "scan_mod_coef");
424  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
425  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
426  scanmod = alloc3d_float(dims[0], dims[1], dims[2]);
427  m = 1;
428  for (i = 0; i < rank; i++)
429  m *= dims[i];
430  r = (float *) malloc(m * sizeof (float *));
431  /* Get epoch start times */
432  attr_index = SDfindattr(sds_id, "k4_epochs");
433  status = SDattrinfo(sds_id, attr_index, attr_name, &num_type, &count);
434  num_k4_epoch = count;
435  k4_epoch = (int32_t *) malloc(count * sizeof (int32_t *));
436  status = SDreadattr(sds_id, attr_index, k4_epoch);
437  /* Get functional type index*/
438  /* attr_index = SDfindattr(sds_id,"ftype index");
439  status = SDattrinfo(sds_id,attr_index,attr_name,&num_type, &count);
440  k4_ftype = (int32_t *) malloc(count * sizeof(int32_t *));
441  status = SDreadattr(sds_id,attr_index,k4_ftype); */
442  /* Get the data...*/
443  status = SDreaddata(sds_id, start, NULL, dims, r);
444  if (status != 0) {
445  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
446  __FILE__, __LINE__, sdsname, cal_path);
447  exit(1);
448  } else {
449  status = SDendaccess(sds_id);
450  for (i = 0; i < dims[0]; i++) {
451  for (j = 0; j < dims[1]; j++) {
452  for (k = 0; k < dims[2]; k++) {
453  m = i * dims[1] * dims[2] + j * dims[2] + k;
454  scanmod[i][j][k] = r[m];
455  }
456  }
457  }
458  }
459 
460  free(r);
461 
462  /* Mirror Side correction*/
463  strcpy(sdsname, "mirror_side_coef");
464  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
465  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
466  msidecor = alloc4d_float(dims[0], dims[1], dims[2], dims[3]);
467  m = 1;
468  for (i = 0; i < rank; i++)
469  m *= dims[i];
470  r = (float *) malloc(m * sizeof (float *));
471  /* Get epoch start times */
472  attr_index = SDfindattr(sds_id, "ms_epochs");
473  status = SDattrinfo(sds_id, attr_index, attr_name, &num_type, &count);
474  num_ms_epoch = count;
475  ms_epoch = (int32_t *) malloc(count * sizeof (int32_t *));
476  status = SDreadattr(sds_id, attr_index, ms_epoch);
477  /* Get functional type index*/
478  /* attr_index = SDfindattr(sds_id,"ftype index");
479  status = SDattrinfo(sds_id,attr_index,attr_name,&num_type, &count);
480  ms_ftype = (int32_t *) malloc(count * sizeof(int32_t *));
481  status = SDreadattr(sds_id,attr_index,ms_ftype); */
482  /* Get the data...*/
483  status = SDreaddata(sds_id, start4, NULL, dims, r);
484  if (status != 0) {
485  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
486  __FILE__, __LINE__, sdsname, cal_path);
487  exit(1);
488  } else {
489  status = SDendaccess(sds_id);
490  for (i = 0; i < dims[0]; i++) {
491  for (j = 0; j < dims[1]; j++) {
492  for (k = 0; k < dims[2]; k++) {
493  for (l = 0; l < dims[3]; l++) {
494  m = i * dims[1] * dims[2] * dims[3] + j * dims[2] * dims[3] + k * dims[3] + l;
495  msidecor[i][j][k][l] = r[m];
496  }
497  }
498  }
499  }
500  }
501 
502  free(r);
503 
504  /* Focal Plane Temperature Array*/
505  strcpy(sdsname, "fp_temp");
506  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
507  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
508  /* Get the data...*/
509  status = SDreaddata(sds_id, start, NULL, dims, fp_temp);
510  if (status != 0) {
511  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
512  __FILE__, __LINE__, sdsname, cal_path);
513  exit(1);
514  } else {
515  status = SDendaccess(sds_id);
516  }
517 
518  /* TDI List Array*/
519  strcpy(sdsname, "TDI_list");
520  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
521  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
522  /* Get the data...*/
523  status = SDreaddata(sds_id, start, NULL, dims, tdi_list);
524  if (status != 0) {
525  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
526  __FILE__, __LINE__, sdsname, cal_path);
527  exit(1);
528  } else {
529  status = SDendaccess(sds_id);
530  }
531 
532  /* Get Reference Date info */
533  attr_index = SDfindattr(sd_id, "Reference Year");
534  status = SDreadattr(sd_id, attr_index, &ref_year);
535  if (status != 0) {
536  printf("-E- %s Line %d: Error reading Reference Year from %s.\n",
537  __FILE__, __LINE__, cal_path);
538  exit(1);
539  }
540  attr_index = SDfindattr(sd_id, "Reference Day");
541  status = SDreadattr(sd_id, attr_index, &ref_day);
542  if (status != 0) {
543  printf("-E- %s Line %d: Error reading Reference Day from %s.\n",
544  __FILE__, __LINE__, cal_path);
545  exit(1);
546  }
547  attr_index = SDfindattr(sd_id, "Reference Minute");
548  status = SDreadattr(sd_id, attr_index, &ref_min);
549  if (status != 0) {
550  printf("-E- %s Line %d: Error reading Reference Minute from %s.\n",
551  __FILE__, __LINE__, cal_path);
552  exit(1);
553  }
554 
555  ref_jsec = yds2unix((int) ref_year, (int) ref_day, ((double) (ref_min)) / 1440.0);
556 
557  /* terminate access to the SD interface and close the cal_path */
558  status = SDend(sd_id);
559 
560  firstCall = 0;
561 }
562 
563 /* --------------------------------------------------------------------------*/
564 /* apply_time_coef - given the functional form index, coefficients */
565 /* and "x" parameter, returns "y" */
566 
567 /* --------------------------------------------------------------------------*/
568 float apply_time_coef(short functype, float coef[NCOEF], double param) {
569  float correction_factor;
570 
571  switch (functype) {
572  // default functype (0) - double exponential
573  case 0:
574  correction_factor = 1.0 / (coef[0] -
575  coef[1] * (1.0 - exp(-coef[2] * param)) -
576  coef[3] * (1.0 - exp(-coef[4] * param)));
577  break;
578  // functype (1) - quadratic
579  case 1:
580  correction_factor = 1.0 / (coef[0] + coef[1] * param + coef[2] * pow(param, 2));
581  break;
582  // functype (2) - exponential + linear
583  case 2:
584  correction_factor = 1.0 / (coef[0] -
585  coef[1] * (1.0 - exp(-coef[2] * param)) -
586  coef[3] * param);
587  break;
588  }
589 
590  return correction_factor;
591 }
592 
593 /* --------------------------------------------------------------------------*/
594 /* get_epoch_idx - returns the index for the epoch */
595 
596 /* --------------------------------------------------------------------------*/
597 int get_epoch_idx(int param, double jsec) {
598  int epidx = 0;
599  int32_t *epochs;
600  int cnt, i;
601 
602  switch (param) {
603  case K1:
604  epochs = k1_epoch;
605  cnt = num_k1_epoch;
606  break;
607  case GR:
608  epochs = gr_epoch;
609  cnt = num_gr_epoch;
610  break;
611  case K3:
612  epochs = k3_epoch;
613  cnt = num_k3_epoch;
614  break;
615  case K4:
616  epochs = k4_epoch;
617  cnt = num_k4_epoch;
618  break;
619  case MSIDE:
620  epochs = ms_epoch;
621  cnt = num_ms_epoch;
622  break;
623  case DARK:
624  epochs = dkrest_epoch;
625  cnt = num_dkrest_epoch;
626  break;
627  }
628  for (i = 0; i < cnt; i++) {
629  if (jsec > (double) epochs[i]) {
630  epidx = i;
631  }
632  }
633  return epidx;
634 }
635 
636 /* --------------------------------------------------------------------------*/
637 /* calc_knees_two - slightly modified version of the original... */
638 
639 /* --------------------------------------------------------------------------*/
640 void calc_knees_two(int16 *tdi, float32 counts[8][4][5], float32 rads[8][4][5]) {
641 
642  int16 dets[4];
643  int32 i, j, k;
644  int32 scnts[4]; /* saturation counts */
645  float32 srads[4]; /* saturation radiance */
646  float32 loc_slopes[4];
647  float32 slopes[NBAND][4][4];
648  int32 cnts[NBAND][4][4];
649  int32 oindex[NDET];
650 
651 
652  for (i = 0; i < NBAND; i++)
653  for (j = 0; j < 4; j++)
654  for (k = 0; k < 4; k++) {
655  slopes[i][j][k] = cnts2rad[i][j][k];
656  cnts[i][j][k] = det_off[i][j][k];
657  }
658 
659  for (i = 0; i < NBAND; i++) {
660  for (j = 0; j < 4; j++)
661  dets[j] = tdi_list[j][tdi[i]] - 1;
662  for (j = 0; j < NGAIN; j++) {
663  for (k = 0; k < NDET; k++) {
664  scnts[k] = 1023 - cnts[i][j][dets[k]];
665  srads[k] = scnts[k] * slopes[i][j][dets[k]];
666  loc_slopes[k] = slopes[i][j][dets[k]];
667  }
668 
669  sort_srads(srads, oindex);
670 
671  rads[i][j][0] = 0;
672  for (k = 1; k < 5; k++)
673  rads[i][j][k] = srads[oindex[k - 1]];
674 
675  counts[i][j][0] = 0;
676  counts[i][j][1] = (scnts[oindex[0]] +
677  srads[oindex[0]] / loc_slopes[oindex[1]] +
678  srads[oindex[0]] / loc_slopes[oindex[2]] +
679  srads[oindex[0]] / loc_slopes[oindex[3]]) / 4.0;
680 
681  counts[i][j][2] = (scnts[oindex[0]] + scnts[oindex[1]] +
682  srads[oindex[1]] / loc_slopes[oindex[2]] +
683  srads[oindex[1]] / loc_slopes[oindex[3]]) / 4.0;
684 
685 
686  counts[i][j][3] = (scnts[oindex[0]] + scnts[oindex[1]] +
687  scnts[oindex[2]] +
688  srads[oindex[2]] / loc_slopes[oindex[3]]) / 4.0;
689 
690  counts[i][j][4] = (scnts[oindex[0]] + scnts[oindex[1]] +
691  scnts[oindex[2]] + scnts[oindex[3]]) / 4.0;
692 
693  }
694  }
695 }
696 /* -------------------------------------------------------------------------- */
697 /* l1b_rad - returns calibrated L1B radiances */
698 
699 /* -------------------------------------------------------------------------- */
700 int32_t l1b_rad(int syear, int sday, int32_t smsec, int32_t msec,
701  char *dtype, int32_t nsta, int32_t ninc, int32_t npix,
702  float *dark_mean, short *gain, short *tdi,
703  short *scan_temp, float *inst_temp, int mside,
704  short *l1a_data, float *l1b_data, cal_mod_struc *cal_mod) {
705  int i, pixel, knee, n, count;
706  int band;
707  int l1_data_lo, l1_data_hi;
708  int dark;
709  float dark_ratio;
710  double delta_t;
711  int gn[NBAND];
712  float k1, gcorr, k3, k4, mirror;
713  int16 k1_ftype, gr_ftype, k4_ftype, ms_ftype;
714  int epoch_idx[NCOEF];
715  double jsec;
716  float coef[NCOEF];
717  float cal_offset = 0.;
718  int gain_flag, tdi_flag;
719  float slope;
720  short count1, count2;
721  short scanpix;
722  int called_calc_knees = 0;
723 
724  jsec = yds2unix(syear, sday, ((double) (msec)) / 1000.0);
725  delta_t = (jsec - (double) ref_jsec) / 86400.;
726 
727  // TDI, gain and scan_temp funny business...Gene E. got some 'spalin' to do ;)
728  for (band = 0; band < NBAND; band++) {
729  if (tdi[band] < 0)
730  tdi[band] = 0;
731  if (tdi[band] > 255)
732  tdi[band] = 255;
733  if (scan_temp[band] < 0)
734  scan_temp[band] = 0;
735  if (scan_temp[band] > 255)
736  scan_temp[band] = 255;
737  // Fix gain 2 / gain 3 telemetry bit flip
738  switch (gain[band]) {
739  case 0:
740  gn[band] = 0;
741  break;
742  case 1:
743  gn[band] = 2;
744  break;
745  case 2:
746  gn[band] = 1;
747  break;
748  case 3:
749  gn[band] = 3;
750  break;
751  }
752 
753  }
754 
755  for (band = 0; band < NBAND && tdi[band] == prev_tdi[band]; band++);
756  if (band < NBAND)
757  tdi_flag = 1;
758  else
759  tdi_flag = 0;
760 
761  // Calc knees ...
762 
763  if (syear != prev_syear || sday != prev_sday
764  || smsec != prev_smsec || tdi_flag) {
765  prev_syear = syear;
766  prev_sday = sday;
767  prev_smsec = smsec;
768 
769  for (band = 0; band < NBAND; band++)
770  prev_tdi[band] = tdi[band];
771 
772  calc_knees_two(tdi, cal_counts, cal_rads);
773  called_calc_knees = 1;
774  }
775 
776 
777  for (band = 0; band < NBAND && gn[band] == prev_gain[band]; band++);
778 
779  if (band < NBAND)
780  gain_flag = 1;
781  else
782  gain_flag = 0;
783  // ... and then generate the 'radiance' LUT...
784  if (called_calc_knees || gain_flag) {
785  called_calc_knees = 0;
786  for (band = 0; band < NBAND; band++)
787  prev_gain[band] = gn[band];
788  for (band = 0; band < NBAND; band++) {
789  for (knee = 1; knee <= 4; knee++) {
790  n = 1;
791  while (((int16) cal_counts[band][gn[band]][knee] ==
792  (int16) cal_counts[band][gn[band]][knee - n]) && n <= knee)
793  n++;
794  count1 = (int16) cal_counts[band][gn[band]][knee - n] + 1;
795  count2 = (int16) cal_counts[band][gn[band]][knee];
796  if (knee == 1)
797  count1 = 0;
798  if (knee == 4)
799  count2 = 1023;
800  slope = (cal_rads[band][gn[band]][knee] -
801  cal_rads[band][gn[band]][knee - n]) /
802  (cal_counts[band][gn[band]][knee] -
803  cal_counts[band][gn[band]][knee - n]);
804  for (count = count1; count <= count2; count++)
805  g_f[band][count] =
806  slope * (count - cal_counts[band][gn[band]][knee - n]) +
807  cal_rads[band][gn[band]][knee - n];
808  }
809  }
810  }
811 
812  // ...and now for the real heavy lifting...
813  epoch_idx[0] = get_epoch_idx(K1, jsec);
814  epoch_idx[1] = get_epoch_idx(K3, jsec);
815  epoch_idx[2] = get_epoch_idx(MSIDE, jsec);
816  epoch_idx[3] = get_epoch_idx(K4, jsec);
817  epoch_idx[4] = get_epoch_idx(GR, jsec);
818 
819  for (band = 0; band < NBAND; band++) {
820 
821  // Radiometric decay correction
822  for (i = 0; i < NCOEF; i++) {
823  coef[i] = radcor[epoch_idx[0]][i][band];
824  }
825  k1_ftype = (int16) coef[NCOEF - 1];
826  k1 = apply_time_coef(k1_ftype, coef, delta_t);
827 
828  gcorr = 1.;
829  if (gain3corr) {
830  for (i = 0; i < NCOEF; i++) {
831  coef[i] = g3corr[epoch_idx[4]][i][band];
832  }
833  // gain3corr is NOT applied in the inverse...so invert it :)
834  gr_ftype = (int16) coef[NCOEF - 1];
835  gcorr = 1. / apply_time_coef(gr_ftype, coef, delta_t);
836  }
837  /* Check for override on system gain and offset */
838  /* and pass back the gain, offset used */
839  if (cal_mod->flag == 1) {
840  k1 = cal_mod->gain[band] * k1;
841  } else if (cal_mod->flag == 2) {
842  cal_offset = cal_mod->offset[band];
843  }
844  else if (cal_mod->flag == 3) {
845  k1 = cal_mod->gain[band] * k1;
846  cal_offset = cal_mod->offset[band];
847  } else {
848  cal_mod->gain[band] = k1;
849  cal_mod->offset[band] = cal_offset;
850  }
851 
852  // Temperature correction
853  for (i = 0; i < NCOEF; i++)
854  coef[i] = fptempcor[epoch_idx[1]][i][band];
855  // k3_ftype = (int16)fptempcor[epoch_idx[1]][NCOEF][band];
856  // k3 = apply_time_coef(k3_ftype,coef,delta_t);
857  k3 = (1.0 + coef[0]*(fp_temp[band][scan_temp[band]] - coef[1]));
858 
859  // Mirror side correction
860  for (i = 0; i < NCOEF; i++)
861  coef[i] = msidecor[epoch_idx[2]][i][mside][band];
862  ms_ftype = (int16) coef[NCOEF - 1];
863  mirror = apply_time_coef(ms_ftype, coef, delta_t);
864  if (mside == 1) mirror = 1.0 / mirror;
865 
866  // Scan modulation (a.k.a. RVS) correction
867  for (i = 0; i < NCOEF; i++)
868  coef[i] = scanmod[epoch_idx[3]][i][band];
869  k4_ftype = (int16) coef[NCOEF - 1];
870  for (pixel = 0; pixel < npix; pixel++) {
871 
872  scanpix = nsta + ninc*pixel;
873  // apply scan modulation correction to earth view data only...
874  if ((strcmp(dtype, "SOL") != 0) && (strcmp(dtype, "TDI") != 0) &&
875  (strcmp(dtype, "IGC") != 0))
876  k4 = apply_time_coef(k4_ftype, coef, scanpix - 643);
877  else k4 = 1.0;
878 
879  // If the cal table has the dark_count array, use it.
880  // If not, use the mean of the scene
881  if (haveDarkRestore) {
882  int dkidx = get_epoch_idx(DARK, jsec);
883 
884  dark = ceil(dark_restore[dkidx][gn[band]][band]);
885  dark_ratio = (float) dark - dark_restore[dkidx][gn[band]][band];
886 
887  } else {
888  dark = ceil(dark_mean[band]);
889  dark_ratio = (float) dark - dark_mean[band];
890  }
891 
892  l1_data_lo = l1a_data[band * npix + pixel] - dark;
893 
894  if (l1_data_lo < 0)
895  l1_data_lo = 0;
896  if (l1_data_lo > 1022)
897  l1_data_lo = 1022;
898 
899  l1_data_hi = l1_data_lo + 1;
900 
901  l1b_data[band * npix + pixel] = (k1 * gcorr * k3 * k4 * mirror * (g_f[band][l1_data_lo]*(1 - dark_ratio) + g_f[band][l1_data_hi] * dark_ratio) + cal_offset);
902 
903  }
904 
905  }
906 
907  return SUCCEED;
908 
909 }
910 
NDET
#define NDET
Definition: get_cal_swf.c:85
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
rads
float rads[BANDS_DIMS_1A][GAINS_DIMS_1A][KNEES_DIMS_1A]
Definition: l1a_seawifs.c:47
r
int r
Definition: decode_rs.h:73
sort_srads
void sort_srads(float32 *srads, int32 *oindex)
Definition: calib_get_cal_misc.c:501
cal_rads
float cal_rads[BANDS_DIMS_1A][GAINS_DIMS_1A][KNEES_DIMS_1A]
Definition: calibrate_l1a.c:75
dark_mean
float dark_mean[8]
Definition: l1a_seawifs.c:34
j
int j
Definition: decode_rs.h:73
status
int status
Definition: l1_czcs_hdf.c:32
yds2unix
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
gain
int16 * gain
Definition: l1_czcs_hdf.c:33
FAIL
#define FAIL
Definition: ObpgReadGrid.h:18
read_caltable
void read_caltable(char *cal_path)
Definition: get_cal_swf.c:152
NULL
#define NULL
Definition: decode_rs.h:63
alloc3d_float
float *** alloc3d_float(int m, int n, int o)
Definition: get_cal_swf.c:30
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
alloc3d_long
int32_t *** alloc3d_long(int m, int n, int o)
Definition: get_cal_swf.c:65
l1a_data
int16_t * l1a_data
Definition: l1a_seawifs.c:84
msec
int32 * msec
Definition: l1_czcs_hdf.c:31
counts
uint8 * counts
Definition: l1_czcs_hdf.c:30
syear
int syear
Definition: l1_czcs_hdf.c:15
NCOEF
#define NCOEF
Definition: get_cal_swf.c:84
smsec
int32 smsec
Definition: l1_czcs_hdf.c:16
getcal_proto.h
K1
#define K1
Definition: get_cal_swf.c:89
call1a_proto.h
sday
int sday
Definition: l1_czcs_hdf.c:15
nsta
int32 nsta
Definition: l1_czcs_hdf.c:27
GR
#define GR
Definition: get_cal_swf.c:94
MSIDE
#define MSIDE
Definition: get_cal_swf.c:92
mside
int32_t mside
Definition: l1_hmodis_hdf.c:752
NTDI
#define NTDI
Definition: get_cal_swf.c:87
cal_l1a.h
slope
float32 slope[]
Definition: l2lists.h:30
alloc4d_float
float **** alloc4d_float(int l, int m, int n, int o)
Definition: get_cal_swf.c:45
ninc
int32 ninc
Definition: l1_czcs_hdf.c:28
NBAND
#define NBAND
Definition: get_cal_swf.c:81
NGAIN
#define NGAIN
Definition: get_cal_swf.c:86
tdi
int16_t tdi[BANDS_DIMS_1A]
Definition: l1a_seawifs.c:37
want_verbose
int want_verbose
Definition: genutils_globals.c:3
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
NK
#define NK
Definition: get_cal_swf.c:95
pixel
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution pixel
Definition: HISTORY.txt:192
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
ref_day
int16_t ref_day
Definition: l1a_seawifs.c:42
l1b_rad
int32_t l1b_rad(int syear, int sday, int32_t smsec, int32_t msec, char *dtype, int32_t nsta, int32_t ninc, int32_t npix, float *dark_mean, short *gain, short *tdi, short *scan_temp, float *inst_temp, int mside, short *l1a_data, float *l1b_data, cal_mod_struc *cal_mod)
Definition: get_cal_swf.c:700
dtype
dtype
Definition: DDataset.hpp:31
calc_knees_two
void calc_knees_two(int16 *tdi, float32 counts[8][4][5], float32 rads[8][4][5])
Definition: get_cal_swf.c:640
rank
Extra metadata that will be written to the HDF4 file l2prod rank
Definition: HOWTO_Add_a_product.txt:80
NTEMP
#define NTEMP
Definition: get_cal_swf.c:82
DARK
#define DARK
Definition: get_cal_swf.c:93
apply_time_coef
float apply_time_coef(short functype, float coef[NCOEF], double param)
Definition: get_cal_swf.c:568
K4
#define K4
Definition: get_cal_swf.c:91
get_epoch_idx
int get_epoch_idx(int param, double jsec)
Definition: get_cal_swf.c:597
ref_year
int16_t ref_year
Definition: l1a_seawifs.c:41
i
int i
Definition: decode_rs.h:71
K3
#define K3
Definition: get_cal_swf.c:90
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
cal_counts
float cal_counts[BANDS_DIMS_1A][GAINS_DIMS_1A][KNEES_DIMS_1A]
Definition: calibrate_l1a.c:74
k
int k
Definition: decode_rs.h:73
npix
int npix
Definition: get_cmp.c:28
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
count
int count
Definition: decode_rs.h:79


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