OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
rayleigh.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 
3 #define NSOL 45
4 #define NSEN 41
5 #define NORDER 3
6 #define NWIND 8
7 
8 static int firstCall = 1;
9 
10 static float ray_sen [NSEN];
11 static float ray_sol [NSOL];
12 static float ray_sigma[NWIND];
13 typedef float ray_array[NWIND][NSOL][NORDER][NSEN];
14 static ray_array *ray_for_i;
15 static ray_array *ray_for_q;
16 static ray_array *ray_for_u;
17 
18 
19 /* ----------------------------------------------------------------- */
20 /* ray_press_wang() - Wang (2004) pressure correction */
21 
22 /* ----------------------------------------------------------------- */
23 float ray_press_wang(float taur, float airmass, float pr) {
24  static float p0 = STDPR;
25  float x = (-(0.6543 - 1.608 * taur) + (0.8192 - 1.2541 * taur) * log(airmass)) * taur*airmass;
26  return ( (1.0 - exp(-x * pr / p0)) / (1.0 - exp(-x)));
27 }
28 
29 
30 static void check_dimension_size(const char* file_name, int nc_id, const char* dim_name, size_t length) {
31  int dim_id;
32  size_t dim_length;
33  int status;
34 
35  status = nc_inq_dimid (nc_id, dim_name, &dim_id);
36  if (status != NC_NOERR) {
37  printf("-E- %s: Error looking for dimension %s from %s.\n", __FILE__, dim_name, file_name);
38  exit(1);
39  }
40  status = nc_inq_dimlen(nc_id, dim_id, &dim_length);
41  if (status != NC_NOERR) {
42  printf("-E- %s: Error reading dimension %s from %s.\n", __FILE__, dim_name, file_name);
43  exit(1);
44  }
45  if(dim_length != length) {
46  printf("-E- %s: Error with size of dimension %s from %s.\n", __FILE__, dim_name, file_name);
47  exit(1);
48  }
49 }
50 
51 
52 static void read_lut_variable(const char* file, int nc_id, const char* var_name, float* data) {
53  int var_id;
54  int status;
55 
56  status = nc_inq_varid(nc_id, var_name, &var_id);
57  if (status != NC_NOERR) {
58  printf("-E- %s: Error looking for variable %s from %s.\n", __FILE__, var_name, file);
59  exit(1);
60  }
61  status = nc_get_var_float(nc_id, var_id, data);
62  if (status != NC_NOERR) {
63  printf("-E- %s: Error reading variable %s from %s.\n", __FILE__, var_name, file);
64  exit(1);
65  }
66 }
67 
68 
69 void read_rayleigh_lut(char* file, int32_t iw, int pol_opt) {
70  int nc_id;
71  int status;
72 
73  /* Open the file and initiate the SD interface */
74  status = nc_open(file, NC_NOWRITE, &nc_id);
75  if (status != NC_NOERR) {
76  printf("-E- %s: Error opening file %s.\n", __FILE__, file);
77  exit(1);
78  } else
79  printf("Loading Rayleigh LUT %s\n", file);
80 
81  // check to make sure the dimensions are correct, so we don't overwrite memory
82  check_dimension_size(file, nc_id, "nrad_ray", NSEN);
83  check_dimension_size(file, nc_id, "nsun_ray", NSOL);
84  check_dimension_size(file, nc_id, "nwind_ray", NWIND);
85  check_dimension_size(file, nc_id, "norder_ray", NORDER);
86 
87  // read the variables
88  read_lut_variable(file, nc_id, "senz", ray_sen);
89  read_lut_variable(file, nc_id, "solz", ray_sol);
90  read_lut_variable(file, nc_id, "sigma", ray_sigma);
91  read_lut_variable(file, nc_id, "i_ray", (float*)ray_for_i[iw]);
92  if (pol_opt > 0) {
93  read_lut_variable(file, nc_id, "q_ray", (float*)ray_for_q[iw]);
94  read_lut_variable(file, nc_id, "u_ray", (float*)ray_for_u[iw]);
95  }
96 
97  nc_close(nc_id);
98 }
99 
100 
101 /* -------------------------------------------------------------------------------
102  * rayleigh() - compute Rayleigh radiances with polarization, per band.
103  *
104  * Description:
105  *
106  * This computes the Rayleigh scattering radiance (including
107  * polarization) for a given solar, viewing, and relative azimuthal
108  * angles, as well as the wind speed for the ocean surface roughness
109  * effects. This is a modification of M. Wangs rayleigh_i subroutine,
110  * incorporating the code to read the formatted hdf files.
111  *
112  * Inputs:
113  *
114  * l1rec l1 data line
115  * ip pixel # to process
116  *
117  * Outputs:
118  *
119  * l1rec l1 record with rayleigh radiances
120  *
121  * Hacked by: B. Franz, December 2002.
122  * C Version: B. Franz, November 2004
123  * W. Robinson, SAIC 20 May 2015 prevent solar zonith table index from
124  * exceeding the table size
125  * W. Robinson, SAIC 7 Feb 2017 Adapt for use of band-dependent viewing
126  * geometry and deal with l1rec itself
127  *
128  * ------------------------------------------------------------------------------- */
129 
130 void rayleigh(l1str *l1rec, int32_t ip) {
131 
132  /* indices */
133  int m;
134  int iwind;
135  int isigma;
136  int isigma1;
137  int isigma2;
138  int isol1;
139  int isol2;
140  int isen1;
141  int isen2;
142 
143  char *tmp_str;
144  char file [FILENAME_MAX] = "";
145  char path [FILENAME_MAX] = "";
146  char wavestr[10] = "";
147  char sensorstr[32] = "";
148 
149  /* working variables */
150  float sigma_m;
151  float p, q, h;
152  float f00, f10, f01, f11;
153  float cosd_phi[NORDER];
154  float sind_phi[NORDER];
155  float ray_i_sig[2];
156  float ray_q_sig[2];
157  float ray_u_sig[2];
158  float ray_i, ray_q, ray_u;
159  float airmass, fac;
160 
161  float *r_solz;
162  float *r_senz;
163  float *r_phi;
164  float *r_csolz;
165  float *r_csenz;
166 
167  int32_t iw, ipw, gmult, ix;
168 
169  int32_t sensorID = l1rec->l1file->sensorID;
170  int32_t evalmask = l1_input->evalmask;
171  int32_t nwave = l1rec->l1file->nbands;
172  int pol_opt = input->pol_opt;
173  float *taur = l1rec->l1file->Tau_r;
174  float *Fo = l1rec->Fo;
175  float pr = l1rec->pr[ip];
176  float ws = l1rec->ws[ip];
177  /*
178  * Note that if band-dependent viewing geometry is enabled for this
179  * instrument/run (geom_per_band != NULL), use that value in computing
180  * the Rayleigh radiance. Also, set the index to the angle for that
181  * band - ix, based on whether the angles are band-dependent
182  */
183  ipw = ip * nwave;
184  if (l1rec->geom_per_band != NULL) {
185  r_solz = &l1rec->geom_per_band->solz[ipw];
186  r_senz = &l1rec->geom_per_band->senz[ipw];
187  r_phi = &l1rec->geom_per_band->delphi[ipw];
188  r_csolz = &l1rec->geom_per_band->csolz[ipw];
189  r_csenz = &l1rec->geom_per_band->csenz[ipw];
190  } else {
191  r_solz = &l1rec->solz[ip];
192  r_senz = &l1rec->senz[ip];
193  r_phi = &l1rec->delphi[ip];
194  r_csolz = &l1rec->csolz[ip];
195  r_csenz = &l1rec->csenz[ip];
196  }
197 
198  setvbuf(stdout, NULL, _IOLBF, 0);
199  setvbuf(stderr, NULL, _IOLBF, 0);
200 
201  if (firstCall) {
202 
203  int32_t *wave;
204 
205  nwave = rdsensorinfo(sensorID, evalmask, "Lambda", (void **) &wave);
206 
207  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
208  printf("OCDATAROOT environment variable is not defined.\n");
209  exit(1);
210  }
211  if ((evalmask & NEWRAYTAB) != 0) {
212  strcpy(path, tmp_str);
213  strcat(path, "/eval/");
214  strcat(path, sensorId2SensorDir(sensorID));
215  strcat(path, "/");
216  } else {
217  strcpy(path, tmp_str);
218  strcat(path, "/");
219  strcat(path, sensorId2SensorDir(sensorID));
220  strcat(path, "/");
221  }
222  printf("\n");
223  if ((ray_for_i = (ray_array *) calloc(nwave, sizeof (ray_array))) == NULL) {
224  printf("-E- : Error allocating memory to ray_for_i\n");
225  exit(FATAL_ERROR);
226  }
227 
228  if ((ray_for_q = (ray_array *) calloc(nwave, sizeof (ray_array))) == NULL) {
229  printf("-E- : Error allocating memory to ray_for_q\n");
230  exit(FATAL_ERROR);
231  }
232  if ((ray_for_u = (ray_array *) calloc(nwave, sizeof (ray_array))) == NULL) {
233  printf("-E- : Error allocating memory to ray_for_u\n");
234  exit(FATAL_ERROR);
235  }
236 
237  strcpy(sensorstr, sensorId2SensorName(sensorID));
238  lowcase(sensorstr);
239 
240  for (iw = 0; iw < nwave; iw++) {
241 
242  sprintf(wavestr, "%d", (int) wave[iw]);
243  file[0] = 0;
244 
245  // try subsensor dir netCDF first
246  if(l1rec->l1file->subsensorID >= 0) {
247  strcpy(file, path);
248  strcat(file, subsensorId2SubsensorDir(l1rec->l1file->subsensorID));
249  strcat(file, "/rayleigh/rayleigh_");
250  strcat(file, sensorstr);
251  strcat(file, "_");
252  strcat(file, wavestr);
253  strcat(file, "_iqu.nc");
254 
255  // then try subsensor dir HDF
256  if(access(file, R_OK) == -1) {
257  strcpy(file, path);
258  strcat(file, subsensorId2SubsensorDir(l1rec->l1file->subsensorID));
259  strcat(file, "/rayleigh/rayleigh_");
260  strcat(file, sensorstr);
261  strcat(file, "_");
262  strcat(file, wavestr);
263  strcat(file, "_iqu.hdf");
264  }
265  }
266 
267  // then try sensor dir netCDF
268  if(access(file, R_OK) == -1) {
269  strcpy(file, path);
270  strcat(file, "rayleigh/rayleigh_");
271  strcat(file, sensorstr);
272  strcat(file, "_");
273  strcat(file, wavestr);
274  strcat(file, "_iqu.nc");
275 
276  // then try sensor dir HDF
277  if(access(file, R_OK) == -1) {
278  strcpy(file, path);
279  strcat(file, "rayleigh/rayleigh_");
280  strcat(file, sensorstr);
281  strcat(file, "_");
282  strcat(file, wavestr);
283  strcat(file, "_iqu.hdf");
284  }
285  }
286  read_rayleigh_lut(file, iw, pol_opt);
287  }
288  firstCall = 0;
289  }
290 
291 
292  /* windspeed index into tables */
293 
294  sigma_m = 0.0731 * sqrt(ws);
295  isigma2 = 1;
296  while (sigma_m > ray_sigma[isigma2] && isigma2 < NWIND)
297  isigma2++;
298  isigma1 = isigma2 - 1;
299 
300 
301 
302  /* interpolate Rayleigh Stokes vectors for each band */
303  gmult = 0;
304  if (l1rec->geom_per_band != NULL) {
305  gmult = 1;
306  }
307  for (iw = 0; iw < nwave; iw++) {
308 
309  ray_i = 0.0;
310  ray_q = 0.0;
311  ray_u = 0.0;
312 
313  /* geometry indices into tables */
314 
315  if ((iw == 0) || (gmult == 1)) {
316  ix = iw * gmult;
317 
318  for (m = 0; m < NORDER; m++) {
319  cosd_phi[m] = cos(r_phi[ix] * m / RADEG);
320  sind_phi[m] = sin(r_phi[ix] * m / RADEG);
321  }
322 
323  /* solar zenith indices */
324 
325  isol1 = ((int) r_solz[ix]) / 2;
326  isol2 = isol1 + 1;
327 
328  /* sensor zenith indices */
329 
330  for (isen2 = 0; isen2 < NSEN; isen2++) {
331  if (r_senz[ix] < ray_sen[isen2])
332  break;
333  }
334  isen1 = isen2 - 1;
335 
336  /* interpolation coefficients */
337  /* for solz > 88 or isol1 >= nsol - 1, use coeffs at end index */
338  if (isol1 >= NSOL - 1) {
339  isol1 = NSOL - 1;
340  isol2 = isol1;
341  p = 1.;
342  } else
343  p = (r_solz[ix] - ray_sol[isol1]) / (ray_sol[isol2] - ray_sol[isol1]);
344  q = (r_senz[ix] - ray_sen[isen1]) / (ray_sen[isen2] - ray_sen[isen1]);
345  airmass = 1. / r_csolz[ix] + 1. / r_csenz[ix];
346  }
347 
348 
349  /* interpolate the Rayleigh coefficients for each windspeed */
350 
351  for (isigma = isigma1; isigma <= isigma2; isigma++) {
352 
353  iwind = isigma - isigma1; /* 0 or 1 */
354 
355  ray_i_sig [iwind] = 0.;
356  ray_q_sig [iwind] = 0.;
357  ray_u_sig [iwind] = 0.;
358 
359  /* I component */
360 
361  for (m = 0; m < NORDER; m++) {
362 
363  f00 = ray_for_i[iw][isigma][isol1][m][isen1];
364  f10 = ray_for_i[iw][isigma][isol2][m][isen1];
365  f01 = ray_for_i[iw][isigma][isol1][m][isen2];
366  f11 = ray_for_i[iw][isigma][isol2][m][isen2];
367 
368  ray_i_sig[iwind] = ray_i_sig[iwind] +
369  ((1. - p)*(1. - q) * f00 + p * q * f11 + p * (1. - q) * f10 + q * (1. - p) * f01) * cosd_phi[m];
370  }
371 
372  if (pol_opt <= 0)
373  continue;
374 
375  /* Q component */
376 
377  for (m = 0; m < NORDER; m++) {
378 
379  f00 = ray_for_q[iw][isigma][isol1][m][isen1];
380  f10 = ray_for_q[iw][isigma][isol2][m][isen1];
381  f01 = ray_for_q[iw][isigma][isol1][m][isen2];
382  f11 = ray_for_q[iw][isigma][isol2][m][isen2];
383 
384  ray_q_sig[iwind] = ray_q_sig[iwind] +
385  ((1. - p)*(1. - q) * f00 + p * q * f11 + p * (1. - q) * f10 + q * (1. - p) * f01) * cosd_phi[m];
386  }
387 
388 
389  /* U component */
390 
391  for (m = 0; m < NORDER; m++) {
392 
393  f00 = ray_for_u[iw][isigma][isol1][m][isen1];
394  f10 = ray_for_u[iw][isigma][isol2][m][isen1];
395  f01 = ray_for_u[iw][isigma][isol1][m][isen2];
396  f11 = ray_for_u[iw][isigma][isol2][m][isen2];
397 
398  ray_u_sig[iwind] = ray_u_sig[iwind] +
399  ((1. - p)*(1. - q) * f00 + p * q * f11 + p * (1. - q) * f10 + q * (1. - p) * f01) * sind_phi[m];
400  }
401 
402  }
403 
404 
405  /* do linear interpolation between wind speeds */
406 
407  if (isigma1 == isigma2) {
408 
409  ray_i = ray_i_sig[0];
410 
411  if (pol_opt > 0) {
412  ray_q = ray_q_sig[0];
413  ray_u = ray_u_sig[0];
414  }
415 
416  } else {
417 
418  h = (sigma_m - ray_sigma[isigma1]) / (ray_sigma[isigma2] - ray_sigma[isigma1]);
419 
420  ray_i = ray_i_sig[0] + (ray_i_sig[1] - ray_i_sig[0]) * h;
421 
422  if (pol_opt > 0) {
423  ray_q = ray_q_sig[0] + (ray_q_sig[1] - ray_q_sig[0]) * h;
424  ray_u = ray_u_sig[0] + (ray_u_sig[1] - ray_u_sig[0]) * h;
425  }
426  }
427  /* Compute the rayleigh radiance */
428  fac = ray_press_wang(taur[iw], airmass, pr);
429  ipw = ip * nwave + iw;
430  l1rec->Lr[ipw] = Fo[iw] * ray_i*fac;
431  l1rec->L_q[ipw] = Fo[iw] * ray_q*fac;
432  l1rec->L_u[ipw] = Fo[iw] * ray_u*fac;
433  }
434 }
data_t q
Definition: decode_rs.h:74
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:240
#define NORDER
Definition: rayleigh.c:5
int status
Definition: l1_czcs_hdf.c:32
char * lowcase(char *instr)
Definition: lowcase.c:10
#define NWIND
Definition: rayleigh.c:6
#define NULL
Definition: decode_rs.h:63
read l1rec
float h[MODELMAX]
Definition: atrem_corl1.h:131
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed file
Definition: HISTORY.txt:413
instr * input
#define fac
string path
Definition: color_dtdb.py:221
l1_input_t * l1_input
Definition: l1_options.c:9
#define FATAL_ERROR
Definition: swl0_parms.h:5
#define NSEN
Definition: rayleigh.c:4
#define RADEG
Definition: czcs_ctl_pt.c:5
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
void read_rayleigh_lut(char *file, int32_t iw, int pol_opt)
Definition: rayleigh.c:69
#define STDPR
Definition: l12_parms.h:85
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:198
void rayleigh(l1str *l1rec, int32_t ip)
Definition: rayleigh.c:130
#define NSOL
Definition: rayleigh.c:3
#define NEWRAYTAB
Definition: l1.h:76
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
const char * subsensorId2SubsensorDir(int subsensorId)
Definition: sensorInfo.c:254
float ray_array[NWIND][NSOL][NORDER][NSEN]
Definition: rayleigh.c:13
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:97
float p[MODELMAX]
Definition: atrem_corl1.h:131
float ray_press_wang(float taur, float airmass, float pr)
Definition: rayleigh.c:23