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