ocssw V2022
aerosol.c
Go to the documentation of this file.
1 /* ======================================================================================== */
2 /* module aerosol.c - functions to facilitate aerosol model selection and application */
3 /* */
4 /* Description: */
5 /* */
6 /* This code replaces the original set of fortran subroutines developed by M.Wang, H.Gordon,*/
7 /* and others (e.g., rho_a_sub_quad, linear_abc, funct_eps, load_aer, load_ss11) as well as */
8 /* the original get_aerosol() developed for MSl12. */
9 /* */
10 /* The functions herein read and interpolate the aerosol model tables, which are now stored */
11 /* as individual HDF files per model. Whfere sensor wavelengths differ from tabulated model */
12 /* wavelengths, interpolation is performed. Efficiencies are gained by remembering and */
13 /* reusing computational results when applicable. */
14 /* */
15 /* Primary Function: */
16 /* ---------------- */
17 /* aerosol() - compute aerosol reflectance using specified algorithm (main interface) */
18 /* */
19 /* Secondary Functions: */
20 /* ------------------- */
21 /* wangaer() - compute aerosol reflectance using Gordon & Wang 1994 algorithm */
22 /* fixedaot() - compute aerosol reflectance for fixed aot(lambda) */
23 /* fixedaer() - compute aerosol reflectance for fixed aerosol model */
24 /* get_angstrom() - compute angstrom coefficient (interface to l2_hdf_generic) */
25 /* diff_tran() - compute Rayleigh-aerosol diffuse trans for selected model pair, both paths */
26 /* */
27 /* Supporting Functions: */
28 /* -------------------- */
29 /* load_aermod() - loads the entire aerosol model table for the specified model list */
30 /* ss_to_ms_coef() - for one model, return coefficients of function relating single */
31 /* scattering to multiple scattering at the input geometry. */
32 /* rhoas_to_rhoa() - SS aerosol reflectance to MS aerosol reflectance */
33 /* rhoa_to_rhoas() - MS aerosol reflectance to SS aerosol reflectance */
34 /* model_select_wang() - M. Wang aerosol model selection process . */
35 /* model_select_angst() - Select model pair based on input Angstrom coefficient */
36 /* model_phase() - return phase function at model wavelengths at the input geometry. */
37 /* model_epsilon() - return model epsilon at input wavelengths for input geometry. */
38 /* model_transmittance() - compute path Rayleigh-aerosol diffuse trans for specified model */
39 /* model_taua() - compute AOT at input wavelengths using specified model */
40 /* aeroob - out-of-band water-vapor correction */
41 /* */
42 /* */
43 /* Written By: B. Franz, SAIC, NASA/GSFC Ocean Biology Processing Group, Summer 2004. */
44 /* W. Robinson, SAIC, 24 MAr 2017, modularized and enhanced for band-dependent */
45 /* geometry */
46 /* */
47 /* ======================================================================================== */
48 
49 #include <float.h>
50 
51 #include "l12_proto.h"
52 
53 #define MAXMODEL MAXAERMOD
54 //#define MAXSOLZ 33
55 //#define MAXSENZ 35
56 //#define MAXPHI 19
57 //#define MAXSCATT 75
58 //#define DTNTHETA 33
59 
60 static int32_t first_solz;
61 static int32_t first_senz;
62 static int32_t first_phi;
63 static int32_t first_scatt;
64 static int32_t first_dtnwave;
65 static int32_t first_dtntheta;
66 
67 static float pi = PI;
68 static double radeg = RADEG;
69 static float p0 = STDPR;
70 
71 static int have_ms = 0;
72 static int have_rh = 0;
73 static int have_sd = 0;
74 static int use_rh = 0;
75 static int use_netcdf = 0;
76 
77 static int32_t Nbands;
78 static int32_t Maxband; /* must be >= NBANDS */
79 
80 typedef struct aermod_struct {
81  char name[32];
82  float rh;
83  int sd;
84 
85  /* angstrom exponent (nbands+1)*/
86  float *angstrom;
87 
88  /* single-scattering albedo(nbands+1), extinction coefficient(nbands+1), phase function */
89  float *albedo;
90  float *extc;
91  float **phase;
92 
93  /* quadratic coefficients for SS to MS relationship */
94  float *acost;
95  float *bcost;
96  float *ccost;
97 
98  /* cubic coefficients for ms_epsilon atmospheric correction ..ZA */
99  float *ams_all;
100  float *bms_all;
101  float *cms_all;
102  float *dms_all;
103  float *ems_all;
104 
105 
106  /* Rayleigh-aerosol diffuse transmittance coeffs */
107  float **dtran_a;
108  float **dtran_b;
109 
110  /* derived quantities */
111  float **lnphase;
112  float **d2phase;
113 
114 } aermodstr;
115 
116 aermodstr* alloc_aermodstr(int nbands, int nscatt, int nphi, int nsolz, int nsenz, int ntheta) {
117  aermodstr *model;
118 
119  model = (aermodstr *) malloc(sizeof (aermodstr));
120  model->angstrom = (float*) malloc(nbands * sizeof (float));
121  model->extc = (float*) malloc(nbands * sizeof (float));
122  if(nscatt > 0) {
123  model->albedo = (float*) malloc(nbands * sizeof (float));
124  model->phase = allocate2d_float(nbands, nscatt);
125  model->lnphase = allocate2d_float(nbands, nscatt);
126  model->d2phase = allocate2d_float(nbands, nscatt);
127  model->acost = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
128  model->bcost = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
129  model->ccost = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
130  } else {
131  model->albedo = NULL;
132  model->phase = NULL;
133  model->lnphase = NULL;
134  model->d2phase = NULL;
135  model->acost = NULL;
136  model->bcost = NULL;
137  model->ccost = NULL;
138  }
139  model->dtran_a = allocate2d_float(nbands, ntheta);
140  model->dtran_b = allocate2d_float(nbands, ntheta);
141  model->ams_all = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
142  model->bms_all = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
143  model->cms_all = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
144  model->dms_all = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
145  model->ems_all = (float*) malloc(nbands * nsolz * nphi * nsenz * sizeof (float));
146 
147  return model;
148 }
149 
150 typedef struct aermodtab_struct {
151  int32_t sensorID;
152 
153  /* table dimensions */
154  int32_t nwave;
155  int32_t nmodel;
156  int32_t nsolz;
157  int32_t nsenz;
158  int32_t nphi;
159  int32_t nscatt;
160  int32_t dtran_nwave;
161  int32_t dtran_ntheta;
162 
163  /* table spectral bands and angles */
164  float *wave;
165  float *solz;
166  float *senz;
167  float *phi;
168  float *scatt;
169 
170  /* diffuse transmittance spectral bands and angles */
171  float *dtran_wave;
172  float *dtran_theta;
173  float *dtran_airmass;
174 
175  aermodstr **model;
176 
177 } aermodtabstr;
178 
179 typedef struct alphaT_struct {
180  int32_t modnum;
181  float angstrom;
182 } alphaTstr;
183 
184 typedef struct epsilonT_struct {
185  int32_t modnum;
186  float eps_obs;
187 } epsilonTstr;
188 
189 /* aerosol table */
190 static aermodtabstr *aertab = NULL;
191 
192 /* structure for carrying the geometry information */
193 typedef struct geom_strdef {
194  int gmult; /* band offset multiplier: 0 for nominal geometry
195  1 for band-dependent geometry */
196  float *senz;
197  float *solz;
198  float *phi;
199  float *csolz; /* cosine of solz */
200  float *csenz; /* cosine of senz */
201  float *airmass;
202  float *airmass_plp;
203  float *airmass_sph;
204 } geom_str;
205 
206 geom_str geom;
207 
208 /* global variable declarations */
209 static int loaded = 0;
210 static int interpol = 0;
211 static int32_t *iwatab;
212 static int32_t *iwdtab;
213 
214 static int32_t iwnir_s = -1;
215 static int32_t iwnir_l = -1;
216 
217 static int imm50 = -1;
218 static int imc50 = -1;
219 static int imc70 = -1;
220 static int imt90 = -1;
221 static int imt99 = -1;
222 //static int wang_modx = 0;
223 
224 static float mu0;
225 static float mu;
226 static float airmass;
227 
228 static int32_t evalmask = 0;
229 static int32_t aer_opt = 0;
230 static float airmass_plp;
231 static float airmass_sph;
232 
233 int cmpfunc(const void * a, const void * b) {
234  if (*(double*) a > *(double*) b) return 1;
235  else if (*(double*) a < *(double*) b) return -1;
236  else return 0;
237 }
238 
239 
240 
241 /* ---------------------------------------------------------------------------------------- */
242 /* first_deriv() - returns first derivative (dy/dx) of 1st or last array indices using a */
243 /* 4-pt Lagrangian interpolation. NOTE: It is assumed that 4 points exist. */
244 
245 /* ---------------------------------------------------------------------------------------- */
246 float first_deriv(float x[], float y[], int n) {
247  float a1, a2, a3, a4, a5, a6, d1;
248 
249  if (n == 0) {
250 
251  a1 = x[0] - x[1];
252  a2 = x[0] - x[2];
253  a3 = x[0] - x[3];
254  a4 = x[1] - x[2];
255  a5 = x[1] - x[3];
256  a6 = x[2] - x[3];
257 
258  d1 = y[0]*(1.0 / a1 + 1.0 / a2 + 1.0 / a3)
259  - a2 * a3 * y[1] / (a1 * a4 * a5)
260  + a1 * a3 * y[2] / (a2 * a4 * a6)
261  - a1 * a2 * y[3] / (a3 * a5 * a6);
262 
263  } else {
264 
265  a1 = x[n - 1] - x[n - 4];
266  a2 = x[n - 1] - x[n - 3];
267  a3 = x[n - 1] - x[n - 2];
268  a4 = x[n - 2] - x[n - 4];
269  a5 = x[n - 2] - x[n - 3];
270  a6 = x[n - 3] - x[n - 4];
271 
272  d1 = -a2 * a3 * y[n - 4] / (a6 * a4 * a1)
273  + a1 * a3 * y[n - 3] / (a6 * a5 * a2)
274  - a1 * a2 * y[n - 2] / (a4 * a5 * a3)
275  + y[n - 1]*(1.0 / a1 + 1.0 / a2 + 1.0 / a3);
276  }
277 
278  return (d1);
279 }
280 
281 
282 static void read_dimension_size(const char* file_name, int nc_id, const char* dim_name, int32_t *length) {
283  int dim_id;
284  size_t dim_length;
285  int status;
286 
287  status = nc_inq_dimid (nc_id, dim_name, &dim_id);
288  if (status != NC_NOERR) {
289  printf("-E- %s: Error looking for dimension %s from %s.\n", __FILE__, dim_name, file_name);
290  exit(1);
291  }
292  status = nc_inq_dimlen(nc_id, dim_id, &dim_length);
293  if (status != NC_NOERR) {
294  printf("-E- %s: Error reading dimension %s from %s.\n", __FILE__, dim_name, file_name);
295  exit(1);
296  }
297  *length = (int32_t) dim_length;
298 }
299 
300 
301 static void read_lut_variable(const char* file, int nc_id, const char* var_name, float* data) {
302  int var_id;
303  int status;
304 
305  status = nc_inq_varid(nc_id, var_name, &var_id);
306  if (status != NC_NOERR) {
307  printf("-E- %s: Error looking for variable %s from %s.\n", __FILE__, var_name, file);
308  exit(1);
309  }
310  status = nc_get_var_float(nc_id, var_id, data);
311  if (status != NC_NOERR) {
312  printf("-E- %s: Error reading variable %s from %s.\n", __FILE__, var_name, file);
313  exit(1);
314  }
315 }
316 
317 
318 /* ---------------------------------------------------------------------------------------- */
319 /* load_aermod() - loads the entire aerosol model table for the specified model list */
320 
321 /* ---------------------------------------------------------------------------------------- */
322 int load_aermod(int32_t sensorID, float wave[], int32_t nwave, char *aermodfile, char models[MAXAERMOD][32], int32_t nmodels) {
323  int nc_id;
324  int status;
325 
326  int32_t aer_nwave, nsolz, nsenz, nphi, nscatt, dtran_nwave, dtran_ntheta;
327 
328  float d1phase1;
329  float d1phaseN;
330  float rh;
331  int16_t sd;
332 
333  char file [FILENAME_MAX] = "";
334  char path [FILENAME_MAX] = "";
335 
336  int iw, im, is, iwbase, i;
337  static int firstCall = 1;
338 
339  if (firstCall == 1) {
340  if ((iwatab = (int32_t *) calloc(nwave, sizeof (int32_t))) == NULL) {
341  printf("Unable to allocate space for iwatab.\n");
342  exit(1);
343  }
344  if ((iwdtab = (int32_t *) calloc(nwave, sizeof (int32_t))) == NULL) {
345  printf("Unable to allocate space for iwdtab.\n");
346  exit(1);
347  }
348  firstCall = 0;
349  }
350 
351  printf("Loading aerosol models from %s\n", aermodfile);
352 
353  for (im = 0; im < nmodels + 1; im++) {
354 
355  if (im < nmodels) {
356 
357  // try netCDF first
358  strcpy(file, path);
359  strcat(file, aermodfile);
360  strcat(file, "_");
361  strcat(file, models[im]);
362  strcat(file, ".nc");
363  use_netcdf = 1;
364 
365  // if not try HDF4
366  if(access(file, R_OK) == -1) {
367  strcpy(file, path);
368  strcat(file, aermodfile);
369  strcat(file, "_");
370  strcat(file, models[im]);
371  strcat(file, ".hdf");
372  use_netcdf = 0;
373  }
374 
375  } else {
376 
377  // first try netCDF
378  strcpy(file, path);
379  strcat(file, aermodfile);
380  strcat(file, "_default.nc");
381  use_netcdf = 1;
382 
383  // if not try HDF4
384  if(access(file, R_OK) == -1) {
385  strcpy(file, path);
386  strcat(file, aermodfile);
387  strcat(file, "_default.hdf");
388  use_netcdf = 0;
389  }
390  }
391 
392  status = nc_open(file, NC_NOWRITE, &nc_id);
393  if (status != NC_NOERR) {
394  printf("-E- %s: Error opening file %s.\n", __FILE__, file);
395  exit(1);
396  }
397 
398  /* read dimensions which should be constant between models */
399  if(use_netcdf) {
400  read_dimension_size(file, nc_id, "wave", &aer_nwave);
401  read_dimension_size(file, nc_id, "solz", &nsolz);
402  read_dimension_size(file, nc_id, "senz", &nsenz);
403  read_dimension_size(file, nc_id, "phi", &nphi);
404  nscatt = 0;
405  read_dimension_size(file, nc_id, "dtran_wave", &dtran_nwave);
406  read_dimension_size(file, nc_id, "dtran_theta", &dtran_ntheta);
407  } else {
408  read_dimension_size(file, nc_id, "nwave", &aer_nwave);
409  read_dimension_size(file, nc_id, "nsolz", &nsolz);
410  read_dimension_size(file, nc_id, "nsenz", &nsenz);
411  read_dimension_size(file, nc_id, "nphi", &nphi);
412  read_dimension_size(file, nc_id, "nscatt", &nscatt);
413  read_dimension_size(file, nc_id, "dtran_nwave", &dtran_nwave);
414  read_dimension_size(file, nc_id, "dtran_ntheta", &dtran_ntheta);
415  }
416  // if(aer_nwave != nwave) {
417  // printf("-E- %s:%d - Error nwave dimension = %d not equal to nwave = %d in file = %s.\n",
418  // __FILE__, __LINE__, aer_nwave, nwave, file);
419  // exit(1);
420  // }
421 
422  if (im == 0) {
423  printf("Number of Wavelengths %d\n", aer_nwave);
424  printf("Number of Solar Zenith Angles %d\n", nsolz);
425  printf("Number of View Zenith Angles %d\n", nsenz);
426  printf("Number of Relative Azimuth Angles %d\n", nphi);
427  printf("Number of Scattering Angles %d\n", nscatt);
428  printf("Number of Diffuse Transmittance Wavelengths %d\n", dtran_nwave);
429  printf("Number of Diffuse Transmittance Zenith Angles %d\n", dtran_ntheta);
430 
431  first_solz = nsolz;
432  first_senz = nsenz;
433  first_phi = nphi;
434  first_scatt = nscatt;
435  first_dtnwave = dtran_nwave;
436  first_dtntheta = dtran_ntheta;
437 
438  // allocate the aerosol table
439  if ((aertab = (aermodtabstr *) calloc(1, sizeof (aermodtabstr))) == NULL) {
440  printf("Unable to allocate space for aerosol table.\n");
441  exit(1);
442  }
443 
444  aertab->nmodel = nmodels;
445  aertab->nwave = aer_nwave;
446  aertab->nsolz = nsolz;
447  aertab->nsenz = nsenz;
448  aertab->nphi = nphi;
449  aertab->nscatt = nscatt;
450 
451  aertab->wave = (float *) malloc(aer_nwave * sizeof (float));
452  aertab->solz = (float *) malloc(nsolz * sizeof (float));
453  aertab->senz = (float *) malloc(nsenz * sizeof (float));
454  aertab->phi = (float *) malloc(nphi * sizeof (float));
455  if(use_netcdf)
456  aertab->scatt = NULL;
457  else
458  aertab->scatt = (float *) malloc(nscatt * sizeof (float));
459 
460  aertab->dtran_nwave = dtran_nwave;
461  aertab->dtran_ntheta = dtran_ntheta;
462 
463  aertab->dtran_wave = (float *) malloc(dtran_nwave * sizeof (float));
464  aertab->dtran_theta = (float *) malloc(dtran_ntheta * sizeof (float));
465  aertab->dtran_airmass = (float *) malloc(dtran_ntheta * sizeof (float));
466 
467  // allocate the model tables
468  if ((aertab->model = (aermodstr **) calloc(1, (nmodels + 1) * sizeof (aermodstr*))) == NULL) {
469  printf("Unable to allocate space for %d aerosol models.\n", nmodels + 1);
470  exit(1);
471  }
472  for (i = 0; i < nmodels + 1; i++) {
473  if ((aertab->model[i] = alloc_aermodstr(aer_nwave + 1, nscatt, nphi, nsolz, nsenz, dtran_ntheta)) == NULL) {
474  printf("Unable to allocate space for aerosol model %d.\n", im);
475  exit(1);
476  }
477  }
478 
479  /* read SDSes which are constant between models */
480 
481  if(!use_netcdf)
482  read_lut_variable(file, nc_id, "scatt", aertab->scatt);
483  read_lut_variable(file, nc_id, "wave", aertab->wave);
484  read_lut_variable(file, nc_id, "solz", aertab->solz);
485  read_lut_variable(file, nc_id, "senz", aertab->senz);
486  read_lut_variable(file, nc_id, "phi", aertab->phi);
487  read_lut_variable(file, nc_id, "dtran_wave", aertab->dtran_wave);
488  read_lut_variable(file, nc_id, "dtran_theta", aertab->dtran_theta);
489 
490  } else {
491  /* check that all the aerosol files at least have the same
492  main dimensions */
493  if ((aertab->nsolz != first_solz) || (aertab->nsenz != first_senz) ||
494  (aertab->nphi != first_phi) || (aertab->nscatt != first_scatt) ||
495  (aertab->dtran_nwave != first_dtnwave) || (aertab->dtran_ntheta != first_dtntheta)) {
496  printf("-E- %s, %d: Error, Aerosol table %s\n",
497  __FILE__, __LINE__, file);
498  printf(" has different dimensions from previous tables\n");
499  exit(1);
500  }
501  }
502 
503  if (im < nmodels)
504  strncpy(aertab->model[im]->name, models[im], 32);
505  else
506  strncpy(aertab->model[im]->name, "default", 32);
507 
508  status = nc_get_att_float(nc_id, NC_GLOBAL, "RelativeHumidity", &rh);
509  if(status != NC_NOERR) {
510  status = nc_get_att_float(nc_id, NC_GLOBAL, "Relative Humidity", &rh);
511  }
512  if (status == NC_NOERR) {
513  aertab->model[im]->rh = rh;
514  have_rh = 1;
515  } else {
516  aertab->model[im]->rh = -1.0;
517  have_rh = 0;
518  }
519 
520  if(use_netcdf) {
521  int fmf;
522  status = nc_get_att_int(nc_id, NC_GLOBAL, "AerosolFMF", &fmf);
523  if(status != NC_NOERR) {
524  printf("-E- %s, %d: Error, Aerosol table %s does not have AerosolFMF\n",
525  __FILE__, __LINE__, file);
526  exit(1);
527  }
528  sd = 100 - fmf;
529  aertab->model[im]->sd = sd;
530  } else {
531  status = nc_get_att_short(nc_id, NC_GLOBAL, "Size Distribution", &sd);
532  if (status == NC_NOERR) {
533  aertab->model[im]->sd = sd;
534  have_sd = 1;
535  } else {
536  aertab->model[im]->sd = -1;
537  have_sd = 0;
538  }
539  }
540 
541  if(!use_netcdf) {
542  read_lut_variable(file, nc_id, "albedo", aertab->model[im]->albedo);
543  read_lut_variable(file, nc_id, "phase", aertab->model[im]->phase[0]);
544  read_lut_variable(file, nc_id, "acost", aertab->model[im]->acost);
545  read_lut_variable(file, nc_id, "bcost", aertab->model[im]->bcost);
546  read_lut_variable(file, nc_id, "ccost", aertab->model[im]->ccost);
547  }
548  read_lut_variable(file, nc_id, "extc", aertab->model[im]->extc);
549 
550  // check for multi-scattering epsilon tables
551  if(nc_inq_varid(nc_id, "ams_all", &status) == NC_NOERR) {
552  have_ms = 1;
553  read_lut_variable(file, nc_id, "ams_all", aertab->model[im]->ams_all);
554  read_lut_variable(file, nc_id, "bms_all", aertab->model[im]->bms_all);
555  read_lut_variable(file, nc_id, "cms_all", aertab->model[im]->cms_all);
556  read_lut_variable(file, nc_id, "dms_all", aertab->model[im]->dms_all);
557  read_lut_variable(file, nc_id, "ems_all", aertab->model[im]->ems_all);
558  }
559 
560  read_lut_variable(file, nc_id, "dtran_a", aertab->model[im]->dtran_a[0]);
561  read_lut_variable(file, nc_id, "dtran_b", aertab->model[im]->dtran_b[0]);
562 
563  nc_close(nc_id);
564 
565  /* compute angstrom exponent for each model wavelength relative to max wavelength */
566  iwbase = windex(865, aertab->wave, aertab->nwave);
567  for (iw = 0; iw < aertab->nwave; iw++) {
568  if (iw != iwbase)
569  aertab->model[im]->angstrom[iw] = -log(aertab->model[im]->extc[iw] / aertab->model[im]->extc[iwbase]) /
570  log(aertab->wave[iw] / aertab->wave[iwbase]);
571  }
572  aertab->model[im]->angstrom[iwbase] = aertab->model[im]->angstrom[iwbase - 1];
573 
574  /* precompute log of phase function and 2nd derivative (for cubic interp) */
575  if(aertab->scatt) {
576  for (iw = 0; iw < aertab->nwave; iw++) {
577  for (is = 0; is < aertab->nscatt; is++) {
578  aertab->model[im]->lnphase[iw][is] = log(aertab->model[im]->phase[iw][is]);
579  }
580  d1phase1 = first_deriv(aertab->scatt, &aertab->model[im]->lnphase[iw][0], 0);
581  d1phaseN = first_deriv(aertab->scatt, &aertab->model[im]->lnphase[iw][0], aertab->nscatt);
582  spline(aertab->scatt,
583  &aertab->model[im]->lnphase[iw][0],
584  aertab->nscatt,
585  d1phase1,
586  d1phaseN,
587  &aertab->model[im]->d2phase[iw][0]);
588  }
589  }
590 
591  /* precompute airmass for diffuse transmittance */
592  for (iw = 0; iw < aertab->dtran_nwave; iw++) {
593  for (is = 0; is < aertab->dtran_ntheta; is++) {
594  aertab->dtran_airmass[is] = 1.0 / cos(aertab->dtran_theta[is] / radeg);
595  }
596  }
597  }
598 
599  aertab->nmodel = nmodels;
600  aertab->sensorID = sensorID;
601 
602  /* Require number of table wavelengths to equal number of sensor wavelengths */
603  if (aertab->nwave != nwave) {
604  printf("Number of aerosol LUT wavelengths (%d) is not equal to number of sensor wavelengths (%d).\n",
605  aertab->nwave,nwave);
606  exit(1);
607  }
608  if (aertab->dtran_nwave != nwave) {
609  printf("Number of aerosol diffue trans LUT wavelengths (%d) is not equal to number of sensor wavelengths (%d).\n",
610  aertab->dtran_nwave,nwave);
611  exit(1);
612  }
613 
614  /* map sensor wavelengths to table wavelengths */
615  printf("Wavelengths - Sensor\tAerosol model\tDiffuse transmittance\n");
616  for (iw = 0; iw < nwave; iw++) {
617  iwatab[iw] = windex(wave[iw], aertab->wave, aertab->nwave);
618  iwdtab[iw] = windex(wave[iw], aertab->dtran_wave, aertab->dtran_nwave);
619  printf("\t %6.2f\t %6.2f\t %6.2f\n", wave[iw],aertab->wave[iwatab[iw]], aertab->dtran_wave[iwdtab[iw]]);
620  }
621 
622 // for (iw = 0; iw < nwave; iw++) {
623 // iwatab[iw] = windex(wave[iw], aertab->wave, aertab->nwave);
624 // iwdtab[iw] = windex(wave[iw], aertab->dtran_wave, aertab->dtran_nwave);
625 // if (fabsf(wave[iw] - aertab->wave[iwatab[iw]]) > 0.5) {
626 // printf("Aerosol model coefficients will be interpolated for %5.1f nm band.\n", wave[iw]);
627 // interpol = 1;
628 // }
629 // }
630 //
631 // /* in case of more aertab wavelengths then sensor wavelengths */
632 // if (aertab->nwave != nwave)
633 // interpol = 1;
634 
635  loaded = 1;
636 
637  return (0);
638 }
639 
640 
641 
642 #define INDEX(iw,isol,iphi,isen) (iw*aertab->nsolz*aertab->nphi*aertab->nsenz + isol*aertab->nphi*aertab->nsenz + iphi*aertab->nsenz + isen)
643 
644 /* ---------------------------------------------------------------------------------------- */
645 /* ss_to_ms_coef() - for one model, return coefficients of function relating single */
646 /* scattering to multiple scattering at the input geometry. */
647 /* */
648 /* This is effectively a C version of M. Wangs linear_a_b_c.f. The program optimizes for */
649 /* multiple calls at the same geometry by computing for all models on the first call with a */
650 /* new geometry. It returns pointers to the internal static arrays of coefficients. */
651 /* */
652 /* B. Franz, 1 June 2004. */
653 /* W. Robinson, SAIC adapt to band-ependent geometry */
654 
655 /* ---------------------------------------------------------------------------------------- */
656 void ss_to_ms_coef(int modnum, geom_str *geom, float **a, float **b, float **c) {
657  static float lastsolz = -999.;
658  static float lastsenz = -999.;
659  static float lastphi = -999.;
660 
661  static int computed[MAXMODEL];
662 
663  static float *a_coef[MAXMODEL];
664  static float *b_coef[MAXMODEL];
665  static float *c_coef[MAXMODEL];
666 
667  static float *p, *q, *r;
668  static float as000, as100, as010, as110, as001, as011, as101, as111;
669  static float ai000, ai100, ai010, ai110, ai001, ai011, ai101, ai111;
670  static float ac000, ac100, ac010, ac110, ac001, ac011, ac101, ac111;
671 
672  static int *isolz1, *isolz2;
673  static int *isenz1, *isenz2;
674  static int *iphi1, *iphi2;
675 
676  static float *p_ar, *q_ar, *r_ar;
677  static float p_cnst, q_cnst, r_cnst;
678  static int *isolz1_ar, *isolz2_ar, isolz1_cnst, isolz2_cnst;
679  static int *isenz1_ar, *isenz2_ar, isenz1_cnst, isenz2_cnst;
680  static int *iphi1_ar, *iphi2_ar, iphi1_cnst, iphi2_cnst;
681 
682  float aphi;
683  float px, qx, rx;
684  int im, iw, i, ig;
685  static int firstCall = 1;
686  static int gmult;
687 
688  if (firstCall == 1) {
689  firstCall = 0;
690  for (i = 0; i < MAXMODEL; i++) {
691  if ((a_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
692  printf("Unable to allocate space for a_coef.\n");
693  exit(1);
694  }
695  if ((b_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
696  printf("Unable to allocate space for rhoa.\n");
697  exit(1);
698  }
699  if ((c_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
700  printf("Unable to allocate space for rhoa.\n");
701  exit(1);
702  }
703  }
704  /* set up indicies, weights for band-dependent or nominal geometry */
705  if ((geom->gmult == 0) || (interpol == 1)) {
706  gmult = 0;
707  p = &p_cnst;
708  q = &q_cnst;
709  r = &r_cnst;
710  isolz1 = &isolz1_cnst;
711  isolz2 = &isolz2_cnst;
712  isenz1 = &isenz1_cnst;
713  isenz2 = &isenz2_cnst;
714  iphi1 = &iphi1_cnst;
715  iphi2 = &iphi2_cnst;
716  } else
717  gmult = 1;
718  {
719  if (((p_ar = (float *) malloc(aertab->nwave * sizeof (float)))
720  == NULL) ||
721  ((q_ar = (float *) malloc(aertab->nwave * sizeof (float)))
722  == NULL) ||
723  ((r_ar = (float *) malloc(aertab->nwave * sizeof (float)))
724  == NULL)) {
725  printf("Unable to allocate space for p, q, r weights.\n");
726  exit(1);
727  }
728  if (((isolz1_ar = (int *) malloc(aertab->nwave * sizeof (int)))
729  == NULL) ||
730  ((isenz1_ar = (int *) malloc(aertab->nwave * sizeof (int)))
731  == NULL) ||
732  ((iphi1_ar = (int *) malloc(aertab->nwave * sizeof (int)))
733  == NULL)) {
734  printf("Unable to allocate space for interp indicies 1.\n");
735  exit(1);
736  }
737  if (((isolz2_ar = (int *) malloc(aertab->nwave * sizeof (int)))
738  == NULL) ||
739  ((isenz2_ar = (int *) malloc(aertab->nwave * sizeof (int)))
740  == NULL) ||
741  ((iphi2_ar = (int *) malloc(aertab->nwave * sizeof (int)))
742  == NULL)) {
743  printf("Unable to allocate space for interp indicies 2.\n");
744  exit(1);
745  }
746  p = p_ar;
747  q = q_ar;
748  r = r_ar;
749  isolz1 = isolz1_ar;
750  isolz2 = isolz2_ar;
751  isenz1 = isenz1_ar;
752  isenz2 = isenz2_ar;
753  iphi1 = iphi1_ar;
754  iphi2 = iphi2_ar;
755  }
756  }
757 
758  if ((geom->solz[0] != lastsolz) || (geom->senz[0] != lastsenz) ||
759  (geom->phi[0] != lastphi)) {
760  for (im = 0; im < aertab->nmodel; im++)
761  computed[im] = 0;
762 
763  for (iw = 0; iw < aertab->nwave; iw++) {
764  ig = iw * gmult;
765  /* find bracketing solar indices */
766  for (i = 0; i < aertab->nsolz; i++) {
767  if (geom->solz[ig] < aertab->solz[i])
768  break;
769  }
770  isolz1[iw] = MAX(i - 1, 0);
771  isolz2[iw] = MIN(i, aertab->nsolz - 1);
772  if (isolz2[iw] != isolz1[iw])
773  r[iw] = (geom->solz[ig] - aertab->solz[isolz1[iw]]) /
774  (aertab->solz[isolz2[iw]] - aertab->solz[isolz1[iw]]);
775  else
776  r[iw] = 0.0;
777 
778  /* find bracketing view indices */
779  for (i = 0; i < aertab->nsenz; i++) {
780  if (geom->senz[ig] < aertab->senz[i])
781  break;
782  }
783  isenz1[iw] = MAX(i - 1, 0);
784  isenz2[iw] = MIN(i, aertab->nsenz - 1);
785  if (isenz2[iw] != isenz1[iw])
786  p[iw] = (geom->senz[ig] - aertab->senz[isenz1[iw]]) /
787  (aertab->senz[isenz2[iw]] - aertab->senz[isenz1[iw]]);
788  else
789  p[iw] = 0.0;
790 
791  /* find bracketing azimuth indices */
792  aphi = fabs(geom->phi[ig]);
793  for (i = 0; i < aertab->nphi; i++) {
794  if (aphi < aertab->phi[i])
795  break;
796  }
797  iphi1[iw] = MAX(i - 1, 0);
798  iphi2[iw] = MIN(i, aertab->nphi - 1);
799  if (iphi2[iw] != iphi1[iw])
800  q[iw] = (aphi - aertab->phi[iphi1[iw]]) /
801  (aertab->phi[iphi2[iw]] - aertab->phi[iphi1[iw]]);
802  else
803  q[iw] = 0.0;
804  if (gmult == 0) break;
805  }
806 
807  /* remember last geometry */
808  lastsolz = geom->solz[0];
809  lastsenz = geom->senz[0];
810  lastphi = geom->phi[0];
811  }
812 
813  if (!computed[modnum]) {
814  im = modnum;
815  computed[modnum] = 1;
816 
817  for (iw = 0; iw < aertab->nwave; iw++) {
818  ig = iw * gmult;
819  px = p[ig];
820  qx = q[ig];
821  rx = r[ig];
822  if (isolz2[ig] == 0) {
823  as000 = aertab->model[im]->acost[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
824  as100 = aertab->model[im]->acost[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
825  as001 = aertab->model[im]->acost[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
826  as101 = aertab->model[im]->acost[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
827 
828  ai000 = aertab->model[im]->bcost[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
829  ai100 = aertab->model[im]->bcost[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
830  ai001 = aertab->model[im]->bcost[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
831  ai101 = aertab->model[im]->bcost[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
832 
833  ac000 = aertab->model[im]->ccost[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
834  ac100 = aertab->model[im]->ccost[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
835  ac001 = aertab->model[im]->ccost[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
836  ac101 = aertab->model[im]->ccost[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
837 
838  a_coef[im][iw] = (1. - px)*(1. - rx) * as000 + px * rx * as101
839  + (1. - px) * rx * as001 + px * (1. - rx) * as100;
840 
841  b_coef[im][iw] = (1. - px)*(1. - rx) * ai000 + px * rx * ai101
842  + (1. - px) * rx * ai001 + px * (1. - qx)*(1. - rx) * ai100;
843 
844  c_coef[im][iw] = (1. - px)*(1. - rx) * ac000 + px * rx * ac101
845  + (1. - px) * rx * ac001 + px * (1. - qx)*(1. - rx) * ac100;
846  } else {
847  as000 = aertab->model[im]->acost[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
848  as100 = aertab->model[im]->acost[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
849  as010 = aertab->model[im]->acost[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
850  as110 = aertab->model[im]->acost[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
851  as001 = aertab->model[im]->acost[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
852  as011 = aertab->model[im]->acost[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
853  as101 = aertab->model[im]->acost[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
854  as111 = aertab->model[im]->acost[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
855 
856  ai000 = aertab->model[im]->bcost[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
857  ai100 = aertab->model[im]->bcost[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
858  ai010 = aertab->model[im]->bcost[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
859  ai110 = aertab->model[im]->bcost[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
860  ai001 = aertab->model[im]->bcost[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
861  ai011 = aertab->model[im]->bcost[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
862  ai101 = aertab->model[im]->bcost[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
863  ai111 = aertab->model[im]->bcost[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
864 
865  ac000 = aertab->model[im]->ccost[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
866  ac100 = aertab->model[im]->ccost[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
867  ac010 = aertab->model[im]->ccost[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
868  ac110 = aertab->model[im]->ccost[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
869  ac001 = aertab->model[im]->ccost[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
870  ac011 = aertab->model[im]->ccost[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
871  ac101 = aertab->model[im]->ccost[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
872  ac111 = aertab->model[im]->ccost[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
873 
874  a_coef[im][iw] = (1. - px)*(1. - qx)*(1. - rx) * as000 + px * qx * rx * as111
875  + px * (1. - qx) * rx * as101 + (1. - px) * qx * (1. - rx) * as010
876  + px * qx * (1. - rx) * as110 + (1. - px)*(1. - qx) * rx * as001
877  + (1. - px) * qx * rx * as011 + px * (1. - qx)*(1. - rx) * as100;
878 
879  b_coef[im][iw] = (1. - px)*(1. - qx)*(1. - rx) * ai000 + px * qx * rx * ai111
880  + px * (1. - qx) * rx * ai101 + (1. - px) * qx * (1. - rx) * ai010
881  + px * qx * (1. - rx) * ai110 + (1. - px)*(1. - qx) * rx * ai001
882  + (1. - px) * qx * rx * ai011 + px * (1. - qx)*(1. - rx) * ai100;
883 
884  c_coef[im][iw] = (1. - px)*(1. - qx)*(1. - rx) * ac000 + px * qx * rx * ac111
885  + px * (1. - qx) * rx * ac101 + (1. - px) * qx * (1. - rx) * ac010
886  + px * qx * (1. - rx) * ac110 + (1. - px)*(1. - qx) * rx * ac001
887  + (1. - px) * qx * rx * ac011 + px * (1. - qx)*(1. - rx) * ac100;
888  }
889  }
890  }
891 
892  /* return pointers to coeffs for this geometry */
893  *a = &a_coef[modnum][0];
894  *b = &b_coef[modnum][0];
895  *c = &c_coef[modnum][0];
896 
897  return;
898 }
899 
900 /* ---------------------------------------------------------------------------------------- */
901 /* fresnel_coef() - computes Fresnel reflectance coefficient for specified index of refr. */
902 
903 /* ---------------------------------------------------------------------------------------- */
904 float fresnel_coef(float mu, float index) {
905  float sq, r2, q1;
906 
907  sq = sqrt(pow(index, 2.0) - 1.0 + pow(mu, 2.0));
908  r2 = pow((mu - sq) / (mu + sq), 2.0);
909  q1 = (1.0 - pow(mu, 2.0) - mu * sq) / (1.0 - pow(mu, 2.0) + mu * sq);
910 
911  return (r2 * (q1 * q1 + 1.0) / 2.0);
912 }
913 
914 
915 
916 /* ---------------------------------------------------------------------------------------- */
917 /* ms_eps_coef() - for a given model, returns ams, bms, cms, dms and ems coefficients to */
918 /* compute ms_reflectance at the input geometry. */
919 /* Also, the program optimizes for multiple calls at the same geometry */
920 /* by computing for all models on the first call with a new geometry. */
921 /* It returns pointers to the internal static arrays of coefficients. */
922 /* */
923 /* Z. Ahmad July 08, 2014 */
924 /* W. Robinson, SAIC 23 Mar 2017 add band-dependent geometry to code */
925 /* M. Zhang 06/19/2019 fixed the interpolation part */
926 
927 /* ---------------------------------------------------------------------------------------- */
928 
929 
930 
931 void ms_eps_coef(int modnum, int32_t iwnir_l, float wave[], geom_str *geom,
932  float **a, float **b, float **c, float **d, float **e)
933  {
934  static float lastsolz = -999.;
935  static float lastsenz = -999.;
936  static float lastphi = -999.;
937 
938  static int computed[MAXMODEL];
939 
940  static float *a_coef[MAXMODEL];
941  static float *b_coef[MAXMODEL];
942  static float *c_coef[MAXMODEL];
943  static float *d_coef[MAXMODEL];
944  static float *e_coef[MAXMODEL];
945 
946  static float *p, *q, *r;
947  static float as000, as100, as010, as110, as001, as011, as101, as111;
948  static float ai000, ai100, ai010, ai110, ai001, ai011, ai101, ai111;
949  static float ac000, ac100, ac010, ac110, ac001, ac011, ac101, ac111;
950  static float ad000, ad100, ad010, ad110, ad001, ad011, ad101, ad111;
951  static float ae000, ae100, ae010, ae110, ae001, ae011, ae101, ae111;
952 
953  static int *isolz1, *isolz2;
954  static int *isenz1, *isenz2;
955  static int *iphi1, *iphi2;
956 
957  static float *p_ar, *q_ar, *r_ar;
958  static float p_cnst, q_cnst, r_cnst;
959  static int *isolz1_ar, *isolz2_ar, isolz1_cnst, isolz2_cnst;
960  static int *isenz1_ar, *isenz2_ar, isenz1_cnst, isenz2_cnst;
961  static int *iphi1_ar, *iphi2_ar, iphi1_cnst, iphi2_cnst;
962  static int gmult;
963 
964  float aphi;
965  float px, qx, rx;
966  int im, iw, i, ig;
967  static int firstCall = 1;
968 
969  if (firstCall == 1) {
970  firstCall = 0;
971  for (i = 0; i < MAXMODEL; i++) {
972  if ((a_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
973  printf("Unable to allocate space for a_coef.\n");
974  exit(1);
975  }
976  if ((b_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
977  printf("Unable to allocate space for b_coef.\n");
978  exit(1);
979  }
980  if ((c_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
981  printf("Unable to allocate space for c_coef.\n");
982  exit(1);
983  }
984  if ((d_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
985  printf("Unable to allocate space for d_coef.\n");
986  exit(1);
987  }
988 
989  if ((e_coef[i] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
990  printf("Unable to allocate space for e_coef.\n");
991  exit(1);
992  }
993 
994  }
995  /* set up indicies, weights for band-dependent or nominal geometry */
996  if ((geom->gmult == 0) || (interpol == 1)) {
997  gmult = 0;
998  p = &p_cnst;
999  q = &q_cnst;
1000  r = &r_cnst;
1001  isolz1 = &isolz1_cnst;
1002  isolz2 = &isolz2_cnst;
1003  isenz1 = &isenz1_cnst;
1004  isenz2 = &isenz2_cnst;
1005  iphi1 = &iphi1_cnst;
1006  iphi2 = &iphi2_cnst;
1007  } else {
1008  gmult = 1;
1009  if (((p_ar = (float *) malloc(aertab->nwave * sizeof (float)))
1010  == NULL) ||
1011  ((q_ar = (float *) malloc(aertab->nwave * sizeof (float)))
1012  == NULL) ||
1013  ((r_ar = (float *) malloc(aertab->nwave * sizeof (float)))
1014  == NULL)) {
1015  printf("Unable to allocate space for p, q, r weights.\n");
1016  exit(1);
1017  }
1018  if (((isolz1_ar = (int *) malloc(aertab->nwave * sizeof (int)))
1019  == NULL) ||
1020  ((isenz1_ar = (int *) malloc(aertab->nwave * sizeof (int)))
1021  == NULL) ||
1022  ((iphi1_ar = (int *) malloc(aertab->nwave * sizeof (int)))
1023  == NULL)) {
1024  printf("Unable to allocate space for interp indicies 1.\n");
1025  exit(1);
1026  }
1027  if (((isolz2_ar = (int *) malloc(aertab->nwave * sizeof (int)))
1028  == NULL) ||
1029  ((isenz2_ar = (int *) malloc(aertab->nwave * sizeof (int)))
1030  == NULL) ||
1031  ((iphi2_ar = (int *) malloc(aertab->nwave * sizeof (int)))
1032  == NULL)) {
1033  printf("Unable to allocate space for interp indicies 2.\n");
1034  exit(1);
1035  }
1036  p = p_ar;
1037  q = q_ar;
1038  r = r_ar;
1039  isolz1 = isolz1_ar;
1040  isolz2 = isolz2_ar;
1041  isenz1 = isenz1_ar;
1042  isenz2 = isenz2_ar;
1043  iphi1 = iphi1_ar;
1044  iphi2 = iphi2_ar;
1045  }
1046  }
1047 
1048  if (geom->solz[0] != lastsolz || geom->senz[0] != lastsenz ||
1049  geom->phi[0] != lastphi) {
1050  for (im = 0; im < aertab->nmodel; im++)
1051  computed[im] = 0;
1052 
1053  for (iw = 0; iw < aertab->nwave; iw++) {
1054  ig = iw * gmult;
1055  /* find bracketing solar indices */
1056  for (i = 0; i < aertab->nsolz; i++) {
1057  if (geom->solz[ig] < aertab->solz[i])
1058  break;
1059  }
1060  isolz1[iw] = MAX(i - 1, 0);
1061  isolz2[iw] = MIN(i, aertab->nsolz - 1);
1062  if (isolz2[iw] != isolz1[iw])
1063  r[iw] = (geom->solz[ig] - aertab->solz[isolz1[iw]]) /
1064  (aertab->solz[isolz2[iw]] - aertab->solz[isolz1[iw]]);
1065  else
1066  r[iw] = 0.0;
1067 
1068  /* find bracketing view indices */
1069  for (i = 0; i < aertab->nsenz; i++) {
1070  if (geom->senz[ig] < aertab->senz[i])
1071  break;
1072  }
1073  isenz1[iw] = MAX(i - 1, 0);
1074  isenz2[iw] = MIN(i, aertab->nsenz - 1);
1075  if (isenz2[iw] != isenz1[iw])
1076  p[iw] = (geom->senz[ig] - aertab->senz[isenz1[iw]]) /
1077  (aertab->senz[isenz2[iw]] - aertab->senz[isenz1[iw]]);
1078  else
1079  p[iw] = 0.0;
1080 
1081  /* find bracketing azimuth indices */
1082  aphi = fabs(geom->phi[ig]);
1083  for (i = 0; i < aertab->nphi; i++) {
1084  if (aphi < aertab->phi[i])
1085  break;
1086  }
1087  iphi1[iw] = MAX(i - 1, 0);
1088  iphi2[iw] = MIN(i, aertab->nphi - 1);
1089  if (iphi2[iw] != iphi1[iw])
1090  q[iw] = (aphi - aertab->phi[iphi1[iw]]) /
1091  (aertab->phi[iphi2[iw]] - aertab->phi[iphi1[iw]]);
1092  else
1093  q[iw] = 0.0;
1094  if (gmult == 0) break;
1095  }
1096 
1097  /* remember last geometry */
1098  lastsolz = geom->solz[0];
1099  lastsenz = geom->senz[0];
1100  lastphi = geom->phi[0];
1101 
1102  }
1103 
1104  im = modnum;
1105 
1106  if (!computed[modnum]) {
1107  im = modnum;
1108  computed[modnum] = 1;
1109 
1110  for (iw = 0; iw < aertab->nwave; iw++) {
1111  ig = iw * gmult;
1112  px = p[ig];
1113  qx = q[ig];
1114  rx = r[ig];
1115  if (isolz2[ig] == 0) {
1116  as000 = aertab->model[im]->ams_all[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1117  as100 = aertab->model[im]->ams_all[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1118  as001 = aertab->model[im]->ams_all[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1119  as101 = aertab->model[im]->ams_all[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1120 
1121  ai000 = aertab->model[im]->bms_all[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1122  ai100 = aertab->model[im]->bms_all[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1123  ai001 = aertab->model[im]->bms_all[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1124  ai101 = aertab->model[im]->bms_all[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1125 
1126  ac000 = aertab->model[im]->cms_all[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1127  ac100 = aertab->model[im]->cms_all[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1128  ac001 = aertab->model[im]->cms_all[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1129  ac101 = aertab->model[im]->cms_all[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1130 
1131  ad000 = aertab->model[im]->dms_all[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1132  ad100 = aertab->model[im]->dms_all[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1133  ad001 = aertab->model[im]->dms_all[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1134  ad101 = aertab->model[im]->dms_all[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1135 
1136  ae000 = aertab->model[im]->ems_all[INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1137  ae100 = aertab->model[im]->ems_all[INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1138  ae001 = aertab->model[im]->ems_all[INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1139  ae101 = aertab->model[im]->ems_all[INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1140 
1141  a_coef[im][iw] = (1. - px)*(1. - rx) * as000 + px * rx * as101
1142  + (1. - px) * rx * as001 + px * (1. - rx) * as100;
1143 
1144  b_coef[im][iw] = (1. - px)*(1. - rx) * ai000 + px * rx * ai101
1145  + (1. - px) * rx * ai001 + px * (1. - qx)*(1. - rx) * ai100;
1146 
1147  c_coef[im][iw] = (1. - px)*(1. - rx) * ac000 + px * rx * ac101
1148  + (1. - px) * rx * ac001 + px * (1. - qx)*(1. - rx) * ac100;
1149 
1150  d_coef[im][iw] = (1. - px)*(1. - rx) * ad000 + px * rx * ad101
1151  + (1. - px) * rx * ad001 + px * (1. - qx)*(1. - rx) * ad100;
1152 
1153  e_coef[im][iw] = (1. - px)*(1. - rx) * ae000 + px * rx * ae101
1154  + (1. - px) * rx * ae001 + px * (1. - qx)*(1. - rx) * ae100;
1155  } else {
1156  /* printf("coeffs: as000,ai000,ac000,ad000,ae000\n"); */
1157 
1158  as000 = aertab->model[im]->ams_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1159  as100 = aertab->model[im]->ams_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1160  as010 = aertab->model[im]->ams_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1161  as110 = aertab->model[im]->ams_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1162  as001 = aertab->model[im]->ams_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1163  as011 = aertab->model[im]->ams_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1164  as101 = aertab->model[im]->ams_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1165  as111 = aertab->model[im]->ams_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1166 
1167  ai000 = aertab->model[im]->bms_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1168  ai100 = aertab->model[im]->bms_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1169  ai010 = aertab->model[im]->bms_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1170  ai110 = aertab->model[im]->bms_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1171  ai001 = aertab->model[im]->bms_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1172  ai011 = aertab->model[im]->bms_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1173  ai101 = aertab->model[im]->bms_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1174  ai111 = aertab->model[im]->bms_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1175 
1176  ac000 = aertab->model[im]->cms_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1177  ac100 = aertab->model[im]->cms_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1178  ac010 = aertab->model[im]->cms_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1179  ac110 = aertab->model[im]->cms_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1180  ac001 = aertab->model[im]->cms_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1181  ac011 = aertab->model[im]->cms_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1182  ac101 = aertab->model[im]->cms_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1183  ac111 = aertab->model[im]->cms_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1184 
1185  ad000 = aertab->model[im]->dms_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1186  ad100 = aertab->model[im]->dms_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1187  ad010 = aertab->model[im]->dms_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1188  ad110 = aertab->model[im]->dms_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1189  ad001 = aertab->model[im]->dms_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1190  ad011 = aertab->model[im]->dms_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1191  ad101 = aertab->model[im]->dms_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1192  ad111 = aertab->model[im]->dms_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1193 
1194  ae000 = aertab->model[im]->ems_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1195  ae100 = aertab->model[im]->ems_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1196  ae010 = aertab->model[im]->ems_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1197  ae110 = aertab->model[im]->ems_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1198  ae001 = aertab->model[im]->ems_all[INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1199  ae011 = aertab->model[im]->ems_all[INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1200  ae101 = aertab->model[im]->ems_all[INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1201  ae111 = aertab->model[im]->ems_all[INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1202 
1203 
1204  a_coef[im][iw] = (1. - rx)*(1. - qx)*(1. - px) * as000 + rx * qx * px * as111
1205  + rx * (1. - qx) * px * as101 + (1. - rx) * qx * (1. - px) * as010
1206  + rx * qx * (1. - px) * as110 + (1. - rx)*(1. - qx) * px * as001
1207  + (1. - rx) * qx * px * as011 + rx * (1. - qx)*(1. - px) * as100;
1208 
1209  b_coef[im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ai000 + rx * qx * px * ai111
1210  + rx * (1. - qx) * px * ai101 + (1. - rx) * qx * (1. - px) * ai010
1211  + rx * qx * (1. - px) * ai110 + (1. - rx)*(1. - qx) * px * ai001
1212  + (1. - rx) * qx * px * ai011 + rx * (1. - qx)*(1. - px) * ai100;
1213 
1214  c_coef[im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ac000 + rx * qx * px * ac111
1215  + rx * (1. - qx) * px * ac101 + (1. - rx) * qx * (1. - px) * ac010
1216  + rx * qx * (1. - px) * ac110 + (1. - rx)*(1. - qx) * px * ac001
1217  + (1. - rx) * qx * px * ac011 + rx * (1. - qx)*(1. - px) * ac100;
1218 
1219  d_coef[im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ad000 + rx * qx * px * ad111
1220  + rx * (1. - qx) * px * ad101 + (1. - rx) * qx * (1. - px) * ad010
1221  + rx * qx * (1. - px) * ad110 + (1. - rx)*(1. - qx) * px * ad001
1222  + (1. - rx) * qx * px * ad011 + rx * (1. - qx)*(1. - px) * ad100;
1223 
1224  e_coef[im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ae000 + rx * qx * px * ae111
1225  + rx * (1. - qx) * px * ae101 + (1. - rx) * qx * (1. - px) * ae010
1226  + rx * qx * (1. - px) * ae110 + (1. - rx)*(1. - qx) * px * ae001
1227  + (1. - rx) * qx * px * ae011 + rx * (1. - qx)*(1. - px) * ae100;
1228 
1229  }
1230  }
1231  }
1232 
1233  // root finding is quadratic, but table includes cubic term, make sure it's zero
1234  if (fabs(d_coef[modnum][iwnir_l]) > 1e-9) {
1235  printf("non zero cubic term found in longest NIR wavelength of aerosol table. Zia!!\n");
1236  exit(1);
1237  }
1238 
1239  /* return pointers to coeffs for this geometry */
1240  *a = &a_coef[modnum][0];
1241  *b = &b_coef[modnum][0];
1242  *c = &c_coef[modnum][0];
1243  *d = &d_coef[modnum][0];
1244  *e = &e_coef[modnum][0];
1245 
1246 
1247  return;
1248 }
1249 
1250 /* ---------------------------------------------------------------------------------------- */
1251 /* model_select_ahmad() - select two aerosol models whose epsilon values bracket the */
1252 /* the observed ms epsilon, eps_obs */
1253 /* */
1254 /* Z Ahmad July 2014. */
1255 
1256 /* ---------------------------------------------------------------------------------------- */
1257 int comp_epsilonT(epsilonTstr *x, epsilonTstr *y) {
1258  return (x->eps_obs < y->eps_obs ? -1 : 1);
1259 }
1260 
1261 void model_select_ahmad(int32_t nmodels, int32_t *mindx, float eps_pred[], float eps_obs, int32_t *modmin,
1262  int32_t *modmax, float *modrat) {
1263  static epsilonTstr epsilonT[MAXAERMOD];
1264 
1265  int im, im1, im2;
1266 
1267  // set-up table to keep model epsilon and model number pairs together
1268 
1269  for (im = 0; im < nmodels; im++) {
1270  epsilonT[im].eps_obs = eps_pred[im];
1271  epsilonT[im].modnum = im;
1272  /* printf("%d %7.4f %7.4f\n",im,eps_pred[im],eps_obs); */
1273  }
1274 
1275  // sort into ascending model epsilon order
1276 
1277  qsort(epsilonT, nmodels, sizeof (epsilonTstr),
1278  (int (*)(const void *, const void *)) comp_epsilonT);
1279 
1280  // find bounding epsilon indices in table
1281 
1282  for (im = 0; im < nmodels; im++) {
1283  if (eps_obs < epsilonT[im].eps_obs)
1284  break;
1285  }
1286 
1287  //if(im==0) //no lower bounding by M. Zhang
1288  //{
1289  // *modmin=-1;
1290  // return;
1291  //}
1292 
1293 
1294  im1 = MAX(MIN(im - 1, nmodels - 1), 0);
1295  im2 = MAX(MIN(im, nmodels - 1), 0);
1296 
1297 
1298  // convert table indices to model indices of the input order
1299  // compute model weighting
1300 
1301  *modmin = epsilonT[im1].modnum;
1302  *modmax = epsilonT[im2].modnum;
1303  *modrat = (eps_obs - epsilonT[im1].eps_obs) / (epsilonT[im2].eps_obs - epsilonT[im1].eps_obs);
1304  /* printf("%d %d %7.4f %7.4f %7.4f\n",im1,im2,eps_obs,epsilonT[im1].eps_obs,epsilonT[im2].eps_obs); */
1305 
1306  // If eps_obs is higer or lower than epsilon from table, then set the weight to 1
1307  if (*modmin == *modmax)
1308  *modrat = 1;
1309 
1310  return;
1311 }
1312 
1313 
1314 /*------------------------------------------------------------------------------------------ */
1315 /* comp_rhoa_ms_eps() - for a given model, computes the AOT and aerosol reflectance */
1316 /* */
1317 /* Z ahmad July2014 */
1318 
1319 /*------------------------------------------------------------------------------------------ */
1320 int comp_rhoa_ms_eps(int32_t nwave, float wave[], geom_str *geom,
1321  float tau_iwnir_l, int32_t modl, float tau_pred[], float rho_pred[])
1322  {
1323  float *ac, *bc, *cc, *dc, *ec;
1324  float ext_modl[nwave];
1325  float lg_tau_pred[nwave];
1326  float lg_rho_pred[nwave];
1327  int iw, iwtab;
1328 
1329  /* get the coefficients for lg_rho vs lg_aot */
1330 
1331  // Zia's function ---> something is wrong
1332  ms_eps_coef(modl, iwnir_l, wave, geom, &ac, &bc, &cc, &dc, &ec);
1333  // ms_eps_coef_cal(modl,nwave,solz,senz,phi,&ac,&bc,&cc,&dc,&ec);
1334 
1335 
1336  /* get the extinction coefficients and compute AOT at all wavelengths */
1337 
1338  for (iw = 0; iw < nwave; iw++) {
1339  iwtab = iwatab[iw];
1340  ext_modl[iw] = aertab->model[modl]->extc[iwtab];
1341  }
1342 
1343  /* printf("tau_pred[iw],tau_iwnir_l\n"); */
1344  for (iw = 0; iw < nwave; iw++) {
1345  tau_pred[iw] = (ext_modl[iw] / ext_modl[iwnir_l]) * tau_iwnir_l;
1346  lg_tau_pred[iw] = log(tau_pred[iw]);
1347  }
1348 
1349  /* compute rho_pred */
1350 
1351  for (iw = 0; iw < nwave; iw++) {
1352  iwtab = iwatab[iw];
1353  lg_rho_pred[iw] = ac[iwtab] +
1354  bc[iwtab] * lg_tau_pred[iw] +
1355  cc[iwtab] * pow(lg_tau_pred[iw], 2) +
1356  dc[iwtab] * pow(lg_tau_pred[iw], 3) +
1357  ec[iwtab] * pow(lg_tau_pred[iw], 4);
1358  rho_pred[iw] = exp(lg_rho_pred[iw]);
1359  }
1360 
1361 
1362  return (0);
1363 }
1364 
1365 
1366 /*------------------------------------------------------------------------------------------ */
1367 /* comp_rhoa_ms_eps() - for a given model, computes the AOT and aerosol reflectance */
1368 /* */
1369 /* Z ahmad July2014 */
1370 /* Modified by Amir Ibrahim January 2015 for linear space coefficients */
1371 
1372 /*------------------------------------------------------------------------------------------ */
1373 int comp_rhoa_ms_eps_lin(int32_t nwave, float wave[], geom_str *geom,
1374  float tau_iwnir_l, int32_t modl, float tau_pred[], float rho_pred[])
1375  {
1376  float *ac, *bc, *cc, *dc, *ec;
1377  float ext_modl[nwave];
1378  float lg_tau_pred[nwave];
1379  float lg_rho_pred[nwave];
1380  int iw, iwtab;
1381 
1382  /* get the coefficients for lg_rho vs lg_aot */
1383 
1384  // Zia's function ---> something is wrong
1385  ms_eps_coef(modl, iwnir_l, wave, geom, &ac, &bc, &cc, &dc, &ec);
1386  // ms_eps_coef_cal(modl,nwave,geom,&ac,&bc,&cc,&dc,&ec);
1387 
1388 
1389  /* get the extinction coefficients and compute AOT at all wavelengths */
1390 
1391  for (iw = 0; iw < nwave; iw++) {
1392  iwtab = iwatab[iw];
1393  ext_modl[iw] = aertab->model[modl]->extc[iwtab];
1394  }
1395 
1396  /* printf("tau_pred[iw],tau_iwnir_l\n"); */
1397  for (iw = 0; iw < nwave; iw++) {
1398  tau_pred[iw] = (ext_modl[iw] / ext_modl[iwnir_l]) * tau_iwnir_l;
1399  lg_tau_pred[iw] = (tau_pred[iw]);
1400  }
1401 
1402  /* compute rho_pred */
1403 
1404  for (iw = 0; iw < nwave; iw++) {
1405  iwtab = iwatab[iw];
1406  lg_rho_pred[iw] = ac[iwtab] +
1407  bc[iwtab] * lg_tau_pred[iw] +
1408  cc[iwtab] * pow(lg_tau_pred[iw], 2) +
1409  dc[iwtab] * pow(lg_tau_pred[iw], 3) +
1410  ec[iwtab] * pow(lg_tau_pred[iw], 4);
1411  rho_pred[iw] = (lg_rho_pred[iw]);
1412  }
1413 
1414 
1415  return (0);
1416 }
1417 
1418 
1419 /*----------------------------------------------------------------------------------------------*/
1420 /* spectral_matching() - calculates the best aerosol model that fits aerosol reflectance from */
1421 /* Red to SWIR channels. The spectral matching is done for 3 rh sets, and it assumes 1 band in */
1422 /* the NIR to be perfectly calibrated (i.e. 869 nm) */
1423 /* Amir Ibrahim - July 2016 */
1424 /*----------------------------------------------------------------------------------------------*/
1425 int spectral_matching(int32_t sensorID, float wave[], int32_t nwave,
1426  int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx1[],
1427  int32_t mindx2[], int32_t mindx3[], geom_str *geom,
1428  float wv, float rhoa[], float rho_aer[], int32_t *mod1_indx,
1429  int32_t *mod2_indx, float *weight, float tau_pred_min[],
1430  float tau_pred_max[], float tg_sol_sm[], float tg_sen_sm[],
1431  float Lt_sm[], int32_t ip) {
1432 
1433 
1434  float *ac, *bc, *cc, *dc, *ec;
1435  float ext_iwnir_cal;
1436  double ax, bx, cx, fx;
1437 
1438  int iw, i;
1439 
1440  int status = 0;
1441 
1442 
1443  float tau_iwnir_cal;
1444  float *lg_tau_all_wav = (float*) malloc(nwave * sizeof (float));
1445  float *lg_rho_all_wav_pred = (float*) malloc(nwave * sizeof (float));
1446  float **tau_all_wav = (float**) malloc(nwave * sizeof (float*));
1447  float **rho_all_wav_pred = (float**) malloc(nwave * sizeof (float*));
1448  ;
1449 
1450 
1451  for (i = 0; i < nwave; i++) {
1452  tau_all_wav[i] = (float*) malloc((nmodels) * sizeof (float));
1453  rho_all_wav_pred[i] = (float*) malloc((nmodels) * sizeof (float));
1454  }
1455 
1456  float *ext_all_wav = (float*) malloc(nwave * sizeof (float));
1457 
1458  int iwtab, im, modl;
1459  float lg_tau_iwnir_cal;
1460 
1461  double *diff = malloc(nwave * sizeof (double));
1462  double *chi = malloc(nmodels * sizeof (double));
1463  double *chi_temp = malloc(nmodels * sizeof (double));
1464 
1465  double diff_2;
1466  double min1, min2;
1467  int min1_indx, min2_indx;
1468 
1469  // float *SNR = malloc(nwave * sizeof(float));
1470  float *noise = malloc(nwave * sizeof (float));
1471  static float scaled_lt;
1472  int iwtab_l, iwtab_cal;
1473 
1474 
1475  if (sensorID == MODISA) {
1476  /* These are the noise coefficient for hmodisa from Erdem */
1477 
1478  float ltSnrFitCoef[16][2] = {
1479  {/*412:C0,C1*/0.05499859, 0.00008340},
1480  {/*443:*/0.02939470, 0.00009380},
1481  {/*469:*/0.11931482, 0.00008195},
1482  {/*488:*/0.01927545, 0.00009450},
1483  {/*531:*/0.01397522, 0.00010040},
1484  {/*547:*/0.01139088, 0.00016480},
1485  {/*555*/0.08769538, 0.00007000},
1486  {/*645:*/0.10406925, 0.00008533},
1487  {/*667L*/0.00496291, 0.00014050},
1488  {/*678L*/0.00427147, 0.00013160},
1489  {/*748*/0.00416994, 0.00021250},
1490  {/*859*/0.04055895, 0.000197550},
1491  {/*869*/0.00312263, 0.00018600},
1492  {/*1240*/0.07877732, 0.00049940},
1493  {/*1640*/100000.1, 100000000.1},
1494  {/*2130*/0.00628912, 0.00021160}
1495  };
1496 
1497  /* scaled noise based on MODISA spatial resolution */
1498  float snr_scale[16] = {
1499  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 2, 2, 2
1500  };
1501  /* Noise model */
1502  for (i = 0; i < nwave; i++) {
1503  scaled_lt = Lt_sm[ip * nwave + i];
1504  noise[i] = (ltSnrFitCoef[i][0] + ltSnrFitCoef[i][1] * scaled_lt * snr_scale[i]);
1505  //SNR[i] = scaled_lt * snr_mult / noise[i]; //noise model based on
1506  }
1507  } else {
1508  for (i = 0; i < nwave; i++)
1509  //SNR[i] = 1;
1510  noise[i] = 1;
1511  }
1512 
1513 
1514 
1515  int32_t *mindx = (int32_t*) malloc(nmodels * sizeof (int32_t));
1516  if (nmodels == 1) {
1517  *mindx = 0;
1518  } else {
1519  if (mindx1[0] == mindx2[0]) {
1520  nmodels = 10;
1521  for (i = 0; i < 10; i++)
1522  mindx[i] = mindx1[i];
1523  } else {
1524  if (nmodels < 30) {
1525  for (i = 0; i < 10; i++) {
1526  mindx[i] = mindx1[i];
1527  mindx[i + 10] = mindx2[i];
1528  }
1529  } else if (mindx3[0] >= mindx1[0]) {
1530  for (i = 0; i < 10; i++) {
1531  mindx[i] = mindx1[i];
1532  mindx[i + 10] = mindx2[i];
1533  mindx[i + 20] = mindx3[i];
1534  }
1535  } else {
1536  for (i = 0; i < 10; i++) {
1537  mindx[i] = mindx3[i];
1538  mindx[i + 10] = mindx1[i];
1539  mindx[i + 20] = mindx2[i];
1540  }
1541  }
1542  }
1543  }
1544 
1545  // compute MS epsilon for all nmodels. note that nmodels may be only a subset of
1546  // the full model suite. mindx maps from subset index to full index in suite
1547  /*
1548  float SNR[nwave];
1549  for (i=0;i<9;i++)
1550  SNR[i] = 0;
1551  SNR[10] = 1;
1552  SNR[11] = 0.3;
1553  SNR[12] = 1;
1554  SNR[13] = 0.5;
1555  SNR[14] = 0;
1556  SNR[15] = 0.2;
1557  */
1558  // printf("nmodels=%d\n",nmodels);
1559  for (im = 0; im < nmodels; im++) {
1560 
1561  modl = mindx[im];
1562  /* compute AOT at longest aerosol wavelength (iwnir_l) */
1563  int iwnir_cal;
1564 
1565  if (sensorID == AVIRIS) {
1566  iwnir_cal = 54; // this is hard coded (need to change) --> Amir
1567  iwtab_cal = iwatab[iwnir_cal];
1568  } else if (sensorID == MODISA) {
1569  iwnir_cal = 12;
1570  iwtab_cal = iwatab[iwnir_cal];
1571 
1572  } else {
1573  printf("Error: Sensor unkown %d\n", sensorID);
1574  return (1);
1575  }
1576  //printf("Using %f nm to estimate optical thickness\n",wave[iwnir_cal]);
1577 
1578  //ms_eps_coef_cal(modl,nwave,geom,&ac,&bc,&cc,&dc,&ec);
1579  ms_eps_coef(modl, iwnir_l, wave, geom, &ac, &bc, &cc, &dc, &ec);
1580 
1581  iwtab_l = iwatab[iwnir_l];
1582  // iwtab_s = iwatab[iwnir_s];
1583 
1584  ax = (double) ac[iwtab_cal] - log((double) rhoa[iwnir_cal]);
1585  bx = (double) bc[iwtab_cal];
1586  cx = (double) cc[iwtab_cal];
1587 
1588  fx = bx * bx - 4.0 * ax*cx;
1589  if (fx > 0.0 && cx != 0.0) {
1590  lg_tau_iwnir_cal = 0.5 * (-bx + sqrt(fx)) / cx;
1591  tau_iwnir_cal = exp(lg_tau_iwnir_cal);
1592  status = 0;
1593  } else {
1594  status = 1;
1595  break;
1596  }
1597 
1598  /* compute AOT at shortest aerosol wavelength (iwnir_s) */
1599 
1600  /* compute reflectance at all wavelength */
1601 
1602  ext_iwnir_cal = aertab->model[modl]->extc[iwtab_cal];
1603  for (iw = 0; iw <= iwtab_l - 1; iw++) {
1604  iwtab = iwatab[iw];
1605  ext_all_wav[iw] = aertab->model[modl]->extc[iwtab];
1606  tau_all_wav[iw][im] = (ext_all_wav[iw] / ext_iwnir_cal) * tau_iwnir_cal;
1607  lg_tau_all_wav[iw] = log(tau_all_wav[iw][im]);
1608  lg_rho_all_wav_pred[iw] = ac[iwtab] + bc[iwtab] * lg_tau_all_wav[iw] +
1609  cc[iwtab] * pow(lg_tau_all_wav[iw], 2) +
1610  dc[iwtab] * pow(lg_tau_all_wav[iw], 3) +
1611  ec[iwtab] * pow(lg_tau_all_wav[iw], 4);
1612  rho_all_wav_pred[iw][im] = exp(lg_rho_all_wav_pred[iw]);
1613 
1614 
1615  }
1616  for (iw = iwnir_s; iw <= iwnir_l; iw++) {
1617  // calculate the rms error
1618  diff[iw] = pow((rhoa[iw] - rho_all_wav_pred[iw][im]), 2);
1619  diff_2 += diff[iw]* (1 / noise[iw]) *(tg_sol_sm[ip * nwave + iw] * tg_sol_sm[ip * nwave + iw]);
1620 
1621  }
1622  chi[im] = sqrt(diff_2 / (iwnir_l - iwnir_s));
1623  diff_2 = 0.0;
1624  // printf("Model number = %d",im);
1625  }
1626 
1627  min1 = 1000;
1628  min2 = 1000;
1629  min1_indx = 0;
1630  min2_indx = 0;
1631 
1632  for (i = 0; i < nmodels; i++)
1633  chi_temp[i] = chi[i];
1634 
1635  qsort(chi_temp, nmodels, sizeof (chi_temp[0]), cmpfunc);
1636  min1 = chi_temp[0];
1637  min2 = chi_temp[1];
1638 
1639  for (i = 0; i < nmodels; i++) {
1640  if (chi[i] == min1)
1641  min1_indx = i;
1642  else if (chi[i] == min2)
1643  min2_indx = i;
1644  }
1645 
1646 
1647  *mod1_indx = 0;
1648  *mod2_indx = 0;
1649 
1650  if (min2 == 0)
1651  min2 = 0.001;
1652  // double weight1 = sqrt(pow(aertab->model[min2_indx]->extc[iwtab_cal] - aertab->model[min1_indx]->extc[iwtab_cal],2))/(aertab->model[min1_indx]->extc[iwtab_cal]);
1653  double weight1 = ((rhoa[iwnir_s] / rhoa[iwnir_l]) - (rho_all_wav_pred[iwnir_s][min1_indx] / rho_all_wav_pred[iwnir_l][min1_indx])) / ((rho_all_wav_pred[iwnir_s][min2_indx] / rho_all_wav_pred[iwnir_l][min2_indx])-(rho_all_wav_pred[iwnir_s][min1_indx] / rho_all_wav_pred[iwnir_l][min1_indx]));
1654  //double weight1 = (min2-min1)/min1 ;
1655  *weight = weight1;
1656  if (nmodels == 1) {
1657  *weight = 1;
1658  weight1 = 1;
1659  }
1660 
1661  *mod1_indx = mindx[min1_indx];
1662  *mod2_indx = mindx[min2_indx];
1663  for (iw = 0; iw <= iwnir_l; iw++) {
1664  rho_aer[iw] = (1 - weight1) * rho_all_wav_pred[iw][min1_indx] + (weight1) * rho_all_wav_pred[iw][min2_indx];
1665  tau_pred_min[iw] = tau_all_wav[iw][min1_indx];
1666  tau_pred_max[iw] = tau_all_wav[iw][min2_indx];
1667  }
1668 
1669  free(lg_tau_all_wav);
1670  free(lg_rho_all_wav_pred);
1671  for (i = 0; i < nwave; i++) {
1672  free(tau_all_wav[i]);
1673  free(rho_all_wav_pred[i]);
1674  }
1675  free(tau_all_wav);
1676  free(rho_all_wav_pred);
1677  free(ext_all_wav);
1678  free(chi_temp);
1679  free(chi);
1680  free(diff);
1681  free(mindx);
1682  // free(SNR);
1683  free(noise);
1684  //free(tg_sol_sm);
1685  // free(tg_sen_sm);
1686 
1687 
1688  return (status);
1689 }
1690 /*------------------------------------------------------------------------------------------*/
1691 /* ahmad_atm_corr() - compute aerosol reflectance at all wavelengths */
1692 /* */
1693 /* Z Ahmad. July 2014 */
1694 
1695 /*---------------------------------------------------------------------------------------- */
1696 int ahmad_atm_corr(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
1697  int32_t nmodels, int32_t mindx[],
1698  geom_str *geom, float wv, float rhoa[],
1699  int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir,
1700  float tau_pred_max[], float tau_pred_min[], float rho_pred_max[], float rho_pred_min[],
1701  float tau_aer[], float rho_aer[]) {
1702 
1703 
1704  float *ac, *bc, *cc, *dc, *ec;
1705  float ext_iwnir_l, ext_iwnir_s;
1706  double ax, bx, cx, fx;
1707  float tau_iwnir_l[nmodels];
1708  float lg_tau_iwnir_s[nmodels], tau_iwnir_s[nmodels];
1709  float lg_rho_iwnir_s_pred[nmodels], rho_iwnir_s_pred[nmodels];
1710  float eps_pred[nmodels];
1711  int im, modl;
1712  int im1, im2;
1713  float mwt;
1714  float lg_tau_iwnir_l;
1715  int iwtab_l, iwtab_s;
1716 
1717 
1718  static float eps_obs;
1719 
1720  int iw;
1721 
1722  int status = 0.0;
1723 
1724 
1725  /* compute the observed epsilon */
1726 
1727  eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
1728 
1729  /* printf("rho_869,rho_748,eps_obs\n"); */
1730  /* printf("%10.5f %10.5f %10.5f\n",rhoa[iwnir_l],rhoa[iwnir_s],eps_obs); */
1731 
1732 
1733  // compute MS epsilon for all nmodels. note that nmodels may be only a subset of
1734  // the full model suite. mindx maps from subset index to full index in suite
1735 
1736 
1737  for (im = 0; im < nmodels; im++) {
1738 
1739  modl = mindx[im]; // index in full model suite, needed for coefficient look-up
1740 
1741  /* compute AOT at longest aerosol wavelength (iwnir_l) */
1742 
1743  // Zia's function ---> something is wrong
1744  //ms_eps_coef(modl,iwnir_l,wave,solz,senz,phi,&ac,&bc,&cc,&dc,&ec);
1745 
1746  // ms_eps_coef_cal(modl,nwave,solz,senz,phi,&ac,&bc,&cc,&dc,&ec);
1747  ms_eps_coef(modl, iwnir_l, wave, geom, &ac, &bc, &cc, &dc, &ec);
1748  iwtab_l = iwatab[iwnir_l];
1749  iwtab_s = iwatab[iwnir_s];
1750 
1751  ax = (double) ac[iwtab_l] - log((double) rhoa[iwnir_l]);
1752  bx = (double) bc[iwtab_l];
1753  cx = (double) cc[iwtab_l];
1754 
1755  fx = bx * bx - 4.0 * ax*cx;
1756  if (fx > 0.0 && cx != 0.0) {
1757  lg_tau_iwnir_l = 0.5 * (-bx + sqrt(fx)) / cx;
1758  tau_iwnir_l[im] = exp(lg_tau_iwnir_l);
1759  } else {
1760  status = 1;
1761  break; // TODO: need to think about it
1762  }
1763 
1764  /* compute AOT at shortest aerosol wavelength (iwnir_s) */
1765 
1766  ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
1767  ext_iwnir_s = aertab->model[modl]->extc[iwtab_s];
1768 
1769  tau_iwnir_s[im] = (ext_iwnir_s / ext_iwnir_l) * tau_iwnir_l[im];
1770  lg_tau_iwnir_s[im] = log(tau_iwnir_s[im]);
1771 
1772 
1773  /* compute reflectance at (iwnir_s) */
1774 
1775  lg_rho_iwnir_s_pred[im] = ac[iwtab_s] + bc[iwtab_s] * lg_tau_iwnir_s[im] +
1776  cc[iwtab_s] * pow(lg_tau_iwnir_s[im], 2) +
1777  dc[iwtab_s] * pow(lg_tau_iwnir_s[im], 3) +
1778  ec[iwtab_s] * pow(lg_tau_iwnir_s[im], 4);
1779 
1780  rho_iwnir_s_pred[im] = exp(lg_rho_iwnir_s_pred[im]);
1781 
1782  /* compute model epsilon */
1783 
1784  eps_pred[im] = rho_iwnir_s_pred[im] / rhoa[iwnir_l];
1785 
1786  }
1787 
1788 
1789  *epsnir = eps_obs;
1790 
1791 
1792  /* now, find the bounding models, but skip this if we only have one */
1793  /* model (as when vicariously calibrating) */
1794 
1795  if (nmodels > 1) {
1796 
1797  /* locate two model_epsilons that bracket the observed epsilon */
1798  /* this will return the model numbers for the subset of models */
1799 
1800  model_select_ahmad(nmodels, mindx, eps_pred, eps_obs, &im1, &im2, &mwt);
1801 
1802  } else {
1803 
1804  /* single model case (allows fixed model by limiting model suite) */
1805 
1806  im1 = 0;
1807  im2 = 0;
1808  mwt = 0.0;
1809  }
1810 
1811  /* map selection to full suite for subsequent processing and return */
1812 
1813  //if(im1<0) // no lower bounding, M. Zhang
1814  //{
1815  // return 1;
1816  // status=1;
1817  //}
1818 
1819 
1820  *modmin = mindx[im1];
1821  *modmax = mindx[im2];
1822  *modrat = mwt;
1823 
1824  /* compute tau_pred and rho_predicted */
1825 
1826  comp_rhoa_ms_eps(nwave, wave, geom, tau_iwnir_l[im1], *modmin, tau_pred_min, rho_pred_min);
1827  comp_rhoa_ms_eps(nwave, wave, geom, tau_iwnir_l[im2], *modmax, tau_pred_max, rho_pred_max);
1828 
1829  /* compute weighted tau_aer and rho_aer */
1830 
1831  for (iw = 0; iw < nwave; iw++) {
1832  tau_aer[iw] = (1.0 - mwt) * tau_pred_min[iw] + mwt * tau_pred_max[iw];
1833  rho_aer[iw] = (1.0 - mwt) * rho_pred_min[iw] + mwt * rho_pred_max[iw];
1834  }
1835 
1836 
1837  return (status);
1838 }
1839 
1840 /*------------------------------------------------------------------------------------------*/
1841 /* ahmad_atm_corr_lin() - compute aerosol reflectance at all wavelengths in linear space */
1842 /* */
1843 /* Z Ahmad. July 2014 */
1844 /* Modified to linear spaced coefficient by Amir Ibrahim January 2017 */
1845 
1846 /*---------------------------------------------------------------------------------------- */
1847 int ahmad_atm_corr_lin(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
1848  int32_t nmodels, int32_t mindx[],
1849  geom_str *geom, float wv, float rhoa[],
1850  int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir,
1851  float tau_pred_max[], float tau_pred_min[], float rho_pred_max[], float rho_pred_min[],
1852  float tau_aer[], float rho_aer[]) {
1853 
1854 
1855  float *ac, *bc, *cc, *dc, *ec;
1856  float ext_iwnir_l, ext_iwnir_s;
1857  double ax, bx, cx, fx;
1858  float tau_iwnir_l[nmodels];
1859  float lg_tau_iwnir_s[nmodels], tau_iwnir_s[nmodels];
1860  float lg_rho_iwnir_s_pred[nmodels], rho_iwnir_s_pred[nmodels];
1861  float eps_pred[nmodels];
1862  int im, modl;
1863  int im1, im2;
1864  float mwt;
1865  float lg_tau_iwnir_l;
1866  int iwtab_l, iwtab_s;
1867 
1868 
1869  static float eps_obs;
1870 
1871  int iw;
1872 
1873  int status = 0.0;
1874 
1875 
1876  /* compute the observed epsilon */
1877 
1878  eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
1879 
1880  /* printf("rho_869,rho_748,eps_obs\n"); */
1881  /* printf("%10.5f %10.5f %10.5f\n",rhoa[iwnir_l],rhoa[iwnir_s],eps_obs); */
1882 
1883 
1884  // compute MS epsilon for all nmodels. note that nmodels may be only a subset of
1885  // the full model suite. mindx maps from subset index to full index in suite
1886 
1887 
1888  for (im = 0; im < nmodels; im++) {
1889 
1890  modl = mindx[im]; // index in full model suite, needed for coefficient look-up
1891 
1892  /* compute AOT at longest aerosol wavelength (iwnir_l) */
1893 
1894  // Zia's function ---> something is wrong
1895  //ms_eps_coef(modl,iwnir_l,wave,solz,senz,phi,&ac,&bc,&cc,&dc,&ec);
1896 
1897  // ms_eps_coef_cal(modl,nwave,solz,senz,phi,&ac,&bc,&cc,&dc,&ec);
1898  ms_eps_coef(modl, iwnir_l, wave, geom, &ac, &bc, &cc, &dc, &ec);
1899  iwtab_l = iwatab[iwnir_l];
1900  iwtab_s = iwatab[iwnir_s];
1901 
1902  ax = (double) ac[iwtab_l]-(double) rhoa[iwnir_l];
1903  bx = (double) bc[iwtab_l];
1904  cx = (double) cc[iwtab_l];
1905 
1906  fx = bx * bx - 4.0 * ax*cx;
1907  if (fx > 0.0 && cx != 0.0) {
1908  lg_tau_iwnir_l = 0.5 * (-bx + sqrt(fx)) / cx;
1909  tau_iwnir_l[im] = (lg_tau_iwnir_l);
1910  } else {
1911  status = 1;
1912  break;
1913  }
1914 
1915  /* compute AOT at shortest aerosol wavelength (iwnir_s) */
1916 
1917  ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
1918  ext_iwnir_s = aertab->model[modl]->extc[iwtab_s];
1919 
1920  tau_iwnir_s[im] = (ext_iwnir_s / ext_iwnir_l) * tau_iwnir_l[im];
1921  lg_tau_iwnir_s[im] = (tau_iwnir_s[im]);
1922 
1923  if (ax > 1e5) {
1924  printf("\nErroneous aerosol option, %d\n", aer_opt);
1925  printf("You are using a log-space LUT. The aerosol LUT coefficients need to be in linear-space. Use aer_opt=-18 instead. \n");
1926  exit(FATAL_ERROR);
1927  }
1928 
1929 
1930  /* compute reflectance at (iwnir_s) */
1931 
1932  lg_rho_iwnir_s_pred[im] = ac[iwtab_s] + bc[iwtab_s] * lg_tau_iwnir_s[im] +
1933  cc[iwtab_s] * pow(lg_tau_iwnir_s[im], 2) +
1934  dc[iwtab_s] * pow(lg_tau_iwnir_s[im], 3) +
1935  ec[iwtab_s] * pow(lg_tau_iwnir_s[im], 4);
1936 
1937  rho_iwnir_s_pred[im] = (lg_rho_iwnir_s_pred[im]);
1938 
1939  /* compute model epsilon */
1940 
1941  eps_pred[im] = rho_iwnir_s_pred[im] / rhoa[iwnir_l];
1942 
1943  }
1944 
1945 
1946  *epsnir = eps_obs;
1947 
1948 
1949  /* now, find the bounding models, but skip this if we only have one */
1950  /* model (as when vicariously calibrating) */
1951 
1952  if (nmodels > 1) {
1953 
1954  /* locate two model_epsilons that bracket the observed epsilon */
1955  /* this will return the model numbers for the subset of models */
1956 
1957  model_select_ahmad(nmodels, mindx, eps_pred, eps_obs, &im1, &im2, &mwt);
1958 
1959  } else {
1960 
1961  /* single model case (allows fixed model by limiting model suite) */
1962 
1963  im1 = 0;
1964  im2 = 0;
1965  mwt = 0.0;
1966  }
1967 
1968  /* map selection to full suite for subsequent processing and return */
1969 
1970  *modmin = mindx[im1];
1971  *modmax = mindx[im2];
1972  *modrat = mwt;
1973 
1974  /* compute tau_pred and rho_predicted */
1975 
1976  comp_rhoa_ms_eps_lin(nwave, wave, geom, tau_iwnir_l[im1], *modmin, tau_pred_min, rho_pred_min);
1977  comp_rhoa_ms_eps_lin(nwave, wave, geom, tau_iwnir_l[im2], *modmax, tau_pred_max, rho_pred_max);
1978 
1979  /* compute weighted tau_aer and rho_aer */
1980 
1981  for (iw = 0; iw < nwave; iw++) {
1982  tau_aer[iw] = (1.0 - mwt) * tau_pred_min[iw] + mwt * tau_pred_max[iw];
1983  rho_aer[iw] = (1.0 - mwt) * rho_pred_min[iw] + mwt * rho_pred_max[iw];
1984  }
1985 
1986 
1987  return (status);
1988 }
1989 
1990 /* ---------------------------------------------------------------------------------------- */
1991 /* model_phase() - return phase function at model wavelengths at the input geometry. */
1992 /* */
1993 /* This is effectively a C version of load_ss.f by M. Wang. The program optimizes for */
1994 /* multiple calls at the same geometry by computing for all models on the first call with a */
1995 /* new geometry. It returns a pointer to the internal static array for the requested model.*/
1996 /* */
1997 /* B. Franz, 1 June 2004. */
1998 /* W. Robinson, SAIC modified to use band-dependent viewing geometry, if avail */
1999 
2000 /* ---------------------------------------------------------------------------------------- */
2001 float *model_phase(int modnum, geom_str *geom) {
2002  static float nw = 1.334;
2003 
2004  static int computed[MAXMODEL];
2005  static float lastsolz = -999.;
2006  static float lastsenz = -999.;
2007  static float lastphi = -999.;
2008  static float *phase[MAXMODEL];
2009  static float *fres1, *fres2;
2010  static float *scatt1, *scatt2;
2011  static float scatt1_cnst, scatt2_cnst, fres1_cnst, fres2_cnst;
2012  static float *scatt1_ar, *scatt2_ar, *fres1_ar, *fres2_ar;
2013  static int firstCall = 1, gmult = 1;
2014 
2015  float phase1, phase2;
2016  int im, iw, ig;
2017 
2018  if (firstCall == 1) {
2019  firstCall = 0;
2020  for (im = 0; im < MAXMODEL; im++) {
2021  if ((phase[im] = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
2022  printf("Unable to allocate space for phase.\n");
2023  exit(1);
2024  }
2025  }
2026  /* set up geometry based scattering and fresnel */
2027  if ((geom->gmult == 0) || (interpol == 1)) {
2028  gmult = 0;
2029  scatt1 = &scatt1_cnst;
2030  scatt2 = &scatt2_cnst;
2031  fres1 = &fres1_cnst;
2032  fres2 = &fres2_cnst;
2033  } else {
2034  gmult = 1;
2035  if ((scatt1_ar =
2036  (float *) malloc(aertab->nwave * sizeof (float))) == NULL) {
2037  printf("Unable to allocate space for scatt1.\n");
2038  exit(1);
2039  }
2040  if ((scatt2_ar =
2041  (float *) malloc(aertab->nwave * sizeof (float))) == NULL) {
2042  printf("Unable to allocate space for scatt2.\n");
2043  exit(1);
2044  }
2045  if ((fres1_ar =
2046  (float *) malloc(aertab->nwave * sizeof (float))) == NULL) {
2047  printf("Unable to allocate space for fres1.\n");
2048  exit(1);
2049  }
2050  if ((fres2_ar =
2051  (float *) malloc(aertab->nwave * sizeof (float))) == NULL) {
2052  printf("Unable to allocate space for fres2.\n");
2053  exit(1);
2054  }
2055  scatt1 = scatt1_ar;
2056  scatt2 = scatt2_ar;
2057  fres1 = fres1_ar;
2058  fres2 = fres2_ar;
2059  }
2060  }
2061 
2062  /* recalculate only if geometry changes */
2063 
2064  if ((geom->solz[0] != lastsolz) || (geom->senz[0] != lastsenz) ||
2065  (geom->phi[0] != lastphi)) {
2066 
2067  float temp;
2068 
2069  for (im = 0; im < aertab->nmodel; im++)
2070  computed[im] = 0;
2071 
2072  /* determine scattering angles (direct and surface reflected) */
2073  for (iw = 0; iw < aertab->nwave; iw++) {
2074  ig = gmult * iw;
2075  temp = sqrt((1.0 - geom->csenz[ig] * geom->csenz[ig]) *
2076  (1.0 - geom->csolz[ig] * geom->csolz[ig])) *
2077  cos(geom->phi[ig] / radeg);
2078  scatt1[iw] = acos(
2079  MAX(-geom->csenz[ig] * geom->csolz[ig] + temp, -1.0)) * radeg;
2080  scatt2[iw] = acos(
2081  MIN(geom->csenz[ig] * geom->csolz[ig] + temp, 1.0)) * radeg;
2082 
2083  /* compute Fresnel coefficients */
2084  fres1[iw] = fresnel_coef(geom->csenz[ig], nw);
2085  fres2[iw] = fresnel_coef(geom->csolz[ig], nw);
2086  if (gmult == 0) break;
2087  }
2088 
2089  lastsolz = geom->solz[0];
2090  lastsenz = geom->senz[0];
2091  lastphi = geom->phi[0];
2092  }
2093 
2094  if (!computed[modnum]) {
2095 
2096  im = modnum;
2097  computed[modnum] = 1;
2098 
2099  /* compute phase function for this geometry, all models */
2100  for (iw = 0; iw < aertab->nwave; iw++) {
2101  ig = gmult * iw;
2102  splint(aertab->scatt,
2103  &aertab->model[im]->lnphase[iw][0],
2104  &aertab->model[im]->d2phase[iw][0],
2105  aertab->nscatt, scatt1[ig], &phase1);
2106  splint(aertab->scatt,
2107  &aertab->model[im]->lnphase[iw][0],
2108  &aertab->model[im]->d2phase[iw][0],
2109  aertab->nscatt, scatt2[ig], &phase2);
2110  // incident diffuse reflected diff dir
2111  phase[im][iw] = exp(phase1) +
2112  exp(phase2)*(fres1[ig] + fres2[ig]);
2113  }
2114  }
2115 
2116  return (&phase[modnum][0]);
2117 }
2118 
2119 
2120 /* ---------------------------------------------------------------------------------------- */
2121 /* aeroob_cf() - out-of-band water-vapor scale factor */
2122 
2123 /* ---------------------------------------------------------------------------------------- */
2124 float aeroob_cf(int modnum, geom_str *geom) {
2125  static int firstCall = 1;
2126  static int iw1;
2127  static int iw2;
2128 
2129  float *phase;
2130  float rhoas1, rhoas2;
2131  float eps;
2132  float cf;
2133 
2134  if (firstCall) {
2135  iw1 = windex(765, aertab->wave, aertab->nwave);
2136  iw2 = windex(865, aertab->wave, aertab->nwave);
2137  if (iw1 == iw2) iw1--;
2138  firstCall = 0;
2139  }
2140 
2141  phase = model_phase(modnum, geom);
2142  rhoas1 = aertab->model[modnum]->albedo[iw1] * phase[iw1] * aertab->model[modnum]->extc[iw1];
2143  rhoas2 = aertab->model[modnum]->albedo[iw2] * phase[iw2] * aertab->model[modnum]->extc[iw2];
2144  eps = rhoas1 / rhoas2;
2145  cf = log(eps) / (aertab->wave[iw2] - aertab->wave[iw1]);
2146 
2147  return (cf);
2148 }
2149 
2150 
2151 /* ---------------------------------------------------------------------------------------- */
2152 /* aeroob - out-of-band water-vapor correction */
2153 
2154 /* ---------------------------------------------------------------------------------------- */
2155 float aeroob(int32_t sensorID, int32_t iw, float airmass, float cf, float wv) {
2156  static float *a01;
2157  static float *a02;
2158  static float *a03;
2159  static float *a04;
2160  static float *a05;
2161  static float *a06;
2162  static float *a07;
2163  static float *a08;
2164  static float *a09;
2165  static float *a10;
2166  static float *a11;
2167  static float *a12;
2168 
2169  static int firstCall = 1;
2170 
2171  float f;
2172  f = 1.0;
2173  if ((input->gas_opt & GAS_TRANS_TBL_BIT) == 0) {
2174  if (firstCall) {
2175  firstCall = 0;
2176  rdsensorinfo(sensorID, evalmask, "oobwv01", (void **) &a01); /* coeff #1 per sensor wave */
2177  rdsensorinfo(sensorID, evalmask, "oobwv02", (void **) &a02);
2178  rdsensorinfo(sensorID, evalmask, "oobwv03", (void **) &a03);
2179  rdsensorinfo(sensorID, evalmask, "oobwv04", (void **) &a04);
2180  rdsensorinfo(sensorID, evalmask, "oobwv05", (void **) &a05);
2181  rdsensorinfo(sensorID, evalmask, "oobwv06", (void **) &a06);
2182  rdsensorinfo(sensorID, evalmask, "oobwv07", (void **) &a07);
2183  rdsensorinfo(sensorID, evalmask, "oobwv08", (void **) &a08);
2184  rdsensorinfo(sensorID, evalmask, "oobwv09", (void **) &a09);
2185  rdsensorinfo(sensorID, evalmask, "oobwv10", (void **) &a10);
2186  rdsensorinfo(sensorID, evalmask, "oobwv11", (void **) &a11);
2187  rdsensorinfo(sensorID, evalmask, "oobwv12", (void **) &a12);
2188  printf("\nLoading water-vapor correction coefficients.\n");
2189  }
2190 
2191  f = (a01[iw] + a02[iw] * airmass + cf * (a03[iw] + a04[iw] * airmass))
2192  + (a05[iw] + a06[iw] * airmass + cf * (a07[iw] + a08[iw] * airmass)) * wv
2193  + (a09[iw] + a10[iw] * airmass + cf * (a11[iw] + a12[iw] * airmass)) * wv*wv;
2194  }
2195  return (f);
2196 }
2197 
2198 
2199 /* ---------------------------------------------------------------------------------------- */
2200 /* rhoa_to_rhoas() - MS aerosol reflectance to SS aerosol reflectance */
2201 /* */
2202 /* B. Franz, 1 June 2004. */
2203 
2204 /* ---------------------------------------------------------------------------------------- */
2205 int rhoa_to_rhoas(int32_t sensorID, int modnum, geom_str *geom, float wv,
2206  float rhoa[], float wave[], int32_t nwave, int iw1, int iw2, float rhoas[]) {
2207  float *ac, *bc, *cc;
2208  double a, b, c;
2209  double f;
2210  int iw, ig;
2211  int iwtab;
2212  float cf;
2213  int status = 0;
2214 
2215  ss_to_ms_coef(modnum, geom, &ac, &bc, &cc);
2216 
2217  cf = aeroob_cf(modnum, geom);
2218 
2219  for (iw = iw1; iw <= iw2; iw++) {
2220  ig = iw * geom->gmult;
2221  if (rhoa[iw] < 1.e-20)
2222  rhoas[iw] = rhoa[iw];
2223  else {
2224  iwtab = iwatab[iw];
2225  a = (double) ac[iwtab];
2226  b = (double) bc[iwtab];
2227  c = (double) cc[iwtab];
2228  f = b * b - 4 * c * (a - log((double) rhoa[iw]));
2229  if (f > 0.00001) { // this was 0.0, but small values caused segfault (BAF, 9/2014)
2230  if (fabs(c) > 1.e-20) {
2231  rhoas[iw] = exp(0.5 * (-b + sqrt(f)) / c);
2232  } else if (fabs(a) > 1.e-20 && fabs(b) > 1.e-20) {
2233  rhoas[iw] = pow(rhoa[iw] / a, 1. / b);
2234  } else {
2235  status = 1;
2236  break;
2237  }
2238  rhoas[iw] = rhoas[iw] / aeroob(sensorID, iw, geom->airmass[ig], cf, wv);
2239  if (!isfinite(rhoas[iw]) || rhoas[iw] < 1.e-20) {
2240  status = 1;
2241  break;
2242  }
2243  } else {
2244  status = 1;
2245  break;
2246  }
2247  }
2248  }
2249 
2250  // return input values and failure status if any wavelengths failed
2251 
2252  if (status != 0) {
2253  for (iw = iw1; iw <= iw2; iw++)
2254  rhoas[iw] = rhoa[iw];
2255  }
2256 
2257  return (status);
2258 }
2259 
2260 
2261 
2262 /* ---------------------------------------------------------------------------------------- */
2263 /* rhoas_to_rhoa() - SS aerosol reflectance to MS aerosol reflectance */
2264 /* */
2265 /* B. Franz, 1 June 2004. */
2266 
2267 /* ---------------------------------------------------------------------------------------- */
2268 void rhoas_to_rhoa(int32_t sensorID, int modnum, geom_str *geom, float wv,
2269  float rhoas[], float wave[], int32_t nwave, int iw1, int iw2, float rhoa[]) {
2270  float *ac, *bc, *cc;
2271  float a, b, c;
2272  float lnrhoas;
2273  int iw, ig;
2274  int iwtab;
2275  float cf;
2276 
2277  ss_to_ms_coef(modnum, geom, &ac, &bc, &cc);
2278 
2279  cf = aeroob_cf(modnum, geom);
2280 
2281  for (iw = iw1; iw <= iw2; iw++) {
2282  ig = iw * geom->gmult;
2283 
2284  /* these changes ensure that rhoa1 -> rhoas -> rhoa2 == rhoa1 */
2285  /* but, that changes everything, slightly (tau) */
2286 
2287  if (rhoas[iw] < 1.e-20)
2288  rhoa[iw] = rhoas[iw];
2289  else {
2290  iwtab = iwatab[iw];
2291  a = ac[iwtab];
2292  b = bc[iwtab];
2293  c = cc[iwtab];
2294  lnrhoas = log(rhoas[iw] * aeroob(sensorID, iw, geom->airmass[ig], cf, wv));
2295  rhoa[iw] = exp(a + b * lnrhoas + c * lnrhoas * lnrhoas);
2296  }
2297 
2298  /*
2299  iwtab = iwatab[iw];
2300  a = ac[iwtab];
2301  b = bc[iwtab];
2302  c = cc[iwtab];
2303  lnrhoas = log(rhoas[iw]);
2304  rhoa[iw] = exp(a + b*lnrhoas + c*lnrhoas*lnrhoas)
2305  * aeroob(sensorID,iw, geom->airmass[ig],cf,wv);
2306  */
2307  }
2308 
2309  return;
2310 }
2311 
2312 
2313 /* ---------------------------------------------------------------------------------------- */
2314 /* model_epsilon() - return model epsilon at input wavelengths for input geometry. */
2315 /* */
2316 /* If the input wavelengths are not equal to the model wavelengths, the model epsilon will */
2317 /* be interpolated to the input wavelengths. It is assumed that the longest (last) input */
2318 /* wavelength is equivalent to the longest wavelength of the model table. Hence, */
2319 /* the function should always be called with the full array of input sensor wavelengths. */
2320 /* */
2321 /* The program optimizes for multiple calls at the same geometry by computing for all */
2322 /* models on the first call with a new geometry. It returns a pointer to the internal */
2323 /* static arrays of epsilon for the requested model. . */
2324 /* */
2325 /* B. Franz, 1 June 2004. */
2326 
2327 /* ---------------------------------------------------------------------------------------- */
2328 float *model_epsilon(int modnum, int32_t iwnir_l, float wave[], int32_t nwave, geom_str *geom) {
2329  static float lastsolz = -999.;
2330  static float lastsenz = -999.;
2331  static float lastphi = -999.;
2332  static int32_t lastiwl = -999;
2333  static int lastmod = -999;
2334  static float *epsilon[MAXMODEL];
2335  static int firstCall = 1;
2336  int i;
2337  float maxwave;
2338  /* recalculate only if geometry changes */
2339 
2340  if (firstCall == 1) {
2341  firstCall = 0;
2342  maxwave = MAX(aertab->nwave, nwave);
2343  for (i = 0; i < MAXMODEL; i++) {
2344  if ((epsilon[i] = (float *) calloc(maxwave, sizeof (float))) == NULL) {
2345  printf("Unable to allocate space for epsilon.\n");
2346  exit(1);
2347  }
2348 
2349  }
2350  }
2351 
2352  /* recalculate only if geometry changes */
2353 
2354  if (modnum != lastmod || geom->solz[0] != lastsolz ||
2355  geom->senz[0] != lastsenz || geom->phi[0] != lastphi ||
2356  iwnir_l != lastiwl) {
2357 
2358  int iwnir = iwatab[iwnir_l];
2359  float *phase;
2360  float *lneps;
2361  float rhoas1, rhoas2;
2362  int im, iw, iwtab;
2363 
2364  if ((lneps = (float *) calloc(aertab->nwave, sizeof (float))) == NULL) {
2365  printf("Unable to allocate space for lneps.\n");
2366  exit(1);
2367  }
2368 
2369  im = modnum;
2370  phase = model_phase(im, geom);
2371  for (iw = 0; iw < aertab->nwave; iw++) {
2372  rhoas1 = aertab->model[im]->albedo[iw] * phase[iw] * aertab->model[im]->extc[iw];
2373  rhoas2 = aertab->model[im]->albedo[iwnir] * phase[iwnir] * aertab->model[im]->extc[iwnir];
2374  epsilon[im][iw] = rhoas1 / rhoas2;
2375  if (interpol)
2376  lneps[iw] = log(epsilon[im][iw]);
2377  }
2378  if (interpol) {
2379  for (iw = 0; iw < nwave; iw++) {
2380  iwtab = iwatab[iw];
2381  if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
2382  epsilon[im][iw] = exp(linterp(aertab->wave, lneps, aertab->nwave, wave[iw]));
2383  else
2384  epsilon[im][iw] = exp(lneps[iwtab]);
2385  }
2386  }
2387 
2388  lastsolz = geom->solz[0];
2389  lastsenz = geom->senz[0];
2390  lastphi = geom->phi[0];
2391  lastiwl = iwnir_l;
2392  lastmod = modnum;
2393  free(lneps);
2394 
2395  }
2396 
2397  return (&epsilon[modnum][0]);
2398 }
2399 
2400 
2401 /* ---------------------------------------------------------------------------------------- */
2402 /* model_select_wang() - M. Wang aerosol model selection process . */
2403 /* */
2404 /* B. Franz, 1 June 2004. */
2405 
2406 /* ---------------------------------------------------------------------------------------- */
2407 int model_select_wang(int32_t sensorID, float wave[], int32_t nwave,
2408  int32_t nmodel, int32_t mindx[], geom_str *geom, float wv,
2409  float rhoa[], int32_t iwnir_s, int32_t iwnir_l, int32_t *modmin,
2410  int32_t *modmax, float *modrat, float *epsnir) {
2411  float *rhoas;
2412  float eps_ret [MAXMODEL];
2413  float eps_mod [MAXMODEL];
2414  float eps_err [MAXMODEL];
2415  int imod [MAXMODEL];
2416  int nmod = nmodel;
2417  float eps_ave;
2418  float *eps;
2419  float err_m;
2420  float err_p;
2421  int jm, im, iim;
2422  int eps_flg = 0;
2423  float wt;
2424  float tot_err;
2425  int itmp;
2426 
2427  *modmin = -1;
2428  *modmax = -1;
2429  *modrat = 0.0;
2430  *epsnir = 0.0;
2431 
2432  if ((rhoas = (float *) calloc(nwave, sizeof (float))) == NULL) {
2433  printf("Unable to allocate space for rhoas.\n");
2434  exit(1);
2435  }
2436 
2437  /* get theoretical and retrieved epsilon for each model, and */
2438  /* compute average retrieved epsilon */
2439 
2440  eps_ave = 0.0;
2441  for (jm = 0; jm < nmod; jm++) {
2442 
2443  im = mindx[jm];
2444 
2445  /* retrieve epsilon in NIR assuming each model */
2446  rhoa_to_rhoas(sensorID, im, geom, wv, rhoa, wave, nwave, iwnir_s, iwnir_l, rhoas);
2447 
2448  if (rhoas[iwnir_l] > 0.0000001)
2449  eps_ret[jm] = rhoas[iwnir_s] / rhoas[iwnir_l];
2450  else
2451  eps_ret[jm] = 0;
2452 
2453  /* get model epsilon for each model at this geometry */
2454  eps = model_epsilon(im, iwnir_l, wave, nwave, geom);
2455  eps_mod[jm] = eps[iwnir_s];
2456 
2457  eps_ave += eps_ret[jm];
2458  }
2459  if (isfinite(eps_ave))
2460  eps_ave /= nmod;
2461  else
2462  eps_ave = 1.0;
2463 
2464 
2465  /* determine average retrieved epsilon for the four retrievals which most */
2466  /* closely match the theoretical model epsilons. the model set is reduced */
2467  /* from the full suite to the final four by iterative outlier rejection. */
2468 
2469  while (nmod > 4) {
2470 
2471  /* compute differences between retrieved and model epsilon */
2472  for (im = 0; im < nmodel; im++) {
2473  imod[im] = im;
2474  eps_err[im] = eps_ave - eps_mod[im];
2475  }
2476 
2477  /* sort model indices by smallest to largest absolute differences */
2478  for (im = 0; im < nmodel - 1; im++) {
2479  for (iim = im + 1; iim < nmodel; iim++)
2480  if (fabs(eps_err[imod[im]]) > fabs(eps_err[imod[iim]])) {
2481  itmp = imod[im ];
2482  imod[im ] = imod[iim];
2483  imod[iim] = itmp;
2484  }
2485  }
2486 
2487  /* recompute average retrieved epsilon over the n-2 closest models */
2488  /* averaging is done as a weighted mean with wt=1-|eps_err|/tot_err */
2489  /* note that the sum of the weights is equal to n-1 */
2490 
2491  nmod = nmod - 2;
2492 
2493  tot_err = 0.0;
2494  for (iim = 0; iim < nmod; iim++) {
2495  im = imod[iim];
2496  tot_err += fabs(eps_err[im]);
2497  }
2498 
2499  eps_ave = 0.0;
2500  for (iim = 0; iim < nmod; iim++) {
2501  im = imod[iim];
2502  wt = 1.0 - fabs(eps_err[im]) / tot_err;
2503  eps_ave += eps_ret[im] * wt;
2504  }
2505  eps_ave /= (nmod - 1);
2506  }
2507 
2508  /* now select the two models which bracket eps_ave */
2509  err_m = 0 - FLT_MAX;
2510  err_p = FLT_MAX;
2511  for (im = 0; im < nmodel; im++) {
2512  eps_err[im] = eps_ave - eps_mod[im];
2513  if (eps_err[im] >= 0.0) {
2514  if (eps_err[im] < err_p) {
2515  err_p = eps_err[im];
2516  *modmin = im;
2517  }
2518  } else {
2519  if (eps_err[im] > err_m) {
2520  err_m = eps_err[im];
2521  *modmax = im;
2522  }
2523  }
2524  }
2525 
2526  /* compute interpolation ratio */
2527  if (*modmin < 0) {
2528  /* no lower-bounding model */
2529  *modmin = *modmax;
2530  *modrat = 0.0;
2531  eps_flg = -1;
2532  } else if (*modmax < 0) {
2533  /* no upper-bounding model */
2534  *modmax = *modmin;
2535  *modrat = 0.0;
2536  eps_flg = 1;
2537  } else
2538  *modrat = (eps_ave - eps_mod[*modmin]) / (eps_mod[*modmax] - eps_mod[*modmin]);
2539 
2540  *modmin = mindx[*modmin];
2541  *modmax = mindx[*modmax];
2542 
2543  /* return retrieved epsilon */
2544  *epsnir = eps_ave;
2545 
2546  free(rhoas);
2547 
2548  return (eps_flg);
2549 }
2550 
2551 
2552 /* ---------------------------------------------------------------------------------------- */
2553 /* model_select_angst() - Select model pair based on input Angstrom coefficient */
2554 /* */
2555 /* B. Franz, 1 August 2004. */
2556 /* ---------------------------------------------------------------------------------------- */
2557 int compalphaT(alphaTstr *x, alphaTstr *y) {
2558  return (x->angstrom < y->angstrom ? -1 : 1);
2559 }
2560 
2561 void model_select_angstrom(float angstrom, int32_t *modmin, int32_t *modmax, float *modrat) {
2562  static alphaTstr alphaT[MAXAERMOD];
2563  static int firstCall = 1;
2564 
2565  int im, im1, im2;
2566 
2567  if (firstCall) {
2568 
2569  int ib = windex(520, aertab->wave, aertab->nwave);
2570 
2571  /* grab angstrom coefficients and sort in ascending order */
2572  for (im = 0; im < aertab->nmodel; im++) {
2573  alphaT[im].angstrom = aertab->model[im]->angstrom[ib];
2574  alphaT[im].modnum = im;
2575  }
2576  qsort(alphaT, aertab->nmodel, sizeof (alphaTstr),
2577  (int (*)(const void *, const void *)) compalphaT);
2578 
2579  firstCall = 0;
2580  }
2581 
2582  for (im = 0; im < aertab->nmodel; im++) {
2583  if (angstrom < alphaT[im].angstrom)
2584  break;
2585  }
2586  im1 = MAX(MIN(im - 1, aertab->nmodel - 1), 0);
2587  im2 = MAX(MIN(im, aertab->nmodel - 1), 0);
2588 
2589  *modmin = alphaT[im1].modnum;
2590  *modmax = alphaT[im2].modnum;
2591 
2592  if (im1 == im2)
2593  *modrat = 1.0;
2594  else
2595  *modrat = (angstrom - alphaT[im1].angstrom) /
2596  (alphaT[im2].angstrom - alphaT[im1].angstrom);
2597 
2598  return;
2599 
2600 }
2601 
2602 /* ---------------------------------------------------------------------------------------- */
2603 /* model_taua() - compute AOT at input wavelengths using specified model */
2604 /* Note by W. Robinson - the aot is determined for SENSOR wavelength nir_l
2605  using the mu, mu0 for the geometry for that band. To get the aot at other
2606  TABLE bands, it starts with the aot(nir_l). So, how would the geometry
2607  change get properly translated to the other table bands? For now, I'm
2608  not accounting for that translation. Possible method would be to (for
2609  interpol=0 only) compute the aot(nir_l, geom(sensor band)) and use
2610  that aot to get the aot(geom(sensor band)) */
2611 
2612 /* ---------------------------------------------------------------------------------------- */
2613 void model_taua(int32_t sensorID, int modnum, float wave[], int32_t nwave,
2614  int32_t iwnir_l, float rhoa[], geom_str *geom, float wv, float taua[]) {
2615  float *aot;
2616  float *lnaot;
2617  float *rhoas;
2618 
2619  int iwnir = iwatab[iwnir_l];
2620  float *phase = model_phase(modnum, geom);
2621  int iw, iwg, iwtab;
2622  float maxwave;
2623 
2624  iwg = iwnir * geom->gmult; /* index for geometry */
2625 
2626  maxwave = MAX(aertab->nwave, nwave);
2627 
2628  if ((rhoas = (float *) calloc(nwave, sizeof (float))) == NULL) {
2629  printf("Unable to allocate space for rhoas.\n");
2630  exit(1);
2631  }
2632  if ((aot = (float *) calloc(maxwave, sizeof (float))) == NULL) {
2633  printf("Unable to allocate space for aot.\n");
2634  exit(1);
2635  }
2636  if ((lnaot = (float *) calloc(maxwave, sizeof (float))) == NULL) {
2637  printf("Unable to allocate space for lnaot.\n");
2638  exit(1);
2639  }
2640 
2641  /* get SS aerosol reflectance at longest sensor wavelength */
2642  rhoa_to_rhoas(sensorID, modnum, geom, wv, rhoa, wave, nwave, iwnir_l, iwnir_l, rhoas);
2643 
2644  /* get aerosol optical thickness at longest sensor wavelength */
2645  aot[iwnir] = rhoas[iwnir_l]*(4.0 * geom->csolz[iwg] * geom->csenz[iwg])
2646  / (phase[iwnir] * aertab->model[modnum]->albedo[iwnir]);
2647 
2648  /* get aerosol optical thickness at all other table wavelengths */
2649  for (iw = 0; iw < aertab->nwave; iw++) {
2650  /* note to self: i actually have aot(865) for czcs at this point */
2651  aot[iw] = aot[iwnir] * aertab->model[modnum]->extc[iw] / aertab->model[modnum]->extc[iwnir];
2652  if (interpol)
2653  lnaot[iw] = log(aot[iw]);
2654  }
2655 
2656  /* interpolate aot from table to sensor wavelengths */
2657  if (interpol) {
2658  for (iw = 0; iw < nwave; iw++) {
2659  iwtab = iwatab[iw];
2660  if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
2661  taua[iw] = exp(linterp(aertab->wave, lnaot, aertab->nwave, wave[iw]));
2662  else
2663  taua[iw] = aot[iwtab];
2664  }
2665  } else
2666  for (iw = 0; iw < nwave; iw++)
2667  taua[iw] = aot[iw];
2668 
2669  free(rhoas);
2670  free(aot);
2671  free(lnaot);
2672  return;
2673 }
2674 
2675 /*------------------------------------------------------------------------------------------ */
2676 /* model_taua_mseps() - for a given model and rhoa[iwnir_l], computes the AOT using MSEPS */
2677 /* coefficients relating rhoa to taua (in log space) */
2678 /*------------------------------------------------------------------------------------------ */
2679 int model_taua_mseps(int32_t modl, float wave[], int32_t nwave, int32_t iwnir_l,
2680  float rhoa[], geom_str *geom, float tau_pred[])
2681  {
2682  float *ac, *bc, *cc, *dc, *ec;
2683  double ax, bx, cx, fx;
2684  float ext_modl[nwave];
2685  float tau_iwnir_l;
2686  int iw, iwtab;
2687  int iwtab_l = iwatab[iwnir_l];
2688 
2689  /* get the coefficients for lg_rho vs lg_aot */
2690  ms_eps_coef(modl, iwnir_l, wave, geom, &ac, &bc, &cc, &dc, &ec);
2691 
2692 
2693  ax = (double) ac[iwtab_l] - log((double) rhoa[iwnir_l]);
2694  bx = (double) bc[iwtab_l];
2695  cx = (double) cc[iwtab_l];
2696 
2697  fx = bx * bx - 4.0 * ax*cx;
2698  if (fx > 0.0 && cx != 0.0) {
2699  tau_iwnir_l = exp(0.5*(-bx + sqrt(fx)) / cx);
2700  } else {
2701  tau_iwnir_l = 0.0;
2702  }
2703 
2704  /* get the extinction coefficients and compute AOT at all wavelengths */
2705  for (iw = 0; iw < nwave; iw++) {
2706  iwtab = iwatab[iw];
2707  ext_modl[iw] = aertab->model[modl]->extc[iwtab];
2708  }
2709  for (iw = 0; iw < nwave; iw++) {
2710  tau_pred[iw] = (ext_modl[iw] / ext_modl[iwnir_l]) * tau_iwnir_l;
2711  }
2712 
2713  return (0);
2714 }
2715 
2716 /* ---------------------------------------------------------------------------------------- */
2717 /* model_transmittance() - compute path Rayleigh-aerosol diffuse trans for specified model */
2718 /* */
2719 /* W. Robinson, SAIC, 24 Mar 2017, modify for band-dependent geometry */
2720 
2721 /* ---------------------------------------------------------------------------------------- */
2722 void model_transmittance(int modnum, float wave[], int32_t nwave,
2723  float *theta, int gmult, float taua[], float dtran[]) {
2724  static int firstCall = 1;
2725  static float *intexp;
2726  static int *inttst;
2727 
2728  int i1, i2, i, ig;
2729  int iw, iwtab;
2730  float a1, b1;
2731  float a2, b2;
2732  float x1, x2;
2733  float y1, y2;
2734  float xbar;
2735  float wt;
2736 
2737  /* disable interpolation, since we now allow LUT wavelengths to differ from sensor nominal wavelengths
2738  BAF, 6/2022
2739 
2740  if (firstCall) {
2741  float taur1, um1;
2742  float taur2, um2;
2743  firstCall = 0;
2744 
2745  intexp = (float *) malloc(nwave * sizeof (float));
2746  inttst = (int *) malloc(nwave * sizeof (int));
2747 
2748  for (iw = 0; iw < nwave; iw++) {
2749  intexp[iw] = 1.0;
2750  inttst[iw] = 0;
2751  iwtab = iwdtab[iw];
2752  if (fabs(wave[iw] - aertab->dtran_wave[iwtab]) > 0.51) {
2753  um1 = aertab->dtran_wave[iwtab] / 1000.0;
2754  um2 = wave[iw] / 1000.0;
2755  taur1 = 0.008569 * pow(um1, -4)*(1.0 + (0.0113 * pow(um1, -2))+(0.00013 * pow(um1, -4)));
2756  taur2 = 0.008569 * pow(um2, -4)*(1.0 + (0.0113 * pow(um2, -2))+(0.00013 * pow(um2, -4)));
2757  intexp[iw] = taur2 / taur1;
2758  inttst[iw] = 1;
2759  printf("Interpolating diffuse transmittance for %d from %f by %f\n",
2760  (int) wave[iw], aertab->dtran_wave[iwtab], intexp[iw]);
2761  }
2762  }
2763  }
2764  */
2765 
2766  /* use coefficients of nearest wavelength */
2767  for (iw = 0; iw < nwave; iw++) {
2768  if ((iw == 0) || (gmult != 0)) {
2769  ig = iw * gmult;
2770  /* find bracketing zenith angle indices */
2771  for (i = 0; i < aertab->dtran_ntheta; i++) {
2772  if (theta[ig] < aertab->dtran_theta[i])
2773  break;
2774  }
2775  if (i == aertab->dtran_ntheta) {
2776  i1 = i - 1;
2777  i2 = i - 1;
2778  wt = 0.0;
2779  } else {
2780  i1 = MIN(MAX(i - 1, 0), aertab->dtran_ntheta - 2);
2781  i2 = i1 + 1;
2782  x1 = aertab->dtran_airmass[i1];
2783  x2 = aertab->dtran_airmass[i2];
2784  xbar = 1.0 / cos(theta[ig] / radeg);
2785  wt = (xbar - x1) / (x2 - x1);
2786  }
2787  }
2788  iwtab = iwdtab[iw];
2789 
2790  a1 = aertab->model[modnum]->dtran_a[iwtab][i1];
2791  b1 = aertab->model[modnum]->dtran_b[iwtab][i1];
2792 
2793  a2 = aertab->model[modnum]->dtran_a[iwtab][i2];
2794  b2 = aertab->model[modnum]->dtran_b[iwtab][i2];
2795 
2796  /* disable interpolation, BAF, 6/2022
2797  if (inttst[iw]) {
2798  a1 = pow(a1, intexp[iw]);
2799  a2 = pow(a2, intexp[iw]);
2800  }
2801  */
2802 
2803  y1 = a1 * exp(-b1 * taua[iw]);
2804  y2 = a2 * exp(-b2 * taua[iw]);
2805 
2806  dtran[iw] = MAX(MIN((1.0 - wt) * y1 + wt*y2, 1.0), 1e-5);
2807  }
2808 
2809  return;
2810 }
2811 
2812 
2813 /* ---------------------------------------------------------------------------------------- */
2814 /* diff_tran() - compute Rayleigh-aerosol diffuse trans for selected model pair, both paths */
2815 /* ---------------------------------------------------------------------------------------- */
2816 void diff_tran(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_l,
2817  geom_str *geom, float wv, float pr, float taur[], int32_t modmin,
2818  int32_t modmax, float modrat, float rhoa[], float taua[], float tsol[],
2819  float tsen[], float tauamin[], float tauamax[], int taua_opt) {
2820  int iw, gmult, ig;
2821  float *tsolmin;
2822  float *tsolmax;
2823  float *tsenmin;
2824  float *tsenmax;
2825 
2826  if ((tsolmin = (float *) calloc(nwave, sizeof (float))) == NULL) {
2827  printf("Unable to allocate space for tsolmin.\n");
2828  exit(1);
2829  }
2830  if ((tsolmax = (float *) calloc(nwave, sizeof (float))) == NULL) {
2831  printf("Unable to allocate space for tsolmax.\n");
2832  exit(1);
2833  }
2834  if ((tsenmin = (float *) calloc(nwave, sizeof (float))) == NULL) {
2835  printf("Unable to allocate space for tsenmin.\n");
2836  exit(1);
2837  }
2838  if ((tsenmax = (float *) calloc(nwave, sizeof (float))) == NULL) {
2839  printf("Unable to allocate space for tsenmax.\n");
2840  exit(1);
2841  }
2842  gmult = geom->gmult;
2843  if (interpol == 1) gmult = 0; /* geom used for model wavelengths */
2844 
2845  /* get AOT per band for each model, if not already computed */
2846  if (taua_opt == 0) {
2847  // using single scattering approximation space (Gordan and Wang)
2848  model_taua(sensorID, modmin, wave, nwave, iwnir_l, rhoa, geom, wv, tauamin);
2849  model_taua(sensorID, modmax, wave, nwave, iwnir_l, rhoa, geom, wv, tauamax);
2850  //printf("%d %d %d %f %f %f\n",taua_opt,modmin, iwnir_l, rhoa[iwnir_l], tauamin[0], tauamin[iwnir_l]);
2851  } else if (taua_opt == 2) {
2852  // using multi-scattering relationship between rhoa and taua (Ahmad)
2853  model_taua_mseps(modmin, wave, nwave, iwnir_l, rhoa, geom, tauamin);
2854  model_taua_mseps(modmax, wave, nwave, iwnir_l, rhoa, geom, tauamax);
2855  //printf("%d %d %d %f %f %f\n",taua_opt,modmin, iwnir_l, rhoa[iwnir_l], tauamin[0], tauamin[iwnir_l]);
2856  }
2857 
2858  /* get diff trans sun to ground, per band for each model */
2859  model_transmittance(modmin, wave, nwave, geom->solz, gmult, tauamin, tsolmin);
2860  model_transmittance(modmax, wave, nwave, geom->solz, gmult, tauamax, tsolmax);
2861 
2862  /* get diff trans ground to sensor, per band for each model */
2863  model_transmittance(modmin, wave, nwave, geom->senz, gmult, tauamin, tsenmin);
2864  model_transmittance(modmax, wave, nwave, geom->senz, gmult, tauamax, tsenmax);
2865 
2866  /* interpolate and pressure correct */
2867  for (iw = 0; iw < nwave; iw++) {
2868  ig = iw * geom->gmult; /* geom used for sensor wavelengths */
2869  taua[iw] = tauamin[iw]*(1.0 - modrat) + tauamax[iw] * modrat;
2870  tsol[iw] = tsolmin[iw]*(1.0 - modrat) + tsolmax[iw] * modrat;
2871  tsen[iw] = tsenmin[iw]*(1.0 - modrat) + tsenmax[iw] * modrat;
2872 
2873  /* correct for pressure difference from standard pressure */
2874  tsol[iw] = tsol[iw] * exp(-0.5 * taur[iw] / geom->csolz[ig]*(pr / p0 - 1));
2875  tsen[iw] = tsen[iw] * exp(-0.5 * taur[iw] / geom->csenz[ig] *(pr / p0 - 1));
2876 
2877  if ((evalmask & TRANSSPHER) != 0) {
2878  /* correct for airmass difference, plane-parallel to spherical atmosphere */
2879  tsol[iw] = pow(tsol[iw], geom->airmass_sph[ig] / geom->airmass_plp[ig]);
2880  tsen[iw] = pow(tsen[iw], geom->airmass_sph[ig] / geom->airmass_plp[ig]);
2881  }
2882  }
2883 
2884  free(tsolmin);
2885  free(tsolmax);
2886  free(tsenmin);
2887  free(tsenmax);
2888 
2889  return;
2890 }
2891 
2892 
2893 
2894 /* ---------------------------------------------------------------------------------------- */
2895 /* smaer() - compute aerosol reflectance using MSEPS approach of Ahmad
2896  * perform spectral matching in reflectance space */
2897 /* output the reflectance of the two bracketing aerosol models */
2898 /*
2899  * Amir Ibrahim, Jul 2016. */
2900 
2901 /* ---------------------------------------------------------------------------------------- */
2902 
2903 int smaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s,
2904  int32_t iwnir_l, int32_t nmodels, int32_t mindx1[], int32_t mindx2[],
2905  int32_t mindx3[], geom_str *geom, float wv, float rhoa[],
2906  float rho_aer[], int32_t *modmin, int32_t *modmax, float *modrat,
2907  float tau_pred_min[], float tau_pred_max[], float tg_sol_sm[],
2908  float tg_sen_sm[], float Lt_sm[], int32_t ip) {
2909  int iw;
2910 
2911  if (!have_ms) {
2912  printf("\nThe multi-scattering spectral matching atmospheric correction method requires\n");
2913  printf("ams_all, bms_all, cms_all, dms_all, ems in the aerosol model tables.\n");
2914  exit(1);
2915  }
2916 
2917  // calculate the weighted aerosol reflectance from Red to SWIR
2918  if (spectral_matching(sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx1, mindx2, mindx3, geom, wv, rhoa,
2919  rho_aer, modmin, modmax, modrat, tau_pred_min, tau_pred_max, tg_sol_sm, tg_sen_sm, Lt_sm, ip) == 0) {
2920 
2921  // extrapolate to VIS
2922  for (iw = 0; iw < nwave; iw++) {
2923  rhoa[iw] = rho_aer[iw];
2924  }
2925  } else
2926  return (1);
2927 
2928  return (0);
2929 }
2930 
2931 /* ---------------------------------------------------------------------------------------- */
2932 /* ahmadaer() - compute aerosol reflectance using MSEPS approach of Ahmad */
2933 /* */
2934 /* Z. Ahmad, August 2014. */
2935 /* ---------------------------------------------------------------------------------------- */
2936 int ahmadaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
2937  int32_t nmodels, int32_t mindx[],
2938  geom_str *geom, float wv, float rhoa[],
2939  int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir,
2940  float tau_pred_min[], float tau_pred_max[]) {
2941  int iw;
2942 
2943  float rho_pred_min[nwave], rho_pred_max[nwave];
2944  float rho_aer[nwave], tau_aer[nwave];
2945 
2946  if (!have_ms) {
2947  printf("\nThe multi-scattering epsilon atmospheric correction method requires\n");
2948  printf("ams_all, bms_all, cms_all, dms_all, ems in the aerosol model tables.\n");
2949  exit(1);
2950  }
2951  if (aer_opt == AERRHMSEPS_lin) {
2952  if (ahmad_atm_corr_lin(sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx, geom, wv,
2953  rhoa, modmin, modmax, modrat, epsnir,
2954  tau_pred_max, tau_pred_min, rho_pred_max, rho_pred_min, tau_aer, rho_aer) != 0)
2955  return (1);
2956  } else {
2957  /* use the ms_epsilon method to get rhoa */
2958  if (ahmad_atm_corr(sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx, geom, wv,
2959  rhoa, modmin, modmax, modrat, epsnir,
2960  tau_pred_max, tau_pred_min, rho_pred_max, rho_pred_min, tau_aer, rho_aer) != 0)
2961  return (1);
2962  }
2963 
2964  for (iw = 0; iw < nwave; iw++) {
2965  rhoa[iw] = rho_aer[iw];
2966  }
2967 
2968  return (0);
2969 }
2970 
2971 
2972 /* ---------------------------------------------------------------------------------------- */
2973 /* wangaer() - compute aerosol reflectance using Gordon & Wang 1994 algorithm */
2974 /* */
2975 /* B. Franz, 1 June 2004. */
2976 
2977 /* ---------------------------------------------------------------------------------------- */
2978 
2979 int wangaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s,
2980  int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *geom,
2981  float wv, float rhoa[], int32_t *modmin, int32_t *modmax,
2982  float *modrat, float *epsnir, float tauamin[], float tauamax[]) {
2983  int modflg;
2984  float *epsmin;
2985  float *epsmax;
2986 
2987  float epsmin1;
2988  float epsmax1;
2989 
2990  float *rhoasmin;
2991  float *rhoasmax;
2992  float *rhoamin;
2993  float *rhoamax;
2994 
2995  float cc = 0.0;
2996  int iw;
2997 
2998  if ((rhoasmin = (float *) calloc(nwave, sizeof (float))) == NULL) {
2999  printf("Unable to allocate space for rhoasmin.\n");
3000  exit(1);
3001  }
3002  if ((rhoasmax = (float *) calloc(nwave, sizeof (float))) == NULL) {
3003  printf("Unable to allocate space for rhoasmax.\n");
3004  exit(1);
3005  }
3006  if ((rhoamin = (float *) calloc(nwave, sizeof (float))) == NULL) {
3007  printf("Unable to allocate space for rhoamin.\n");
3008  exit(1);
3009  }
3010  if ((rhoamax = (float *) calloc(nwave, sizeof (float))) == NULL) {
3011  printf("Unable to allocate space for rhoasmax.\n");
3012  exit(1);
3013  }
3014 
3015  /* find upper and lower-bounding models */
3016  modflg = model_select_wang(sensorID, wave, nwave, nmodels, mindx, geom, wv,
3017  rhoa, iwnir_s, iwnir_l, modmin, modmax, modrat, epsnir);
3018 
3019  /* if no lower-bounding aerosol model, set-up for extrapolation */
3020  if (modflg < 0)
3021  cc = log(*epsnir) / (wave[iwnir_l] - wave[iwnir_s]);
3022 
3023  /* get model epsilon for each bounding model, all wavelengths */
3024  epsmin = model_epsilon(*modmin, iwnir_l, wave, nwave, geom);
3025  epsmax = model_epsilon(*modmax, iwnir_l, wave, nwave, geom);
3026 
3027  /* get SS aerosol reflectance at longest wavelength for the two models */
3028  if (rhoa_to_rhoas(sensorID, *modmin, geom, wv, rhoa, wave, nwave, iwnir_l, iwnir_l, rhoasmin) != 0) {
3029  free(rhoamin);
3030  free(rhoasmin);
3031  free(rhoamax);
3032  free(rhoasmax);
3033  return (1);
3034  }
3035  if (rhoa_to_rhoas(sensorID, *modmax, geom, wv, rhoa, wave, nwave, iwnir_l, iwnir_l, rhoasmax) != 0) {
3036  free(rhoamin);
3037  free(rhoasmin);
3038  free(rhoamax);
3039  free(rhoasmax);
3040  return (1);
3041  }
3042  /* compute SS aerosol reflectance in all bands */
3043  for (iw = 0; iw < nwave; iw++) {
3044 
3045  epsmin1 = epsmin[iw];
3046  epsmax1 = epsmax[iw];
3047 
3048  if (modflg < 0) {
3049  epsmin1 = exp(cc * (wave[iwnir_l] - wave[iw]));
3050  epsmax1 = epsmin1;
3051  }
3052 
3053  rhoasmin[iw] = rhoasmin[iwnir_l] * epsmin1;
3054  rhoasmax[iw] = rhoasmax[iwnir_l] * epsmax1;
3055  }
3056 
3057  /* compute MS aerosol reflectance in visible bands */
3058  rhoas_to_rhoa(sensorID, *modmin, geom, wv, rhoasmin, wave, nwave, 0, nwave - 1, rhoamin);
3059  rhoas_to_rhoa(sensorID, *modmax, geom, wv, rhoasmax, wave, nwave, 0, nwave - 1, rhoamax);
3060 
3061  /* interpolate between upper and lower-bounding models */
3062  for (iw = 0; iw < nwave; iw++) {
3063  rhoa[iw] = rhoamin[iw]*(1.0 - (*modrat)) + rhoamax[iw]*(*modrat);
3064  }
3065 
3066  model_taua(sensorID, *modmin, wave, nwave, iwnir_l, rhoa, geom, wv, tauamin);
3067  model_taua(sensorID, *modmax, wave, nwave, iwnir_l, rhoa, geom, wv, tauamax);
3068 
3069  free(rhoamin);
3070  free(rhoasmin);
3071  free(rhoamax);
3072  free(rhoasmax);
3073 
3074  return (0);
3075 }
3076 
3077 /* ---------------------------------------------------------------------------------------- */
3078 /* model_select_franz() - Franz aerosol model selection process. */
3079 /* */
3080 /* B. Franz, 1 February 2009. */
3081 
3082 /* ---------------------------------------------------------------------------------------- */
3083 
3084 typedef struct rhoaT_struct {
3085  int32_t modnum;
3086  float rhoa;
3087  float eps;
3088 } rhoaTstr;
3089 
3090 int comp_rhoaT(rhoaTstr *x, rhoaTstr *y) {
3091  return (x->rhoa < y->rhoa ? -1 : 1);
3092 }
3093 
3094 int model_select_franz(int32_t sensorID, float wave[], int32_t nwave,
3095  int32_t nmodel, int32_t mindx[], geom_str *geom, float wv,
3096  float rhoa[], int32_t iwnir_s, int32_t iwnir_l, int32_t *modmin,
3097  int32_t *modmax, float *modrat, float *epsnir) {
3098 
3099  float *rhoas;
3100  float *rhoa_tmp;
3101  rhoaTstr rhoa_tab[MAXMODEL];
3102 
3103  float *eps;
3104  int jm, im;
3105  int jm1, jm2;
3106  float wt;
3107 
3108  *modmin = -1;
3109  *modmax = -1;
3110  *modrat = 0.0;
3111  *epsnir = 0.0;
3112 
3113  if ((rhoas = (float *) calloc(nwave, sizeof (float))) == NULL) {
3114  printf("Unable to allocate space for rhoas.\n");
3115  exit(1);
3116  }
3117  if ((rhoa_tmp = (float *) calloc(nwave, sizeof (float))) == NULL) {
3118  printf("Unable to allocate space for rhoas_tmp.\n");
3119  exit(1);
3120  }
3121 
3122  // predict MS aerosol reflectance assuming each model
3123 
3124  for (jm = 0; jm < nmodel; jm++) {
3125 
3126  im = mindx[jm];
3127 
3128  // get model epsilon at this geometry
3129  eps = model_epsilon(im, iwnir_l, wave, nwave, geom);
3130 
3131  // get SS aerosol reflectance at iwnir_l
3132  if (rhoa_to_rhoas(sensorID, im, geom, wv, rhoa, wave, nwave, iwnir_l, iwnir_l, rhoas) != 0) {
3133  free(rhoas);
3134  free(rhoa_tmp);
3135  return (1);
3136  }
3137  // get SS aerosol reflectance at iwnir_s
3138  rhoas[iwnir_s] = rhoas[iwnir_l] * eps[iwnir_s];
3139 
3140  // get MS aerosol reflectance at iwnir_s and save
3141  rhoas_to_rhoa(sensorID, im, geom, wv, rhoas, wave, nwave, iwnir_s, iwnir_s, rhoa_tmp);
3142 
3143  rhoa_tab[jm].modnum = im;
3144  rhoa_tab[jm].rhoa = rhoa_tmp[iwnir_s];
3145  rhoa_tab[jm].eps = eps[iwnir_s];
3146  }
3147 
3148  // put results in ascending order of predicted rhoa[iwnir_s]
3149  qsort(rhoa_tab, nmodel, sizeof (rhoaTstr), (int (*)(const void *, const void *)) comp_rhoaT);
3150 
3151  // compare observed rhoa with model predictions at iwnir_s to select models
3152  for (jm = 0; jm < nmodel; jm++) {
3153  if (rhoa_tab[jm].rhoa > rhoa[iwnir_s])
3154  break;
3155  }
3156  if (jm == 0) {
3157  jm1 = 0;
3158  jm2 = 1;
3159  } else if (jm == nmodel) {
3160  jm1 = nmodel - 2;
3161  jm2 = nmodel - 1;
3162  } else {
3163  jm1 = jm - 1;
3164  jm2 = jm1 + 1;
3165  }
3166  wt = (rhoa[iwnir_s] - rhoa_tab[jm1].rhoa) / (rhoa_tab[jm2].rhoa - rhoa_tab[jm1].rhoa);
3167 
3168  *modmin = rhoa_tab[jm1].modnum;
3169  *modmax = rhoa_tab[jm2].modnum;
3170  *modrat = wt;
3171  *epsnir = rhoa_tab[jm1].eps * (1.0 - wt) + rhoa_tab[jm2].eps*wt;
3172 
3173  free(rhoas);
3174  free(rhoa_tmp);
3175 
3176  return (0);
3177 }
3178 
3179 
3180 /* ---------------------------------------------------------------------------------------- */
3181 /* order_models is the sorting function used in rhaer */
3182 
3183 /* ---------------------------------------------------------------------------------------- */
3184 static int order_models(const void *p1, const void *p2) {
3185  aermodstr *x = *(aermodstr **) p1;
3186  aermodstr *y = *(aermodstr **) p2;
3187 
3188  if (x->rh == y->rh) {
3189  if (x->sd > y->sd)
3190  return ( 1);
3191  else
3192  return (-1);
3193  } else {
3194  if (x->rh > y->rh)
3195  return ( 1);
3196  else
3197  return (-1);
3198  }
3199 }
3200 
3201 
3202 /* ---------------------------------------------------------------------------------------- */
3203 /* rhaer() - compute aerosol reflectance using RH descrimination + desired selection scheme */
3204 
3205 /* ---------------------------------------------------------------------------------------- */
3206 int rhaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
3207  geom_str *geom, float wv, float rh, float pr, float taur[], float rhoa[],
3208  int32_t *modmin1, int32_t *modmax1, float *modrat1, int32_t *modmin2, int32_t *modmax2, float *modrat2,
3209  float *eps, float taua[], float tsol[], float tsen[], float tg_sol_sm[], float tg_sen_sm[], float Lt_sm[], int32_t ip) {
3210  static int firstCall = 1;
3211  static int nrh;
3212  static float rhtab[MAXAERMOD];
3213  static int nsd;
3214  static int sdtab[MAXAERMOD];
3215 
3216  float *rhoa1;
3217  float *rhoa2;
3218  float *taua1;
3219  float *taua2;
3220  float *tsol1;
3221  float *tsol2;
3222  float *tsen1;
3223  float *tsen2;
3224  float eps1=1.0;
3225  float eps2=1.0;
3226  //float modrat1;
3227  //float modrat2;
3228 
3229  float *tau_pred_min1;
3230  float *tau_pred_max1;
3231  float *tau_pred_min2;
3232  float *tau_pred_max2;
3233 
3234  //int32_t modmin1;
3235  //int32_t modmin2;
3236  //int32_t modmax1;
3237  //int32_t modmax2;
3238  int32_t mindx1[MAXAERMOD];
3239  int32_t mindx2[MAXAERMOD];
3240  int32_t mindx3[MAXAERMOD]; // Third model set defined --> Amir
3241  int irh1, irh2, irh;
3242  int irh3; // Third RH index --> Amir
3243  int isd;
3244  float wt;
3245 
3246 
3247  int iw, im;
3248 
3249  if (firstCall) {
3250  firstCall = 0;
3251  float lastrh = -1.0;
3252  int lastsd = -1;
3253  if (!have_rh) {
3254  printf("-E- %s line %d: This aerosol selection method requires models with a Relative Humidity attribute and size distribution.\n",
3255  __FILE__, __LINE__);
3256  exit(1);
3257  }
3258 
3259  // need in order of rh and sd within rh
3260  qsort(aertab->model, aertab->nmodel, sizeof (aermodstr*), (int (*)(const void *, const void *)) order_models);
3261 
3262  // count the number of model humidities and the number of model size distributions
3263  // note that use of a single model suite will yield nrh=1, which inherently avoids RH weighting that case
3264 
3265  nsd = 0;
3266  nrh = 0;
3267 
3268  for (im = 0; im < aertab->nmodel; im++) {
3269  if (aertab->model[im]->rh != lastrh) {
3270  rhtab[nrh] = aertab->model[im]->rh;
3271  lastrh = rhtab[nrh];
3272  nrh++;
3273  }
3274  if (nrh == 1 && aertab->model[im]->sd != lastsd) {
3275  sdtab[nsd] = aertab->model[im]->sd;
3276  lastsd = sdtab[nsd];
3277  nsd++;
3278  }
3279  }
3280  if (nrh * nsd != aertab->nmodel) {
3281  printf("-E- %s line %d: number of humidities (%d) x number of size distributions (%d) must equal number of models (%d).\n",
3282  __FILE__, __LINE__, nrh, nsd, aertab->nmodel);
3283  exit(1);
3284  } else {
3285  printf("%d aerosol models: %d humidities x %d size fractions\n", aertab->nmodel, nrh, nsd);
3286  for (irh = 0; irh < nrh; irh++) {
3287  for (isd = 0; isd < nsd; isd++) {
3288  im = irh * nsd + isd;
3289  printf("model %d, rh=%f, sd=%d, alpha=%f, name=%s\n",
3290  im, aertab->model[im]->rh, aertab->model[im]->sd, aertab->model[im]->angstrom[1], aertab->model[im]->name);
3291  }
3292  }
3293  }
3294  }
3295 
3296  // initialize
3297  if ((taua1 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3298  printf("Unable to allocate space for taua1.\n");
3299  exit(1);
3300  }
3301  if ((taua2 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3302  printf("Unable to allocate space for taua2.\n");
3303  exit(1);
3304  }
3305  if ((tsol1 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3306  printf("Unable to allocate space for tsol1.\n");
3307  exit(1);
3308  }
3309  if ((tsol2 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3310  printf("Unable to allocate space for tsol2.\n");
3311  exit(1);
3312  }
3313  if ((rhoa1 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3314  printf("Unable to allocate space for rhoa1.\n");
3315  exit(1);
3316  }
3317  if ((rhoa2 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3318  printf("Unable to allocate space for rhoa2.\n");
3319  exit(1);
3320  }
3321  if ((tsen1 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3322  printf("Unable to allocate space for tsen1.\n");
3323  exit(1);
3324  }
3325  if ((tsen2 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3326  printf("Unable to allocate space for tsen2.\n");
3327  exit(1);
3328  }
3329  if ((tau_pred_min1 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3330  printf("Unable to allocate space for tau_pred_min1.\n");
3331  exit(1);
3332  }
3333  if ((tau_pred_min2 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3334  printf("Unable to allocate space for tau_pred_min2.\n");
3335  exit(1);
3336  }
3337  if ((tau_pred_max1 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3338  printf("Unable to allocate space for tau_pred_max1.\n");
3339  exit(1);
3340  }
3341  if ((tau_pred_max2 = (float *) malloc(nwave * sizeof (float))) == NULL) {
3342  printf("Unable to allocate space for tau_pred_max2.\n");
3343  exit(1);
3344  }
3345 
3346  //*modmin = -1;
3347  //*modmax = -1;
3348  //*modrat = 0;
3349  //*modmin2 = -1;
3350  //*modmax2 = -1;
3351  //*modrat2 = 0;
3352 
3353  for (iw = 0; iw < nwave; iw++) {
3354  taua[iw] = -1.0;
3355  tsol[iw] = -1.0;
3356  tsen[iw] = -1.0;
3357  rhoa1[iw] = rhoa[iw];
3358  rhoa2[iw] = rhoa[iw];
3359  rhoa [iw] = BAD_FLT;
3360  }
3361 
3362 
3363  // adjust rh for spectral matching
3364  if (aer_opt == AERRHSM) {
3365  if (rh >= 95) {
3366  printf("Warning rh is greater than 95%%. Reset to 94%% rh=%f\n", rh);
3367  rh = 94;
3368  }
3369  }
3370 
3371  // find RH index and wts
3372  if (nrh == 1 || rhtab[0] > rh) { // actual RH < smallest model RH or only one model RH
3373  irh1 = 0;
3374  irh2 = 0;
3375  wt = 0.0;
3376  } else if (rhtab[nrh - 1] < rh) { // actual RH > smallestlargest model RH
3377  irh1 = nrh - 1;
3378  irh2 = nrh - 1;
3379  wt = 0.0;
3380  } else {
3381  for (irh = 0; irh < nrh; irh++) {
3382  if (rhtab[irh] > rh)
3383  break;
3384  }
3385  irh1 = MIN(MAX(0, irh - 1), nrh - 2);
3386  irh2 = irh1 + 1;
3387  wt = (rh - rhtab[irh1]) / (rhtab[irh2] - rhtab[irh1]);
3388  }
3389 
3390  // calculate the third closest RH index for the models suite -->Amir
3391  if (rh - rhtab[irh1] <= (rhtab[irh2] - rhtab[irh1]) / 2)
3392  if (irh1 > 1)
3393  irh3 = irh1 - 1;
3394  else
3395  irh3 = irh1;
3396  else
3397  if (irh2 < 10)
3398  irh3 = irh2 + 1;
3399  else
3400  irh3 = irh2;
3401 
3402  for (im = 0; im < nsd; im++) {
3403  mindx3[im] = irh3 * nsd + im;
3404  }
3405 
3406  // set indices of active model sets
3407 
3408  for (im = 0; im < nsd; im++) {
3409  mindx1[im] = irh1 * nsd + im;
3410  mindx2[im] = irh2 * nsd + im;
3411  }
3412 
3413  // compute aerosol reflectances, aot, diffuse trans, eps from first model set
3414 
3415  /* perform spectral matching from Red to SWIR in radiance space, based on Ahmad
3416  * Multi-scattering coeff method --> Amir*/
3417  if (aer_opt == AERRHSM) {
3418 
3419  int nmodels = 30;
3420 
3421  if (mindx3[0] == mindx1[0] || mindx3[0] == mindx2[0] || mindx3[0] > 79 || mindx3[0] < 0)
3422  nmodels = 20;
3423 
3424  else if (mindx1[0] == mindx2[0])
3425  nmodels = 10;
3426 
3427  int32_t vcal_flag = 1;
3428  if ((aertab->nmodel) == vcal_flag)
3429  nmodels = 1;
3430 
3431  int modmin, modmax;
3432  float modrat;
3433 
3434  float rho_aer[nwave];
3435  // printf("%zu %zu %zu\n",mindx1[0],mindx2[0],mindx3[0]);
3436 
3437  if (smaer(sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx1, mindx1, mindx1,
3438  geom, wv, rhoa1, rho_aer, &modmin, &modmax, &modrat, tau_pred_min1, tau_pred_max1, tg_sol_sm, tg_sen_sm, Lt_sm, ip) == 0) {
3439 
3440  *modmin1 = modmin;
3441  *modmax1 = modmax;
3442  *modrat1 = modrat;
3443  *modmin2 = modmin;
3444  *modmax2 = modmax;
3445  *modrat2 = modrat;
3446  //printf("modmin=%f\n",modrat);
3447  //tau_pred_min2 = tau_pred_min1;
3448  //tau_pred_max2 = tau_pred_max1;
3449  diff_tran(sensorID, wave, nwave, iwnir_l, geom, wv, pr, taur,
3450  *modmin1, *modmax1, *modrat1, rhoa1, taua1, tsol1, tsen1, tau_pred_min1, tau_pred_max1, 1);
3451  //return (0);
3452  } //second run
3453  if (smaer(sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx2, mindx2, mindx2,
3454  geom, wv, rhoa2, rho_aer, &modmin, &modmax, &modrat, tau_pred_min2, tau_pred_max2, tg_sol_sm, tg_sen_sm, Lt_sm, ip) == 0) {
3455 
3456  *modmin1 = modmin;
3457  *modmax1 = modmax;
3458  *modrat1 = modrat;
3459  *modmin2 = modmin;
3460  *modmax2 = modmax;
3461  *modrat2 = modrat;
3462  //printf("modmin=%f\n",modrat);
3463  //tau_pred_min2 = tau_pred_min1;
3464  //tau_pred_max2 = tau_pred_max1;
3465  diff_tran(sensorID, wave, nwave, iwnir_l, geom, wv, pr, taur,
3466  *modmin2, *modmax2, *modrat2, rhoa2, taua2, tsol2, tsen2, tau_pred_min2, tau_pred_max2, 1);
3467  //return (0);
3468  } else {
3469  free(taua1);
3470  free(taua2);
3471  free(tsol1);
3472  free(tsol2);
3473  free(tsen1);
3474  free(tsen2);
3475  free(rhoa1);
3476  free(rhoa2);
3477  free(tau_pred_min1);
3478  free(tau_pred_max1);
3479  free(tau_pred_min2);
3480  free(tau_pred_max2);
3481  return (1);
3482  }
3483  } else if (aer_opt == AERRHMSEPS || aer_opt == AERRHMSEPS_lin) {
3484  if (ahmadaer(sensorID, wave, nwave, iwnir_s, iwnir_l, nsd, mindx1,
3485  geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1) != 0) {
3486  free(taua1);
3487  free(taua2);
3488  free(tsol1);
3489  free(tsol2);
3490  free(tsen1);
3491  free(tsen2);
3492  free(rhoa1);
3493  free(rhoa2);
3494  free(tau_pred_min1);
3495  free(tau_pred_max1);
3496  free(tau_pred_min2);
3497  free(tau_pred_max2);
3498  return (1);
3499 
3500  }
3501  } else {
3502  if (wangaer(sensorID, wave, nwave, iwnir_s, iwnir_l, nsd, mindx1,
3503  geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1) != 0) {
3504  free(taua1);
3505  free(taua2);
3506  free(tsol1);
3507  free(tsol2);
3508  free(tsen1);
3509  free(tsen2);
3510  free(rhoa1);
3511  free(rhoa2);
3512  free(tau_pred_min1);
3513  free(tau_pred_max1);
3514  free(tau_pred_min2);
3515  free(tau_pred_max2);
3516  return (1);
3517  }
3518  }
3519 
3520  if (aer_opt != AERRHSM)
3521  diff_tran(sensorID, wave, nwave, iwnir_l, geom, wv, pr, taur,
3522  *modmin1, *modmax1, *modrat1, rhoa1, taua1, tsol1, tsen1, tau_pred_min1, tau_pred_max1, 1);
3523 
3524  // compute aerosol reflectances, aot, diffuse trans, eps from second model set (if needed)
3525 
3526  if (irh2 != irh1) {
3527 
3528  if (aer_opt == AERRHSM) {
3529  for (iw = 0; iw < nwave; iw++) {
3530  rhoa[iw] = rhoa1[iw]*(1 - wt) + rhoa2[iw] * wt;
3531  taua[iw] = taua1[iw]*(1 - wt) + taua2[iw] * wt;
3532  tsol[iw] = tsol1[iw]*(1 - wt) + tsol2[iw] * wt;
3533  tsen[iw] = tsen1[iw]*(1 - wt) + tsen2[iw] * wt;
3534  }
3535  *eps = eps1;
3536  return (1); // A second model set interpolation is not needed if the spectral matching is active
3537  }
3538 
3539  if (aer_opt == AERRHMSEPS || aer_opt == AERRHMSEPS_lin) {
3540  if (ahmadaer(sensorID, wave, nwave, iwnir_s, iwnir_l, nsd, mindx2,
3541  geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2) != 0) {
3542  free(taua1);
3543  free(taua2);
3544  free(tsol1);
3545  free(tsol2);
3546  free(tsen1);
3547  free(tsen2);
3548  free(rhoa1);
3549  free(rhoa2);
3550  free(tau_pred_min1);
3551  free(tau_pred_max1);
3552  free(tau_pred_min2);
3553  free(tau_pred_max2);
3554  return (1);
3555  }
3556  } else {
3557  if (wangaer(sensorID, wave, nwave, iwnir_s, iwnir_l, nsd, mindx2,
3558  geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2) != 0) {
3559  free(taua1);
3560  free(taua2);
3561  free(tsol1);
3562  free(tsol2);
3563  free(tsen1);
3564  free(tsen2);
3565  free(rhoa1);
3566  free(rhoa2);
3567  free(tau_pred_min1);
3568  free(tau_pred_max1);
3569  free(tau_pred_min2);
3570  free(tau_pred_max2);
3571  return (1);
3572  }
3573  }
3574 
3575  diff_tran(sensorID, wave, nwave, iwnir_l, geom, wv, pr, taur,
3576  *modmin2, *modmax2, *modrat2, rhoa2, taua2, tsol2, tsen2, tau_pred_min2, tau_pred_max2, 1);
3577 
3578  for (iw = 0; iw < nwave; iw++) {
3579  rhoa[iw] = rhoa1[iw]*(1 - wt) + rhoa2[iw] * wt;
3580  taua[iw] = taua1[iw]*(1 - wt) + taua2[iw] * wt;
3581  tsol[iw] = tsol1[iw]*(1 - wt) + tsol2[iw] * wt;
3582  tsen[iw] = tsen1[iw]*(1 - wt) + tsen2[iw] * wt;
3583  }
3584  *eps = eps1 * (1 - wt) + eps2*wt;
3585 
3586  } else {
3587 
3588  for (iw = 0; iw < nwave; iw++) {
3589  rhoa[iw] = rhoa1[iw];
3590  taua[iw] = taua1[iw];
3591  tsol[iw] = tsol1[iw];
3592  tsen[iw] = tsen1[iw];
3593  }
3594  *eps = eps1;
3595  }
3596 
3597  // store model info for the dominant RH only
3598 
3599  //if (wt < 0.5) {
3600  // *modmin = modmin1;
3601  // *modmax = modmax1;
3602  // *modrat = modrat1;
3603  //} else {
3604  // *modmin = modmin2;
3605  // *modmax = modmax2;
3606  // *modrat = modrat2;
3607  //}
3608 
3609  free(taua1);
3610  free(taua2);
3611  free(tsol1);
3612  free(tsol2);
3613  free(tsen1);
3614  free(tsen2);
3615  free(rhoa1);
3616  free(rhoa2);
3617  free(tau_pred_min1);
3618  free(tau_pred_max1);
3619  free(tau_pred_min2);
3620  free(tau_pred_max2);
3621 
3622  return (0);
3623 }
3624 
3625 
3626 /* ---------------------------------------------------------------------------------------- */
3627 /* fixedaer() - compute aerosol reflectance for fixed aerosol model */
3628 /* */
3629 /* B. Franz, August 2004. */
3630 
3631 /* ---------------------------------------------------------------------------------------- */
3632 int fixedaer(int32_t sensorID, int32_t modnum, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
3633  char models[MAXAERMOD][32], int32_t nmodels,
3634  geom_str *geom, float wv, float rhoa[], float *epsnir) {
3635  float *eps;
3636  float *rhoas;
3637  int iw;
3638 
3639  if ((rhoas = (float *) calloc(nwave, sizeof (float))) == NULL) {
3640  printf("Unable to allocate space for rhoas.\n");
3641  exit(1);
3642  }
3643 
3644  if (rhoa[iwnir_l] < 0.0) {
3645  *epsnir = BAD_FLT;
3646  // for (iw=0; iw<nwave; iw++)
3647  // rhoas[iw] = BAD_FLT;
3648  free(rhoas);
3649  return (1);
3650  }
3651 
3652  /* get model epsilon for all wavelengths at this geometry */
3653  eps = model_epsilon(modnum, iwnir_l, wave, nwave, geom);
3654 
3655  /* get SS aerosol reflectance at longest wavelength */
3656  if (rhoa_to_rhoas(sensorID, modnum, geom, wv, rhoa, wave, nwave, iwnir_l, iwnir_l, rhoas) != 0) {
3657  printf("Error getting rhoas\n");
3658  free(rhoas);
3659  return (1);
3660  }
3661 
3662  /* compute SS aerosol reflectance in visible bands */
3663  for (iw = 0; iw < nwave; iw++) {
3664  rhoas[iw] = rhoas[iwnir_l] * eps[iw];
3665  }
3666 
3667  /* compute MS aerosol reflectance in visible bands */
3668  rhoas_to_rhoa(sensorID, modnum, geom, wv, rhoas, wave, nwave, 0, nwave - 1, rhoa);
3669 
3670  if (iwnir_s == iwnir_l)
3671  *epsnir = eps[iwnir_l - 1];
3672  else
3673  *epsnir = eps[iwnir_s];
3674 
3675  free(rhoas);
3676 
3677  return (0);
3678 }
3679 
3680 
3681 /* ---------------------------------------------------------------------------------------- */
3682 /* fixedmodpair() - compute aerosol reflectance for fixed model pair */
3683 /* */
3684 /* B. Franz, August 2004. */
3685 
3686 /* ---------------------------------------------------------------------------------------- */
3687 int fixedmodpair(int32_t sensorID, float wave[], int32_t nwave,
3688  int32_t iwnir_s, int32_t iwnir_l, geom_str *geom, float wv,
3689  int32_t modmin, int32_t modmax, float modrat, float rhoa[], float *eps) {
3690  float *rhoa1;
3691  float *rhoa2;
3692  float eps1;
3693  float eps2;
3694  int iw;
3695  int status;
3696 
3697  if ((rhoa1 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3698  printf("Unable to allocate space for rhoa1.\n");
3699  exit(1);
3700  }
3701  if ((rhoa2 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3702  printf("Unable to allocate space for rhoa2.\n");
3703  exit(1);
3704  }
3705 
3706  if (modmin < 0 || modmin >= input->naermodels ||
3707  modmax < 0 || modmax >= input->naermodels ||
3708  modrat < 0.0 || modrat > 1.0) {
3709  printf("Bogus input for fixed model pair: %d %d %f\n", modmin + 1, modmax + 1, modrat);
3710  exit(1);
3711  }
3712 
3713 
3714  if (rhoa[iwnir_l] > input->rhoamin) {
3715 
3716  rhoa2[iwnir_l] = rhoa1[iwnir_l] = rhoa[iwnir_l];
3717 
3718  /* get aerosol reflectance in visible for fixed model, and epsilon */
3719  status = fixedaer(sensorID, modmin, wave, nwave, iwnir_s, iwnir_l,
3720  input->aermodels, input->naermodels, geom, wv, rhoa1, &eps1);
3721  status = fixedaer(sensorID, modmax, wave, nwave, iwnir_s, iwnir_l,
3722  input->aermodels, input->naermodels, geom, wv, rhoa2, &eps2);
3723 
3724  /* convert aerosol relectance to radiance */
3725  if (status == 0) {
3726  for (iw = 0; iw < nwave; iw++) {
3727  if (iw != iwnir_l) // without this check tLw-La may go slight negative
3728  rhoa[iw] = MAX((1.0 - modrat) * rhoa1[iw] + modrat * rhoa2[iw], 0.0);
3729  }
3730  *eps = (1.0 - modrat) * eps1 + modrat*eps2;
3731  }
3732 
3733  } else if (rhoa[iwnir_l] > -(input->rhoamin)) {
3734 
3735  /* if input NIR is near zero, assume a small white aerosol signal */
3736  *eps = 1.0;
3737  for (iw = 0; iw < nwave; iw++) {
3738  rhoa[iw] = MAX(rhoa[iwnir_l], 1e-6);
3739  }
3740 
3741  status = 0;
3742 
3743  } else {
3744 
3745  /* if input NIR is very negative, fail the case */
3746  status = 1;
3747  }
3748 
3749  free(rhoa1);
3750  free(rhoa2);
3751 
3752  return (status);
3753 }
3754 
3755 
3756 /* ---------------------------------------------------------------------------------------- */
3757 /* fixedaot() - compute aerosol reflectance for fixed aot(lambda) */
3758 /* */
3759 /* B. Franz, August 2004. */
3760 
3761 /* ---------------------------------------------------------------------------------------- */
3762 
3763 int fixedaot(int32_t sensorID, float aot[], float wave[], int32_t nwave,
3764  int32_t iwnir_s, int32_t iwnir_l, geom_str *geom, float wv,
3765  int32_t *modmin, int32_t *modmax, float *modrat, float rhoa[],
3766  float *epsnir) {
3767  static int firstCall = 1;
3768  static int angst_band1 = -1;
3769  static int angst_band2 = -1;
3770 
3771  float *phase1;
3772  float *phase2;
3773  float *f1;
3774  float *f2;
3775  float *lnf1;
3776  float *lnf2;
3777  float *rhoas1;
3778  float *rhoas2;
3779  float *rhoa1;
3780  float *rhoa2;
3781  float eps1;
3782  float eps2;
3783  float angstrom;
3784  int ig, gmult, iw, iwtab;
3785  float maxwave;
3786 
3787  maxwave = MAX(aertab->nwave, nwave);
3788 
3789  if ((rhoa1 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3790  printf("Unable to allocate space for rhoa1.\n");
3791  exit(1);
3792  }
3793  if ((rhoa2 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3794  printf("Unable to allocate space for rhoa2.\n");
3795  exit(1);
3796  }
3797  if ((rhoas1 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3798  printf("Unable to allocate space for rhoas1.\n");
3799  exit(1);
3800  }
3801  if ((rhoas2 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3802  printf("Unable to allocate space for rhoas2.\n");
3803  exit(1);
3804  }
3805  if ((f1 = (float *) calloc(maxwave, sizeof (float))) == NULL) {
3806  printf("Unable to allocate space for f1.\n");
3807  exit(1);
3808  }
3809  if ((f2 = (float *) calloc(maxwave, sizeof (float))) == NULL) {
3810  printf("Unable to allocate space for f2.\n");
3811  exit(1);
3812  }
3813  if ((lnf1 = (float *) calloc(maxwave, sizeof (float))) == NULL) {
3814  printf("Unable to allocate space for lnf1.\n");
3815  exit(1);
3816  }
3817  if ((lnf2 = (float *) calloc(maxwave, sizeof (float))) == NULL) {
3818  printf("Unable to allocate space for lnf2.\n");
3819  exit(1);
3820  }
3821 
3822  if (firstCall) {
3823  angst_band1 = windex(520, wave, nwave);
3824  angst_band2 = windex(865, wave, nwave);
3825  firstCall = 0;
3826  }
3827 
3828  /* bail on negative input aot */
3829  for (iw = 0; iw < nwave; iw++) {
3830  if (aot[iw] < 0.0) {
3831  free(rhoa1);
3832  free(rhoa2);
3833  free(rhoas1);
3834  free(rhoas2);
3835  free(f1);
3836  free(f2);
3837  free(lnf1);
3838  free(lnf2);
3839  return (1);
3840  }
3841  }
3842 
3843  /* compute angstrom and use to select bounding models */
3844  if (aot[iwnir_l] > 0.0)
3845  angstrom = -log(aot[angst_band1] / aot[angst_band2]) /
3846  log(wave[angst_band1] / wave[angst_band2]);
3847  else
3848  angstrom = 0.0;
3849 
3850  model_select_angstrom(angstrom, modmin, modmax, modrat);
3851 
3852 
3853  /* get model phase function for all wavelengths at this geometry for the two models */
3854  phase1 = model_phase(*modmin, geom);
3855  phase2 = model_phase(*modmax, geom);
3856 
3857  gmult = (interpol == 1) ? 0 : geom->gmult;
3858 
3859  /* compute factor for SS approximation, set-up for interpolation */
3860  for (iw = 0; iw < aertab->nwave; iw++) {
3861  ig = gmult * iw;
3862  f1[iw] = aertab->model[*modmin]->albedo[iw] *
3863  phase1[iw] / 4.0 / geom->csolz[ig] / geom->csenz[ig];
3864  f2[iw] = aertab->model[*modmax]->albedo[iw] *
3865  phase2[iw] / 4.0 / geom->csolz[ig] / geom->csenz[ig];
3866  if (interpol) {
3867  lnf1[iw] = log(f1[iw]);
3868  lnf2[iw] = log(f2[iw]);
3869  }
3870  }
3871 
3872  /* compute SS aerosol reflectance */
3873  if (interpol) {
3874  for (iw = 0; iw < nwave; iw++) {
3875  iwtab = iwatab[iw];
3876  if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0) {
3877  rhoas1[iw] = aot[iw] * exp(linterp(aertab->wave, lnf1, aertab->nwave, wave[iw]));
3878  rhoas2[iw] = aot[iw] * exp(linterp(aertab->wave, lnf2, aertab->nwave, wave[iw]));
3879  } else {
3880  rhoas1[iw] = aot[iw] * f1[iwtab];
3881  rhoas2[iw] = aot[iw] * f2[iwtab];
3882  }
3883  }
3884  } else {
3885  for (iw = 0; iw < nwave; iw++) {
3886  iwtab = iwatab[iw];
3887  rhoas1[iw] = aot[iw] * f1[iwtab];
3888  rhoas2[iw] = aot[iw] * f2[iwtab];
3889  }
3890  }
3891  eps1 = rhoas1[iwnir_s] / rhoas1[iwnir_l];
3892  eps2 = rhoas2[iwnir_s] / rhoas2[iwnir_l];
3893 
3894  /* compute MS aerosol reflectance */
3895  rhoas_to_rhoa(sensorID, *modmin, geom, wv, rhoas1, wave, nwave, 0, nwave - 1, rhoa1);
3896  rhoas_to_rhoa(sensorID, *modmax, geom, wv, rhoas2, wave, nwave, 0, nwave - 1, rhoa2);
3897 
3898  for (iw = 0; iw < nwave; iw++) {
3899  rhoa[iw] = (1.0 - *modrat) * rhoa1[iw] + *modrat * rhoa2[iw];
3900  }
3901  *epsnir = (1.0 - *modrat) * eps1 + *modrat * eps2;
3902 
3903  free(rhoa1);
3904  free(rhoa2);
3905  free(rhoas1);
3906  free(rhoas2);
3907  free(f1);
3908  free(f2);
3909  free(lnf1);
3910  free(lnf2);
3911 
3912  return (0);
3913 }
3914 
3915 
3916 /* ---------------------------------------------------------------------------------------- */
3917 /* aerosol() - compute aerosol reflectance using specified algorithm */
3918 /* */
3919 /* B. Franz, 1 June 2004. */
3920 
3921 /* ---------------------------------------------------------------------------------------- */
3922 int aerosol(l2str *l2rec, int32_t aer_opt_in, aestr *aerec, int32_t ip,
3923  float wave[], int32_t nwave, int32_t iwnir_s_in, int32_t iwnir_l_in,
3924  float F0_in[], float La1_in[], float La2_out[],
3925  float t_sol_out[], float t_sen_out[], float *eps, float taua_out[],
3926  int32_t *modmin, int32_t *modmax, float *modrat,
3927  int32_t *modmin2, int32_t *modmax2, float *modrat2) {
3928  static int firstCall = 1;
3929  static int32_t mindx[MAXAERMOD];
3930 
3931  int status = 1;
3932  float *rhoa;
3933  float *radref;
3934  float temp;
3935  float *aot;
3936  float angstrom;
3937 
3938  int iw, ipb;
3939 
3940  float *F0;
3941  float *taur;
3942  float *La1;
3943  float *La2;
3944  float *t_sol;
3945  float *t_sen;
3946  float *taua;
3947  float *taua_pred_min;
3948  float *taua_pred_max;
3949 
3950  l1str *l1rec = l2rec->l1rec;
3951 
3952  float *tg_sol_sm = l2rec->l1rec->tg_sol;
3953  float *tg_sen_sm = l2rec->l1rec->tg_sen;
3954  float *Lt_sm = l2rec->l1rec->Lt;
3955 
3956  int32_t sensorID = l1rec->l1file->sensorID;
3957  float solz = l1rec->solz [ip];
3958  float senz = l1rec->senz [ip];
3959  float wv = l1rec->wv [ip];
3960  float rh = l1rec->rh [ip];
3961  float pr = l1rec->pr [ip];
3962  int taua_opt;
3963 
3964  if (firstCall == 1) {
3965  Nbands = nwave;
3966  Maxband = nwave + 1; /* Must be >= Nbands */
3967  }
3968 
3969  if ((rhoa = (float *) calloc(nwave, sizeof (float))) == NULL) {
3970  printf("Unable to allocate space for rhoa.\n");
3971  exit(1);
3972  }
3973  if ((radref = (float *) calloc(nwave, sizeof (float))) == NULL) {
3974  printf("Unable to allocate space for raderef.\n");
3975  exit(1);
3976  }
3977  if ((F0 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3978  printf("Unable to allocate space for F0.\n");
3979  exit(1);
3980  }
3981  if ((taur = (float *) calloc(nwave, sizeof (float))) == NULL) {
3982  printf("Unable to allocate space for taur.\n");
3983  exit(1);
3984  }
3985  if ((La1 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3986  printf("Unable to allocate space for La1.\n");
3987  exit(1);
3988  }
3989  if ((La2 = (float *) calloc(nwave, sizeof (float))) == NULL) {
3990  printf("Unable to allocate space for La2.\n");
3991  exit(1);
3992  }
3993  if ((t_sol = (float *) calloc(nwave, sizeof (float))) == NULL) {
3994  printf("Unable to allocate space for t_sol.\n");
3995  exit(1);
3996  }
3997  if ((t_sen = (float *) calloc(nwave, sizeof (float))) == NULL) {
3998  printf("Unable to allocate space for t_sen.\n");
3999  exit(1);
4000  }
4001  if ((taua = (float *) calloc(nwave, sizeof (float))) == NULL) {
4002  printf("Unable to allocate space for taua.\n");
4003  exit(1);
4004  }
4005  if ((taua_pred_min = (float *) calloc(nwave, sizeof (float))) == NULL) {
4006  printf("Unable to allocate space for taua_pred_min.\n");
4007  exit(1);
4008  }
4009  if ((taua_pred_max = (float *) calloc(nwave, sizeof (float))) == NULL) {
4010  printf("Unable to allocate space for taua_pred_max.\n");
4011  exit(1);
4012  }
4013 
4014  /* set up the geometry structure to be passed */
4015  if (l1rec->geom_per_band == NULL) {
4016  /* nominal geometry setup */
4017  geom.senz = &l1rec->senz[ip];
4018  geom.solz = &l1rec->solz[ip];
4019  geom.phi = &l1rec->delphi[ip];
4020  geom.csolz = &l1rec->csolz[ip];
4021  geom.csenz = &l1rec->csenz[ip];
4022  geom.gmult = 0;
4023 
4024  geom.airmass = (float *) malloc(sizeof (float));
4025  *geom.airmass = 1. / geom.csolz[0] + 1. / geom.csenz[0];
4026  if ((evalmask & TRANSSPHER) != 0) {
4027  geom.airmass_sph = (float *) malloc(sizeof (float));
4028  geom.airmass_plp = (float *) malloc(sizeof (float));
4029  *geom.airmass_sph = ky_airmass(geom.solz[0]) +
4030  ky_airmass(geom.senz[0]);
4031  *geom.airmass_plp = pp_airmass(geom.solz[0]) +
4032  pp_airmass(geom.senz[0]);
4033  }
4034  } else {
4035  ipb = ip * Nbands;
4036  geom.senz = &l1rec->geom_per_band->senz[ipb];
4037  geom.solz = &l1rec->geom_per_band->solz[ipb];
4038  geom.phi = &l1rec->geom_per_band->delphi[ipb];
4039  geom.csolz = &l1rec->geom_per_band->csolz[ipb];
4040  geom.csenz = &l1rec->geom_per_band->csenz[ipb];
4041  geom.gmult = 1;
4042 
4043  geom.airmass = (float *) malloc(Nbands * sizeof (float));
4044  for (iw = 0; iw < Nbands; iw++) {
4045  geom.airmass[iw] = 1. / geom.csolz[iw] + 1. / geom.csenz[iw];
4046  }
4047  if ((evalmask & TRANSSPHER) != 0) {
4048  geom.airmass_plp = (float *) malloc(Nbands * sizeof (float));
4049  geom.airmass_sph = (float *) malloc(Nbands * sizeof (float));
4050  for (iw = 0; iw < Nbands; iw++) {
4051  geom.airmass_plp[iw] = pp_airmass(geom.solz[iw]) +
4052  pp_airmass(geom.senz[iw]);
4053  geom.airmass_sph[iw] = ky_airmass(geom.solz[iw]) +
4054  ky_airmass(geom.senz[iw]);
4055  }
4056  }
4057  }
4058 
4059  /* set static global evaluation level */
4060  evalmask = l1_input->evalmask;
4061  aer_opt = aer_opt_in;
4062 
4063  /* transfer inputs per band to inputs per wavelength */
4064  for (iw = 0; iw < nwave; iw++) {
4065  F0 [iw] = F0_in [iw];
4066  taur[iw] = l1rec->l1file->Tau_r[iw];
4067  La1 [iw] = La1_in[iw];
4068  if (iw == iwnir_s_in) iwnir_s = iw;
4069  if (iw == iwnir_l_in) iwnir_l = iw;
4070  }
4071 
4072  /* compute total airmass (static global) */
4073  mu0 = cos(solz / radeg);
4074  mu = cos(senz / radeg);
4075  airmass = 1.0 / mu0 + 1.0 / mu;
4076  if ((evalmask & TRANSSPHER) != 0) {
4077  airmass_plp = pp_airmass(solz) + pp_airmass(senz);
4078  airmass_sph = ky_airmass(solz) + ky_airmass(senz);
4079  }
4080 
4081  /* initialize epsilon and diffuse transmittance to defaults */
4082  *eps = 1.0;
4083  *modmin = BAD_INT;
4084  *modmax = BAD_INT;
4085  *modrat = BAD_FLT;
4086  *modmin2 = BAD_INT;
4087  *modmax2 = BAD_INT;
4088  *modrat2 = BAD_FLT;
4089  for (iw = 0; iw < nwave; iw++) {
4090  taua [iw] = 0.0;
4091  t_sol [iw] = 1.0;
4092  t_sen [iw] = 1.0;
4093  La2 [iw] = 0.0;
4094  radref[iw] = pi / F0[iw] / geom.csolz[iw * geom.gmult];
4095  }
4096 
4097  /* load aerosol model tables */
4098  if (!loaded) {
4099  int32_t im;
4100  load_aermod(sensorID, wave, nwave, input->aermodfile, input->aermodels, input->naermodels);
4101  for (im = 0; im < aertab->nmodel; im++) mindx[im] = im;
4102  if (have_rh && aertab->nmodel >= 30) {
4103  printf("\nLimiting aerosol models based on RH.\n");
4104  use_rh = 1;
4105  }
4106  }
4107  /* do not permit use of geom_per_band if interpolation is done */
4108  /* WDR note that it CAN work, but currently, we don't want users to
4109  use this combination */
4110  if ((interpol == 1) && (geom.gmult == 1)) {
4111  fprintf(stderr, "-E- %s line %d: Interpolated aerosol tables are\n",
4112  __FILE__, __LINE__);
4113  fprintf(stderr, " not permitted for use with band-dependent geometry, set geom_per_band=0\n");
4114  exit(1);
4115  }
4116 
4117  /* Change the aerosol option if need be */
4118  if (use_rh)
4119  switch (aer_opt) {
4120  case AERWANG: aer_opt = AERRH;
4121  break;
4122  case AERWANGNIR: aer_opt = AERRHNIR;
4123  break;
4124  case AERWANGSWIR: aer_opt = AERRHSWIR;
4125  break;
4126  case AERMUMM: aer_opt = AERRHMUMM;
4127  break;
4128  }
4129 
4130 
4131  if (firstCall) {
4132  if (aer_opt > 0 && aer_opt <= MAXAERMOD) {
4133  printf("\nUsing fixed aerosol model #%d (%s)\n", aer_opt, input->aermodels[aer_opt - 1]);
4134  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4135  } else {
4136  switch (aer_opt) {
4137  case AERRHSM: // Spectral Matching --> Amir
4138  printf("\nUsing Spectral Matching of aerosols reflectance for\n");
4139  printf("wavelength from %4.1f nm to %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
4140  break;
4141  case AERWHITE:
4142  printf("\nUsing white-aerosol approximation\n");
4143  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4144  break;
4145  case AERWANG:
4146  case AERRH:
4147  printf("\nUsing Gordon & Wang aerosol model selection\n");
4148  printf("Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
4149  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4150  break;
4151  case AERWANGNIR:
4152  case AERRHNIR:
4153  printf("\nUsing Gordon & Wang aerosol model selection\n");
4154  printf(" and NIR correction with up to %d iterations\n", input->aer_iter_max);
4155  printf("Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
4156  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4157  break;
4158  case AERWANGSWIR:
4159  case AERRHSWIR:
4160  printf("\nUsing Gordon & Wang aerosol model selection with NIR/SWIR switching.\n");
4161  printf("NIR correction with up to %d iterations\n", input->aer_iter_max);
4162  printf("NIR bands at %d and %d nm\n", input->aer_wave_short, input->aer_wave_long);
4163  printf("SWIR bands at %d and %d nm\n\n", input->aer_swir_short, input->aer_swir_long);
4164  break;
4165  case AERMUMM:
4166  case AERRHMUMM:
4167  printf("\nUsing Gordon & Wang aerosol model selection\n");
4168  printf(" and MUMM correction\n");
4169  printf("Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
4170  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4171  break;
4172  case AERRHMSEPS:
4173  printf("\nUsing multi-scattering aerosol model selection.\n");
4174  printf(" and NIR correction with up to %d iterations\n", input->aer_iter_max);
4175  printf("Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
4176  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4177  break;
4178  case AERRHMSEPS_lin:
4179  printf("\nUsing multi-scattering aerosol model selection in linear space.\n");
4180  printf(" and NIR correction with up to %d iterations\n", input->aer_iter_max);
4181  printf("Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
4182  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4183  break;
4184  case FIXAOT:
4185  printf("\nUsing fixed, input aerosol optical thicknesses for aerosol selection.\n");
4186  break;
4187  case FIXMODPAIR:
4188  printf("\nUsing multiple scattering aerosols with fixed model pair\n");
4189  break;
4190  case FIXMODPAIRNIR:
4191  printf("\nUsing multiple scattering aerosols with fixed model pair\n");
4192  printf(" and NIR iteration with up to %d iterations\n", input->aer_iter_max);
4193  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4194  break;
4195  case FIXANGSTROM:
4196  printf("\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
4197  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4198  break;
4199  case FIXANGSTROMNIR:
4200  printf("\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
4201  printf(" and NIR iteration with up to %d iterations\n", input->aer_iter_max);
4202  printf("Extrapolating from %4.1f nm\n", wave[iwnir_l]);
4203  break;
4204  default:
4205  printf("\nErroneous aerosol option, %d\n", aer_opt);
4206  exit(FATAL_ERROR);
4207  break;
4208  }
4209  }
4210  firstCall = 0;
4211  }
4212 
4213 
4214  /* convert input aerosol radiances to relectance */
4215  for (iw = iwnir_s; iw <= iwnir_l; iw += MAX(iwnir_l - iwnir_s, 1))
4216  rhoa[iw] = La1[iw] * radref[iw];
4217 
4218 
4219  /* compute aerosol using requested method */
4220  /* -------------------------------------- */
4221 
4222  switch (aer_opt) {
4223 
4224  case AERWANG: case AERWANGNIR: case AERWANGSWIR: case AERMUMM:
4225 
4226  if (iwnir_l <= iwnir_s || wave[iwnir_s] < 600) {
4227  printf("Aerosol selection bands must be greater than 600nm with short wave less than long wave (%d,%d)\n", iwnir_l, iwnir_s);
4228  exit(1);
4229  }
4230 
4231  /* Multi-Scattering with Gordon & Wang Aerosol Selection */
4232  /* ----------------------------------------------------- */
4233 
4234  /* convert input NIR aerosol radiances to reflectance */
4235  for (iw = iwnir_s; iw <= iwnir_l; iw++) {
4236  rhoa[iw] = La1[iw] * radref[iw];
4237  }
4238 
4239  /* require sufficient signal in two NIR bands */
4240  if (rhoa[iwnir_s] > input->rhoamin && rhoa[iwnir_l] > input->rhoamin) {
4241 
4242  /* require MS epsilon to be reasonable */
4243  if (La1[iwnir_s] / La1[iwnir_l] > 0.1) {
4244 
4245  status = wangaer(sensorID, wave, nwave, iwnir_s, iwnir_l,
4246  aertab->nmodel, mindx,
4247  &geom, wv, rhoa, modmin, modmax, modrat, eps, taua, taua); // this taua is not used //
4248 
4249  if (status == 0)
4250  for (iw = 0; iw < nwave; iw++)
4251  La2[iw] = rhoa[iw] / radref[iw];
4252  }
4253 
4254  } else if (rhoa[iwnir_s] > -(input->rhoamin) && rhoa[iwnir_l] > -(input->rhoamin)) {
4255 
4256  /* if input NIR is near zero, assume a small white aerosol signal */
4257  *eps = 1.0;
4258  *modmin = aertab->nmodel;
4259  *modmax = aertab->nmodel;
4260  *modrat = 0.0;
4261  temp = MAX(rhoa[iwnir_l], 1e-6);
4262  for (iw = 0; iw < nwave; iw++) {
4263  rhoa[iw] = temp;
4264  La2 [iw] = rhoa[iw] / radref[iw];
4265  }
4266 
4267  status = 0;
4268 
4269  } else {
4270 
4271  /* if input NIR is very negative, fail the case */
4272  status = 1;
4273  }
4274 
4275  break;
4276 
4277  case AERRH:
4278  case AERRHNIR:
4279  case AERRHSWIR:
4280  case AERRHMUMM:
4281  case AERRHMSEPS:
4282  case AERRHMSEPS_lin:
4283  case AERRHSM: // spectral matching --> Amir
4284 
4285  if (iwnir_l <= iwnir_s || wave[iwnir_s] < 600) {
4286  printf("Aerosol selection bands must be greater than 600nm with short wave less than long wave");
4287  exit(1);
4288  }
4289 
4290  /* Multi-Scattering with RH-based selection */
4291  /* ----------------------------------------------------- */
4292 
4293  /* convert input NIR aerosol radiances to relectance */
4294  for (iw = iwnir_s; iw <= iwnir_l; iw++) {
4295  rhoa[iw] = La1[iw] * radref[iw];
4296  }
4297 
4298  /* require sufficient signal in two NIR bands */
4299  if (rhoa[iwnir_s] > input->rhoamin && rhoa[iwnir_l] > input->rhoamin) {
4300 
4301  /* require MS epsilon to be reasonable */
4302  if (rhoa[iwnir_s] / rhoa[iwnir_l] > 0.1 && rhoa[iwnir_s] / rhoa[iwnir_l] < 10.0) {
4303 
4304  status = rhaer(sensorID, wave, nwave, iwnir_s, iwnir_l,
4305  &geom, wv, rh, pr, taur, rhoa,
4306  modmin, modmax, modrat, modmin2, modmax2, modrat2, eps, taua, t_sol, t_sen, tg_sol_sm, tg_sen_sm, Lt_sm, ip);
4307  status = 0;
4308  if (status == 0)
4309  for (iw = 0; iw < nwave; iw++)
4310  La2[iw] = rhoa[iw] / radref[iw];
4311  }
4312 
4313  } else if (rhoa[iwnir_s] > -(input->rhoamin) && rhoa[iwnir_l] > -(input->rhoamin)) {
4314 
4315  /* if input NIR is near zero, assume a small white aerosol signal */
4316  *eps = 1.0;
4317  *modmin = aertab->nmodel;
4318  *modmax = aertab->nmodel;
4319  *modrat = 0.0;
4320  *modmin2 = aertab->nmodel;
4321  *modmax2 = aertab->nmodel;
4322  *modrat2 = 0.0;
4323  temp = MAX(rhoa[iwnir_l], 1e-6);
4324  for (iw = 0; iw < nwave; iw++) {
4325  rhoa[iw] = temp;
4326  La2 [iw] = rhoa[iw] / radref[iw];
4327  }
4328 
4329  // The diff_tran function needs taua. In this exception case, taua was not yet computed. The taua_opt
4330  // tells diff_tran what method to use for the computation (0=SS, 2=MS). I used the aer_opt to decide,
4331  // but we could use the global have_ms flag. For now I want to limit the impact, as some sensors are
4332  // run with the G&W SS aerosol option, but their tables "have_ms". BAF, 6/2022.
4333 
4334  if (aer_opt == AERRHMSEPS || aer_opt == AERRHMSEPS_lin)
4335  taua_opt = 2;
4336  else
4337  taua_opt = 0;
4338 
4339  diff_tran(sensorID, wave, nwave, iwnir_l, &geom, wv, pr, taur,
4340  *modmin, *modmax, *modrat, rhoa, taua, t_sol, t_sen, taua_pred_min, taua_pred_max, taua_opt);
4341 
4342  status = 0;
4343 
4344  } else {
4345 
4346  /* if input NIR is very negative, fail the case */
4347  status = 1;
4348  }
4349 
4350  break;
4351 
4352  case AERWHITE:
4353 
4354  /* White Aerosol */
4355  /* ------------- */
4356 
4357  if (La1[iwnir_l] > 0.0) {
4358 
4359  *eps = 1.0;
4360  *modmin = 0;
4361  *modmax = 0;
4362  *modrat = 0.0;
4363 
4364  for (iw = 0; iw < nwave; iw++) {
4365  La2 [iw] = *eps * F0[iw] / F0[iwnir_l] * La1[iwnir_l];
4366  rhoa[iw] = La2[iw] * radref[iw];
4367  }
4368 
4369 
4370  status = 0;
4371  }
4372  break;
4373 
4374  case FIXMODPAIR: case FIXMODPAIRNIR:
4375 
4376  /* Multi-Scattering with Fixed Model Pair */
4377  /* -------------------------------------- */
4378 
4379  if (aerec != NULL && aerec->mode == ON) {
4380  *modmin = aerec->mod_min[ip] - 1;
4381  *modmax = aerec->mod_max[ip] - 1;
4382  *modrat = aerec->mod_rat[ip];
4383  } else {
4384  *modmin = input->aermodmin - 1;
4385  *modmax = input->aermodmax - 1;
4386  *modrat = input->aermodrat;
4387  }
4388 
4389  status = fixedmodpair(sensorID, wave, nwave, iwnir_s, iwnir_l, &geom, wv,
4390  *modmin, *modmax, *modrat, rhoa, eps);
4391  if (status == 0) {
4392  for (iw = 0; iw < nwave; iw++)
4393  La2[iw] = rhoa[iw] / radref[iw];
4394  }
4395 
4396  break;
4397 
4398  case FIXANGSTROM: case FIXANGSTROMNIR:
4399 
4400  if (input->aer_angstrom > -2) {
4401  angstrom = input->aer_angstrom;
4402  } else {
4403  int16_t year, day;
4404  double sec;
4405  unix2yds(l2rec->l1rec->scantime, &year, &day, &sec);
4406  angstrom = bin_climatology(input->aerbinfile, day, l1rec->lon[ip], l1rec->lat[ip], "angstrom");
4407  }
4408 
4409  if (angstrom > -2) {
4410 
4411  model_select_angstrom(angstrom, modmin, modmax, modrat);
4412 
4413  status = fixedmodpair(sensorID, wave, nwave, iwnir_s, iwnir_l, &geom, wv,
4414  *modmin, *modmax, *modrat, rhoa, eps);
4415  } else
4416  status = 1;
4417 
4418  if (status == 0) {
4419  for (iw = 0; iw < nwave; iw++)
4420  La2[iw] = rhoa[iw] / radref[iw];
4421  }
4422 
4423  break;
4424 
4425  case FIXAOT:
4426 
4427  /* Multi-Scattering with fixed AOT */
4428  /* ------------------------------- */
4429  if (aerec != NULL && aerec->mode == ON)
4430  aot = &aerec->taua[ip * Nbands];
4431  else
4432  aot = input->taua;
4433 
4434  status = fixedaot(sensorID, aot, wave, nwave, iwnir_s, iwnir_l, &geom, wv,
4435  modmin, modmax, modrat, rhoa, eps);
4436  if (status == 0) {
4437  for (iw = 0; iw < nwave; iw++)
4438  La2[iw] = rhoa[iw] / radref[iw];
4439  }
4440 
4441  break;
4442 
4443  default:
4444 
4445  /* Multi-Scattering with Fixed Model */
4446  /* --------------------------------- */
4447 
4448  *modmin = aer_opt - 1;
4449  *modmax = aer_opt - 1;
4450  *modrat = 0.0;
4451 
4452  if (*modmin < 0 || *modmin > input->naermodels - 1) {
4453  printf("Invalid aerosol option: %d\n", *modmin);
4454  exit(1);
4455  }
4456 
4457  /* convert input NIR aerosol radiance to relectance */
4458  rhoa[iwnir_l] = La1[iwnir_l] * radref[iwnir_l];
4459 
4460  /* get aerosol reflectance in visible for fixed model, and epsilon */
4461  status = fixedaer(sensorID, *modmin, wave, nwave, iwnir_s, iwnir_l,
4462  input->aermodels, input->naermodels,
4463  &geom, wv, rhoa, eps);
4464 
4465  /* convert aerosol relectance to radiance */
4466  if (status == 0)
4467  for (iw = 0; iw < nwave; iw++)
4468  La2[iw] = rhoa[iw] / radref[iw];
4469  break;
4470  }
4471 
4472  /* compute diffuse transmittance through aerosol and Rayleigh, if not yet computed */
4473  if (status == 0 && aer_opt != AERRHNIR && aer_opt != AERRHSWIR && aer_opt != AERRH && aer_opt != AERRHMSEPS && aer_opt != AERRHSM) {
4474 
4475  diff_tran(sensorID, wave, nwave, iwnir_l, &geom, wv, pr, taur,
4476  *modmin, *modmax, *modrat, rhoa, taua, t_sol, t_sen, taua_pred_min, taua_pred_max, 0);
4477  }
4478 
4479  /* transfer outputs per wavelength to outputs per band */
4480  for (iw = 0; iw < nwave; iw++) {
4481  La2_out [iw] = La2 [iw];
4482  taua_out [iw] = taua [iw];
4483  t_sol_out[iw] = t_sol[iw];
4484  t_sen_out[iw] = t_sen[iw];
4485  }
4486 
4487  /* model numbers are reported as 1-based */
4488  *modmin = *modmin + 1;
4489  *modmax = *modmax + 1;
4490  *modmin2 = *modmin2 + 1;
4491  *modmax2 = *modmax2 + 1;
4492 
4493  free(rhoa);
4494  free(radref);
4495  free(F0);
4496  free(taur);
4497  free(La1);
4498  free(La2);
4499  free(t_sol);
4500  free(t_sen);
4501  free(taua);
4502  free(taua_pred_min);
4503  free(taua_pred_max);
4504  free(geom.airmass);
4505  if ((evalmask & TRANSSPHER) != 0) {
4506  free(geom.airmass_plp);
4507  free(geom.airmass_sph);
4508  }
4509 
4510  return (status);
4511 }
4512 
4513 
4514 
4515 /* --------------------------------------------------------------- */
4516 /* get_angstrom.c - compute angstrom coefficient */
4517 /* */
4518 /* Inputs: */
4519 /* l2rec - level-2 structure containing one complete scan */
4520 /* after atmospheric correction. */
4521 /* band = band number 0-7 */
4522 /* Outputs: */
4523 /* angst - angstrom coefficient */
4524 /* */
4525 /* Written By: B. A. Franz, SAIC, August 2004 */
4526 /* */
4527 
4528 /* --------------------------------------------------------------- */
4529 void get_angstrom(l2str *l2rec, int band, float angst[]) {
4530  static int firstCall = 1;
4531  static int32_t ib2;
4532  static float wave2;
4533 
4534  int32_t ip;
4535  int32_t ib1;
4536  float wave1;
4537  float aot1;
4538  float aot2;
4539 
4540  l1str *l1rec = l2rec->l1rec;
4541  filehandle *l1file = l1rec->l1file;
4542 
4543  if (firstCall) {
4544  ib2 = windex(865.0, l1file->fwave, l1file->nbands);
4545  wave2 = l1file->fwave[ib2];
4546  firstCall = 0;
4547  }
4548 
4549  if (band < 0)
4550  ib1 = windex(443.0, l1file->fwave, l1file->nbands);
4551  else
4552  ib1 = band;
4553  wave1 = l1file->fwave[ib1];
4554 
4555  for (ip = 0; ip < l1rec->npix; ip++) {
4556  aot1 = l2rec->taua[ip * l1file->nbands + ib1];
4557  aot2 = l2rec->taua[ip * l1file->nbands + ib2];
4558  if (aot1 > 0.0 && aot2 > 0.0)
4559  angst[ip] = -log(aot1 / aot2) / log(wave1 / wave2);
4560  else if (aot1 == 0.0 && aot2 == 0.0)
4561  angst[ip] = 0.0;
4562  else {
4563  angst[ip] = BAD_FLT;
4564  l1rec->flags[ip] |= PRODFAIL;
4565  }
4566  // printf("angst = %f aot1= %f aot2= %f ib1=%d ib2=%d\n",angst[ip],aot1,aot2,ib1,ib2);
4567  }
4568 
4569  return;
4570 }
4571 
4572 
4573 /* --------------------------------------------------------------- */
4574 /* get_ms_epsilon.c - */
4575 
4576 /* --------------------------------------------------------------- */
4577 void get_ms_epsilon(l2str *l2rec, float eps[]) {
4578  int32_t ip, ipb1, ipb2;
4579 
4580  l1str *l1rec = l2rec->l1rec;
4581  filehandle *l1file = l1rec->l1file;
4582 
4583  for (ip = 0; ip < l1rec->npix; ip++) {
4584  ipb1 = ip * l1file->nbands + iwnir_s;
4585  ipb2 = ip * l1file->nbands + iwnir_l;
4586  if (l2rec->La[ipb2] > 0.0) {
4587  eps[ip] = l2rec->La[ipb1] / l2rec->La[ipb2]
4588  * l1rec->Fo[iwnir_l] / l1rec->Fo[iwnir_s];
4589  } else {
4590  /* "should" indicate ATMFAIL, but just to be safe */
4591  eps[ip] = BAD_FLT;
4592  l1rec->flags[ip] |= PRODFAIL;
4593  }
4594  }
4595 
4596  return;
4597 }
rhoaT_struct::rhoa
float rhoa
Definition: aerosol.c:3086
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
index
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 the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
MAX
#define MAX(A, B)
Definition: swl0_utils.h:26
aermodtab_struct::senz
float * senz
Definition: aerosol.c:166
comp_rhoa_ms_eps
int comp_rhoa_ms_eps(int32_t nwave, float wave[], geom_str *geom, float tau_iwnir_l, int32_t modl, float tau_pred[], float rho_pred[])
Definition: aerosol.c:1320
q
data_t q
Definition: decode_rs.h:74
ss_to_ms_coef
void ss_to_ms_coef(int modnum, geom_str *geom, float **a, float **b, float **c)
Definition: aerosol.c:656
MIN
#define MIN(x, y)
Definition: rice.h:169
python.test.test_tmatrix.epsilon
int epsilon
Definition: test_tmatrix.py:12
r
int r
Definition: decode_rs.h:73
pp_airmass
float pp_airmass(float theta)
Definition: airmass.c:14
ms_eps_coef
void ms_eps_coef(int modnum, int32_t iwnir_l, float wave[], geom_str *geom, float **a, float **b, float **c, float **d, float **e)
Definition: aerosol.c:931
AERRH
#define AERRH
Definition: l12_parms.h:34
ahmad_atm_corr
int ahmad_atm_corr(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *geom, float wv, float rhoa[], int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir, float tau_pred_max[], float tau_pred_min[], float rho_pred_max[], float rho_pred_min[], float tau_aer[], float rho_aer[])
Definition: aerosol.c:1696
aermod_struct::d2phase
float ** d2phase
Definition: aerosol.c:112
geom_strdef::senz
float * senz
Definition: aerosol.c:196
day
int32_t day
Definition: atrem_corl1v2.h:308
aermodtab_struct::sensorID
int32_t sensorID
Definition: aerosol.c:151
status
int status
Definition: l1_czcs_hdf.c:32
model_select_wang
int model_select_wang(int32_t sensorID, float wave[], int32_t nwave, int32_t nmodel, int32_t mindx[], geom_str *geom, float wv, float rhoa[], int32_t iwnir_s, int32_t iwnir_l, int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
Definition: aerosol.c:2407
rhoaT_struct::eps
float eps
Definition: aerosol.c:3087
FIXMODPAIRNIR
#define FIXMODPAIRNIR
Definition: l12_parms.h:27
geom_strdef::airmass_sph
float * airmass_sph
Definition: aerosol.c:203
f1
float f1(float x)
aeroob_cf
float aeroob_cf(int modnum, geom_str *geom)
Definition: aerosol.c:2124
INDEX
#define INDEX(iw, isol, iphi, isen)
Definition: aerosol.c:642
mapgen_overlay.ax
ax
Definition: mapgen_overlay.py:179
rhoas_to_rhoa
void rhoas_to_rhoa(int32_t sensorID, int modnum, geom_str *geom, float wv, float rhoas[], float wave[], int32_t nwave, int iw1, int iw2, float rhoa[])
Definition: aerosol.c:2268
AERRHSM
#define AERRHSM
Definition: l12_parms.h:37
aermod_struct::ams_all
float * ams_all
Definition: aerosol.c:99
bin_climatology
float bin_climatology(char *l3file, int32_t day, float lon, float lat, char *prodname)
Definition: bin_climatology.c:287
compalphaT
int compalphaT(alphaTstr *x, alphaTstr *y)
Definition: aerosol.c:2557
viirs_ler_luts::pr
real pr
Definition: viirs_ler_luts.f95:38
NULL
#define NULL
Definition: decode_rs.h:63
F0
map< string, float > F0
Definition: DDSensor.cpp:39
alloc_aermodstr
aermodstr * alloc_aermodstr(int nbands, int nscatt, int nphi, int nsolz, int nsenz, int ntheta)
Definition: aerosol.c:116
aermod_struct::bcost
float * bcost
Definition: aerosol.c:95
FIXANGSTROMNIR
#define FIXANGSTROMNIR
Definition: l12_parms.h:29
aermod_struct::name
char name[32]
Definition: aerosol.c:81
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
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
aermodtab_struct::nmodel
int32_t nmodel
Definition: aerosol.c:155
aermod_struct::ems_all
float * ems_all
Definition: aerosol.c:103
pi
const double pi
Definition: VcstCmnConsts.h:473
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
aermod_struct::dtran_a
float ** dtran_a
Definition: aerosol.c:107
ahmadaer
int ahmadaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *geom, float wv, float rhoa[], int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir, float tau_pred_min[], float tau_pred_max[])
Definition: aerosol.c:2936
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
ON
#define ON
Definition: l1.h:43
is
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 offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat rather resolving GSFcd01518 Added a LogReport Changed to create the Average Temperatures vdata even if the number of scans is
Definition: HISTORY.txt:301
AERWHITE
#define AERWHITE
Definition: l12_parms.h:22
TRANSSPHER
#define TRANSSPHER
Definition: l1.h:83
Nbands
@ Nbands
Definition: hybrid.c:24
computed
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor and incalculable values of SV or BB averages was moved outside the loop over frames in Emissive_Cal c since none of these quantities are frame dependent Initialization of b1 and XMS values in Preprocess c routine Process_OBCENG_Emiss was moved inside the detector loops The code was altered so that if up to five scans are dropped between the leading middle or middle trailing the leading or trailing granule will still be used in emissive calibration to form a cross granule average QA bits and are set for a gap between the leading middle and middle trailing granules respectively This may in rare instances lead to a change in emissive calibration coefficients for scans at the beginning or end of a granule A small bug in the Band correction algorithm was corrected an uncertainty value was being checked against an upper bound whereas the proper quantity to be checked was the corresponding which is the product of the Band radiance times the ratio of the Band to Band scaling factors times the LUT correction value for that detector In addition a new LUT which allows for a frame offset with regard to the Band radiance was added A LUT which switches the correction off or on was also added Changes which do not affect scientific output of the the pixel is flagged with the newly created flag and the number of pixels for which this occurs is counted in the QA_common table The array of b1s in Preprocess c was being initialized to outside the loop over which meant that if b1 could not be computed
Definition: HISTORY.txt:627
aermodtab_struct
Definition: aerosol.c:150
model_taua
void model_taua(int32_t sensorID, int modnum, float wave[], int32_t nwave, int32_t iwnir_l, float rhoa[], geom_str *geom, float wv, float taua[])
Definition: aerosol.c:2613
aerosol
int aerosol(l2str *l2rec, int32_t aer_opt_in, aestr *aerec, int32_t ip, float wave[], int32_t nwave, int32_t iwnir_s_in, int32_t iwnir_l_in, float F0_in[], float La1_in[], float La2_out[], float t_sol_out[], float t_sen_out[], float *eps, float taua_out[], int32_t *modmin, int32_t *modmax, float *modrat, int32_t *modmin2, int32_t *modmax2, float *modrat2)
Definition: aerosol.c:3922
spacetrack_tles.year
year
Definition: spacetrack_tles.py:170
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
AERRHNIR
#define AERRHNIR
Definition: l12_parms.h:24
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
aermod_struct::albedo
float * albedo
Definition: aerosol.c:89
model_phase
float * model_phase(int modnum, geom_str *geom)
Definition: aerosol.c:2001
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
im
int32_t im
Definition: atrem_corl1.h:161
input
instr * input
Definition: main_l2binmatch.cpp:27
geom_strdef::gmult
int gmult
Definition: aerosol.c:194
aermodtab_struct::nsolz
int32_t nsolz
Definition: aerosol.c:156
l12_proto.h
PI
#define PI
Definition: l3_get_org.c:6
f2
float f2(float y)
aermod_struct::lnphase
float ** lnphase
Definition: aerosol.c:111
aermod_struct::phase
float ** phase
Definition: aerosol.c:91
aermodtab_struct::dtran_nwave
int32_t dtran_nwave
Definition: aerosol.c:160
geom
geom_str geom
Definition: aerosol.c:206
FIXMODPAIR
#define FIXMODPAIR
Definition: l12_parms.h:26
ahmad_atm_corr_lin
int ahmad_atm_corr_lin(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *geom, float wv, float rhoa[], int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir, float tau_pred_max[], float tau_pred_min[], float rho_pred_max[], float rho_pred_min[], float tau_aer[], float rho_aer[])
Definition: aerosol.c:1847
diff_tran
void diff_tran(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_l, geom_str *geom, float wv, float pr, float taur[], int32_t modmin, int32_t modmax, float modrat, float rhoa[], float taua[], float tsol[], float tsen[], float tauamin[], float tauamax[], int taua_opt)
Definition: aerosol.c:2816
f
double precision function f(R1)
Definition: tmd.lp.f:1454
epsilonT_struct::eps_obs
float eps_obs
Definition: aerosol.c:186
AERWANG
#define AERWANG
Definition: l12_parms.h:23
rhaer
int rhaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, geom_str *geom, float wv, float rh, float pr, float taur[], float rhoa[], int32_t *modmin1, int32_t *modmax1, float *modrat1, int32_t *modmin2, int32_t *modmax2, float *modrat2, float *eps, float taua[], float tsol[], float tsen[], float tg_sol_sm[], float tg_sen_sm[], float Lt_sm[], int32_t ip)
Definition: aerosol.c:3206
linterp
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
aermod_struct::angstrom
float * angstrom
Definition: aerosol.c:86
color_dtdb.path
string path
Definition: color_dtdb.py:221
hico.value_locate.x
list x
Definition: value_locate.py:64
aermodtab_struct::dtran_wave
float * dtran_wave
Definition: aerosol.c:171
alphaT_struct
Definition: aerosol.c:179
geom_strdef
Definition: aerosol.c:193
color_dtdb.y
y
Definition: color_dtdb.py:430
color_dtdb.bc
int bc
Definition: color_dtdb.py:185
aermod_struct::extc
float * extc
Definition: aerosol.c:90
aermodtab_struct::model
aermodstr ** model
Definition: aerosol.c:175
aermodtab_struct::dtran_airmass
float * dtran_airmass
Definition: aerosol.c:173
epsilonT_struct::modnum
int32_t modnum
Definition: aerosol.c:185
aermodtab_struct::dtran_theta
float * dtran_theta
Definition: aerosol.c:172
AERRHSWIR
#define AERRHSWIR
Definition: l12_parms.h:33
geom_strdef::solz
float * solz
Definition: aerosol.c:197
comp_epsilonT
int comp_epsilonT(epsilonTstr *x, epsilonTstr *y)
Definition: aerosol.c:1257
aermodtab_struct::nphi
int32_t nphi
Definition: aerosol.c:158
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
MAXAERMOD
#define MAXAERMOD
Definition: l12_parms.h:21
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
get_angstrom
void get_angstrom(l2str *l2rec, int band, float angst[])
Definition: aerosol.c:4529
aermodtab_struct::nsenz
int32_t nsenz
Definition: aerosol.c:157
unix2yds
void unix2yds(double usec, short *year, short *day, double *secs)
MAXMODEL
#define MAXMODEL
Definition: aerosol.c:53
aermod_struct::rh
float rh
Definition: aerosol.c:82
geom_strdef::csolz
float * csolz
Definition: aerosol.c:199
model_transmittance
void model_transmittance(int modnum, float wave[], int32_t nwave, float *theta, int gmult, float taua[], float dtran[])
Definition: aerosol.c:2722
ac
float * ac[MAX_BANDS]
Definition: pml_iop_tables.c:26
AERMUMM
#define AERMUMM
Definition: l12_parms.h:32
diff
subroutine diff(x, conec, n, dconecno, dn, dconecmk, units, u, inno, i, outno, o, input, deriv)
Definition: ffnet.f:205
ky_airmass
float ky_airmass(float theta)
Definition: airmass.c:4
aermodtab_struct::dtran_ntheta
int32_t dtran_ntheta
Definition: aerosol.c:161
geom_strdef::csenz
float * csenz
Definition: aerosol.c:200
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
geom_strdef::airmass_plp
float * airmass_plp
Definition: aerosol.c:202
geom_strdef::phi
float * phi
Definition: aerosol.c:198
aermod_struct::ccost
float * ccost
Definition: aerosol.c:96
aermodtab_struct::scatt
float * scatt
Definition: aerosol.c:168
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
AVIRIS
#define AVIRIS
Definition: sensorDefs.h:30
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
b1
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor b1
Definition: HISTORY.txt:576
aermod_struct::sd
int sd
Definition: aerosol.c:83
AERRHMSEPS
#define AERRHMSEPS
Definition: l12_parms.h:36
fresnel_coef
float fresnel_coef(float mu, float index)
Definition: aerosol.c:904
FIXANGSTROM
#define FIXANGSTROM
Definition: l12_parms.h:28
aermodtab_struct::wave
float * wave
Definition: aerosol.c:164
rhoa_to_rhoas
int rhoa_to_rhoas(int32_t sensorID, int modnum, geom_str *geom, float wv, float rhoa[], float wave[], int32_t nwave, int iw1, int iw2, float rhoas[])
Definition: aerosol.c:2205
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
STDPR
#define STDPR
Definition: l12_parms.h:85
rhoaT_struct
Definition: aerosol.c:3084
wangaer
int wangaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *geom, float wv, float rhoa[], int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir, float tauamin[], float tauamax[])
Definition: aerosol.c:2979
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
phase
subroutine phase
Definition: phase.f:2
model_select_franz
int model_select_franz(int32_t sensorID, float wave[], int32_t nwave, int32_t nmodel, int32_t mindx[], geom_str *geom, float wv, float rhoa[], int32_t iwnir_s, int32_t iwnir_l, int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
Definition: aerosol.c:3094
aermod_struct::acost
float * acost
Definition: aerosol.c:94
aermodtab_struct::nscatt
int32_t nscatt
Definition: aerosol.c:159
aermod_struct::dtran_b
float ** dtran_b
Definition: aerosol.c:108
AERRHMUMM
#define AERRHMUMM
Definition: l12_parms.h:35
rhoaT_struct::modnum
int32_t modnum
Definition: aerosol.c:3085
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
aermod_struct::dms_all
float * dms_all
Definition: aerosol.c:102
cmpfunc
int cmpfunc(const void *a, const void *b)
Definition: aerosol.c:233
model_epsilon
float * model_epsilon(int modnum, int32_t iwnir_l, float wave[], int32_t nwave, geom_str *geom)
Definition: aerosol.c:2328
fabs
#define fabs(a)
Definition: misc.h:93
fixedaot
int fixedaot(int32_t sensorID, float aot[], float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, geom_str *geom, float wv, int32_t *modmin, int32_t *modmax, float *modrat, float rhoa[], float *epsnir)
Definition: aerosol.c:3763
aermod_struct::bms_all
float * bms_all
Definition: aerosol.c:100
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
GAS_TRANS_TBL_BIT
#define GAS_TRANS_TBL_BIT
Definition: l12_parms.h:62
splint
subroutine splint(xa, ya, y2a, n, x, y)
Definition: atrem_app_refl_f90_cubeio.f:3270
BAD_INT
#define BAD_INT
Definition: genutils.h:23
strncpy
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 and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)
spectral_matching
int spectral_matching(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx1[], int32_t mindx2[], int32_t mindx3[], geom_str *geom, float wv, float rhoa[], float rho_aer[], int32_t *mod1_indx, int32_t *mod2_indx, float *weight, float tau_pred_min[], float tau_pred_max[], float tg_sol_sm[], float tg_sen_sm[], float Lt_sm[], int32_t ip)
Definition: aerosol.c:1425
comp_rhoa_ms_eps_lin
int comp_rhoa_ms_eps_lin(int32_t nwave, float wave[], geom_str *geom, float tau_iwnir_l, int32_t modl, float tau_pred[], float rho_pred[])
Definition: aerosol.c:1373
comp_rhoaT
int comp_rhoaT(rhoaTstr *x, rhoaTstr *y)
Definition: aerosol.c:3090
aermod_struct
Definition: aerosol.c:80
model_select_ahmad
void model_select_ahmad(int32_t nmodels, int32_t *mindx, float eps_pred[], float eps_obs, int32_t *modmin, int32_t *modmax, float *modrat)
Definition: aerosol.c:1261
first_deriv
float first_deriv(float x[], float y[], int n)
Definition: aerosol.c:246
alphaT_struct::angstrom
float angstrom
Definition: aerosol.c:181
rdsensorinfo
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
allocate2d_float
float ** allocate2d_float(size_t h, size_t w)
Allocate a two-dimensional array of type float of a given size.
Definition: allocate2d.c:125
model
int32_t model
Definition: atrem_corl1.h:132
for
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
aermodtab_struct::solz
float * solz
Definition: aerosol.c:165
model_select_angstrom
void model_select_angstrom(float angstrom, int32_t *modmin, int32_t *modmax, float *modrat)
Definition: aerosol.c:2561
AERRHMSEPS_lin
#define AERRHMSEPS_lin
Definition: l12_parms.h:38
c
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 as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
geom_strdef::airmass
float * airmass
Definition: aerosol.c:201
AERWANGNIR
#define AERWANGNIR
Definition: l12_parms.h:25
aermod_struct::cms_all
float * cms_all
Definition: aerosol.c:101
aermodtab_struct::phi
float * phi
Definition: aerosol.c:167
smaer
int smaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, int32_t nmodels, int32_t mindx1[], int32_t mindx2[], int32_t mindx3[], geom_str *geom, float wv, float rhoa[], float rho_aer[], int32_t *modmin, int32_t *modmax, float *modrat, float tau_pred_min[], float tau_pred_max[], float tg_sol_sm[], float tg_sen_sm[], float Lt_sm[], int32_t ip)
Definition: aerosol.c:2903
i
int i
Definition: decode_rs.h:71
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
model_taua_mseps
int model_taua_mseps(int32_t modl, float wave[], int32_t nwave, int32_t iwnir_l, float rhoa[], geom_str *geom, float tau_pred[])
Definition: aerosol.c:2679
aermodtab_struct::nwave
int32_t nwave
Definition: aerosol.c:154
a
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
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:97
aeroob
float aeroob(int32_t sensorID, int32_t iw, float airmass, float cf, float wv)
Definition: aerosol.c:2155
MODISA
#define MODISA
Definition: sensorDefs.h:19
FIXAOT
#define FIXAOT
Definition: l12_parms.h:30
viirs_ler_luts::p1
real p1
Definition: viirs_ler_luts.f95:37
get_ms_epsilon
void get_ms_epsilon(l2str *l2rec, float eps[])
Definition: aerosol.c:4577
AERWANGSWIR
#define AERWANGSWIR
Definition: l12_parms.h:31
load_aermod
int load_aermod(int32_t sensorID, float wave[], int32_t nwave, char *aermodfile, char models[MAXAERMOD][32], int32_t nmodels)
Definition: aerosol.c:322
epsilonT_struct
Definition: aerosol.c:184
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
fixedmodpair
int fixedmodpair(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, geom_str *geom, float wv, int32_t modmin, int32_t modmax, float modrat, float rhoa[], float *eps)
Definition: aerosol.c:3687
alphaT_struct::modnum
int32_t modnum
Definition: aerosol.c:180
fixedaer
int fixedaer(int32_t sensorID, int32_t modnum, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l, char models[MAXAERMOD][32], int32_t nmodels, geom_str *geom, float wv, float rhoa[], float *epsnir)
Definition: aerosol.c:3632
color_dtdb.min1
float min1
Definition: color_dtdb.py:283