8 static int firstCall = 1;
10 static float ray_sen [
NSEN];
11 static float ray_sol [
NSOL];
12 static float ray_sigma[
NWIND];
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)));
30 static void check_dimension_size(
const char* file_name,
int nc_id,
const char* dim_name,
size_t length) {
35 status = nc_inq_dimid (nc_id, dim_name, &dim_id);
37 printf(
"-E- %s: Error looking for dimension %s from %s.\n", __FILE__, dim_name, file_name);
40 status = nc_inq_dimlen(nc_id, dim_id, &dim_length);
42 printf(
"-E- %s: Error reading dimension %s from %s.\n", __FILE__, dim_name, file_name);
45 if(dim_length != length) {
46 printf(
"-E- %s: Error with size of dimension %s from %s.\n", __FILE__, dim_name, file_name);
52 static void read_lut_variable(
const char*
file,
int nc_id,
const char* var_name,
float*
data) {
56 status = nc_inq_varid(nc_id, var_name, &var_id);
58 printf(
"-E- %s: Error looking for variable %s from %s.\n", __FILE__, var_name,
file);
61 status = nc_get_var_float(nc_id, var_id,
data);
63 printf(
"-E- %s: Error reading variable %s from %s.\n", __FILE__, var_name,
file);
76 printf(
"-E- %s: Error opening file %s.\n", __FILE__,
file);
79 printf(
"Loading Rayleigh LUT %s\n",
file);
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);
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]);
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]);
144 char file [FILENAME_MAX] =
"";
145 char path [FILENAME_MAX] =
"";
146 char wavestr[10] =
"";
147 char sensorstr[32] =
"";
152 float f00, f10, f01, f11;
158 float ray_i, ray_q, ray_u;
167 int32_t iw, ipw, gmult, ix;
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;
176 float ws =
l1rec->ws[ip];
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];
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];
207 if ((tmp_str = getenv(
"OCDATAROOT")) ==
NULL) {
208 printf(
"OCDATAROOT environment variable is not defined.\n");
213 strcat(
path,
"/eval/");
224 printf(
"-E- : Error allocating memory to ray_for_i\n");
229 printf(
"-E- : Error allocating memory to ray_for_q\n");
233 printf(
"-E- : Error allocating memory to ray_for_u\n");
240 for (iw = 0; iw < nwave; iw++) {
242 sprintf(wavestr,
"%d", (
int) wave[iw]);
246 if(
l1rec->l1file->subsensorID >= 0) {
249 strcat(
file,
"/rayleigh/rayleigh_");
250 strcat(
file, sensorstr);
252 strcat(
file, wavestr);
253 strcat(
file,
"_iqu.nc");
256 if(access(
file, R_OK) == -1) {
259 strcat(
file,
"/rayleigh/rayleigh_");
260 strcat(
file, sensorstr);
262 strcat(
file, wavestr);
263 strcat(
file,
"_iqu.hdf");
268 if(access(
file, R_OK) == -1) {
270 strcat(
file,
"rayleigh/rayleigh_");
271 strcat(
file, sensorstr);
273 strcat(
file, wavestr);
274 strcat(
file,
"_iqu.nc");
277 if(access(
file, R_OK) == -1) {
279 strcat(
file,
"rayleigh/rayleigh_");
280 strcat(
file, sensorstr);
282 strcat(
file, wavestr);
283 strcat(
file,
"_iqu.hdf");
294 sigma_m = 0.0731 * sqrt(ws);
296 while (sigma_m > ray_sigma[isigma2] && isigma2 <
NWIND)
298 isigma1 = isigma2 - 1;
307 for (iw = 0; iw < nwave; iw++) {
315 if ((iw == 0) || (gmult == 1)) {
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);
325 isol1 = ((
int) r_solz[ix]) / 2;
330 for (isen2 = 0; isen2 <
NSEN; isen2++) {
331 if (r_senz[ix] < ray_sen[isen2])
338 if (isol1 >=
NSOL - 1) {
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];
351 for (isigma = isigma1; isigma <= isigma2; isigma++) {
353 iwind = isigma - isigma1;
355 ray_i_sig [iwind] = 0.;
356 ray_q_sig [iwind] = 0.;
357 ray_u_sig [iwind] = 0.;
361 for (m = 0; m <
NORDER; m++) {
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];
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];
377 for (m = 0; m <
NORDER; m++) {
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];
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];
391 for (m = 0; m <
NORDER; m++) {
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];
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];
407 if (isigma1 == isigma2) {
409 ray_i = ray_i_sig[0];
412 ray_q = ray_q_sig[0];
413 ray_u = ray_u_sig[0];
418 h = (sigma_m - ray_sigma[isigma1]) / (ray_sigma[isigma2] - ray_sigma[isigma1]);
420 ray_i = ray_i_sig[0] + (ray_i_sig[1] - ray_i_sig[0]) *
h;
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;
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;