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
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  int var_id;
73  int use_netcdf=0;
74 
75  /* Open the file and initiate the SD interface */
76  status = nc_open(file, NC_NOWRITE, &nc_id);
77  if (status != NC_NOERR) {
78  printf("-E- %s: Error opening file %s.\n", __FILE__, file);
79  exit(1);
80  } else
81  printf("Loading Rayleigh LUT %s\n", file);
82 
83  // check to make sure the dimensions are correct, so we don't overwrite memory
84  if(strstr(file,"_iqu.nc"))
85  use_netcdf=1;
86 
87  if(use_netcdf){
88  check_dimension_size(file, nc_id, "sensor_zenith", NSEN);
89  check_dimension_size(file, nc_id, "solar_zenith", NSOL);
90  check_dimension_size(file, nc_id, "wave_mean_square_slope", NWIND);
91  check_dimension_size(file, nc_id, "fourier_component", NORDER);
92 
93  // read the variables
94  read_lut_variable(file, nc_id, "sensor_zenith", ray_sen);
95  read_lut_variable(file, nc_id, "solar_zenith", ray_sol);
96  }
97  else{
98  check_dimension_size(file, nc_id, "nrad_ray", NSEN);
99  check_dimension_size(file, nc_id, "nsun_ray", NSOL);
100  check_dimension_size(file, nc_id, "nwind_ray", NWIND);
101  check_dimension_size(file, nc_id, "norder_ray", NORDER);
102 
103  // read the variables
104  read_lut_variable(file, nc_id, "senz", ray_sen);
105  read_lut_variable(file, nc_id, "solz", ray_sol);
106  }
107 
108  status = nc_inq_varid(nc_id, "sigma", &var_id);
109  if(status==NC_NOERR)
110  read_lut_variable(file, nc_id, "sigma", ray_sigma);
111  else{
112  status = nc_inq_varid(nc_id, "wave_mean_square_slope", &var_id);
113 
114  if(status!=NC_NOERR){
115  printf("-E- %s: Error looking for variable sigma or wave_mean_square_slope from %s.\n", __FILE__, file);
116  exit(1);
117  }
118  read_lut_variable(file, nc_id, "wave_mean_square_slope", ray_sigma);
119  }
120 
121  if(use_netcdf){
122  read_lut_variable(file, nc_id, "stokes_i_rayleigh", (float*)ray_for_i[iw]);
123  if (pol_opt > 0) {
124  read_lut_variable(file, nc_id, "stokes_q_rayleigh", (float*)ray_for_q[iw]);
125  read_lut_variable(file, nc_id, "stokes_u_rayleigh", (float*)ray_for_u[iw]);
126  }
127  }
128  else{
129  read_lut_variable(file, nc_id, "i_ray", (float*)ray_for_i[iw]);
130  if (pol_opt > 0) {
131  read_lut_variable(file, nc_id, "q_ray", (float*)ray_for_q[iw]);
132  read_lut_variable(file, nc_id, "u_ray", (float*)ray_for_u[iw]);
133  }
134  }
135 
136  nc_close(nc_id);
137 }
138 
139 
140 /* -------------------------------------------------------------------------------
141  * rayleigh() - compute Rayleigh radiances with polarization, per band.
142  *
143  * Description:
144  *
145  * This computes the Rayleigh scattering radiance (including
146  * polarization) for a given solar, viewing, and relative azimuthal
147  * angles, as well as the wind speed for the ocean surface roughness
148  * effects. This is a modification of M. Wangs rayleigh_i subroutine,
149  * incorporating the code to read the formatted hdf files.
150  *
151  * Inputs:
152  *
153  * l1rec l1 data line
154  * ip pixel # to process
155  *
156  * Outputs:
157  *
158  * l1rec l1 record with rayleigh radiances
159  *
160  * Hacked by: B. Franz, December 2002.
161  * C Version: B. Franz, November 2004
162  * W. Robinson, SAIC 20 May 2015 prevent solar zonith table index from
163  * exceeding the table size
164  * W. Robinson, SAIC 7 Feb 2017 Adapt for use of band-dependent viewing
165  * geometry and deal with l1rec itself
166  *
167  * ------------------------------------------------------------------------------- */
168 
169 void rayleigh(l1str *l1rec, int32_t ip) {
170 
171  /* indices */
172  int m;
173  int iwind;
174  int isigma;
175  int isigma1;
176  int isigma2;
177  int isol1;
178  int isol2;
179  int isen1;
180  int isen2;
181 
182  char *tmp_str;
183  char file [FILENAME_MAX] = "";
184  char path [FILENAME_MAX] = "";
185  char wavestr[10] = "";
186  char sensorstr[32] = "";
187 
188  /* working variables */
189  float sigma_m;
190  float p, q, h;
191  float f00, f10, f01, f11;
192  float cosd_phi[NORDER];
193  float sind_phi[NORDER];
194  float ray_i_sig[2];
195  float ray_q_sig[2];
196  float ray_u_sig[2];
197  float ray_i, ray_q, ray_u;
198  float airmass, fac;
199 
200  float *r_solz;
201  float *r_senz;
202  float *r_phi;
203  float *r_csolz;
204  float *r_csenz;
205 
206  int32_t iw, ipw, gmult, ix;
207 
208  int32_t sensorID = l1rec->l1file->sensorID;
209  int32_t evalmask = l1_input->evalmask;
210  int32_t nwave = l1rec->l1file->nbands;
211  int pol_opt = input->pol_opt;
212  float *taur = l1rec->l1file->Tau_r;
213  float *Fo = l1rec->Fo;
214  float pr = l1rec->pr[ip];
215  float ws = l1rec->ws[ip];
216  /*
217  * Note that if band-dependent viewing geometry is enabled for this
218  * instrument/run (geom_per_band != NULL), use that value in computing
219  * the Rayleigh radiance. Also, set the index to the angle for that
220  * band - ix, based on whether the angles are band-dependent
221  */
222  ipw = ip * nwave;
223  if (l1rec->geom_per_band != NULL) {
224  r_solz = &l1rec->geom_per_band->solz[ipw];
225  r_senz = &l1rec->geom_per_band->senz[ipw];
226  r_phi = &l1rec->geom_per_band->delphi[ipw];
227  r_csolz = &l1rec->geom_per_band->csolz[ipw];
228  r_csenz = &l1rec->geom_per_band->csenz[ipw];
229  } else {
230  r_solz = &l1rec->solz[ip];
231  r_senz = &l1rec->senz[ip];
232  r_phi = &l1rec->delphi[ip];
233  r_csolz = &l1rec->csolz[ip];
234  r_csenz = &l1rec->csenz[ip];
235  }
236 
237  setvbuf(stdout, NULL, _IOLBF, 0);
238  setvbuf(stderr, NULL, _IOLBF, 0);
239 
240  if (firstCall) {
241 
242  int32_t *wave;
243 
244  nwave = rdsensorinfo(sensorID, evalmask, "Lambda", (void **) &wave);
245 
246  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
247  printf("OCDATAROOT environment variable is not defined.\n");
248  exit(1);
249  }
250  if ((evalmask & NEWRAYTAB) != 0) {
251  strcpy(path, tmp_str);
252  strcat(path, "/eval/");
253  strcat(path, sensorId2SensorDir(sensorID));
254  strcat(path, "/");
255  } else {
256  strcpy(path, tmp_str);
257  strcat(path, "/");
258  strcat(path, sensorId2SensorDir(sensorID));
259  strcat(path, "/");
260  }
261  printf("\n");
262  if ((ray_for_i = (ray_array *) calloc(nwave, sizeof (ray_array))) == NULL) {
263  printf("-E- : Error allocating memory to ray_for_i\n");
264  exit(FATAL_ERROR);
265  }
266 
267  if ((ray_for_q = (ray_array *) calloc(nwave, sizeof (ray_array))) == NULL) {
268  printf("-E- : Error allocating memory to ray_for_q\n");
269  exit(FATAL_ERROR);
270  }
271  if ((ray_for_u = (ray_array *) calloc(nwave, sizeof (ray_array))) == NULL) {
272  printf("-E- : Error allocating memory to ray_for_u\n");
273  exit(FATAL_ERROR);
274  }
275 
276  strcpy(sensorstr, sensorId2SensorName(sensorID));
277  lowcase(sensorstr);
278 
279  for (iw = 0; iw < nwave; iw++) {
280 
281  sprintf(wavestr, "%d", (int) wave[iw]);
282  file[0] = 0;
283 
284  // try subsensor dir netCDF first
285  if(l1rec->l1file->subsensorID >= 0) {
286  strcpy(file, path);
287  strcat(file, subsensorId2SubsensorDir(l1rec->l1file->subsensorID));
288  strcat(file, "/rayleigh/rayleigh_");
289  strcat(file, sensorstr);
290  strcat(file, "_");
291  strcat(file, wavestr);
292  strcat(file, "_iqu.nc");
293 
294  // then try subsensor dir HDF
295  if(access(file, R_OK) == -1) {
296  strcpy(file, path);
297  strcat(file, subsensorId2SubsensorDir(l1rec->l1file->subsensorID));
298  strcat(file, "/rayleigh/rayleigh_");
299  strcat(file, sensorstr);
300  strcat(file, "_");
301  strcat(file, wavestr);
302  strcat(file, "_iqu.hdf");
303  }
304  }
305 
306  // then try sensor dir netCDF
307  if(access(file, R_OK) == -1) {
308  strcpy(file, path);
309  strcat(file, "rayleigh/rayleigh_");
310  strcat(file, sensorstr);
311  strcat(file, "_");
312  strcat(file, wavestr);
313  strcat(file, "_iqu.nc");
314 
315  // then try sensor dir HDF
316  if(access(file, R_OK) == -1) {
317  strcpy(file, path);
318  strcat(file, "rayleigh/rayleigh_");
319  strcat(file, sensorstr);
320  strcat(file, "_");
321  strcat(file, wavestr);
322  strcat(file, "_iqu.hdf");
323  }
324  }
325  read_rayleigh_lut(file, iw, pol_opt);
326  }
327  firstCall = 0;
328  }
329 
330 
331  /* windspeed index into tables */
332 
333  sigma_m = 0.0731 * sqrt(ws);
334  isigma2 = 1;
335  while (sigma_m > ray_sigma[isigma2] && isigma2 < NWIND)
336  isigma2++;
337  isigma1 = isigma2 - 1;
338 
339 
340 
341  /* interpolate Rayleigh Stokes vectors for each band */
342  gmult = 0;
343  if (l1rec->geom_per_band != NULL) {
344  gmult = 1;
345  }
346  for (iw = 0; iw < nwave; iw++) {
347 
348  ray_i = 0.0;
349  ray_q = 0.0;
350  ray_u = 0.0;
351 
352  /* geometry indices into tables */
353 
354  if ((iw == 0) || (gmult == 1)) {
355  ix = iw * gmult;
356 
357  for (m = 0; m < NORDER; m++) {
358  cosd_phi[m] = cos(r_phi[ix] * m / RADEG);
359  sind_phi[m] = sin(r_phi[ix] * m / RADEG);
360  }
361 
362  /* solar zenith indices */
363 
364  isol1 = ((int) r_solz[ix]) / 2;
365  isol2 = isol1 + 1;
366 
367  /* sensor zenith indices */
368 
369  for (isen2 = 0; isen2 < NSEN; isen2++) {
370  if (r_senz[ix] < ray_sen[isen2])
371  break;
372  }
373  isen1 = isen2 - 1;
374 
375  /* interpolation coefficients */
376  /* for solz > 88 or isol1 >= nsol - 1, use coeffs at end index */
377  if (isol1 >= NSOL - 1) {
378  isol1 = NSOL - 1;
379  isol2 = isol1;
380  p = 1.;
381  } else
382  p = (r_solz[ix] - ray_sol[isol1]) / (ray_sol[isol2] - ray_sol[isol1]);
383  q = (r_senz[ix] - ray_sen[isen1]) / (ray_sen[isen2] - ray_sen[isen1]);
384  airmass = 1. / r_csolz[ix] + 1. / r_csenz[ix];
385  }
386 
387 
388  /* interpolate the Rayleigh coefficients for each windspeed */
389 
390  for (isigma = isigma1; isigma <= isigma2; isigma++) {
391 
392  iwind = isigma - isigma1; /* 0 or 1 */
393 
394  ray_i_sig [iwind] = 0.;
395  ray_q_sig [iwind] = 0.;
396  ray_u_sig [iwind] = 0.;
397 
398  /* I component */
399 
400  for (m = 0; m < NORDER; m++) {
401 
402  f00 = ray_for_i[iw][isigma][isol1][m][isen1];
403  f10 = ray_for_i[iw][isigma][isol2][m][isen1];
404  f01 = ray_for_i[iw][isigma][isol1][m][isen2];
405  f11 = ray_for_i[iw][isigma][isol2][m][isen2];
406 
407  ray_i_sig[iwind] = ray_i_sig[iwind] +
408  ((1. - p)*(1. - q) * f00 + p * q * f11 + p * (1. - q) * f10 + q * (1. - p) * f01) * cosd_phi[m];
409  }
410 
411  if (pol_opt <= 0)
412  continue;
413 
414  /* Q component */
415 
416  for (m = 0; m < NORDER; m++) {
417 
418  f00 = ray_for_q[iw][isigma][isol1][m][isen1];
419  f10 = ray_for_q[iw][isigma][isol2][m][isen1];
420  f01 = ray_for_q[iw][isigma][isol1][m][isen2];
421  f11 = ray_for_q[iw][isigma][isol2][m][isen2];
422 
423  ray_q_sig[iwind] = ray_q_sig[iwind] +
424  ((1. - p)*(1. - q) * f00 + p * q * f11 + p * (1. - q) * f10 + q * (1. - p) * f01) * cosd_phi[m];
425  }
426 
427 
428  /* U component */
429 
430  for (m = 0; m < NORDER; m++) {
431 
432  f00 = ray_for_u[iw][isigma][isol1][m][isen1];
433  f10 = ray_for_u[iw][isigma][isol2][m][isen1];
434  f01 = ray_for_u[iw][isigma][isol1][m][isen2];
435  f11 = ray_for_u[iw][isigma][isol2][m][isen2];
436 
437  ray_u_sig[iwind] = ray_u_sig[iwind] +
438  ((1. - p)*(1. - q) * f00 + p * q * f11 + p * (1. - q) * f10 + q * (1. - p) * f01) * sind_phi[m];
439  }
440 
441  }
442 
443 
444  /* do linear interpolation between wind speeds */
445 
446  if (isigma1 == isigma2) {
447 
448  ray_i = ray_i_sig[0];
449 
450  if (pol_opt > 0) {
451  ray_q = ray_q_sig[0];
452  ray_u = ray_u_sig[0];
453  }
454 
455  } else {
456 
457  h = (sigma_m - ray_sigma[isigma1]) / (ray_sigma[isigma2] - ray_sigma[isigma1]);
458 
459  ray_i = ray_i_sig[0] + (ray_i_sig[1] - ray_i_sig[0]) * h;
460 
461  if (pol_opt > 0) {
462  ray_q = ray_q_sig[0] + (ray_q_sig[1] - ray_q_sig[0]) * h;
463  ray_u = ray_u_sig[0] + (ray_u_sig[1] - ray_u_sig[0]) * h;
464  }
465  }
466  /* Compute the rayleigh radiance */
467  fac = ray_press_wang(taur[iw], airmass, pr);
468  ipw = ip * nwave + iw;
469  l1rec->Lr[ipw] = Fo[iw] * ray_i*fac;
470  l1rec->L_q[ipw] = Fo[iw] * ray_q*fac;
471  l1rec->L_u[ipw] = Fo[iw] * ray_u*fac;
472  }
473 }
q
data_t q
Definition: decode_rs.h:74
sensorId2SensorDir
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:315
NORDER
#define NORDER
Definition: rayleigh.c:5
status
int status
Definition: l1_czcs_hdf.c:32
lowcase
char * lowcase(char *instr)
Definition: lowcase.c:10
NWIND
#define NWIND
Definition: rayleigh.c:6
NULL
#define NULL
Definition: decode_rs.h:63
viirs_ler_luts::pr
real pr
Definition: viirs_ler_luts.f95:38
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
h
float h[MODELMAX]
Definition: atrem_corl1.h:131
file
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
input
instr * input
Definition: main_l2binmatch.cpp:27
l12_proto.h
fac
#define fac
Definition: atmocor1_land.c:116
color_dtdb.path
string path
Definition: color_dtdb.py:221
hico.value_locate.x
list x
Definition: value_locate.py:64
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
NSEN
#define NSEN
Definition: rayleigh.c:4
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
read_rayleigh_lut
void read_rayleigh_lut(char *file, int32_t iw, int pol_opt)
Definition: rayleigh.c:69
STDPR
#define STDPR
Definition: l12_parms.h:85
sensorId2SensorName
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:273
rayleigh
void rayleigh(l1str *l1rec, int32_t ip)
Definition: rayleigh.c:169
NSOL
#define NSOL
Definition: rayleigh.c:3
NEWRAYTAB
#define NEWRAYTAB
Definition: l1.h:77
rdsensorinfo
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
subsensorId2SubsensorDir
const char * subsensorId2SubsensorDir(int subsensorId)
Definition: sensorInfo.c:329
ray_array
float ray_array[NWIND][NSOL][NORDER][NSEN]
Definition: rayleigh.c:13
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:91
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
ray_press_wang
float ray_press_wang(float taur, float airmass, float pr)
Definition: rayleigh.c:23


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