53 #define MAXMODEL MAXAERMOD
61 static double radeg =
RADEG;
62 static float p0 =
STDPR;
64 static int have_ms = 0;
65 static int have_rh = 0;
66 static int have_sd = 0;
67 static int use_rh = 0;
70 static int32_t Maxband;
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));
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));
172 static aermodtabstr *aertab =
NULL;
191 static int loaded = 0;
192 static int interpol = 0;
193 static int32_t *iwatab;
194 static int32_t *iwdtab;
196 static int32_t iwnir_s = -1;
197 static int32_t iwnir_l = -1;
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;
208 static float airmass;
210 static int32_t evalmask = 0;
211 static int32_t aer_opt = 0;
212 static float airmass_plp;
213 static float airmass_sph;
216 if (*(
double*)
a > *(
double*)
b)
return 1;
217 else if (*(
double*)
a < *(
double*)
b)
return -1;
229 float a1, a2, a3, a4, a5, a6, d1;
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);
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];
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);
274 int32 dims[H4_MAX_VAR_DIMS];
276 int32
start[4] = {0, 0, 0, 0};
278 int32_t mwave, msolz, msenz, mphi, mscatt, dtran_mwave, dtran_mtheta;
285 char name [H4_MAX_NC_NAME] =
"";
286 char sdsname[H4_MAX_NC_NAME] =
"";
287 char file [FILENAME_MAX] =
"";
288 char path [FILENAME_MAX] =
"";
290 int iw,
im,
is, iwbase,
i;
291 static int firstCall = 1;
293 if (firstCall == 1) {
294 if ((iwatab = (int32_t *) calloc(nwave,
sizeof (int32_t))) ==
NULL) {
295 printf(
"Unable to allocate space for iwatab.\n");
298 if ((iwdtab = (int32_t *) calloc(nwave,
sizeof (int32_t))) ==
NULL) {
299 printf(
"Unable to allocate space for iwdtab.\n");
305 printf(
"Loading aerosol models from %s\n", aermodfile);
307 for (
im = 0;
im < nmodels + 1;
im++) {
312 strcat(
file, aermodfile);
315 strcat(
file,
".hdf");
318 sd_id = SDstart(
file, DFACC_RDONLY);
320 printf(
"-E- %s: Error opening file %s.\n",
328 strcat(
file, aermodfile);
329 strcat(
file,
"_default.hdf");
332 sd_id = SDstart(
file, DFACC_RDONLY);
334 printf(
"-E- %s: Error opening file %s.\n",
342 status = SDreadattr(sd_id, SDfindattr(sd_id,
"Number of Wavelengths"),
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);
359 int32_t nwave = mwave;
360 int32_t nsolz = msolz;
361 int32_t nsenz = msenz;
363 int32_t nscatt = mscatt;
364 int32_t dtran_nwave = dtran_mwave;
365 int32_t dtran_ntheta = dtran_mtheta;
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);
373 printf(
"Number of Diffuse Transmittance Wavelengths %d\n", dtran_nwave);
374 printf(
"Number of Diffuse Transmittance Zenith Angles %d\n", dtran_ntheta);
378 if ((aertab = (aermodtabstr *) calloc(1,
sizeof (aermodtabstr))) ==
NULL) {
379 printf(
"Unable to allocate space for aerosol table.\n");
383 aertab->nmodel = nmodels;
384 aertab->nwave = nwave;
385 aertab->nsolz = nsolz;
386 aertab->nsenz = nsenz;
388 aertab->nscatt = nscatt;
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));
396 aertab->dtran_nwave = dtran_nwave;
397 aertab->dtran_ntheta = dtran_ntheta;
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));
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);
409 for (
i = 0;
i < nmodels + 1;
i++) {
411 printf(
"Unable to allocate space for aerosol model %d.\n",
im);
419 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
423 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
426 status = SDendaccess(sds_id);
430 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
434 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
437 status = SDendaccess(sds_id);
441 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
445 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
448 status = SDendaccess(sds_id);
452 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
456 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
459 status = SDendaccess(sds_id);
463 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
467 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
470 status = SDendaccess(sds_id);
473 strcpy(sdsname,
"dtran_wave");
474 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
478 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
481 status = SDendaccess(sds_id);
484 strcpy(sdsname,
"dtran_theta");
485 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
489 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
492 status = SDendaccess(sds_id);
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");
513 strncpy(aertab->model[
im]->name,
"default", 32);
515 status = SDreadattr(sd_id, SDfindattr(sd_id,
"Relative Humidity"), &rh);
517 aertab->model[
im]->rh = rh;
520 aertab->model[
im]->rh = -1.0;
522 status = SDreadattr(sd_id, SDfindattr(sd_id,
"Size Distribution"), &
sd);
524 aertab->model[
im]->sd =
sd;
527 aertab->model[
im]->sd = -1;
529 strcpy(sdsname,
"albedo");
530 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
534 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
537 status = SDendaccess(sds_id);
541 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
545 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
548 status = SDendaccess(sds_id);
565 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
569 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
572 status = SDendaccess(sds_id);
576 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
580 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
583 status = SDendaccess(sds_id);
587 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
591 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
594 status = SDendaccess(sds_id);
598 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
602 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
605 status = SDendaccess(sds_id);
609 if (SDnametoindex(sd_id,
"ams_all") !=
FAIL) {
613 strcpy(sdsname,
"ams_all");
614 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
618 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
621 status = SDendaccess(sds_id);
624 strcpy(sdsname,
"bms_all");
625 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
629 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
632 status = SDendaccess(sds_id);
635 strcpy(sdsname,
"cms_all");
636 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
640 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
643 status = SDendaccess(sds_id);
646 strcpy(sdsname,
"dms_all");
647 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
651 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
654 status = SDendaccess(sds_id);
657 strcpy(sdsname,
"ems_all");
658 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
662 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
665 status = SDendaccess(sds_id);
670 strcpy(sdsname,
"dtran_a");
671 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
675 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
678 status = SDendaccess(sds_id);
681 strcpy(sdsname,
"dtran_b");
682 sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
686 printf(
"-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
file);
689 status = SDendaccess(sds_id);
697 iwbase =
windex(865, aertab->wave, aertab->nwave);
698 for (iw = 0; iw < aertab->nwave; iw++) {
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]);
703 aertab->model[
im]->angstrom[iwbase] = aertab->model[
im]->angstrom[iwbase - 1];
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]);
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);
713 &aertab->model[
im]->lnphase[iw][0],
717 &aertab->model[
im]->d2phase[iw][0]);
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);
728 aertab->nmodel = nmodels;
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]);
742 if (aertab->nwave != nwave)
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;
755 if (imm50 >= 0 && imc50 >= 0 && imc70 >= 0 && imt90 >= 0 && imt99 >= 0) {
757 printf(
"\nM. Wang model cross-over correction enabled.\n");
765 #define INDEX(iw,isol,iphi,isen) (iw*aertab->nsolz*aertab->nphi*aertab->nsenz + isol*aertab->nphi*aertab->nsenz + iphi*aertab->nsenz + isen)
780 static float lastsolz = -999.;
781 static float lastsenz = -999.;
782 static float lastphi = -999.;
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;
795 static int *isolz1, *isolz2;
796 static int *isenz1, *isenz2;
797 static int *iphi1, *iphi2;
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;
808 static int firstCall = 1;
811 if (firstCall == 1) {
814 if ((a_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
815 printf(
"Unable to allocate space for a_coef.\n");
818 if ((b_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
819 printf(
"Unable to allocate space for rhoa.\n");
822 if ((c_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
823 printf(
"Unable to allocate space for rhoa.\n");
828 if ((
geom->gmult == 0) || (interpol == 1)) {
833 isolz1 = &isolz1_cnst;
834 isolz2 = &isolz2_cnst;
835 isenz1 = &isenz1_cnst;
836 isenz2 = &isenz2_cnst;
842 if (((p_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
844 ((q_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
846 ((r_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
848 printf(
"Unable to allocate space for p, q, r weights.\n");
851 if (((isolz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
853 ((isenz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
855 ((iphi1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
857 printf(
"Unable to allocate space for interp indicies 1.\n");
860 if (((isolz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
862 ((isenz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
864 ((iphi2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
866 printf(
"Unable to allocate space for interp indicies 2.\n");
881 if ((
geom->solz[0] != lastsolz) || (
geom->senz[0] != lastsenz) ||
882 (
geom->phi[0] != lastphi)) {
883 for (
im = 0;
im < aertab->nmodel;
im++)
886 for (iw = 0; iw < aertab->nwave; iw++) {
889 for (
i = 0;
i < aertab->nsolz;
i++) {
890 if (
geom->solz[ig] < aertab->solz[
i])
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]]);
902 for (
i = 0;
i < aertab->nsenz;
i++) {
903 if (
geom->senz[ig] < aertab->senz[
i])
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]]);
916 for (
i = 0;
i < aertab->nphi;
i++) {
917 if (aphi < aertab->phi[
i])
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]]);
927 if (gmult == 0)
break;
931 lastsolz =
geom->solz[0];
932 lastsenz =
geom->senz[0];
933 lastphi =
geom->phi[0];
940 for (iw = 0; iw < aertab->nwave; iw++) {
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])];
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])];
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])];
961 a_coef[
im][iw] = (1. - px)*(1. - rx) * as000 + px * rx * as101
962 + (1. - px) * rx * as001 + px * (1. - rx) * as100;
964 b_coef[
im][iw] = (1. - px)*(1. - rx) * ai000 + px * rx * ai101
965 + (1. - px) * rx * ai001 + px * (1. - qx)*(1. - rx) * ai100;
967 c_coef[
im][iw] = (1. - px)*(1. - rx) * ac000 + px * rx * ac101
968 + (1. - px) * rx * ac001 + px * (1. - qx)*(1. - rx) * ac100;
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])];
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])];
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])];
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;
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;
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;
1016 *
a = &a_coef[modnum][0];
1017 *
b = &b_coef[modnum][0];
1018 *
c = &c_coef[modnum][0];
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);
1034 return (r2 * (q1 * q1 + 1.0) / 2.0);
1055 float **
a,
float **
b,
float **
c,
float **d,
float **e)
1057 static float lastsolz = -999.;
1058 static float lastsenz = -999.;
1059 static float lastphi = -999.;
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;
1076 static int *isolz1, *isolz2;
1077 static int *isenz1, *isenz2;
1078 static int *iphi1, *iphi2;
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;
1090 static int firstCall = 1;
1092 if (firstCall == 1) {
1095 if ((a_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1096 printf(
"Unable to allocate space for a_coef.\n");
1099 if ((b_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1100 printf(
"Unable to allocate space for b_coef.\n");
1103 if ((c_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1104 printf(
"Unable to allocate space for c_coef.\n");
1107 if ((d_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1108 printf(
"Unable to allocate space for d_coef.\n");
1112 if ((e_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1113 printf(
"Unable to allocate space for e_coef.\n");
1119 if ((
geom->gmult == 0) || (interpol == 1)) {
1124 isolz1 = &isolz1_cnst;
1125 isolz2 = &isolz2_cnst;
1126 isenz1 = &isenz1_cnst;
1127 isenz2 = &isenz2_cnst;
1128 iphi1 = &iphi1_cnst;
1129 iphi2 = &iphi2_cnst;
1132 if (((p_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1134 ((q_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1136 ((r_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1138 printf(
"Unable to allocate space for p, q, r weights.\n");
1141 if (((isolz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1143 ((isenz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1145 ((iphi1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1147 printf(
"Unable to allocate space for interp indicies 1.\n");
1150 if (((isolz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1152 ((isenz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1154 ((iphi2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1156 printf(
"Unable to allocate space for interp indicies 2.\n");
1171 if (
geom->solz[0] != lastsolz ||
geom->senz[0] != lastsenz ||
1172 geom->phi[0] != lastphi) {
1173 for (
im = 0;
im < aertab->nmodel;
im++)
1176 for (iw = 0; iw < aertab->nwave; iw++) {
1179 for (
i = 0;
i < aertab->nsolz;
i++) {
1180 if (
geom->solz[ig] < aertab->solz[
i])
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]]);
1192 for (
i = 0;
i < aertab->nsenz;
i++) {
1193 if (
geom->senz[ig] < aertab->senz[
i])
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]]);
1206 for (
i = 0;
i < aertab->nphi;
i++) {
1207 if (aphi < aertab->phi[
i])
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]]);
1217 if (gmult == 0)
break;
1221 lastsolz =
geom->solz[0];
1222 lastsenz =
geom->senz[0];
1223 lastphi =
geom->phi[0];
1233 for (iw = 0; iw < aertab->nwave; iw++) {
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])];
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])];
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])];
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])];
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])];
1264 a_coef[
im][iw] = (1. - px)*(1. - rx) * as000 + px * rx * as101
1265 + (1. - px) * rx * as001 + px * (1. - rx) * as100;
1267 b_coef[
im][iw] = (1. - px)*(1. - rx) * ai000 + px * rx * ai101
1268 + (1. - px) * rx * ai001 + px * (1. - qx)*(1. - rx) * ai100;
1270 c_coef[
im][iw] = (1. - px)*(1. - rx) * ac000 + px * rx * ac101
1271 + (1. - px) * rx * ac001 + px * (1. - qx)*(1. - rx) * ac100;
1273 d_coef[
im][iw] = (1. - px)*(1. - rx) * ad000 + px * rx * ad101
1274 + (1. - px) * rx * ad001 + px * (1. - qx)*(1. - rx) * ad100;
1276 e_coef[
im][iw] = (1. - px)*(1. - rx) * ae000 + px * rx * ae101
1277 + (1. - px) * rx * ae001 + px * (1. - qx)*(1. - rx) * ae100;
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])];
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])];
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])];
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])];
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])];
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;
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;
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;
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;
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;
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");
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];
1381 return (
x->eps_obs <
y->eps_obs ? -1 : 1);
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) {
1392 for (
im = 0;
im < nmodels;
im++) {
1393 epsilonT[
im].eps_obs = eps_pred[
im];
1394 epsilonT[
im].modnum =
im;
1400 qsort(epsilonT, nmodels,
sizeof (epsilonTstr),
1405 for (
im = 0;
im < nmodels;
im++) {
1406 if (eps_obs < epsilonT[
im].eps_obs)
1417 im1 =
MAX(
MIN(
im - 1, nmodels - 1), 0);
1418 im2 =
MAX(
MIN(
im, nmodels - 1), 0);
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);
1430 if (*modmin == *modmax)
1444 float tau_iwnir_l, int32_t modl,
float tau_pred[],
float rho_pred[])
1446 float *
ac, *
bc, *cc, *dc, *ec;
1447 float ext_modl[nwave];
1448 float lg_tau_pred[nwave];
1449 float lg_rho_pred[nwave];
1461 for (iw = 0; iw < nwave; iw++) {
1463 ext_modl[iw] = aertab->model[modl]->extc[iwtab];
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]);
1474 for (iw = 0; iw < nwave; 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]);
1497 float tau_iwnir_l, int32_t modl,
float tau_pred[],
float rho_pred[])
1499 float *
ac, *
bc, *cc, *dc, *ec;
1500 float ext_modl[nwave];
1501 float lg_tau_pred[nwave];
1502 float lg_rho_pred[nwave];
1514 for (iw = 0; iw <= nwave; iw++) {
1516 ext_modl[iw] = aertab->model[modl]->extc[iwtab];
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]);
1527 for (iw = 0; iw <= nwave; 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]);
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) {
1558 float *
ac, *
bc, *cc, *dc, *ec;
1559 float ext_iwnir_cal;
1560 double ax, bx, cx, fx;
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*));
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));
1580 float *ext_all_wav = (
float*) malloc(nwave *
sizeof (
float));
1582 int iwtab,
im, modl;
1583 float lg_tau_iwnir_cal;
1585 double *diff = malloc(nwave *
sizeof (
double));
1586 double *chi = malloc(nmodels *
sizeof (
double));
1587 double *chi_temp = malloc(nmodels *
sizeof (
double));
1591 int min1_indx, min2_indx;
1594 float *noise = malloc(nwave *
sizeof (
float));
1595 static float scaled_lt;
1596 int iwtab_l, iwtab_cal;
1602 float ltSnrFitCoef[16][2] = {
1603 {0.05499859, 0.00008340},
1604 {0.02939470, 0.00009380},
1605 {0.11931482, 0.00008195},
1606 {0.01927545, 0.00009450},
1607 {0.01397522, 0.00010040},
1608 {0.01139088, 0.00016480},
1609 {0.08769538, 0.00007000},
1610 {0.10406925, 0.00008533},
1611 {0.00496291, 0.00014050},
1612 {0.00427147, 0.00013160},
1613 {0.00416994, 0.00021250},
1614 {0.04055895, 0.000197550},
1615 {0.00312263, 0.00018600},
1616 {0.07877732, 0.00049940},
1617 {100000.1, 100000000.1},
1618 {0.00628912, 0.00021160}
1622 float snr_scale[16] = {
1623 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 2, 2, 2
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]);
1632 for (
i = 0;
i < nwave;
i++)
1639 int32_t *mindx = (int32_t*) malloc(nmodels *
sizeof (int32_t));
1643 if (mindx1[0] == mindx2[0]) {
1645 for (
i = 0;
i < 10;
i++)
1646 mindx[
i] = mindx1[
i];
1649 for (
i = 0;
i < 10;
i++) {
1650 mindx[
i] = mindx1[
i];
1651 mindx[
i + 10] = mindx2[
i];
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];
1660 for (
i = 0;
i < 10;
i++) {
1661 mindx[
i] = mindx3[
i];
1662 mindx[
i + 10] = mindx1[
i];
1663 mindx[
i + 20] = mindx2[
i];
1683 for (
im = 0;
im < nmodels;
im++) {
1691 iwtab_cal = iwatab[iwnir_cal];
1694 iwtab_cal = iwatab[iwnir_cal];
1697 printf(
"Error: Sensor unkown %d\n",
sensorID);
1705 iwtab_l = iwatab[iwnir_l];
1708 ax = (double)
ac[iwtab_cal] - log((
double) rhoa[iwnir_cal]);
1709 bx = (double)
bc[iwtab_cal];
1710 cx = (double) cc[iwtab_cal];
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);
1726 ext_iwnir_cal = aertab->model[modl]->extc[iwtab_cal];
1727 for (iw = 0; iw <= iwtab_l - 1; 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]);
1740 for (iw = iwnir_s; iw <= iwnir_l; iw++) {
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]);
1746 chi[
im] = sqrt(diff_2 / (iwnir_l - iwnir_s));
1756 for (
i = 0;
i < nmodels;
i++)
1757 chi_temp[
i] = chi[
i];
1759 qsort(chi_temp, nmodels,
sizeof (chi_temp[0]),
cmpfunc);
1763 for (
i = 0;
i < nmodels;
i++) {
1766 else if (chi[
i] == min2)
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]));
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];
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]);
1800 free(rho_all_wav_pred);
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[]) {
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];
1838 float lg_tau_iwnir_l;
1839 int iwtab_l, iwtab_s;
1842 static float eps_obs;
1851 eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
1861 for (
im = 0;
im < nmodels;
im++) {
1872 iwtab_l = iwatab[iwnir_l];
1873 iwtab_s = iwatab[iwnir_s];
1875 ax = (double)
ac[iwtab_l] - log((
double) rhoa[iwnir_l]);
1876 bx = (double)
bc[iwtab_l];
1877 cx = (double) cc[iwtab_l];
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);
1890 ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
1891 ext_iwnir_s = aertab->model[modl]->extc[iwtab_s];
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]);
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);
1904 rho_iwnir_s_pred[
im] = exp(lg_rho_iwnir_s_pred[
im]);
1908 eps_pred[
im] = rho_iwnir_s_pred[
im] / rhoa[iwnir_l];
1944 *modmin = mindx[im1];
1945 *modmax = mindx[im2];
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];
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[]) {
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];
1989 float lg_tau_iwnir_l;
1990 int iwtab_l, iwtab_s;
1993 static float eps_obs;
2002 eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
2012 for (
im = 0;
im < nmodels;
im++) {
2023 iwtab_l = iwatab[iwnir_l];
2024 iwtab_s = iwatab[iwnir_s];
2026 ax = (double)
ac[iwtab_l]-(
double) rhoa[iwnir_l];
2027 bx = (double)
bc[iwtab_l];
2028 cx = (double) cc[iwtab_l];
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);
2041 ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
2042 ext_iwnir_s = aertab->model[modl]->extc[iwtab_s];
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]);
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");
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);
2061 rho_iwnir_s_pred[
im] = (lg_rho_iwnir_s_pred[
im]);
2065 eps_pred[
im] = rho_iwnir_s_pred[
im] / rhoa[iwnir_l];
2094 *modmin = mindx[im1];
2095 *modmax = mindx[im2];
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];
2126 static float nw = 1.334;
2129 static float lastsolz = -999.;
2130 static float lastsenz = -999.;
2131 static float lastphi = -999.;
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;
2139 float phase1, phase2;
2142 if (firstCall == 1) {
2145 if ((
phase[
im] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
2146 printf(
"Unable to allocate space for phase.\n");
2151 if ((
geom->gmult == 0) || (interpol == 1)) {
2153 scatt1 = &scatt1_cnst;
2154 scatt2 = &scatt2_cnst;
2155 fres1 = &fres1_cnst;
2156 fres2 = &fres2_cnst;
2160 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2161 printf(
"Unable to allocate space for scatt1.\n");
2165 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2166 printf(
"Unable to allocate space for scatt2.\n");
2170 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2171 printf(
"Unable to allocate space for fres1.\n");
2175 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2176 printf(
"Unable to allocate space for fres2.\n");
2188 if ((
geom->solz[0] != lastsolz) || (
geom->senz[0] != lastsenz) ||
2189 (
geom->phi[0] != lastphi)) {
2193 for (
im = 0;
im < aertab->nmodel;
im++)
2197 for (iw = 0; iw < aertab->nwave; 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);
2203 MAX(-
geom->csenz[ig] *
geom->csolz[ig] + temp, -1.0)) * radeg;
2205 MIN(
geom->csenz[ig] *
geom->csolz[ig] + temp, 1.0)) * radeg;
2210 if (gmult == 0)
break;
2213 lastsolz =
geom->solz[0];
2214 lastsenz =
geom->senz[0];
2215 lastphi =
geom->phi[0];
2224 for (iw = 0; iw < aertab->nwave; iw++) {
2227 &aertab->model[
im]->lnphase[iw][0],
2228 &aertab->model[
im]->d2phase[iw][0],
2229 aertab->nscatt, scatt1[ig], &phase1);
2231 &aertab->model[
im]->lnphase[iw][0],
2232 &aertab->model[
im]->d2phase[iw][0],
2233 aertab->nscatt, scatt2[ig], &phase2);
2236 exp(phase2)*(fres1[ig] + fres2[ig]);
2240 return (&
phase[modnum][0]);
2249 static int firstCall = 1;
2254 float rhoas1, rhoas2;
2259 iw1 =
windex(765, aertab->wave, aertab->nwave);
2260 iw2 =
windex(865, aertab->wave, aertab->nwave);
2261 if (iw1 == iw2) iw1--;
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]);
2293 static int firstCall = 1;
2311 printf(
"\nLoading water-vapor correction coefficients.\n");
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;
2329 float rhoa[],
float wave[], int32_t nwave,
int iw1,
int iw2,
float rhoas[]) {
2330 float *
ac, *
bc, *cc;
2342 for (iw = iw1; iw <= iw2; iw++) {
2343 ig = iw *
geom->gmult;
2344 if (rhoa[iw] < 1.e-20)
2345 rhoas[iw] = rhoa[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]));
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);
2362 if (!isfinite(rhoas[iw]) || rhoas[iw] < 1.e-20) {
2376 for (iw = iw1; iw <= iw2; iw++)
2377 rhoas[iw] = rhoa[iw];
2392 float rhoas[],
float wave[], int32_t nwave,
int iw1,
int iw2,
float rhoa[]) {
2393 float *
ac, *
bc, *cc;
2404 for (iw = iw1; iw <= iw2; iw++) {
2405 ig = iw *
geom->gmult;
2410 if (rhoas[iw] < 1.e-20)
2411 rhoa[iw] = rhoas[iw];
2418 rhoa[iw] = exp(
a +
b * lnrhoas +
c * lnrhoas * lnrhoas);
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;
2458 static int firstCall = 1;
2463 if (firstCall == 1) {
2465 maxwave =
MAX(aertab->nwave, nwave);
2467 if ((
epsilon[
i] = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
2468 printf(
"Unable to allocate space for epsilon.\n");
2477 if (modnum != lastmod ||
geom->solz[0] != lastsolz ||
2478 geom->senz[0] != lastsenz ||
geom->phi[0] != lastphi ||
2479 iwnir_l != lastiwl) {
2481 int iwnir = iwatab[iwnir_l];
2484 float rhoas1, rhoas2;
2487 if ((lneps = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
2488 printf(
"Unable to allocate space for lneps.\n");
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];
2502 for (iw = 0; iw < nwave; iw++) {
2504 if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
2505 epsilon[
im][iw] = exp(
linterp(aertab->wave, lneps, aertab->nwave, wave[iw]));
2511 lastsolz =
geom->solz[0];
2512 lastsenz =
geom->senz[0];
2513 lastphi =
geom->phi[0];
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) {
2555 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
2556 printf(
"Unable to allocate space for rhoas.\n");
2564 for (jm = 0; jm < nmod; jm++) {
2571 if (rhoas[iwnir_l] > 0.0000001)
2572 eps_ret[jm] = rhoas[iwnir_s] / rhoas[iwnir_l];
2578 eps_mod[jm] = eps[iwnir_s];
2580 eps_ave += eps_ret[jm];
2582 if (isfinite(eps_ave))
2595 for (
im = 0;
im < nmodel;
im++) {
2597 eps_err[
im] = eps_ave - eps_mod[
im];
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]])) {
2605 imod[
im ] = imod[iim];
2617 for (iim = 0; iim < nmod; iim++) {
2619 tot_err +=
fabs(eps_err[
im]);
2623 for (iim = 0; iim < nmod; iim++) {
2625 wt = 1.0 -
fabs(eps_err[
im]) / tot_err;
2626 eps_ave += eps_ret[
im] * wt;
2628 eps_ave /= (nmod - 1);
2632 err_m = 0 - 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];
2642 if (eps_err[
im] > err_m) {
2643 err_m = eps_err[
im];
2650 if (wang_modx && eps_mod[imc50] > eps_mod[imt99]) {
2651 if (*modmax == imt90)
2653 else if (*modmax == imc50 && *modmin == imt99)
2655 else if (*modmin == imm50)
2665 }
else if (*modmax < 0) {
2671 *modrat = (eps_ave - eps_mod[*modmin]) / (eps_mod[*modmax] - eps_mod[*modmin]);
2673 *modmin = mindx[*modmin];
2674 *modmax = mindx[*modmax];
2692 return (
x->angstrom <
y->angstrom ? -1 : 1);
2697 static int firstCall = 1;
2703 int ib =
windex(520, aertab->wave, aertab->nwave);
2706 for (
im = 0;
im < aertab->nmodel;
im++) {
2707 alphaT[
im].angstrom = aertab->model[
im]->angstrom[ib];
2708 alphaT[
im].modnum =
im;
2710 qsort(alphaT, aertab->nmodel, sizeof (alphaTstr),
2711 (
int (*)(
const void *,
const void *))
compalphaT);
2716 for (
im = 0;
im < aertab->nmodel;
im++) {
2717 if (angstrom < alphaT[
im].angstrom)
2720 im1 =
MAX(
MIN(
im - 1, aertab->nmodel - 1), 0);
2721 im2 =
MAX(
MIN(
im, aertab->nmodel - 1), 0);
2723 *modmin = alphaT[im1].modnum;
2724 *modmax = alphaT[im2].modnum;
2729 *modrat = (angstrom - alphaT[im1].angstrom) /
2730 (alphaT[im2].angstrom - alphaT[im1].angstrom);
2748 int32_t iwnir_l,
float rhoa[], geom_str *
geom,
float wv,
float taua[]) {
2753 int iwnir = iwatab[iwnir_l];
2758 iwg = iwnir *
geom->gmult;
2760 maxwave =
MAX(aertab->nwave, nwave);
2762 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
2763 printf(
"Unable to allocate space for rhoas.\n");
2766 if ((
aot = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
2767 printf(
"Unable to allocate space for aot.\n");
2770 if ((lnaot = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
2771 printf(
"Unable to allocate space for lnaot.\n");
2779 aot[iwnir] = rhoas[iwnir_l]*(4.0 *
geom->csolz[iwg] *
geom->csenz[iwg])
2780 / (
phase[iwnir] * aertab->model[modnum]->albedo[iwnir]);
2783 for (iw = 0; iw < aertab->nwave; iw++) {
2785 aot[iw] =
aot[iwnir] * aertab->model[modnum]->extc[iw] / aertab->model[modnum]->extc[iwnir];
2787 lnaot[iw] = log(
aot[iw]);
2792 for (iw = 0; iw < nwave; iw++) {
2794 if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
2795 taua[iw] = exp(
linterp(aertab->wave, lnaot, aertab->nwave, wave[iw]));
2797 taua[iw] =
aot[iwtab];
2800 for (iw = 0; iw < nwave; iw++)
2817 float *theta,
int gmult,
float taua[],
float dtran[]) {
2818 static int firstCall = 1;
2819 static float *intexp;
2836 intexp = (
float *) malloc(nwave *
sizeof (
float));
2837 inttst = (
int *) malloc(nwave *
sizeof (
int));
2839 for (iw = 0; iw < nwave; 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;
2850 printf(
"Interpolating diffuse transmittance for %d from %f by %f\n",
2851 (
int) wave[iw], aertab->dtran_wave[iwtab], intexp[iw]);
2857 for (iw = 0; iw < nwave; iw++) {
2858 if ((iw == 0) || (gmult != 0)) {
2861 for (
i = 0;
i < aertab->dtran_ntheta;
i++) {
2862 if (theta[ig] < aertab->dtran_theta[
i])
2865 if (
i == aertab->dtran_ntheta) {
2870 i1 =
MIN(
MAX(
i - 1, 0), aertab->dtran_ntheta - 2);
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);
2880 a1 = aertab->model[modnum]->dtran_a[iwtab][i1];
2881 b1 = aertab->model[modnum]->dtran_b[iwtab][i1];
2883 a2 = aertab->model[modnum]->dtran_a[iwtab][i2];
2884 b2 = aertab->model[modnum]->dtran_b[iwtab][i2];
2887 a1 = pow(a1, intexp[iw]);
2888 a2 = pow(a2, intexp[iw]);
2891 y1 = a1 * exp(-
b1 * taua[iw]);
2892 y2 = a2 * exp(-b2 * taua[iw]);
2894 dtran[iw] =
MAX(
MIN((1.0 - wt) * y1 + wt*y2, 1.0), 1e-5);
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) {
2915 if ((tsolmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
2916 printf(
"Unable to allocate space for tsolmin.\n");
2919 if ((tsolmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
2920 printf(
"Unable to allocate space for tsolmax.\n");
2923 if ((tsenmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
2924 printf(
"Unable to allocate space for tsenmin.\n");
2927 if ((tsenmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
2928 printf(
"Unable to allocate space for tsenmax.\n");
2931 gmult =
geom->gmult;
2932 if (interpol == 1) gmult = 0;
2949 for (iw = 0; iw < nwave; iw++) {
2950 ig = iw *
geom->gmult;
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;
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));
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]);
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) {
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");
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) {
3004 for (iw = 0; iw < nwave; iw++) {
3005 rhoa[iw] = rho_aer[iw];
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[]) {
3026 float rho_pred_min[nwave], rho_pred_max[nwave];
3027 float rho_aer[nwave], tau_aer[nwave];
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");
3036 rhoa, modmin, modmax, modrat, epsnir,
3037 tau_pred_max, tau_pred_min, rho_pred_max, rho_pred_min, tau_aer, rho_aer) != 0)
3042 rhoa, modmin, modmax, modrat, epsnir,
3043 tau_pred_max, tau_pred_min, rho_pred_max, rho_pred_min, tau_aer, rho_aer) != 0)
3047 for (iw = 0; iw < nwave; iw++) {
3048 rhoa[iw] = rho_aer[iw];
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[]) {
3081 if ((rhoasmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3082 printf(
"Unable to allocate space for rhoasmin.\n");
3085 if ((rhoasmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3086 printf(
"Unable to allocate space for rhoasmax.\n");
3089 if ((rhoamin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3090 printf(
"Unable to allocate space for rhoamin.\n");
3093 if ((rhoamax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3094 printf(
"Unable to allocate space for rhoasmax.\n");
3100 rhoa, iwnir_s, iwnir_l, modmin, modmax, modrat, epsnir);
3104 cc = log(*epsnir) / (wave[iwnir_l] - wave[iwnir_s]);
3126 for (iw = 0; iw < nwave; iw++) {
3128 epsmin1 = epsmin[iw];
3129 epsmax1 = epsmax[iw];
3132 epsmin1 = exp(cc * (wave[iwnir_l] - wave[iw]));
3136 rhoasmin[iw] = rhoasmin[iwnir_l] * epsmin1;
3137 rhoasmax[iw] = rhoasmax[iwnir_l] * epsmax1;
3145 for (iw = 0; iw < nwave; iw++) {
3146 rhoa[iw] = rhoamin[iw]*(1.0 - (*modrat)) + rhoamax[iw]*(*modrat);
3174 return (
x->rhoa <
y->rhoa ? -1 : 1);
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) {
3196 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3197 printf(
"Unable to allocate space for rhoas.\n");
3200 if ((rhoa_tmp = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3201 printf(
"Unable to allocate space for rhoas_tmp.\n");
3207 for (jm = 0; jm < nmodel; jm++) {
3221 rhoas[iwnir_s] = rhoas[iwnir_l] * eps[iwnir_s];
3226 rhoa_tab[jm].modnum =
im;
3227 rhoa_tab[jm].rhoa = rhoa_tmp[iwnir_s];
3228 rhoa_tab[jm].eps = eps[iwnir_s];
3232 qsort(rhoa_tab, nmodel,
sizeof (rhoaTstr), (
int (*)(
const void *,
const void *))
comp_rhoaT);
3235 for (jm = 0; jm < nmodel; jm++) {
3236 if (rhoa_tab[jm].rhoa > rhoa[iwnir_s])
3242 }
else if (jm == nmodel) {
3249 wt = (rhoa[iwnir_s] - rhoa_tab[jm1].rhoa) / (rhoa_tab[jm2].rhoa - rhoa_tab[jm1].rhoa);
3251 *modmin = rhoa_tab[jm1].modnum;
3252 *modmax = rhoa_tab[jm2].modnum;
3254 *epsnir = rhoa_tab[jm1].eps * (1.0 - wt) + rhoa_tab[jm2].eps*wt;
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[]) {
3286 if ((rhoasmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3287 printf(
"Unable to allocate space for rhoasmin.\n");
3290 if ((rhoasmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3291 printf(
"Unable to allocate space for rhoasmax.\n");
3294 if ((rhoamin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3295 printf(
"Unable to allocate space for rhoamin.\n");
3298 if ((rhoamax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3299 printf(
"Unable to allocate space for rhoamax.\n");
3305 if (
model_select_franz(
sensorID, wave, nwave, nmodels, mindx,
geom, wv, rhoa, iwnir_s, iwnir_l,
3306 modmin, modmax, modrat, epsnir) != 0)
3329 for (iw = 0; iw < nwave; iw++) {
3331 epsmin1 = epsmin[iw];
3332 epsmax1 = epsmax[iw];
3334 rhoasmin[iw] = rhoasmin[iwnir_l] * epsmin1;
3335 rhoasmax[iw] = rhoasmax[iwnir_l] * epsmax1;
3343 for (iw = 0; iw < nwave; iw++) {
3344 rhoa[iw] = rhoamin[iw]*(1.0 - (*modrat)) + rhoamax[iw]*(*modrat);
3362 static int order_models(
const void *
p1,
const void *p2) {
3363 aermodstr *
x = *(aermodstr **)
p1;
3364 aermodstr *
y = *(aermodstr **) p2;
3366 if (
x->rh ==
y->rh) {
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;
3407 float *tau_pred_min1;
3408 float *tau_pred_max1;
3409 float *tau_pred_min2;
3410 float *tau_pred_max2;
3419 int irh1, irh2, irh;
3429 float lastrh = -1.0;
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__);
3438 qsort(aertab->model, aertab->nmodel, sizeof (aermodstr*), (
int (*)(
const void *,
const void *)) order_models);
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];
3452 if (nrh == 1 && aertab->model[
im]->sd != lastsd) {
3453 sdtab[nsd] = aertab->model[
im]->sd;
3454 lastsd = sdtab[nsd];
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);
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);
3475 if ((taua1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3476 printf(
"Unable to allocate space for taua1.\n");
3479 if ((taua2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3480 printf(
"Unable to allocate space for taua2.\n");
3483 if ((tsol1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3484 printf(
"Unable to allocate space for tsol1.\n");
3487 if ((tsol2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3488 printf(
"Unable to allocate space for tsol2.\n");
3491 if ((rhoa1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3492 printf(
"Unable to allocate space for rhoa1.\n");
3495 if ((rhoa2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3496 printf(
"Unable to allocate space for rhoa2.\n");
3499 if ((tsen1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3500 printf(
"Unable to allocate space for tsen1.\n");
3503 if ((tsen2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3504 printf(
"Unable to allocate space for tsen2.\n");
3507 if ((tau_pred_min1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3508 printf(
"Unable to allocate space for tau_pred_min1.\n");
3511 if ((tau_pred_min2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3512 printf(
"Unable to allocate space for tau_pred_min2.\n");
3515 if ((tau_pred_max1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3516 printf(
"Unable to allocate space for tau_pred_max1.\n");
3519 if ((tau_pred_max2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
3520 printf(
"Unable to allocate space for tau_pred_max2.\n");
3531 for (iw = 0; iw < nwave; iw++) {
3535 rhoa1[iw] = rhoa[iw];
3536 rhoa2[iw] = rhoa[iw];
3544 printf(
"Warning rh is greater than 95%%. Reset to 94%% rh=%f\n", rh);
3550 if (nrh == 1 || rhtab[0] > rh) {
3554 }
else if (rhtab[nrh - 1] < rh) {
3559 for (irh = 0; irh < nrh; irh++) {
3560 if (rhtab[irh] > rh)
3563 irh1 =
MIN(
MAX(0, irh - 1), nrh - 2);
3565 wt = (rh - rhtab[irh1]) / (rhtab[irh2] - rhtab[irh1]);
3569 if (rh - rhtab[irh1] <= (rhtab[irh2] - rhtab[irh1]) / 2)
3580 for (
im = 0;
im < nsd;
im++) {
3581 mindx3[
im] = irh3 * nsd +
im;
3586 for (
im = 0;
im < nsd;
im++) {
3587 mindx1[
im] = irh1 * nsd +
im;
3588 mindx2[
im] = irh2 * nsd +
im;
3599 if (mindx3[0] == mindx1[0] || mindx3[0] == mindx2[0] || mindx3[0] > 79 || mindx3[0] < 0)
3602 else if (mindx1[0] == mindx2[0])
3605 int32_t vcal_flag = 1;
3606 if ((aertab->nmodel) == vcal_flag)
3612 float rho_aer[nwave];
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) {
3628 *modmin1, *modmax1, *modrat1, rhoa1, taua1, tsol1, tsen1, tau_pred_min1, tau_pred_max1, 1);
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) {
3644 *modmin2, *modmax2, *modrat2, rhoa2, taua2, tsol2, tsen2, tau_pred_min2, tau_pred_max2, 1);
3655 free(tau_pred_min1);
3656 free(tau_pred_max1);
3657 free(tau_pred_min2);
3658 free(tau_pred_max2);
3663 geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1) != 0) {
3672 free(tau_pred_min1);
3673 free(tau_pred_max1);
3674 free(tau_pred_min2);
3675 free(tau_pred_max2);
3681 geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1) != 0) {
3690 free(tau_pred_min1);
3691 free(tau_pred_max1);
3692 free(tau_pred_min2);
3693 free(tau_pred_max2);
3698 geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1) != 0) {
3707 free(tau_pred_min1);
3708 free(tau_pred_max1);
3709 free(tau_pred_min2);
3710 free(tau_pred_max2);
3717 *modmin1, *modmax1, *modrat1, rhoa1, taua1, tsol1, tsen1, tau_pred_min1, tau_pred_max1, 1);
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;
3736 geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2) != 0) {
3745 free(tau_pred_min1);
3746 free(tau_pred_max1);
3747 free(tau_pred_min2);
3748 free(tau_pred_max2);
3753 geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2) != 0) {
3762 free(tau_pred_min1);
3763 free(tau_pred_max1);
3764 free(tau_pred_min2);
3765 free(tau_pred_max2);
3770 geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2) != 0) {
3779 free(tau_pred_min1);
3780 free(tau_pred_max1);
3781 free(tau_pred_min2);
3782 free(tau_pred_max2);
3788 *modmin2, *modmax2, *modrat2, rhoa2, taua2, tsol2, tsen2, tau_pred_min2, tau_pred_max2, 1);
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;
3796 *eps = eps1 * (1 - wt) + eps2*wt;
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];
3829 free(tau_pred_min1);
3830 free(tau_pred_max1);
3831 free(tau_pred_min2);
3832 free(tau_pred_max2);
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) {
3851 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3852 printf(
"Unable to allocate space for rhoas.\n");
3856 if (rhoa[iwnir_l] < 0.0) {
3869 printf(
"Error getting rhoas\n");
3875 for (iw = 0; iw < nwave; iw++) {
3876 rhoas[iw] = rhoas[iwnir_l] * eps[iw];
3882 if (iwnir_s == iwnir_l)
3883 *epsnir = eps[iwnir_l - 1];
3885 *epsnir = eps[iwnir_s];
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) {
3909 if ((rhoa1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3910 printf(
"Unable to allocate space for rhoa1.\n");
3913 if ((rhoa2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3914 printf(
"Unable to allocate space for rhoa2.\n");
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);
3926 if (rhoa[iwnir_l] >
input->rhoamin) {
3928 rhoa2[iwnir_l] = rhoa1[iwnir_l] = rhoa[iwnir_l];
3938 for (iw = 0; iw < nwave; iw++) {
3940 rhoa[iw] =
MAX((1.0 - modrat) * rhoa1[iw] + modrat * rhoa2[iw], 0.0);
3942 *eps = (1.0 - modrat) * eps1 + modrat*eps2;
3945 }
else if (rhoa[iwnir_l] > -(
input->rhoamin)) {
3949 for (iw = 0; iw < nwave; iw++) {
3950 rhoa[iw] =
MAX(rhoa[iwnir_l], 1e-6);
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[],
3979 static int firstCall = 1;
3980 static int angst_band1 = -1;
3981 static int angst_band2 = -1;
3996 int ig, gmult, iw, iwtab;
3999 maxwave =
MAX(aertab->nwave, nwave);
4001 if ((rhoa1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4002 printf(
"Unable to allocate space for rhoa1.\n");
4005 if ((rhoa2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4006 printf(
"Unable to allocate space for rhoa2.\n");
4009 if ((rhoas1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4010 printf(
"Unable to allocate space for rhoas1.\n");
4013 if ((rhoas2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4014 printf(
"Unable to allocate space for rhoas2.\n");
4017 if ((
f1 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
4018 printf(
"Unable to allocate space for f1.\n");
4021 if ((
f2 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
4022 printf(
"Unable to allocate space for f2.\n");
4025 if ((lnf1 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
4026 printf(
"Unable to allocate space for lnf1.\n");
4029 if ((lnf2 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
4030 printf(
"Unable to allocate space for lnf2.\n");
4035 angst_band1 =
windex(520, wave, nwave);
4036 angst_band2 =
windex(865, wave, nwave);
4041 for (iw = 0; iw < nwave; iw++) {
4042 if (
aot[iw] < 0.0) {
4056 if (
aot[iwnir_l] > 0.0)
4057 angstrom = -log(
aot[angst_band1] /
aot[angst_band2]) /
4058 log(wave[angst_band1] / wave[angst_band2]);
4069 gmult = (interpol == 1) ? 0 :
geom->gmult;
4072 for (iw = 0; iw < aertab->nwave; 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];
4079 lnf1[iw] = log(
f1[iw]);
4080 lnf2[iw] = log(
f2[iw]);
4086 for (iw = 0; iw < nwave; 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]));
4092 rhoas1[iw] =
aot[iw] *
f1[iwtab];
4093 rhoas2[iw] =
aot[iw] *
f2[iwtab];
4097 for (iw = 0; iw < nwave; iw++) {
4099 rhoas1[iw] =
aot[iw] *
f1[iwtab];
4100 rhoas2[iw] =
aot[iw] *
f2[iwtab];
4103 eps1 = rhoas1[iwnir_s] / rhoas1[iwnir_l];
4104 eps2 = rhoas2[iwnir_s] / rhoas2[iwnir_l];
4110 for (iw = 0; iw < nwave; iw++) {
4111 rhoa[iw] = (1.0 - *modrat) * rhoa1[iw] + *modrat * rhoa2[iw];
4113 *epsnir = (1.0 - *modrat) * eps1 + *modrat * eps2;
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;
4159 float *taua_pred_min;
4160 float *taua_pred_max;
4162 l1str *
l1rec = l2rec->l1rec;
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;
4171 float wv =
l1rec->wv [ip];
4172 float rh =
l1rec->rh [ip];
4175 if (firstCall == 1) {
4177 Maxband = nwave + 1;
4180 if ((rhoa = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4181 printf(
"Unable to allocate space for rhoa.\n");
4184 if ((radref = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4185 printf(
"Unable to allocate space for raderef.\n");
4188 if ((F0 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4189 printf(
"Unable to allocate space for F0.\n");
4192 if ((taur = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4193 printf(
"Unable to allocate space for taur.\n");
4196 if ((La1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4197 printf(
"Unable to allocate space for La1.\n");
4200 if ((La2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4201 printf(
"Unable to allocate space for La2.\n");
4204 if ((t_sol = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4205 printf(
"Unable to allocate space for t_sol.\n");
4208 if ((t_sen = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4209 printf(
"Unable to allocate space for t_sen.\n");
4212 if ((taua = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4213 printf(
"Unable to allocate space for taua.\n");
4216 if ((taua_pred_min = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4217 printf(
"Unable to allocate space for taua_pred_min.\n");
4220 if ((taua_pred_max = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4221 printf(
"Unable to allocate space for taua_pred_max.\n");
4235 geom.airmass = (
float *) malloc(
sizeof (
float));
4236 *
geom.airmass = 1. /
geom.csolz[0] + 1. /
geom.csenz[0];
4238 geom.airmass_sph = (
float *) malloc(
sizeof (
float));
4239 geom.airmass_plp = (
float *) malloc(
sizeof (
float));
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];
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];
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++) {
4271 evalmask =
input->evalmask;
4272 aer_opt = aer_opt_in;
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;
4284 mu0 = cos(
solz / radeg);
4285 mu = cos(
senz / radeg);
4286 airmass = 1.0 / mu0 + 1.0 / mu;
4300 for (iw = 0; iw < nwave; iw++) {
4305 radref[iw] =
pi / F0[iw] /
geom.csolz[iw *
geom.gmult];
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");
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");
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]);
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]);
4353 printf(
"\nUsing white-aerosol approximation\n");
4354 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
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]);
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]);
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);
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]);
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]);
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]);
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]);
4402 printf(
"\nUsing fixed, input aerosol optical thicknesses for aerosol selection.\n");
4405 printf(
"\nUsing multiple scattering aerosols with fixed model pair\n");
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]);
4413 printf(
"\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
4414 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
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]);
4422 printf(
"\nErroneous aerosol option, %d\n", aer_opt);
4432 for (iw = iwnir_s; iw <= iwnir_l; iw +=
MAX(iwnir_l - iwnir_s, 1))
4433 rhoa[iw] = La1[iw] * radref[iw];
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);
4452 for (iw = iwnir_s; iw <= iwnir_l; iw++) {
4453 rhoa[iw] = La1[iw] * radref[iw];
4457 if (rhoa[iwnir_s] >
input->rhoamin && rhoa[iwnir_l] >
input->rhoamin) {
4460 if (La1[iwnir_s] / La1[iwnir_l] > 0.1) {
4463 aertab->nmodel, mindx,
4464 &
geom, wv, rhoa, modmin, modmax, modrat, eps, taua, taua);
4467 for (iw = 0; iw < nwave; iw++)
4468 La2[iw] = rhoa[iw] / radref[iw];
4471 }
else if (rhoa[iwnir_s] > -(
input->rhoamin) && rhoa[iwnir_l] > -(
input->rhoamin)) {
4475 *modmin = aertab->nmodel;
4476 *modmax = aertab->nmodel;
4478 temp =
MAX(rhoa[iwnir_l], 1e-6);
4479 for (iw = 0; iw < nwave; iw++) {
4481 La2 [iw] = rhoa[iw] / radref[iw];
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");
4512 for (iw = iwnir_s; iw <= iwnir_l; iw++) {
4513 rhoa[iw] = La1[iw] * radref[iw];
4517 if (rhoa[iwnir_s] >
input->rhoamin && rhoa[iwnir_l] >
input->rhoamin) {
4520 if (rhoa[iwnir_s] / rhoa[iwnir_l] > 0.1 && rhoa[iwnir_s] / rhoa[iwnir_l] < 10.0) {
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);
4527 for (iw = 0; iw < nwave; iw++)
4528 La2[iw] = rhoa[iw] / radref[iw];
4531 }
else if (rhoa[iwnir_s] > -(
input->rhoamin) && rhoa[iwnir_l] > -(
input->rhoamin)) {
4535 *modmin = aertab->nmodel;
4536 *modmax = aertab->nmodel;
4538 *modmin2 = aertab->nmodel;
4539 *modmax2 = aertab->nmodel;
4541 temp =
MAX(rhoa[iwnir_l], 1e-6);
4542 for (iw = 0; iw < nwave; iw++) {
4544 La2 [iw] = rhoa[iw] / radref[iw];
4548 *modmin, *modmax, *modrat, rhoa, taua, t_sol, t_sen, taua_pred_min, taua_pred_max, 0);
4565 if (La1[iwnir_l] > 0.0) {
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];
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];
4592 *modmin =
input->aermodmin - 1;
4593 *modmax =
input->aermodmax - 1;
4594 *modrat =
input->aermodrat;
4598 *modmin, *modmax, *modrat, rhoa, eps);
4600 for (iw = 0; iw < nwave; iw++)
4601 La2[iw] = rhoa[iw] / radref[iw];
4608 if (
input->aer_angstrom > -2) {
4609 angstrom =
input->aer_angstrom;
4613 unix2yds(l2rec->l1rec->scantime, &year, &
day, &sec);
4617 if (angstrom > -2) {
4622 *modmin, *modmax, *modrat, rhoa, eps);
4627 for (iw = 0; iw < nwave; iw++)
4628 La2[iw] = rhoa[iw] / radref[iw];
4637 if (aerec !=
NULL && aerec->mode ==
ON)
4643 modmin, modmax, modrat, rhoa, eps);
4645 for (iw = 0; iw < nwave; iw++)
4646 La2[iw] = rhoa[iw] / radref[iw];
4656 *modmin = aer_opt - 1;
4657 *modmax = aer_opt - 1;
4660 if (*modmin < 0 || *modmin >
input->naermodels - 1) {
4661 printf(
"Invalid aerosol option: %d\n", *modmin);
4666 rhoa[iwnir_l] = La1[iwnir_l] * radref[iwnir_l];
4671 &
geom, wv, rhoa, eps);
4675 for (iw = 0; iw < nwave; iw++)
4676 La2[iw] = rhoa[iw] / radref[iw];
4684 *modmin, *modmax, *modrat, rhoa, taua, t_sol, t_sen, taua_pred_min, taua_pred_max, 0);
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];
4696 *modmin = *modmin + 1;
4697 *modmax = *modmax + 1;
4698 *modmin2 = *modmin2 + 1;
4699 *modmax2 = *modmax2 + 1;
4710 free(taua_pred_min);
4711 free(taua_pred_max);
4714 free(
geom.airmass_plp);
4715 free(
geom.airmass_sph);
4738 static int firstCall = 1;
4748 l1str *
l1rec = l2rec->l1rec;
4753 wave2 =
l1file->fwave[ib2];
4761 wave1 =
l1file->fwave[ib1];
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)
4786 int32_t ip, ipb1, ipb2;
4788 l1str *
l1rec = l2rec->l1rec;
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]