57 #define MAXMODEL MAXAERMOD
64 static int32_t first_solz;
65 static int32_t first_senz;
66 static int32_t first_phi;
67 static int32_t first_scatt;
68 static int32_t first_dtnwave;
69 static int32_t first_dtntheta;
70 static int32_t first_npc;
71 static int32_t first_ntau_870;
74 static double radeg =
RADEG;
75 static float p0 =
STDPR;
77 static int have_ms = 0;
78 static int have_rh = 0;
79 static int have_sd = 0;
80 static int use_rh = 0;
81 static int use_netcdf = 0;
82 static int use_pca_lut = 0;
85 static int32_t Maxband;
138 static int cmp(
const void *
a,
const void *
b) {
139 struct str *a1 = (
struct str *)
a;
140 struct str *a2 = (
struct str*)
b;
141 if((*a1).value<(*a2).value)
return -1;
142 else if((*a1).value>(*a2).value)
return 1;
146 aermodstr*
alloc_aermodstr(
int nbands,
int nscatt,
int nphi,
int nsolz,
int nsenz,
int ntheta,
int npc,
int ntau_870) {
150 model = (aermodstr *) calloc(1,
sizeof(aermodstr));
151 model->angstrom = (
float*) malloc(
nbands *
sizeof (
float));
152 model->extc = (
float*) malloc(
nbands *
sizeof (
float));
161 model->pc_mean_rhoa = (
float*) malloc(
nbands *
sizeof(
float));
170 model->pc_mean_td = (
float*) malloc(
nbands *
sizeof(
float));
173 model->tau_870 = (
float*) malloc(ntau_870 *
sizeof(
float));
176 model->pc = (
float*) malloc(npc *
sizeof(
float));
179 model->albedo = (
float*) malloc(
nbands *
sizeof (
float));
183 model->acost = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
184 model->bcost = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
185 model->ccost = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
189 model->ams_all = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
190 model->bms_all = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
191 model->cms_all = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
192 model->dms_all = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
193 model->ems_all = (
float*) malloc(
nbands * nsolz * nphi * nsenz *
sizeof (
float));
240 static aermodtabstr *aertab =
NULL;
259 static int loaded = 0;
260 static int interpol = 0;
261 static int32_t *iwatab;
262 static int32_t *iwdtab;
264 static int32_t iwnir_s = -1;
265 static int32_t iwnir_l = -1;
269 static float airmass;
271 static int32_t evalmask = 0;
272 static int32_t aer_opt = 0;
273 static float airmass_plp;
274 static float airmass_sph;
277 if (*(
double*)
a > *(
double*)
b)
return 1;
278 else if (*(
double*)
a < *(
double*)
b)
return -1;
290 float a1, a2, a3, a4, a5, a6, d1;
301 d1 =
y[0]*(1.0 / a1 + 1.0 / a2 + 1.0 / a3)
302 - a2 * a3 *
y[1] / (a1 * a4 * a5)
303 + a1 * a3 *
y[2] / (a2 * a4 * a6)
304 - a1 * a2 *
y[3] / (a3 * a5 * a6);
308 a1 =
x[n - 1] -
x[n - 4];
309 a2 =
x[n - 1] -
x[n - 3];
310 a3 =
x[n - 1] -
x[n - 2];
311 a4 =
x[n - 2] -
x[n - 4];
312 a5 =
x[n - 2] -
x[n - 3];
313 a6 =
x[n - 3] -
x[n - 4];
315 d1 = -a2 * a3 *
y[n - 4] / (a6 * a4 * a1)
316 + a1 * a3 *
y[n - 3] / (a6 * a5 * a2)
317 - a1 * a2 *
y[n - 2] / (a4 * a5 * a3)
318 +
y[n - 1]*(1.0 / a1 + 1.0 / a2 + 1.0 / a3);
325 static void read_dimension_size(
const char* file_name,
int nc_id,
const char* dim_name, int32_t *length) {
330 status = nc_inq_dimid (nc_id, dim_name, &dim_id);
332 printf(
"-E- %s: Error looking for dimension %s from %s.\n", __FILE__, dim_name, file_name);
335 status = nc_inq_dimlen(nc_id, dim_id, &dim_length);
337 printf(
"-E- %s: Error reading dimension %s from %s.\n", __FILE__, dim_name, file_name);
340 *length = (int32_t) dim_length;
344 static void read_lut_variable(
const char*
file,
int nc_id,
const char* var_name,
float*
data) {
348 status = nc_inq_varid(nc_id, var_name, &var_id);
350 printf(
"-E- %s: Error looking for variable %s from %s.\n", __FILE__, var_name,
file);
353 status = nc_get_var_float(nc_id, var_id,
data);
355 printf(
"-E- %s: Error reading variable %s from %s.\n", __FILE__, var_name,
file);
384 int32_t aer_nwave, nsolz, nsenz, nphi, nscatt, dtran_nwave, dtran_ntheta, npc, ntau_870;
391 char file [FILENAME_MAX] =
"";
392 char path [FILENAME_MAX] =
"";
394 int iw,
im,
is, iwbase,
i;
395 static int firstCall = 1;
397 if (firstCall == 1) {
398 if ((iwatab = (int32_t *) calloc(nwave,
sizeof (int32_t))) ==
NULL) {
399 printf(
"Unable to allocate space for iwatab.\n");
402 if ((iwdtab = (int32_t *) calloc(nwave,
sizeof (int32_t))) ==
NULL) {
403 printf(
"Unable to allocate space for iwdtab.\n");
409 printf(
"Loading aerosol models from %s\n", aermodfile);
411 for (
im = 0;
im < nmodels + 1;
im++) {
417 strcat(
file, aermodfile);
424 if(access(
file, R_OK) == -1) {
426 strcat(
file, aermodfile);
429 strcat(
file,
".hdf");
437 strcat(
file, aermodfile);
438 strcat(
file,
"_default.nc");
442 if(access(
file, R_OK) == -1) {
444 strcat(
file, aermodfile);
445 strcat(
file,
"_default.hdf");
452 printf(
"-E- %s: Error opening file %s.\n", __FILE__,
file);
467 status = nc_inq_attlen(nc_id, NC_GLOBAL,
"title", &attlen);
469 printf(
"-E- %s: could not find title in %s.\n", __FILE__,
file);
472 char*
title = malloc(
sizeof(
char) * attlen);
473 nc_get_att_text(nc_id, NC_GLOBAL,
"title",
title);
474 if(strcasestr(
title,
"pca aerosol model data")) {
480 read_dimension_size(
file, nc_id,
"wavelength", &aer_nwave);
481 read_dimension_size(
file, nc_id,
"solar_zenith", &nsolz);
482 read_dimension_size(
file, nc_id,
"sensor_zenith", &nsenz);
483 read_dimension_size(
file, nc_id,
"relative_azimuth", &nphi);
484 read_dimension_size(
file, nc_id,
"principal_component", &npc);
485 read_dimension_size(
file, nc_id,
"aerosol_optical_thickness", &ntau_870);
487 read_dimension_size(
file, nc_id,
"nwave", &aer_nwave);
488 read_dimension_size(
file, nc_id,
"nsolz", &nsolz);
489 read_dimension_size(
file, nc_id,
"nsenz", &nsenz);
490 read_dimension_size(
file, nc_id,
"nphi", &nphi);
491 read_dimension_size(
file, nc_id,
"dtran_nwave", &dtran_nwave);
492 read_dimension_size(
file, nc_id,
"dtran_ntheta", &dtran_ntheta);
495 read_dimension_size(
file, nc_id,
"nwave", &aer_nwave);
496 read_dimension_size(
file, nc_id,
"nsolz", &nsolz);
497 read_dimension_size(
file, nc_id,
"nsenz", &nsenz);
498 read_dimension_size(
file, nc_id,
"nphi", &nphi);
499 read_dimension_size(
file, nc_id,
"nscatt", &nscatt);
500 read_dimension_size(
file, nc_id,
"dtran_nwave", &dtran_nwave);
501 read_dimension_size(
file, nc_id,
"dtran_ntheta", &dtran_ntheta);
510 printf(
"Number of Wavelengths %d\n", aer_nwave);
511 printf(
"Number of Solar Zenith Angles %d\n", nsolz);
512 printf(
"Number of View Zenith Angles %d\n", nsenz);
513 printf(
"Number of Relative Azimuth Angles %d\n", nphi);
515 printf(
"Number of Principal Components %d\n", npc);
516 printf(
"Number of tau 870 Angles %d\n", ntau_870);
518 printf(
"Number of Scattering Angles %d\n", nscatt);
519 printf(
"Number of Diffuse Transmittance Wavelengths %d\n", dtran_nwave);
520 printf(
"Number of Diffuse Transmittance Zenith Angles %d\n", dtran_ntheta);
526 first_scatt = nscatt;
527 first_dtnwave = dtran_nwave;
528 first_dtntheta = dtran_ntheta;
530 first_ntau_870 = ntau_870;
533 if ((aertab = (aermodtabstr *) calloc(1,
sizeof (aermodtabstr))) ==
NULL) {
534 printf(
"Unable to allocate space for aerosol table.\n");
538 aertab->nmodel = nmodels;
539 aertab->nwave = aer_nwave;
540 aertab->nsolz = nsolz;
541 aertab->nsenz = nsenz;
543 aertab->nscatt = nscatt;
544 aertab->dtran_nwave = dtran_nwave;
545 aertab->dtran_ntheta = dtran_ntheta;
549 aertab->ntau_870 = ntau_870;
552 aertab->wave = (
float *) malloc(aer_nwave *
sizeof (
float));
553 aertab->solz = (
float *) malloc(nsolz *
sizeof (
float));
554 aertab->senz = (
float *) malloc(nsenz *
sizeof (
float));
555 aertab->phi = (
float *) malloc(nphi *
sizeof (
float));
557 aertab->scatt =
NULL;
559 aertab->scatt = (
float *) malloc(nscatt *
sizeof (
float));
562 aertab->dtran_wave = (
float *) malloc(dtran_nwave *
sizeof (
float));
563 aertab->dtran_theta = (
float *) malloc(dtran_ntheta *
sizeof (
float));
564 aertab->dtran_airmass = (
float *) malloc(dtran_ntheta *
sizeof (
float));
568 if ((aertab->model = (aermodstr **) calloc(1, (nmodels + 1) *
sizeof (aermodstr*))) ==
NULL) {
569 printf(
"Unable to allocate space for %d aerosol models.\n", nmodels + 1);
572 for (
i = 0;
i < nmodels + 1;
i++) {
573 if ((aertab->model[
i] =
alloc_aermodstr(aer_nwave, nscatt, nphi, nsolz, nsenz, dtran_ntheta, npc, ntau_870)) ==
NULL) {
574 printf(
"Unable to allocate space for aerosol model %d.\n",
im);
582 read_lut_variable(
file, nc_id,
"scatt", aertab->scatt);
584 read_lut_variable(
file, nc_id,
"wave", aertab->wave);
585 read_lut_variable(
file, nc_id,
"solz", aertab->solz);
586 read_lut_variable(
file, nc_id,
"senz", aertab->senz);
587 read_lut_variable(
file, nc_id,
"phi", aertab->phi);
591 read_lut_variable(
file, nc_id,
"wavelength", aertab->wave);
592 read_lut_variable(
file, nc_id,
"solar_zenith", aertab->solz);
593 read_lut_variable(
file, nc_id,
"sensor_zenith", aertab->senz);
594 read_lut_variable(
file, nc_id,
"relative_azimuth", aertab->phi);
597 read_lut_variable(
file, nc_id,
"wave", aertab->wave);
598 read_lut_variable(
file, nc_id,
"solz", aertab->solz);
599 read_lut_variable(
file, nc_id,
"senz",aertab->senz);
600 read_lut_variable(
file, nc_id,
"phi", aertab->phi);
604 read_lut_variable(
file, nc_id,
"dtran_wave", aertab->dtran_wave);
605 read_lut_variable(
file, nc_id,
"dtran_theta", aertab->dtran_theta);
610 if ((aertab->nsolz != first_solz) || (aertab->nsenz != first_senz) ||
611 (aertab->nphi != first_phi) || (aertab->nscatt != first_scatt) ||
612 (aertab->dtran_nwave != first_dtnwave) || (aertab->dtran_ntheta != first_dtntheta) ||
613 (use_pca_lut && aertab->ntau_870 != first_ntau_870) || (use_pca_lut && aertab->npc != first_npc)) {
614 printf(
"-E- %s, %d: Error, Aerosol table %s\n",
615 __FILE__, __LINE__,
file);
616 printf(
" has different dimensions from previous tables\n");
624 strncpy(aertab->model[
im]->name,
"default", 32);
626 status = nc_get_att_float(nc_id, NC_GLOBAL,
"RelativeHumidity", &rh);
628 status = nc_get_att_float(nc_id, NC_GLOBAL,
"Relative Humidity", &rh);
631 status = nc_get_att_float(nc_id, NC_GLOBAL,
"relative_humidity", &rh);
636 aertab->model[
im]->rh = rh;
639 aertab->model[
im]->rh = -1.0;
645 status = nc_get_att_float(nc_id, NC_GLOBAL,
"AerosolFMF", &fmf);
647 status = nc_get_att_float(nc_id, NC_GLOBAL,
"fine_mode_fraction", &fmf);
649 printf(
"-E- %s, %d: Error, Aerosol table %s does not have AerosolFMF\n",
650 __FILE__, __LINE__,
file);
655 aertab->model[
im]->sd = sd;
657 status = nc_get_att_short(nc_id, NC_GLOBAL,
"Size Distribution", &sd);
659 aertab->model[
im]->sd = sd;
662 aertab->model[
im]->sd = -1;
668 read_lut_variable(
file, nc_id,
"extc", aertab->model[
im]->extc);
671 read_lut_variable(
file, nc_id,
"extinction_coefficient", aertab->model[
im]->extc);
673 read_lut_variable(
file, nc_id,
"extc", aertab->model[
im]->extc);
677 read_lut_variable(
file, nc_id,
"aerosol_reflectance_score", aertab->model[
im]->pc_rhoa[0][0][0][0]);
678 read_lut_variable(
file, nc_id,
"aerosol_reflectance_component", aertab->model[
im]->pc_components_rhoa[0]);
679 read_lut_variable(
file, nc_id,
"aerosol_reflectance_mean", aertab->model[
im]->pc_mean_rhoa);
680 read_lut_variable(
file, nc_id,
"diffuse_transmittance_score", aertab->model[
im]->pc_td[0][0]);
681 read_lut_variable(
file, nc_id,
"diffuse_transmittance_component", aertab->model[
im]->pc_components_td[0]);
682 read_lut_variable(
file, nc_id,
"diffuse_transmittance_mean", aertab->model[
im]->pc_mean_td);
683 read_lut_variable(
file, nc_id,
"angstrom", aertab->model[
im]->angstrom);
684 read_lut_variable(
file, nc_id,
"aerosol_optical_thickness", aertab->model[
im]->tau_870);
687 read_lut_variable(
file, nc_id,
"albedo", aertab->model[
im]->albedo);
688 read_lut_variable(
file, nc_id,
"phase", aertab->model[
im]->phase[0]);
689 read_lut_variable(
file, nc_id,
"acost", aertab->model[
im]->acost);
690 read_lut_variable(
file, nc_id,
"bcost", aertab->model[
im]->bcost);
691 read_lut_variable(
file, nc_id,
"ccost", aertab->model[
im]->ccost);
695 if(nc_inq_varid(nc_id,
"ams_all", &
status) == NC_NOERR) {
697 read_lut_variable(
file, nc_id,
"ams_all", aertab->model[
im]->ams_all);
698 read_lut_variable(
file, nc_id,
"bms_all", aertab->model[
im]->bms_all);
699 read_lut_variable(
file, nc_id,
"cms_all", aertab->model[
im]->cms_all);
700 read_lut_variable(
file, nc_id,
"dms_all", aertab->model[
im]->dms_all);
701 read_lut_variable(
file, nc_id,
"ems_all", aertab->model[
im]->ems_all);
704 read_lut_variable(
file, nc_id,
"dtran_a", aertab->model[
im]->dtran_a[0]);
705 read_lut_variable(
file, nc_id,
"dtran_b", aertab->model[
im]->dtran_b[0]);
712 iwbase =
windex(865, aertab->wave, aertab->nwave);
713 for (iw = 0; iw < aertab->nwave; iw++) {
715 aertab->model[
im]->angstrom[iw] = -log(aertab->model[
im]->extc[iw] / aertab->model[
im]->extc[iwbase]) /
716 log(aertab->wave[iw] / aertab->wave[iwbase]);
718 aertab->model[
im]->angstrom[iwbase] = aertab->model[
im]->angstrom[iwbase - 1];
722 for (iw = 0; iw < aertab->nwave; iw++) {
723 for (
is = 0;
is < aertab->nscatt;
is++) {
724 aertab->model[
im]->lnphase[iw][
is] = log(aertab->model[
im]->phase[iw][
is]);
726 d1phase1 =
first_deriv(aertab->scatt, &aertab->model[
im]->lnphase[iw][0], 0);
727 d1phaseN =
first_deriv(aertab->scatt, &aertab->model[
im]->lnphase[iw][0], aertab->nscatt);
729 &aertab->model[
im]->lnphase[iw][0],
733 &aertab->model[
im]->d2phase[iw][0]);
740 aertab->nmodel = nmodels;
744 if (aertab->nwave != nwave) {
745 printf(
"Number of aerosol LUT wavelengths (%d) is not equal to number of sensor wavelengths (%d).\n",
746 aertab->nwave,nwave);
752 for (
is = 0;
is < aertab->dtran_ntheta;
is++) {
753 aertab->dtran_airmass[
is] = 1.0 / cos(aertab->dtran_theta[
is] / radeg);
756 if (aertab->dtran_nwave != nwave) {
757 printf(
"Number of aerosol diffue trans LUT wavelengths (%d) is not equal to number of sensor wavelengths (%d).\n",
758 aertab->dtran_nwave,nwave);
763 printf(
"Wavelengths - Sensor\tAerosol model\tDiffuse transmittance\n");
764 for (iw = 0; iw < nwave; iw++) {
765 iwatab[iw] =
windex(wave[iw], aertab->wave, aertab->nwave);
766 iwdtab[iw] =
windex(wave[iw], aertab->dtran_wave, aertab->dtran_nwave);
767 printf(
"\t %6.2f\t %6.2f\t %6.2f\n", wave[iw],aertab->wave[iwatab[iw]], aertab->dtran_wave[iwdtab[iw]]);
779 #define INDEX(iw,isol,iphi,isen) (iw*aertab->nsolz*aertab->nphi*aertab->nsenz + isol*aertab->nphi*aertab->nsenz + iphi*aertab->nsenz + isen)
794 static float lastsolz = -999.;
795 static float lastsenz = -999.;
796 static float lastphi = -999.;
804 static float *
p, *
q, *
r;
805 static float as000, as100, as010, as110, as001, as011, as101, as111;
806 static float ai000, ai100, ai010, ai110, ai001, ai011, ai101, ai111;
807 static float ac000, ac100, ac010, ac110, ac001, ac011, ac101, ac111;
809 static int *isolz1, *isolz2;
810 static int *isenz1, *isenz2;
811 static int *iphi1, *iphi2;
813 static float *p_ar, *q_ar, *r_ar;
814 static float p_cnst, q_cnst, r_cnst;
815 static int *isolz1_ar, *isolz2_ar, isolz1_cnst, isolz2_cnst;
816 static int *isenz1_ar, *isenz2_ar, isenz1_cnst, isenz2_cnst;
817 static int *iphi1_ar, *iphi2_ar, iphi1_cnst, iphi2_cnst;
822 static int firstCall = 1;
825 if (firstCall == 1) {
828 if ((a_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
829 printf(
"Unable to allocate space for a_coef.\n");
832 if ((b_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
833 printf(
"Unable to allocate space for rhoa.\n");
836 if ((c_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
837 printf(
"Unable to allocate space for rhoa.\n");
842 if ((
geom->gmult == 0) || (interpol == 1)) {
847 isolz1 = &isolz1_cnst;
848 isolz2 = &isolz2_cnst;
849 isenz1 = &isenz1_cnst;
850 isenz2 = &isenz2_cnst;
856 if (((p_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
858 ((q_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
860 ((r_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
862 printf(
"Unable to allocate space for p, q, r weights.\n");
865 if (((isolz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
867 ((isenz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
869 ((iphi1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
871 printf(
"Unable to allocate space for interp indicies 1.\n");
874 if (((isolz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
876 ((isenz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
878 ((iphi2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
880 printf(
"Unable to allocate space for interp indicies 2.\n");
895 if ((
geom->solz[0] != lastsolz) || (
geom->senz[0] != lastsenz) ||
896 (
geom->phi[0] != lastphi)) {
897 for (
im = 0;
im < aertab->nmodel;
im++)
900 for (iw = 0; iw < aertab->nwave; iw++) {
903 for (
i = 0;
i < aertab->nsolz;
i++) {
904 if (
geom->solz[ig] < aertab->solz[
i])
907 isolz1[iw] =
MAX(
i - 1, 0);
908 isolz2[iw] =
MIN(
i, aertab->nsolz - 1);
909 if (isolz2[iw] != isolz1[iw])
910 r[iw] = (
geom->solz[ig] - aertab->solz[isolz1[iw]]) /
911 (aertab->solz[isolz2[iw]] - aertab->solz[isolz1[iw]]);
916 for (
i = 0;
i < aertab->nsenz;
i++) {
917 if (
geom->senz[ig] < aertab->senz[
i])
920 isenz1[iw] =
MAX(
i - 1, 0);
921 isenz2[iw] =
MIN(
i, aertab->nsenz - 1);
922 if (isenz2[iw] != isenz1[iw])
923 p[iw] = (
geom->senz[ig] - aertab->senz[isenz1[iw]]) /
924 (aertab->senz[isenz2[iw]] - aertab->senz[isenz1[iw]]);
930 for (
i = 0;
i < aertab->nphi;
i++) {
931 if (aphi < aertab->phi[
i])
934 iphi1[iw] =
MAX(
i - 1, 0);
935 iphi2[iw] =
MIN(
i, aertab->nphi - 1);
936 if (iphi2[iw] != iphi1[iw])
937 q[iw] = (aphi - aertab->phi[iphi1[iw]]) /
938 (aertab->phi[iphi2[iw]] - aertab->phi[iphi1[iw]]);
941 if (gmult == 0)
break;
945 lastsolz =
geom->solz[0];
946 lastsenz =
geom->senz[0];
947 lastphi =
geom->phi[0];
954 for (iw = 0; iw < aertab->nwave; iw++) {
959 if (isolz2[ig] == 0) {
960 as000 = aertab->model[
im]->acost[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
961 as100 = aertab->model[
im]->acost[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
962 as001 = aertab->model[
im]->acost[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
963 as101 = aertab->model[
im]->acost[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
965 ai000 = aertab->model[
im]->bcost[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
966 ai100 = aertab->model[
im]->bcost[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
967 ai001 = aertab->model[
im]->bcost[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
968 ai101 = aertab->model[
im]->bcost[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
970 ac000 = aertab->model[
im]->ccost[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
971 ac100 = aertab->model[
im]->ccost[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
972 ac001 = aertab->model[
im]->ccost[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
973 ac101 = aertab->model[
im]->ccost[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
975 a_coef[
im][iw] = (1. - px)*(1. - rx) * as000 + px * rx * as101
976 + (1. - px) * rx * as001 + px * (1. - rx) * as100;
978 b_coef[
im][iw] = (1. - px)*(1. - rx) * ai000 + px * rx * ai101
979 + (1. - px) * rx * ai001 + px * (1. - qx)*(1. - rx) * ai100;
981 c_coef[
im][iw] = (1. - px)*(1. - rx) * ac000 + px * rx * ac101
982 + (1. - px) * rx * ac001 + px * (1. - qx)*(1. - rx) * ac100;
984 as000 = aertab->model[
im]->acost[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
985 as100 = aertab->model[
im]->acost[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
986 as010 = aertab->model[
im]->acost[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
987 as110 = aertab->model[
im]->acost[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
988 as001 = aertab->model[
im]->acost[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
989 as011 = aertab->model[
im]->acost[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
990 as101 = aertab->model[
im]->acost[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
991 as111 = aertab->model[
im]->acost[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
993 ai000 = aertab->model[
im]->bcost[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
994 ai100 = aertab->model[
im]->bcost[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
995 ai010 = aertab->model[
im]->bcost[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
996 ai110 = aertab->model[
im]->bcost[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
997 ai001 = aertab->model[
im]->bcost[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
998 ai011 = aertab->model[
im]->bcost[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
999 ai101 = aertab->model[
im]->bcost[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1000 ai111 = aertab->model[
im]->bcost[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1002 ac000 = aertab->model[
im]->ccost[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1003 ac100 = aertab->model[
im]->ccost[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1004 ac010 = aertab->model[
im]->ccost[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1005 ac110 = aertab->model[
im]->ccost[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1006 ac001 = aertab->model[
im]->ccost[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1007 ac011 = aertab->model[
im]->ccost[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1008 ac101 = aertab->model[
im]->ccost[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1009 ac111 = aertab->model[
im]->ccost[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1011 a_coef[
im][iw] = (1. - px)*(1. - qx)*(1. - rx) * as000 + px * qx * rx * as111
1012 + px * (1. - qx) * rx * as101 + (1. - px) * qx * (1. - rx) * as010
1013 + px * qx * (1. - rx) * as110 + (1. - px)*(1. - qx) * rx * as001
1014 + (1. - px) * qx * rx * as011 + px * (1. - qx)*(1. - rx) * as100;
1016 b_coef[
im][iw] = (1. - px)*(1. - qx)*(1. - rx) * ai000 + px * qx * rx * ai111
1017 + px * (1. - qx) * rx * ai101 + (1. - px) * qx * (1. - rx) * ai010
1018 + px * qx * (1. - rx) * ai110 + (1. - px)*(1. - qx) * rx * ai001
1019 + (1. - px) * qx * rx * ai011 + px * (1. - qx)*(1. - rx) * ai100;
1021 c_coef[
im][iw] = (1. - px)*(1. - qx)*(1. - rx) * ac000 + px * qx * rx * ac111
1022 + px * (1. - qx) * rx * ac101 + (1. - px) * qx * (1. - rx) * ac010
1023 + px * qx * (1. - rx) * ac110 + (1. - px)*(1. - qx) * rx * ac001
1024 + (1. - px) * qx * rx * ac011 + px * (1. - qx)*(1. - rx) * ac100;
1030 *
a = &a_coef[modnum][0];
1031 *
b = &b_coef[modnum][0];
1032 *
c = &c_coef[modnum][0];
1044 sq = sqrt(pow(
index, 2.0) - 1.0 + pow(mu, 2.0));
1045 r2 = pow((mu - sq) / (mu + sq), 2.0);
1046 q1 = (1.0 - pow(mu, 2.0) - mu * sq) / (1.0 - pow(mu, 2.0) + mu * sq);
1048 return (r2 * (q1 * q1 + 1.0) / 2.0);
1069 float **
a,
float **
b,
float **
c,
float **d,
float **e)
1071 static float lastsolz = -999.;
1072 static float lastsenz = -999.;
1073 static float lastphi = -999.;
1083 static float *
p, *
q, *
r;
1084 static float as000, as100, as010, as110, as001, as011, as101, as111;
1085 static float ai000, ai100, ai010, ai110, ai001, ai011, ai101, ai111;
1086 static float ac000, ac100, ac010, ac110, ac001, ac011, ac101, ac111;
1087 static float ad000, ad100, ad010, ad110, ad001, ad011, ad101, ad111;
1088 static float ae000, ae100, ae010, ae110, ae001, ae011, ae101, ae111;
1090 static int *isolz1, *isolz2;
1091 static int *isenz1, *isenz2;
1092 static int *iphi1, *iphi2;
1094 static float *p_ar, *q_ar, *r_ar;
1095 static float p_cnst, q_cnst, r_cnst;
1096 static int *isolz1_ar, *isolz2_ar, isolz1_cnst, isolz2_cnst;
1097 static int *isenz1_ar, *isenz2_ar, isenz1_cnst, isenz2_cnst;
1098 static int *iphi1_ar, *iphi2_ar, iphi1_cnst, iphi2_cnst;
1104 static int firstCall = 1;
1106 if (firstCall == 1) {
1109 if ((a_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1110 printf(
"Unable to allocate space for a_coef.\n");
1113 if ((b_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1114 printf(
"Unable to allocate space for b_coef.\n");
1117 if ((c_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1118 printf(
"Unable to allocate space for c_coef.\n");
1121 if ((d_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1122 printf(
"Unable to allocate space for d_coef.\n");
1126 if ((e_coef[
i] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
1127 printf(
"Unable to allocate space for e_coef.\n");
1133 if ((
geom->gmult == 0) || (interpol == 1)) {
1138 isolz1 = &isolz1_cnst;
1139 isolz2 = &isolz2_cnst;
1140 isenz1 = &isenz1_cnst;
1141 isenz2 = &isenz2_cnst;
1142 iphi1 = &iphi1_cnst;
1143 iphi2 = &iphi2_cnst;
1146 if (((p_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1148 ((q_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1150 ((r_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1152 printf(
"Unable to allocate space for p, q, r weights.\n");
1155 if (((isolz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1157 ((isenz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1159 ((iphi1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1161 printf(
"Unable to allocate space for interp indicies 1.\n");
1164 if (((isolz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1166 ((isenz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1168 ((iphi2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1170 printf(
"Unable to allocate space for interp indicies 2.\n");
1185 if (
geom->solz[0] != lastsolz ||
geom->senz[0] != lastsenz ||
1186 geom->phi[0] != lastphi) {
1187 for (
im = 0;
im < aertab->nmodel;
im++)
1190 for (iw = 0; iw < aertab->nwave; iw++) {
1193 for (
i = 0;
i < aertab->nsolz;
i++) {
1194 if (
geom->solz[ig] < aertab->solz[
i])
1197 isolz1[iw] =
MAX(
i - 1, 0);
1198 isolz2[iw] =
MIN(
i, aertab->nsolz - 1);
1199 if (isolz2[iw] != isolz1[iw])
1200 r[iw] = (
geom->solz[ig] - aertab->solz[isolz1[iw]]) /
1201 (aertab->solz[isolz2[iw]] - aertab->solz[isolz1[iw]]);
1206 for (
i = 0;
i < aertab->nsenz;
i++) {
1207 if (
geom->senz[ig] < aertab->senz[
i])
1210 isenz1[iw] =
MAX(
i - 1, 0);
1211 isenz2[iw] =
MIN(
i, aertab->nsenz - 1);
1212 if (isenz2[iw] != isenz1[iw])
1213 p[iw] = (
geom->senz[ig] - aertab->senz[isenz1[iw]]) /
1214 (aertab->senz[isenz2[iw]] - aertab->senz[isenz1[iw]]);
1220 for (
i = 0;
i < aertab->nphi;
i++) {
1221 if (aphi < aertab->phi[
i])
1224 iphi1[iw] =
MAX(
i - 1, 0);
1225 iphi2[iw] =
MIN(
i, aertab->nphi - 1);
1226 if (iphi2[iw] != iphi1[iw])
1227 q[iw] = (aphi - aertab->phi[iphi1[iw]]) /
1228 (aertab->phi[iphi2[iw]] - aertab->phi[iphi1[iw]]);
1231 if (gmult == 0)
break;
1235 lastsolz =
geom->solz[0];
1236 lastsenz =
geom->senz[0];
1237 lastphi =
geom->phi[0];
1247 for (iw = 0; iw < aertab->nwave; iw++) {
1252 if (isolz2[ig] == 0) {
1253 as000 = aertab->model[
im]->ams_all[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1254 as100 = aertab->model[
im]->ams_all[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1255 as001 = aertab->model[
im]->ams_all[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1256 as101 = aertab->model[
im]->ams_all[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1258 ai000 = aertab->model[
im]->bms_all[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1259 ai100 = aertab->model[
im]->bms_all[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1260 ai001 = aertab->model[
im]->bms_all[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1261 ai101 = aertab->model[
im]->bms_all[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1263 ac000 = aertab->model[
im]->cms_all[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1264 ac100 = aertab->model[
im]->cms_all[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1265 ac001 = aertab->model[
im]->cms_all[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1266 ac101 = aertab->model[
im]->cms_all[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1268 ad000 = aertab->model[
im]->dms_all[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1269 ad100 = aertab->model[
im]->dms_all[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1270 ad001 = aertab->model[
im]->dms_all[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1271 ad101 = aertab->model[
im]->dms_all[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1273 ae000 = aertab->model[
im]->ems_all[
INDEX(iw, isolz1[ig], 0, isenz1[ig])];
1274 ae100 = aertab->model[
im]->ems_all[
INDEX(iw, isolz1[ig], 0, isenz2[ig])];
1275 ae001 = aertab->model[
im]->ems_all[
INDEX(iw, isolz2[ig], 0, isenz1[ig])];
1276 ae101 = aertab->model[
im]->ems_all[
INDEX(iw, isolz2[ig], 0, isenz2[ig])];
1278 a_coef[
im][iw] = (1. - px)*(1. - rx) * as000 + px * rx * as101
1279 + (1. - px) * rx * as001 + px * (1. - rx) * as100;
1281 b_coef[
im][iw] = (1. - px)*(1. - rx) * ai000 + px * rx * ai101
1282 + (1. - px) * rx * ai001 + px * (1. - qx)*(1. - rx) * ai100;
1284 c_coef[
im][iw] = (1. - px)*(1. - rx) * ac000 + px * rx * ac101
1285 + (1. - px) * rx * ac001 + px * (1. - qx)*(1. - rx) * ac100;
1287 d_coef[
im][iw] = (1. - px)*(1. - rx) * ad000 + px * rx * ad101
1288 + (1. - px) * rx * ad001 + px * (1. - qx)*(1. - rx) * ad100;
1290 e_coef[
im][iw] = (1. - px)*(1. - rx) * ae000 + px * rx * ae101
1291 + (1. - px) * rx * ae001 + px * (1. - qx)*(1. - rx) * ae100;
1295 as000 = aertab->model[
im]->ams_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1296 as100 = aertab->model[
im]->ams_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1297 as010 = aertab->model[
im]->ams_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1298 as110 = aertab->model[
im]->ams_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1299 as001 = aertab->model[
im]->ams_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1300 as011 = aertab->model[
im]->ams_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1301 as101 = aertab->model[
im]->ams_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1302 as111 = aertab->model[
im]->ams_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1304 ai000 = aertab->model[
im]->bms_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1305 ai100 = aertab->model[
im]->bms_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1306 ai010 = aertab->model[
im]->bms_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1307 ai110 = aertab->model[
im]->bms_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1308 ai001 = aertab->model[
im]->bms_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1309 ai011 = aertab->model[
im]->bms_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1310 ai101 = aertab->model[
im]->bms_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1311 ai111 = aertab->model[
im]->bms_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1313 ac000 = aertab->model[
im]->cms_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1314 ac100 = aertab->model[
im]->cms_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1315 ac010 = aertab->model[
im]->cms_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1316 ac110 = aertab->model[
im]->cms_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1317 ac001 = aertab->model[
im]->cms_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1318 ac011 = aertab->model[
im]->cms_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1319 ac101 = aertab->model[
im]->cms_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1320 ac111 = aertab->model[
im]->cms_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1322 ad000 = aertab->model[
im]->dms_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1323 ad100 = aertab->model[
im]->dms_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1324 ad010 = aertab->model[
im]->dms_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1325 ad110 = aertab->model[
im]->dms_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1326 ad001 = aertab->model[
im]->dms_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1327 ad011 = aertab->model[
im]->dms_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1328 ad101 = aertab->model[
im]->dms_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1329 ad111 = aertab->model[
im]->dms_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1331 ae000 = aertab->model[
im]->ems_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz1[ig])];
1332 ae100 = aertab->model[
im]->ems_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz1[ig])];
1333 ae010 = aertab->model[
im]->ems_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz1[ig])];
1334 ae110 = aertab->model[
im]->ems_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz1[ig])];
1335 ae001 = aertab->model[
im]->ems_all[
INDEX(iw, isolz1[ig], iphi1[ig], isenz2[ig])];
1336 ae011 = aertab->model[
im]->ems_all[
INDEX(iw, isolz1[ig], iphi2[ig], isenz2[ig])];
1337 ae101 = aertab->model[
im]->ems_all[
INDEX(iw, isolz2[ig], iphi1[ig], isenz2[ig])];
1338 ae111 = aertab->model[
im]->ems_all[
INDEX(iw, isolz2[ig], iphi2[ig], isenz2[ig])];
1341 a_coef[
im][iw] = (1. - rx)*(1. - qx)*(1. - px) * as000 + rx * qx * px * as111
1342 + rx * (1. - qx) * px * as101 + (1. - rx) * qx * (1. - px) * as010
1343 + rx * qx * (1. - px) * as110 + (1. - rx)*(1. - qx) * px * as001
1344 + (1. - rx) * qx * px * as011 + rx * (1. - qx)*(1. - px) * as100;
1346 b_coef[
im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ai000 + rx * qx * px * ai111
1347 + rx * (1. - qx) * px * ai101 + (1. - rx) * qx * (1. - px) * ai010
1348 + rx * qx * (1. - px) * ai110 + (1. - rx)*(1. - qx) * px * ai001
1349 + (1. - rx) * qx * px * ai011 + rx * (1. - qx)*(1. - px) * ai100;
1351 c_coef[
im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ac000 + rx * qx * px * ac111
1352 + rx * (1. - qx) * px * ac101 + (1. - rx) * qx * (1. - px) * ac010
1353 + rx * qx * (1. - px) * ac110 + (1. - rx)*(1. - qx) * px * ac001
1354 + (1. - rx) * qx * px * ac011 + rx * (1. - qx)*(1. - px) * ac100;
1356 d_coef[
im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ad000 + rx * qx * px * ad111
1357 + rx * (1. - qx) * px * ad101 + (1. - rx) * qx * (1. - px) * ad010
1358 + rx * qx * (1. - px) * ad110 + (1. - rx)*(1. - qx) * px * ad001
1359 + (1. - rx) * qx * px * ad011 + rx * (1. - qx)*(1. - px) * ad100;
1361 e_coef[
im][iw] = (1. - rx)*(1. - qx)*(1. - px) * ae000 + rx * qx * px * ae111
1362 + rx * (1. - qx) * px * ae101 + (1. - rx) * qx * (1. - px) * ae010
1363 + rx * qx * (1. - px) * ae110 + (1. - rx)*(1. - qx) * px * ae001
1364 + (1. - rx) * qx * px * ae011 + rx * (1. - qx)*(1. - px) * ae100;
1371 if (
fabs(d_coef[modnum][iwnir_l]) > 1e-9) {
1372 printf(
"non zero cubic term found in longest NIR wavelength of aerosol table. Zia!!\n");
1377 *
a = &a_coef[modnum][0];
1378 *
b = &b_coef[modnum][0];
1379 *
c = &c_coef[modnum][0];
1380 *d = &d_coef[modnum][0];
1381 *e = &e_coef[modnum][0];
1393 static float lastsolz = -999.;
1394 static float lastsenz = -999.;
1395 static float lastphi = -999.;
1399 static float ***pc_coef;
1400 static int ntau_870,npc;
1402 static float *
p, *
q, *
r;
1403 static float as000, as100, as010, as110, as001, as011, as101, as111;
1405 static int *isolz1, *isolz2;
1406 static int *isenz1, *isenz2;
1407 static int *iphi1, *iphi2;
1409 static float *p_ar, *q_ar, *r_ar;
1410 static int *isolz1_ar, *isolz2_ar;
1411 static int *isenz1_ar, *isenz2_ar;
1412 static int *iphi1_ar, *iphi2_ar;
1417 int im, iw,
i, ig, itau,ipc;
1418 static int firstCall = 1;
1420 if (firstCall == 1) {
1422 ntau_870=aertab->ntau_870;
1425 pc_coef=(
float ***)calloc(
MAXMODEL,
sizeof(
float **));
1428 if ((pc_coef[
i] = (
float **) calloc(ntau_870,
sizeof (
float*))) ==
NULL) {
1429 printf(
"Unable to allocate space for pc_coef.\n");
1432 for(itau=0;itau<ntau_870;itau++){
1433 if ((pc_coef[
i][itau] = (
float *) calloc(aertab->npc, sizeof (
float))) ==
NULL) {
1434 printf(
"Unable to allocate space for pc_coef.\n");
1440 if ((
geom->gmult == 0) || (interpol == 1)) {
1445 if (((p_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1447 ((q_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1449 ((r_ar = (
float *) malloc(aertab->nwave * sizeof (
float)))
1451 printf(
"Unable to allocate space for p, q, r weights.\n");
1454 if (((isolz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1456 ((isenz1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1458 ((iphi1_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1460 printf(
"Unable to allocate space for interp indicies 1.\n");
1463 if (((isolz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1465 ((isenz2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1467 ((iphi2_ar = (
int *) malloc(aertab->nwave * sizeof (
int)))
1469 printf(
"Unable to allocate space for interp indicies 2.\n");
1483 if (
geom->solz[0] != lastsolz ||
geom->senz[0] != lastsenz ||
1484 geom->phi[0] != lastphi) {
1485 for (
im = 0;
im < aertab->nmodel;
im++)
1488 for (iw = 0; iw < aertab->nwave; iw++) {
1491 for (
i = 0;
i < aertab->nsolz;
i++) {
1492 if (
geom->solz[ig] < aertab->solz[
i])
1495 isolz1[iw] =
MAX(
i - 1, 0);
1496 isolz2[iw] =
MIN(
i, aertab->nsolz - 1);
1497 if (isolz2[iw] != isolz1[iw])
1498 r[iw] = (
geom->solz[ig] - aertab->solz[isolz1[iw]]) /
1499 (aertab->solz[isolz2[iw]] - aertab->solz[isolz1[iw]]);
1504 for (
i = 0;
i < aertab->nsenz;
i++) {
1505 if (
geom->senz[ig] < aertab->senz[
i])
1508 isenz1[iw] =
MAX(
i - 1, 0);
1509 isenz2[iw] =
MIN(
i, aertab->nsenz - 1);
1510 if (isenz2[iw] != isenz1[iw])
1511 p[iw] = (
geom->senz[ig] - aertab->senz[isenz1[iw]]) /
1512 (aertab->senz[isenz2[iw]] - aertab->senz[isenz1[iw]]);
1518 for (
i = 0;
i < aertab->nphi;
i++) {
1519 if (aphi < aertab->phi[
i])
1522 iphi1[iw] =
MAX(
i - 1, 0);
1523 iphi2[iw] =
MIN(
i, aertab->nphi - 1);
1524 if (iphi2[iw] != iphi1[iw])
1525 q[iw] = (aphi - aertab->phi[iphi1[iw]]) /
1526 (aertab->phi[iphi2[iw]] - aertab->phi[iphi1[iw]]);
1529 if (gmult == 0)
break;
1533 lastsolz =
geom->solz[0];
1534 lastsenz =
geom->senz[0];
1535 lastphi =
geom->phi[0];
1550 if (isolz2[ig] == 0) {
1552 for(itau=0;itau<ntau_870;itau++)
1553 for(ipc=0;ipc<npc;ipc++){
1554 as000 = aertab->model[
im]->pc_rhoa[itau][isolz1[ig]][0][isenz1[ig]][ipc];
1555 as100 = aertab->model[
im]->pc_rhoa[itau][isolz1[ig]][0][isenz2[ig]][ipc];
1556 as001 = aertab->model[
im]->pc_rhoa[itau][isolz2[ig]][0][isenz1[ig]][ipc];
1557 as101 = aertab->model[
im]->pc_rhoa[itau][isolz2[ig]][0][isenz2[ig]][ipc];
1559 pc_coef[
im][itau][ipc] = (1. - px)*(1. - rx) * as000 + px * rx * as101
1560 + (1. - px) * rx * as001 + px * (1. - rx) * as100;
1565 for(itau=0;itau<ntau_870;itau++)
1566 for(ipc=0;ipc<npc;ipc++){
1567 as000 = aertab->model[
im]->pc_rhoa[itau][isolz1[ig]][iphi1[ig]][isenz1[ig]][ipc];
1568 as100 = aertab->model[
im]->pc_rhoa[itau][isolz2[ig]][iphi1[ig]][isenz1[ig]][ipc];
1569 as010 = aertab->model[
im]->pc_rhoa[itau][isolz1[ig]][iphi2[ig]][isenz1[ig]][ipc];
1570 as110 = aertab->model[
im]->pc_rhoa[itau][isolz2[ig]][iphi2[ig]][isenz1[ig]][ipc];
1571 as001 = aertab->model[
im]->pc_rhoa[itau][isolz1[ig]][iphi1[ig]][isenz2[ig]][ipc];
1572 as011 = aertab->model[
im]->pc_rhoa[itau][isolz1[ig]][iphi2[ig]][isenz2[ig]][ipc];
1573 as101 = aertab->model[
im]->pc_rhoa[itau][isolz2[ig]][iphi1[ig]][isenz2[ig]][ipc];
1574 as111 = aertab->model[
im]->pc_rhoa[itau][isolz2[ig]][iphi2[ig]][isenz2[ig]][ipc];
1576 pc_coef[
im][itau][ipc] = (1. - rx)*(1. - qx)*(1. - px) * as000 + rx * qx * px * as111
1577 + rx * (1. - qx) * px * as101 + (1. - rx) * qx * (1. - px) * as010
1578 + rx * qx * (1. - px) * as110 + (1. - rx)*(1. - qx) * px * as001
1579 + (1. - rx) * qx * px * as011 + rx * (1. - qx)*(1. - px) * as100;
1586 for(itau=0;itau<ntau_870;itau++)
1587 for(ipc=0;ipc<npc;ipc++){
1588 pc_rhoa[itau][ipc]=pc_coef[modnum][itau][ipc];
1609 return (
x->eps_obs <
y->eps_obs ? -1 : 1);
1612 void model_select_ahmad(int32_t nmodels, int32_t *mindx,
float eps_pred[],
float eps_obs, int32_t *modmin,
1613 int32_t *modmax,
float *modrat) {
1620 for (
im = 0;
im < nmodels;
im++) {
1621 epsilonT[
im].eps_obs = eps_pred[
im];
1622 epsilonT[
im].modnum =
im;
1628 qsort(epsilonT, nmodels,
sizeof (epsilonTstr),
1633 for (
im = 0;
im < nmodels;
im++) {
1634 if (eps_obs < epsilonT[
im].eps_obs)
1645 im1 =
MAX(
MIN(
im - 1, nmodels - 1), 0);
1646 im2 =
MAX(
MIN(
im, nmodels - 1), 0);
1652 *modmin = epsilonT[im1].modnum;
1653 *modmax = epsilonT[im2].modnum;
1654 *modrat = (eps_obs - epsilonT[im1].eps_obs) / (epsilonT[im2].eps_obs - epsilonT[im1].eps_obs);
1658 if (*modmin == *modmax)
1672 float tau_iwnir_l, int32_t modl,
float tau_pred[],
float rho_pred[],
float derv_rhoa[],
float derv_taua_tmp[])
1674 float *
ac, *
bc, *cc, *dc, *ec;
1675 float ext_modl[nwave];
1676 float lg_tau_pred[nwave];
1677 float lg_rho_pred[nwave];
1689 for (iw = 0; iw < nwave; iw++) {
1691 ext_modl[iw] = aertab->model[modl]->extc[iwtab];
1695 for (iw = 0; iw < nwave; iw++) {
1696 tau_pred[iw] = (ext_modl[iw] / ext_modl[iwnir_l]) * tau_iwnir_l;
1697 lg_tau_pred[iw] = log(tau_pred[iw]);
1701 derv_taua_tmp[iw] = ext_modl[iw] / ext_modl[iwnir_l];
1702 derv_rhoa[iw] = 1 / tau_pred[iw]
1703 * (ext_modl[iw] / ext_modl[iwnir_l]);
1709 for (iw = 0; iw < nwave; iw++) {
1711 lg_rho_pred[iw] =
ac[iwtab] +
1712 bc[iwtab] * lg_tau_pred[iw] +
1713 cc[iwtab] * pow(lg_tau_pred[iw], 2) +
1714 dc[iwtab] * pow(lg_tau_pred[iw], 3) +
1715 ec[iwtab] * pow(lg_tau_pred[iw], 4);
1716 rho_pred[iw] = exp(lg_rho_pred[iw]);
1719 derv_rhoa[iw] *= (
bc[iwtab] + 2 * cc[iwtab] * lg_tau_pred[iw]
1720 + 3 * dc[iwtab] * pow(lg_tau_pred[iw], 2)
1721 + 4 * ec[iwtab] * pow(lg_tau_pred[iw], 3));
1722 derv_rhoa[iw] *= rho_pred[iw];
1737 float tau_iwnir_l, int32_t modl,
float tau_pred[],
float rho_pred[],
float derv_rhoa[],
float derv_taua_tmp[])
1739 static int firstcall=1;
1740 static float **pc_rhoa;
1741 static int npc, ntau_870;
1742 aermodstr *aermod=aertab->model[modl];
1744 int iw,itau,ipc, itau1,itau2;
1749 ntau_870=aertab->ntau_870;
1751 pc_rhoa=(
float **)malloc(aertab->ntau_870*
sizeof(
float *));
1753 for(itau=0;itau<ntau_870;itau++)
1754 pc_rhoa[itau]=(
float *)malloc(aertab->npc*
sizeof(
float));
1763 for(itau=0;itau<ntau_870;itau++){
1764 if(tau_iwnir_l<aermod->tau_870[itau])
1767 itau1=
MAX(itau-1,0);
1768 itau2=
MIN(itau,ntau_870-1);
1772 if(itau2==ntau_870-1)
1777 for(iw=0;iw<nwave;iw++){
1781 for(ipc=0;ipc<npc;ipc++){
1782 rhoa1+=pc_rhoa[itau1][ipc]*aermod->pc_components_rhoa[ipc][iw];
1783 rhoa2+=pc_rhoa[itau2][ipc]*aermod->pc_components_rhoa[ipc][iw];
1785 rhoa1=rhoa1+aermod->pc_mean_rhoa[iw];
1786 rhoa2=rhoa2+aermod->pc_mean_rhoa[iw];
1787 rho_pred[iw]=rhoa1+(rhoa2-rhoa1)*(tau_iwnir_l-aermod->tau_870[itau1])/(aermod->tau_870[itau2]-aermod->tau_870[itau1]);
1788 tau_pred[iw]=tau_iwnir_l*aermod->extc[iw]/aermod->extc[iwnir_l];
1791 derv_taua_tmp[iw]=aermod->extc[iw]/aermod->extc[iwnir_l];
1792 derv_rhoa [iw]=(rhoa2-rhoa1)/(aermod->tau_870[itau2]-aermod->tau_870[itau1]);
1807 float tau_iwnir_l, int32_t modl,
float tau_pred[],
float rho_pred[])
1809 float *
ac, *
bc, *cc, *dc, *ec;
1810 float ext_modl[nwave];
1811 float lg_tau_pred[nwave];
1812 float lg_rho_pred[nwave];
1824 for (iw = 0; iw < nwave; iw++) {
1826 ext_modl[iw] = aertab->model[modl]->extc[iwtab];
1830 for (iw = 0; iw < nwave; iw++) {
1831 tau_pred[iw] = (ext_modl[iw] / ext_modl[iwnir_l]) * tau_iwnir_l;
1832 lg_tau_pred[iw] = (tau_pred[iw]);
1837 for (iw = 0; iw < nwave; iw++) {
1839 lg_rho_pred[iw] =
ac[iwtab] +
1840 bc[iwtab] * lg_tau_pred[iw] +
1841 cc[iwtab] * pow(lg_tau_pred[iw], 2) +
1842 dc[iwtab] * pow(lg_tau_pred[iw], 3) +
1843 ec[iwtab] * pow(lg_tau_pred[iw], 4);
1844 rho_pred[iw] = (lg_rho_pred[iw]);
1859 int32_t nmodels, int32_t mindx[],
1860 geom_str *
geom,
float wv,
float rhoa[],
1861 int32_t *modmin, int32_t *modmax,
float *modrat,
float *epsnir,
1862 float tau_pred_max[],
float tau_pred_min[],
float rho_pred_max[],
float rho_pred_min[],
1863 float tau_aer[],
float rho_aer[],
int ip,uncertainty_t *uncertainty) {
1865 static int firstcall=1;
1866 static float **pc_rhoa;
1867 static int npc, ntau_870;
1869 float tau_iwnir_l[nmodels];
1870 float rho_iwnir_s_pred[nmodels];
1871 float eps_pred[nmodels];
1876 float derv_taua_rhoa[nmodels];
1877 float derv_eps_rhoa_l[nmodels];
1878 float *derv_rhoa_min=
NULL,*derv_rhoa_max=
NULL;
1879 float *derv_taua_min=
NULL,*derv_taua_max=
NULL;
1880 float derv_mwt_rhoa_s=0., derv_mwt_rhoa_l=0.,derv_mwt_taua_l=0.,derv_mwt_rhow_l=0.;
1884 float derv_eps_Lrc_s,derv_eps_Lrc_l,derv_eps_taua_l,derv_eps_rhow_l,derv_Lg_taua_l;
1885 float derv_eps_mod[4][nmodels];
1888 derv_eps_Lrc_s = uncertainty->derv_eps_Lrc_s;
1889 derv_eps_Lrc_l = uncertainty->derv_eps_Lrc_l;
1890 derv_eps_taua_l = uncertainty->derv_eps_taua_l;
1891 derv_eps_rhow_l = uncertainty->derv_eps_rhow_l;
1892 derv_Lg_taua_l = uncertainty->derv_Lg_taua[iwnir_l];
1894 derv_rhoa_min=(
float *)malloc(nwave*
sizeof(
float));
1895 derv_rhoa_max=(
float *)malloc(nwave*
sizeof(
float));
1896 derv_taua_min=(
float *)malloc(nwave*
sizeof(
float));
1897 derv_taua_max=(
float *)malloc(nwave*
sizeof(
float));
1903 static float *rhoa_modl_s,*rhoa_modl_l;
1907 ntau_870=aertab->ntau_870;
1909 rhoa_modl_s=(
float *)malloc(aertab->ntau_870*
sizeof(
float));
1910 rhoa_modl_l=(
float *)malloc(aertab->ntau_870*
sizeof(
float));
1911 pc_rhoa=(
float **)malloc(aertab->ntau_870*
sizeof(
float *));
1913 for(itau=0;itau<aertab->ntau_870;itau++)
1914 pc_rhoa[itau]=(
float *)malloc(aertab->npc*
sizeof(
float));
1920 eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
1930 for (
im = 0;
im < nmodels;
im++) {
1933 aermod=aertab->model[modl];
1937 for(itau=0;itau<ntau_870;itau++){
1938 rhoa_modl_l[itau]=0;
1939 rhoa_modl_s[itau]=0;
1941 for(ipc=0;ipc<npc;ipc++){
1942 rhoa_modl_s[itau]+=pc_rhoa[itau][ipc]*aermod->pc_components_rhoa[ipc][iwnir_s];
1943 rhoa_modl_l[itau]+=pc_rhoa[itau][ipc]*aermod->pc_components_rhoa[ipc][iwnir_l];
1945 rhoa_modl_s[itau]+=aermod->pc_mean_rhoa[iwnir_s];
1946 rhoa_modl_l[itau]+=aermod->pc_mean_rhoa[iwnir_l];
1950 if(rhoa[iwnir_l]<rhoa_modl_l[0])
1952 else if (rhoa[iwnir_l]>rhoa_modl_l[ntau_870-1])
1955 for(itau=1;itau<ntau_870;itau++){
1956 if(rhoa[iwnir_l]<rhoa_modl_l[itau] && rhoa[iwnir_l]>rhoa_modl_l[itau-1])
1960 tau_iwnir_l[
im]=aermod->tau_870[itau]+(aermod->tau_870[itau]-aermod->tau_870[itau-1])*(rhoa[iwnir_l]-rhoa_modl_l[itau])/(rhoa_modl_l[itau]-rhoa_modl_l[itau-1]);
1961 rho_iwnir_s_pred[
im]=rhoa_modl_s[itau]+(rhoa_modl_s[itau]-rhoa_modl_s[itau-1])*(tau_iwnir_l[
im]-aermod->tau_870[itau])/(aermod->tau_870[itau]-aermod->tau_870[itau-1]);
1964 eps_pred[
im] = rho_iwnir_s_pred[
im] / rhoa[iwnir_l];
1968 derv_taua_rhoa[
im]=(aermod->tau_870[itau]-aermod->tau_870[itau-1])/(rhoa_modl_l[itau]-rhoa_modl_l[itau-1]);
1969 derv_eps_rhoa_l[
im]=derv_taua_rhoa[
im]*((rhoa_modl_s[itau]-rhoa_modl_s[itau-1])/(aermod->tau_870[itau]-aermod->tau_870[itau-1]));
1971 derv_eps_rhoa_l[
im]=derv_eps_rhoa_l[
im]/rhoa[iwnir_l]-rho_iwnir_s_pred[
im]/rhoa[iwnir_l]/rhoa[iwnir_l];
1972 derv_eps_mod[1][
im]=derv_eps_rhoa_l[
im];
1973 derv_eps_mod[2][
im] = -derv_eps_rhoa_l[
im] * uncertainty->derv_Lg_taua[iwnir_l];
1974 derv_eps_mod[3][
im] = -derv_eps_rhoa_l[
im];
1993 for(iw=iwnir_s;iw<=iwnir_l;iw++)
1994 uncertainty->derv_modrat_rhorc[iw-iwnir_s]=0.;
1995 uncertainty->derv_modrat_taua_l=0.;
1998 derv_mwt_rhoa_s = 1 / (eps_pred[im2] - eps_pred[im1])
2003 derv_mwt_rhoa_l = 1 / (eps_pred[im2] - eps_pred[im1])
2005 derv_mwt_rhoa_l -= ((mwt / (eps_pred[im2] - eps_pred[im1])
2006 * derv_eps_mod[1][im2]));
2007 derv_mwt_rhoa_l += ((eps_obs - eps_pred[im2])
2008 / pow(eps_pred[im2] - eps_pred[im1], 2)
2009 * derv_eps_mod[1][im1]);
2011 derv_mwt_taua_l = 1 / (eps_pred[im2] - eps_pred[im1])
2013 derv_mwt_taua_l -= ((mwt / (eps_pred[im2] - eps_pred[im1])
2014 * derv_eps_mod[2][im2]));
2015 derv_mwt_taua_l += ((eps_obs - eps_pred[im2])
2016 / pow(eps_pred[im2] - eps_pred[im1], 2)
2017 * derv_eps_mod[2][im1]);
2019 derv_mwt_rhow_l = 1 / (eps_pred[im2] - eps_pred[im1])
2021 derv_mwt_rhow_l -= ((mwt / (eps_pred[im2] - eps_pred[im1])
2022 * derv_eps_mod[3][im2]));
2023 derv_mwt_rhow_l += ((eps_obs - eps_pred[im2])
2024 / pow(eps_pred[im2] - eps_pred[im1], 2)
2025 * derv_eps_mod[3][im1]);
2027 uncertainty->derv_modrat_rhorc[0] = derv_mwt_rhoa_s;
2028 uncertainty->derv_modrat_rhorc[iwnir_l-iwnir_s] = derv_mwt_rhoa_l;
2029 uncertainty->derv_modrat_taua_l = derv_mwt_taua_l;
2030 uncertainty->derv_modrat_rhow_l = derv_mwt_rhow_l;
2044 *modmin = mindx[im1];
2045 *modmax = mindx[im2];
2050 if(
comp_rhoa_pca(nwave, wave,
geom, tau_iwnir_l[im1], *modmin, tau_pred_min, rho_pred_min,derv_rhoa_min,derv_taua_min) !=0)
2052 if(
comp_rhoa_pca(nwave, wave,
geom, tau_iwnir_l[im2], *modmax, tau_pred_max, rho_pred_max,derv_rhoa_max,derv_taua_max) !=0)
2057 for (iw = 0; iw < nwave; iw++) {
2058 tau_aer[iw] = (1.0 - mwt) * tau_pred_min[iw] + mwt * tau_pred_max[iw];
2059 rho_aer[iw] = (1.0 - mwt) * rho_pred_min[iw] + mwt * rho_pred_max[iw];
2063 uncertainty->derv_La_rhorc[iw][iwnir_l-iwnir_s] = (1.0 - mwt) * derv_rhoa_min[iw]
2064 * derv_taua_rhoa[im1]
2065 + mwt * derv_rhoa_max[iw] * derv_taua_rhoa[im2]
2066 + (rho_pred_max[iw] - rho_pred_min[iw]) * derv_mwt_rhoa_l;
2067 uncertainty->derv_La_taua_l[iw] = -(1.0 - mwt) * derv_rhoa_min[iw]
2068 * derv_taua_rhoa[im1] * derv_Lg_taua_l
2069 - mwt * derv_rhoa_max[iw] * derv_taua_rhoa[im2]
2071 + (rho_pred_max[iw] - rho_pred_min[iw]) * derv_mwt_taua_l;
2072 uncertainty->derv_La_rhorc[iw][0] = (rho_pred_max[iw] - rho_pred_min[iw])
2074 uncertainty->derv_La_rhow_l[iw] = -(1.0 - mwt) * derv_rhoa_min[iw]
2075 * derv_taua_rhoa[im1]
2076 - mwt * derv_rhoa_max[iw] * derv_taua_rhoa[im2]
2077 + (rho_pred_max[iw] - rho_pred_min[iw]) * derv_mwt_rhow_l;
2081 uncertainty->derv_taua_rhorc[iw][iwnir_l-iwnir_s]= (1.0 - mwt) * derv_taua_min[iw]
2082 * derv_taua_rhoa[im1]
2083 + mwt * derv_taua_max[iw] * derv_taua_rhoa[im2]
2084 + (tau_pred_max[iw] - tau_pred_min[iw]) * derv_mwt_rhoa_l;
2085 uncertainty->derv_taua_taua_l[iw] = -(1.0 - mwt) * derv_taua_min[iw]
2086 * derv_taua_rhoa[im1] * derv_Lg_taua_l
2087 - mwt * derv_taua_max[iw] * derv_taua_rhoa[im2]
2089 + (tau_pred_max[iw] - tau_pred_min[iw]) * derv_mwt_taua_l;
2090 uncertainty->derv_taua_rhorc[iw][0] = (tau_pred_max[iw] - tau_pred_min[iw])
2092 uncertainty->derv_taua_rhow_l[iw] = -(1.0 - mwt) * derv_taua_min[iw]
2093 * derv_taua_rhoa[im1]
2094 - mwt * derv_taua_max[iw] * derv_taua_rhoa[im2]
2095 + (tau_pred_max[iw] - tau_pred_min[iw]) * derv_mwt_rhow_l;
2100 uncertainty->derv_taua_min_rhorc_l[iw] = derv_taua_min[iw] * derv_taua_rhoa[im1];
2101 uncertainty->derv_taua_min_taua_l[iw] = -derv_taua_min[iw] * derv_taua_rhoa[im1]
2103 uncertainty->derv_taua_min_rhow_l[iw] = -derv_taua_min[iw] * derv_taua_rhoa[im1];
2107 uncertainty->derv_taua_max_rhorc_l[iw]= derv_taua_max[iw] * derv_taua_rhoa[im2];
2108 uncertainty->derv_taua_max_taua_l [iw]= -derv_taua_max[iw] * derv_taua_rhoa[im2]
2110 uncertainty->derv_taua_max_rhow_l [iw]= -derv_taua_max[iw] * derv_taua_rhoa[im2];
2116 free(derv_rhoa_min);
2117 free(derv_rhoa_max);
2118 free(derv_taua_min);
2119 free(derv_taua_max);
2131 int32_t nmodels, int32_t mindx[],
2132 geom_str *
geom,
float wv,
float rhoa[],
2133 int32_t *modmin, int32_t *modmax,
float *modrat,
float *epsnir,
2134 float tau_pred_max[],
float tau_pred_min[],
float rho_pred_max[],
float rho_pred_min[],
2135 float tau_aer[],
float rho_aer[],
int ip,uncertainty_t *uncertainty) {
2143 eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
2147 uncertainty->derv_eps_Lrc_s = 1 / rhoa[iwnir_l];
2148 uncertainty->derv_eps_Lrc_l = -eps_obs / rhoa[iwnir_l];
2150 if (uncertainty->ratio_rhow[iwnir_s - iwnir_s] != 0.) {
2151 uncertainty->derv_eps_rhow_l =
2152 -1 / (rhoa[iwnir_l]) * uncertainty->ratio_rhow[iwnir_s - iwnir_s];
2153 uncertainty->derv_eps_rhow_l += eps_obs / rhoa[iwnir_l];
2156 status =
ahmad_atm_corr_pca(
sensorID,wave, nwave, iwnir_s, iwnir_l, nmodels, mindx,
geom, wv, rhoa,
2157 modmin, modmax, modrat, epsnir, tau_pred_max, tau_pred_min, rho_pred_max,
2158 rho_pred_min, tau_aer, rho_aer, ip, uncertainty);
2161 float *
ac, *
bc, *cc, *dc, *ec;
2162 float ext_iwnir_l, ext_iwnir_s;
2163 double ax, bx, cx, fx;
2164 float tau_iwnir_l[nmodels];
2165 float lg_tau_iwnir_s[nmodels], tau_iwnir_s[nmodels];
2166 float lg_rho_iwnir_s_pred[nmodels], rho_iwnir_s_pred[nmodels];
2167 float eps_pred[nmodels];
2171 float lg_tau_iwnir_l;
2172 int iwtab_l, iwtab_s;
2174 float derv_taua_rhoa[nmodels];
2175 float derv_eps_rhoa_l[nmodels];
2176 float *derv_rhoa_min=
NULL,*derv_rhoa_max=
NULL;
2177 float *derv_taua_min=
NULL,*derv_taua_max=
NULL;
2178 float derv_mwt_rhoa_s=0., derv_mwt_rhoa_l=0.,derv_mwt_taua_l=0.,derv_mwt_rhow_l=0.;
2180 float derv_eps_Lrc_s,derv_eps_Lrc_l,derv_eps_taua_l,derv_eps_rhow_l,derv_Lg_taua_l;
2181 float derv_eps_mod[4][nmodels];
2191 derv_rhoa_min = (
float *)malloc(nwave *
sizeof(
float));
2192 derv_rhoa_max = (
float *)malloc(nwave *
sizeof(
float));
2193 derv_taua_min = (
float *)malloc(nwave *
sizeof(
float));
2194 derv_taua_max = (
float *)malloc(nwave *
sizeof(
float));
2196 derv_eps_Lrc_s = uncertainty->derv_eps_Lrc_s;
2197 derv_eps_Lrc_l = uncertainty->derv_eps_Lrc_l;
2198 derv_eps_taua_l = uncertainty->derv_eps_taua_l;
2199 derv_eps_rhow_l = uncertainty->derv_eps_rhow_l;
2200 derv_Lg_taua_l = uncertainty->derv_Lg_taua[iwnir_l];
2203 for (
im = 0;
im < nmodels;
im++) {
2214 iwtab_l = iwatab[iwnir_l];
2215 iwtab_s = iwatab[iwnir_s];
2217 ax = (
double)
ac[iwtab_l] - log((
double) rhoa[iwnir_l]);
2219 cx = (
double) cc[iwtab_l];
2221 fx = bx * bx - 4.0 *
ax*cx;
2222 if (fx > 0.0 && cx != 0.0) {
2223 lg_tau_iwnir_l = 0.5 * (-bx + sqrt(fx)) / cx;
2225 derv_taua_rhoa[
im] = 1 / rhoa[iwnir_l] / sqrt(fx);
2227 tau_iwnir_l[
im] = exp(lg_tau_iwnir_l);
2229 derv_taua_rhoa[
im]*=tau_iwnir_l[
im];
2237 ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
2238 ext_iwnir_s = aertab->model[modl]->extc[iwtab_s];
2240 tau_iwnir_s[
im] = (ext_iwnir_s / ext_iwnir_l) * tau_iwnir_l[
im];
2241 lg_tau_iwnir_s[
im] = log(tau_iwnir_s[
im]);
2244 derv_eps_rhoa_l[
im]=1/tau_iwnir_s[
im]*(ext_iwnir_s / ext_iwnir_l)*derv_taua_rhoa[
im];
2249 lg_rho_iwnir_s_pred[
im] =
ac[iwtab_s] +
bc[iwtab_s] * lg_tau_iwnir_s[
im] +
2250 cc[iwtab_s] * pow(lg_tau_iwnir_s[
im], 2) +
2251 dc[iwtab_s] * pow(lg_tau_iwnir_s[
im], 3) +
2252 ec[iwtab_s] * pow(lg_tau_iwnir_s[
im], 4);
2255 derv_eps_rhoa_l[
im]=(
bc[iwtab_s]+ 2*cc[iwtab_s]*lg_tau_iwnir_s[
im]+3*dc[iwtab_s]*pow(lg_tau_iwnir_s[
im], 2)+4*ec[iwtab_s] * pow(lg_tau_iwnir_s[
im], 3) )*derv_eps_rhoa_l[
im];
2257 rho_iwnir_s_pred[
im] = exp(lg_rho_iwnir_s_pred[
im]);
2260 derv_eps_rhoa_l[
im]=rho_iwnir_s_pred[
im]*derv_eps_rhoa_l[
im];
2264 eps_pred[
im] = rho_iwnir_s_pred[
im] / rhoa[iwnir_l];
2268 derv_eps_rhoa_l[
im]=derv_eps_rhoa_l[
im]/rhoa[iwnir_l]-rho_iwnir_s_pred[
im]/rhoa[iwnir_l]/rhoa[iwnir_l];
2269 derv_eps_mod[1][
im]=derv_eps_rhoa_l[
im];
2270 derv_eps_mod[2][
im] = -derv_eps_rhoa_l[
im] * uncertainty->derv_Lg_taua[iwnir_l];
2271 derv_eps_mod[3][
im] = -derv_eps_rhoa_l[
im];
2291 for(iw=iwnir_s;iw<=iwnir_l;iw++)
2292 uncertainty->derv_modrat_rhorc[iw-iwnir_s]=0.;
2293 uncertainty->derv_modrat_taua_l=0.;
2296 derv_mwt_rhoa_s = 1 / (eps_pred[im2] - eps_pred[im1])
2301 derv_mwt_rhoa_l = 1 / (eps_pred[im2] - eps_pred[im1])
2303 derv_mwt_rhoa_l -= ((mwt / (eps_pred[im2] - eps_pred[im1])
2304 * derv_eps_mod[1][im2]));
2305 derv_mwt_rhoa_l += ((eps_obs - eps_pred[im2])
2306 / pow(eps_pred[im2] - eps_pred[im1], 2)
2307 * derv_eps_mod[1][im1]);
2309 derv_mwt_taua_l = 1 / (eps_pred[im2] - eps_pred[im1])
2311 derv_mwt_taua_l -= ((mwt / (eps_pred[im2] - eps_pred[im1])
2312 * derv_eps_mod[2][im2]));
2313 derv_mwt_taua_l += ((eps_obs - eps_pred[im2])
2314 / pow(eps_pred[im2] - eps_pred[im1], 2)
2315 * derv_eps_mod[2][im1]);
2317 derv_mwt_rhow_l = 1 / (eps_pred[im2] - eps_pred[im1])
2319 derv_mwt_rhow_l -= ((mwt / (eps_pred[im2] - eps_pred[im1])
2320 * derv_eps_mod[3][im2]));
2321 derv_mwt_rhow_l += ((eps_obs - eps_pred[im2])
2322 / pow(eps_pred[im2] - eps_pred[im1], 2)
2323 * derv_eps_mod[3][im1]);
2325 uncertainty->derv_modrat_rhorc[0] = derv_mwt_rhoa_s;
2326 uncertainty->derv_modrat_rhorc[iwnir_l-iwnir_s] = derv_mwt_rhoa_l;
2327 uncertainty->derv_modrat_taua_l = derv_mwt_taua_l;
2328 uncertainty->derv_modrat_rhow_l = derv_mwt_rhow_l;
2342 *modmin = mindx[im1];
2343 *modmax = mindx[im2];
2348 comp_rhoa_ms_eps(nwave, wave,
geom, tau_iwnir_l[im1], *modmin, tau_pred_min, rho_pred_min,derv_rhoa_min,derv_taua_min);
2349 comp_rhoa_ms_eps(nwave, wave,
geom, tau_iwnir_l[im2], *modmax, tau_pred_max, rho_pred_max,derv_rhoa_max,derv_taua_max);
2353 for (iw = 0; iw < nwave; iw++) {
2354 tau_aer[iw] = (1.0 - mwt) * tau_pred_min[iw] + mwt * tau_pred_max[iw];
2355 rho_aer[iw] = (1.0 - mwt) * rho_pred_min[iw] + mwt * rho_pred_max[iw];
2359 uncertainty->derv_La_rhorc[iw][iwnir_l-iwnir_s] = (1.0 - mwt) * derv_rhoa_min[iw]
2360 * derv_taua_rhoa[im1]
2361 + mwt * derv_rhoa_max[iw] * derv_taua_rhoa[im2]
2362 + (rho_pred_max[iw] - rho_pred_min[iw]) * derv_mwt_rhoa_l;
2363 uncertainty->derv_La_taua_l[iw] = -(1.0 - mwt) * derv_rhoa_min[iw]
2364 * derv_taua_rhoa[im1] * derv_Lg_taua_l
2365 - mwt * derv_rhoa_max[iw] * derv_taua_rhoa[im2]
2367 + (rho_pred_max[iw] - rho_pred_min[iw]) * derv_mwt_taua_l;
2368 uncertainty->derv_La_rhorc[iw][0] = (rho_pred_max[iw] - rho_pred_min[iw])
2370 uncertainty->derv_La_rhow_l[iw] = -(1.0 - mwt) * derv_rhoa_min[iw]
2371 * derv_taua_rhoa[im1]
2372 - mwt * derv_rhoa_max[iw] * derv_taua_rhoa[im2]
2373 + (rho_pred_max[iw] - rho_pred_min[iw]) * derv_mwt_rhow_l;
2377 uncertainty->derv_taua_rhorc[iw][iwnir_l-iwnir_s]= (1.0 - mwt) * derv_taua_min[iw]
2378 * derv_taua_rhoa[im1]
2379 + mwt * derv_taua_max[iw] * derv_taua_rhoa[im2]
2380 + (tau_pred_max[iw] - tau_pred_min[iw]) * derv_mwt_rhoa_l;
2381 uncertainty->derv_taua_taua_l[iw] = -(1.0 - mwt) * derv_taua_min[iw]
2382 * derv_taua_rhoa[im1] * derv_Lg_taua_l
2383 - mwt * derv_taua_max[iw] * derv_taua_rhoa[im2]
2385 + (tau_pred_max[iw] - tau_pred_min[iw]) * derv_mwt_taua_l;
2386 uncertainty->derv_taua_rhorc[iw][0] = (tau_pred_max[iw] - tau_pred_min[iw])
2388 uncertainty->derv_taua_rhow_l[iw] = -(1.0 - mwt) * derv_taua_min[iw]
2389 * derv_taua_rhoa[im1]
2390 - mwt * derv_taua_max[iw] * derv_taua_rhoa[im2]
2391 + (tau_pred_max[iw] - tau_pred_min[iw]) * derv_mwt_rhow_l;
2396 uncertainty->derv_taua_min_rhorc_l[iw] = derv_taua_min[iw] * derv_taua_rhoa[im1];
2397 uncertainty->derv_taua_min_taua_l[iw] = -derv_taua_min[iw] * derv_taua_rhoa[im1]
2399 uncertainty->derv_taua_min_rhow_l[iw] = -derv_taua_min[iw] * derv_taua_rhoa[im1];
2403 uncertainty->derv_taua_max_rhorc_l[iw]= derv_taua_max[iw] * derv_taua_rhoa[im2];
2404 uncertainty->derv_taua_max_taua_l [iw]= -derv_taua_max[iw] * derv_taua_rhoa[im2]
2406 uncertainty->derv_taua_max_rhow_l [iw]= -derv_taua_max[iw] * derv_taua_rhoa[im2];
2412 free(derv_rhoa_min);
2413 free(derv_rhoa_max);
2414 free(derv_taua_min);
2415 free(derv_taua_max);
2429 int32_t nmodels, int32_t mindx[],
2430 geom_str *
geom,
float wv,
float rhoa[],
2431 int32_t *modmin, int32_t *modmax,
float *modrat,
float *epsnir,
2432 float tau_pred_max[],
float tau_pred_min[],
float rho_pred_max[],
float rho_pred_min[],
2433 float tau_aer[],
float rho_aer[]) {
2436 float *
ac, *
bc, *cc, *dc, *ec;
2437 float ext_iwnir_l, ext_iwnir_s;
2438 double ax, bx, cx, fx;
2439 float tau_iwnir_l[nmodels];
2440 float lg_tau_iwnir_s[nmodels], tau_iwnir_s[nmodels];
2441 float lg_rho_iwnir_s_pred[nmodels], rho_iwnir_s_pred[nmodels];
2442 float eps_pred[nmodels];
2446 float lg_tau_iwnir_l;
2447 int iwtab_l, iwtab_s;
2450 static float eps_obs;
2459 eps_obs = rhoa[iwnir_s] / rhoa[iwnir_l];
2469 for (
im = 0;
im < nmodels;
im++) {
2480 iwtab_l = iwatab[iwnir_l];
2481 iwtab_s = iwatab[iwnir_s];
2483 ax = (
double)
ac[iwtab_l]-(
double) rhoa[iwnir_l];
2485 cx = (
double) cc[iwtab_l];
2487 fx = bx * bx - 4.0 *
ax*cx;
2488 if (fx > 0.0 && cx != 0.0) {
2489 lg_tau_iwnir_l = 0.5 * (-bx + sqrt(fx)) / cx;
2490 tau_iwnir_l[
im] = (lg_tau_iwnir_l);
2498 ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
2499 ext_iwnir_s = aertab->model[modl]->extc[iwtab_s];
2501 tau_iwnir_s[
im] = (ext_iwnir_s / ext_iwnir_l) * tau_iwnir_l[
im];
2502 lg_tau_iwnir_s[
im] = (tau_iwnir_s[
im]);
2505 printf(
"\nErroneous aerosol option, %d\n", aer_opt);
2506 printf(
"You are using a log-space LUT. The aerosol LUT coefficients need to be in linear-space. Use aer_opt=-18 instead. \n");
2513 lg_rho_iwnir_s_pred[
im] =
ac[iwtab_s] +
bc[iwtab_s] * lg_tau_iwnir_s[
im] +
2514 cc[iwtab_s] * pow(lg_tau_iwnir_s[
im], 2) +
2515 dc[iwtab_s] * pow(lg_tau_iwnir_s[
im], 3) +
2516 ec[iwtab_s] * pow(lg_tau_iwnir_s[
im], 4);
2518 rho_iwnir_s_pred[
im] = (lg_rho_iwnir_s_pred[
im]);
2522 eps_pred[
im] = rho_iwnir_s_pred[
im] / rhoa[iwnir_l];
2551 *modmin = mindx[im1];
2552 *modmax = mindx[im2];
2562 for (iw = 0; iw < nwave; iw++) {
2563 tau_aer[iw] = (1.0 - mwt) * tau_pred_min[iw] + mwt * tau_pred_max[iw];
2564 rho_aer[iw] = (1.0 - mwt) * rho_pred_min[iw] + mwt * rho_pred_max[iw];
2583 static float nw = 1.334;
2586 static float lastsolz = -999.;
2587 static float lastsenz = -999.;
2588 static float lastphi = -999.;
2590 static float *fres1, *fres2;
2591 static float *scatt1, *scatt2;
2592 static float scatt1_cnst, scatt2_cnst, fres1_cnst, fres2_cnst;
2593 static float *scatt1_ar, *scatt2_ar, *fres1_ar, *fres2_ar;
2594 static int firstCall = 1, gmult = 1;
2596 float phase1, phase2;
2599 if (firstCall == 1) {
2602 if ((
phase[
im] = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
2603 printf(
"Unable to allocate space for phase.\n");
2608 if ((
geom->gmult == 0) || (interpol == 1)) {
2610 scatt1 = &scatt1_cnst;
2611 scatt2 = &scatt2_cnst;
2612 fres1 = &fres1_cnst;
2613 fres2 = &fres2_cnst;
2617 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2618 printf(
"Unable to allocate space for scatt1.\n");
2622 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2623 printf(
"Unable to allocate space for scatt2.\n");
2627 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2628 printf(
"Unable to allocate space for fres1.\n");
2632 (
float *) malloc(aertab->nwave * sizeof (
float))) ==
NULL) {
2633 printf(
"Unable to allocate space for fres2.\n");
2645 if ((
geom->solz[0] != lastsolz) || (
geom->senz[0] != lastsenz) ||
2646 (
geom->phi[0] != lastphi)) {
2650 for (
im = 0;
im < aertab->nmodel;
im++)
2654 for (iw = 0; iw < aertab->nwave; iw++) {
2656 temp = sqrt((1.0 -
geom->csenz[ig] *
geom->csenz[ig]) *
2657 (1.0 -
geom->csolz[ig] *
geom->csolz[ig])) *
2658 cos(
geom->phi[ig] / radeg);
2660 MAX(-
geom->csenz[ig] *
geom->csolz[ig] + temp, -1.0)) * radeg;
2662 MIN(
geom->csenz[ig] *
geom->csolz[ig] + temp, 1.0)) * radeg;
2667 if (gmult == 0)
break;
2670 lastsolz =
geom->solz[0];
2671 lastsenz =
geom->senz[0];
2672 lastphi =
geom->phi[0];
2681 for (iw = 0; iw < aertab->nwave; iw++) {
2684 &aertab->model[
im]->lnphase[iw][0],
2685 &aertab->model[
im]->d2phase[iw][0],
2686 aertab->nscatt, scatt1[ig], &phase1);
2688 &aertab->model[
im]->lnphase[iw][0],
2689 &aertab->model[
im]->d2phase[iw][0],
2690 aertab->nscatt, scatt2[ig], &phase2);
2693 exp(phase2)*(fres1[ig] + fres2[ig]);
2697 return (&
phase[modnum][0]);
2706 static int firstCall = 1;
2711 float rhoas1, rhoas2;
2716 iw1 =
windex(765, aertab->wave, aertab->nwave);
2717 iw2 =
windex(865, aertab->wave, aertab->nwave);
2718 if (iw1 == iw2) iw1--;
2723 rhoas1 = aertab->model[modnum]->albedo[iw1] *
phase[iw1] * aertab->model[modnum]->extc[iw1];
2724 rhoas2 = aertab->model[modnum]->albedo[iw2] *
phase[iw2] * aertab->model[modnum]->extc[iw2];
2725 eps = rhoas1 / rhoas2;
2726 cf = log(
eps) / (aertab->wave[iw2] - aertab->wave[iw1]);
2750 static int firstCall = 1;
2769 printf(
"\nLoading water-vapor correction coefficients.\n");
2772 f = (a01[iw] + a02[iw] * airmass +
cf * (a03[iw] + a04[iw] * airmass))
2773 + (a05[iw] + a06[iw] * airmass +
cf * (a07[iw] + a08[iw] * airmass)) * wv
2774 + (a09[iw] + a10[iw] * airmass +
cf * (a11[iw] + a12[iw] * airmass)) * wv*wv;
2787 float rhoa[],
float wave[], int32_t nwave,
int iw1,
int iw2,
float rhoas[]) {
2788 float *
ac, *
bc, *cc;
2800 for (iw = iw1; iw <= iw2; iw++) {
2801 ig = iw *
geom->gmult;
2802 if (rhoa[iw] < 1.e-20)
2803 rhoas[iw] = rhoa[iw];
2809 f =
b *
b - 4 *
c * (
a - log((
double) rhoa[iw]));
2811 if (
fabs(
c) > 1.e-20) {
2812 rhoas[iw] = exp(0.5 * (-
b + sqrt(
f)) /
c);
2813 }
else if (
fabs(
a) > 1.e-20 &&
fabs(
b) > 1.e-20) {
2814 rhoas[iw] = pow(rhoa[iw] /
a, 1. /
b);
2820 if (!isfinite(rhoas[iw]) || rhoas[iw] < 1.e-20) {
2834 for (iw = iw1; iw <= iw2; iw++)
2835 rhoas[iw] = rhoa[iw];
2850 float rhoas[],
float wave[], int32_t nwave,
int iw1,
int iw2,
float rhoa[]) {
2851 float *
ac, *
bc, *cc;
2862 for (iw = iw1; iw <= iw2; iw++) {
2863 ig = iw *
geom->gmult;
2868 if (rhoas[iw] < 1.e-20)
2869 rhoa[iw] = rhoas[iw];
2876 rhoa[iw] = exp(
a +
b * lnrhoas +
c * lnrhoas * lnrhoas);
2910 static float lastsolz = -999.;
2911 static float lastsenz = -999.;
2912 static float lastphi = -999.;
2913 static int32_t lastiwl = -999;
2914 static int lastmod = -999;
2916 static int firstCall = 1;
2921 if (firstCall == 1) {
2923 maxwave =
MAX(aertab->nwave, nwave);
2925 if ((
epsilon[
i] = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
2926 printf(
"Unable to allocate space for epsilon.\n");
2935 if (modnum != lastmod ||
geom->solz[0] != lastsolz ||
2936 geom->senz[0] != lastsenz ||
geom->phi[0] != lastphi ||
2937 iwnir_l != lastiwl) {
2939 int iwnir = iwatab[iwnir_l];
2942 float rhoas1, rhoas2;
2945 if ((lneps = (
float *) calloc(aertab->nwave, sizeof (
float))) ==
NULL) {
2946 printf(
"Unable to allocate space for lneps.\n");
2952 for (iw = 0; iw < aertab->nwave; iw++) {
2953 rhoas1 = aertab->model[
im]->albedo[iw] *
phase[iw] * aertab->model[
im]->extc[iw];
2954 rhoas2 = aertab->model[
im]->albedo[iwnir] *
phase[iwnir] * aertab->model[
im]->extc[iwnir];
2960 for (iw = 0; iw < nwave; iw++) {
2962 if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
2963 epsilon[
im][iw] = exp(
linterp(aertab->wave, lneps, aertab->nwave, wave[iw]));
2969 lastsolz =
geom->solz[0];
2970 lastsenz =
geom->senz[0];
2971 lastphi =
geom->phi[0];
2989 int32_t nmodel, int32_t mindx[], geom_str *
geom,
float wv,
2990 float rhoa[], int32_t iwnir_s, int32_t iwnir_l, int32_t *modmin,
2991 int32_t *modmax,
float *modrat,
float *epsnir) {
3013 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3014 printf(
"Unable to allocate space for rhoas.\n");
3022 for (jm = 0; jm < nmod; jm++) {
3029 if (rhoas[iwnir_l] > 0.0000001)
3030 eps_ret[jm] = rhoas[iwnir_s] / rhoas[iwnir_l];
3036 eps_mod[jm] =
eps[iwnir_s];
3038 eps_ave += eps_ret[jm];
3040 if (isfinite(eps_ave))
3053 for (
im = 0;
im < nmodel;
im++) {
3055 eps_err[
im] = eps_ave - eps_mod[
im];
3059 for (
im = 0;
im < nmodel - 1;
im++) {
3060 for (iim =
im + 1; iim < nmodel; iim++)
3061 if (
fabs(eps_err[imod[
im]]) >
fabs(eps_err[imod[iim]])) {
3063 imod[
im ] = imod[iim];
3075 for (iim = 0; iim < nmod; iim++) {
3077 tot_err +=
fabs(eps_err[
im]);
3081 for (iim = 0; iim < nmod; iim++) {
3083 wt = 1.0 -
fabs(eps_err[
im]) / tot_err;
3084 eps_ave += eps_ret[
im] * wt;
3086 eps_ave /= (nmod - 1);
3090 err_m = 0 - FLT_MAX;
3092 for (
im = 0;
im < nmodel;
im++) {
3093 eps_err[
im] = eps_ave - eps_mod[
im];
3094 if (eps_err[
im] >= 0.0) {
3095 if (eps_err[
im] < err_p) {
3096 err_p = eps_err[
im];
3100 if (eps_err[
im] > err_m) {
3101 err_m = eps_err[
im];
3113 }
else if (*modmax < 0) {
3119 *modrat = (eps_ave - eps_mod[*modmin]) / (eps_mod[*modmax] - eps_mod[*modmin]);
3121 *modmin = mindx[*modmin];
3122 *modmax = mindx[*modmax];
3139 return (
x->angstrom <
y->angstrom ? -1 : 1);
3144 static int firstCall = 1;
3150 int ib =
windex(520, aertab->wave, aertab->nwave);
3153 for (
im = 0;
im < aertab->nmodel;
im++) {
3154 alphaT[
im].angstrom = aertab->model[
im]->angstrom[ib];
3155 alphaT[
im].modnum =
im;
3157 qsort(alphaT, aertab->nmodel, sizeof (alphaTstr),
3158 (
int (*)(
const void *,
const void *))
compalphaT);
3163 for (
im = 0;
im < aertab->nmodel;
im++) {
3164 if (angstrom < alphaT[
im].angstrom)
3167 im1 =
MAX(
MIN(
im - 1, aertab->nmodel - 1), 0);
3168 im2 =
MAX(
MIN(
im, aertab->nmodel - 1), 0);
3170 *modmin = alphaT[im1].modnum;
3171 *modmax = alphaT[im2].modnum;
3176 *modrat = (angstrom - alphaT[im1].angstrom) /
3177 (alphaT[im2].angstrom - alphaT[im1].angstrom);
3195 int32_t iwnir_l,
float rhoa[], geom_str *
geom,
float wv,
float taua[]) {
3200 int iwnir = iwatab[iwnir_l];
3205 iwg = iwnir *
geom->gmult;
3207 maxwave =
MAX(aertab->nwave, nwave);
3209 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3210 printf(
"Unable to allocate space for rhoas.\n");
3213 if ((aot = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
3214 printf(
"Unable to allocate space for aot.\n");
3217 if ((lnaot = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
3218 printf(
"Unable to allocate space for lnaot.\n");
3226 aot[iwnir] = rhoas[iwnir_l]*(4.0 *
geom->csolz[iwg] *
geom->csenz[iwg])
3227 / (
phase[iwnir] * aertab->model[modnum]->albedo[iwnir]);
3230 for (iw = 0; iw < aertab->nwave; iw++) {
3232 aot[iw] = aot[iwnir] * aertab->model[modnum]->extc[iw] / aertab->model[modnum]->extc[iwnir];
3234 lnaot[iw] = log(aot[iw]);
3239 for (iw = 0; iw < nwave; iw++) {
3241 if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
3242 taua[iw] = exp(
linterp(aertab->wave, lnaot, aertab->nwave, wave[iw]));
3244 taua[iw] = aot[iwtab];
3247 for (iw = 0; iw < nwave; iw++)
3261 float rhoa[], geom_str *
geom,
float tau_pred[])
3263 static int firstcall=1;
3264 static float **pc_rhoa;
3265 static int npc, ntau_870;
3266 aermodstr *aermod=aertab->model[modl];
3269 static float *rhoa_l;
3273 ntau_870=aertab->ntau_870;
3275 pc_rhoa=(
float **)malloc(aertab->ntau_870*
sizeof(
float *));
3277 for(itau=0;itau<ntau_870;itau++)
3278 pc_rhoa[itau]=(
float *)malloc(aertab->npc*
sizeof(
float));
3279 rhoa_l=(
float *)malloc(ntau_870*
sizeof(
float));
3284 for(itau=0;itau<ntau_870;itau++){
3286 for(ipc=0;ipc<npc;ipc++){
3287 rhoa_l[itau]+=pc_rhoa[itau][ipc]*aermod->pc_components_rhoa[ipc][iwnir_l];
3288 rhoa_l[itau]+=aermod->pc_mean_rhoa[iwnir_l];
3292 for(itau=0;itau<ntau_870;itau++){
3293 if(rhoa[iwnir_l]<rhoa_l[itau] && rhoa[iwnir_l]>rhoa_l[itau-1])
3296 tau_pred[iwnir_l]=aermod->tau_870[itau]+(aermod->tau_870[itau]-aermod->tau_870[itau-1])*(rhoa[iwnir_l]-rhoa_l[itau])/(rhoa_l[itau]-rhoa_l[itau-1]);
3298 for(iw=0;iw<nwave;iw++)
3299 tau_pred[iw]=tau_pred[iwnir_l]*aermod->extc[iw]/aermod->extc[iwnir_l];
3310 float rhoa[], geom_str *
geom,
float tau_pred[])
3312 float *
ac, *
bc, *cc, *dc, *ec;
3313 double ax, bx, cx, fx;
3314 float ext_modl[nwave];
3317 int iwtab_l = iwatab[iwnir_l];
3323 ax = (
double)
ac[iwtab_l] - log((
double) rhoa[iwnir_l]);
3325 cx = (
double) cc[iwtab_l];
3327 fx = bx * bx - 4.0 *
ax*cx;
3328 if (fx > 0.0 && cx != 0.0) {
3329 tau_iwnir_l = exp(0.5*(-bx + sqrt(fx)) / cx);
3335 for (iw = 0; iw < nwave; iw++) {
3337 ext_modl[iw] = aertab->model[modl]->extc[iwtab];
3339 for (iw = 0; iw < nwave; iw++) {
3340 tau_pred[iw] = (ext_modl[iw] / ext_modl[iwnir_l]) * tau_iwnir_l;
3354 float *theta,
int gmult,
float taua[],
float dtran[],
float dt[]) {
3395 for (iw = 0; iw < nwave; iw++) {
3396 if ((iw == 0) || (gmult != 0)) {
3399 for (
i = 0;
i < aertab->dtran_ntheta;
i++) {
3400 if (theta[ig] < aertab->dtran_theta[
i])
3403 if (
i == aertab->dtran_ntheta) {
3408 i1 =
MIN(
MAX(
i - 1, 0), aertab->dtran_ntheta - 2);
3410 x1 = aertab->dtran_airmass[i1];
3411 x2 = aertab->dtran_airmass[i2];
3412 xbar = 1.0 / cos(theta[ig] / radeg);
3413 wt = (xbar - x1) / (x2 - x1);
3418 a1 = aertab->model[modnum]->dtran_a[iwtab][i1];
3419 b1 = aertab->model[modnum]->dtran_b[iwtab][i1];
3421 a2 = aertab->model[modnum]->dtran_a[iwtab][i2];
3422 b2 = aertab->model[modnum]->dtran_b[iwtab][i2];
3431 y1 = a1 * exp(-
b1 * taua[iw]);
3432 y2 = a2 * exp(-b2 * taua[iw]);
3434 dtran[iw] =
MAX(
MIN((1.0 - wt) * y1 + wt*y2, 1.0), 1e-5);
3437 if(
fabs(dtran[iw]-1.0)<0.0000001 ||
fabs(dtran[iw]-1e-5)<0.0000001 )
3440 dt[iw]=-(1.0-wt)*y1*
b1-wt*y2*b2;
3447 float *theta,
int gmult,
float taua[],
float dtran[],
float dt[]) {
3449 static int firstcall=1;
3450 static float *pc_td, *
solz;
3451 static float *deriv_pc;
3452 static int npc, ntau_870,nsolz,nir_base;
3453 aermodstr *aermod=aertab->model[modnum];
3454 float a00,a01,a10,a11;
3456 int iw,itau,itheta,ipc,itau1,itau2,itheta1,itheta2;
3457 float r,
p,deriv_r=0.;
3463 nsolz=aertab->nsenz;
3464 ntau_870=aertab->ntau_870;
3466 pc_td=(
float *)malloc(aertab->npc*
sizeof(
float));
3468 deriv_pc=(
float *)malloc(npc*
sizeof(
float));
3473 if(taua[nir_base]<aermod->tau_870[0]){
3476 }
else if(taua[nir_base]>=aermod->tau_870[ntau_870-1]){
3480 for(itau=0;itau<ntau_870;itau++)
3481 if(taua[nir_base]<aermod->tau_870[itau])
3488 deriv_r=1/(aermod->tau_870[itau2]-aermod->tau_870[itau1]);
3489 r=(taua[nir_base]-aermod->tau_870[itau1])/(aermod->tau_870[itau2]-aermod->tau_870[itau1]);
3494 for(itheta=0;itheta<nsolz;itheta++)
3495 if(*theta<
solz[itheta])
3497 itheta1=
MAX(itheta-1,0);
3498 itheta2=
MIN(itheta,nsolz-1);
3499 if(itheta1!=itheta2)
3504 for(ipc=0;ipc<npc;ipc++){
3505 a00=aermod->pc_td[itau1][itheta1][ipc];
3506 a01=aermod->pc_td[itau1][itheta2][ipc];
3507 a10=aermod->pc_td[itau2][itheta1][ipc];
3508 a11=aermod->pc_td[itau2][itheta2][ipc];
3510 pc_td[ipc]=a00*(1-
r)*(1-
p)+a01*(1-
r)*
p+a10*
r*(1-
p)+a11*
r*
p;
3512 deriv_pc[ipc]=(-a00*(1-
p)-a01*
p+a10*(1-
p)+a11*
p)*deriv_r;
3515 for(iw=0;iw<nwave;iw++){
3519 for(ipc=0;ipc<npc;ipc++){
3520 dtran[iw]+=pc_td[ipc]*aermod->pc_components_td[ipc][iw];
3522 dt[iw]+=aermod->pc_components_td[ipc][iw]*deriv_pc[ipc];
3524 dtran[iw]+=aermod->pc_mean_td[iw];
3525 dtran[iw]=exp(dtran[iw]);
3532 geom_str *
geom,
float wv,
float pr,
float taur[], int32_t modmin,
3533 int32_t modmax,
float modrat,
float rhoa[],
float taua[],
float tsol[],
3534 float tsen[],
float tauamin[],
float tauamax[],
int taua_opt,
int ip,uncertainty_t *uncertainty) {
3536 int iw, gmult, ig,inir,ib;
3543 static int firstcall=1;
3544 static int32_t nbands_ac, *acbands_index;
3549 nbands_ac=
input->nbands_ac;
3550 acbands_index=
input->acbands_index;
3554 if ((tsolmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3555 printf(
"Unable to allocate space for tsolmin.\n");
3558 if ((tsolmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3559 printf(
"Unable to allocate space for tsolmax.\n");
3562 if ((tsenmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3563 printf(
"Unable to allocate space for tsenmin.\n");
3566 if ((tsenmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3567 printf(
"Unable to allocate space for tsenmax.\n");
3574 if ((dtmin = (
float *) calloc(nwave,
sizeof(
float))) ==
NULL) {
3575 printf(
"Unable to allocate space for dtmin.\n");
3578 if ((dtmax = (
float *) calloc(nwave,
sizeof(
float))) ==
NULL) {
3579 printf(
"Unable to allocate space for dtmax.\n");
3585 gmult =
geom->gmult;
3586 if (interpol == 1) gmult = 0;
3589 if (taua_opt == 0) {
3594 }
else if (taua_opt == 2) {
3605 for (iw = 0; iw < nwave; iw++) {
3606 ig = iw *
geom->gmult;
3607 tsol[iw] = tsolmin[iw]*(1.0 - modrat) + tsolmax[iw] * modrat;
3611 for(ib=0;ib<nbands_ac;ib++){
3612 inir=acbands_index[ib]-iwnir_s;
3613 uncertainty->derv_tsol_rhorc[iw][inir]=(tsolmax[iw]-tsolmin[iw])*uncertainty->derv_modrat_rhorc[inir];
3614 if(acbands_index[ib]==iwnir_l){
3615 uncertainty->derv_tsol_rhorc[iw][inir]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhorc_l[iwnir_l];
3616 uncertainty->derv_tsol_rhorc[iw][inir]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_rhorc_l[iwnir_l];
3619 uncertainty->derv_tsol_taua_l[iw]=(tsolmax[iw]-tsolmin[iw])*uncertainty->derv_modrat_taua_l;
3620 uncertainty->derv_tsol_taua_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_taua_l[iw];
3621 uncertainty->derv_tsol_taua_l[iw]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_taua_l[iw];
3623 uncertainty->derv_tsol_rhow_l[iw]=(tsolmax[iw]-tsolmin[iw])*uncertainty->derv_modrat_rhow_l;
3624 uncertainty->derv_tsol_rhow_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhow_l[iw];
3625 uncertainty->derv_tsol_rhow_l[iw]+=modrat *dtmax[iw]*uncertainty->derv_taua_max_rhow_l[iw];
3629 double tmp_pressure_diff = exp(-0.5 * taur[iw] /
geom->csolz[ig]*(
pr / p0 - 1));
3630 tsol[iw] = tsol[iw] * tmp_pressure_diff;
3635 for(ib=0;ib<nbands_ac;ib++){
3636 inir=acbands_index[ib]-iwnir_s;
3637 uncertainty->derv_tsol_rhorc[iw][inir]*=tmp_pressure_diff;
3639 uncertainty->derv_tsol_taua_l[iw]*=tmp_pressure_diff;
3640 uncertainty->derv_tsol_rhow_l[iw]*=tmp_pressure_diff;
3646 tsol[iw] = pow(tsol[iw],
geom->airmass_sph[ig] /
geom->airmass_plp[ig]);
3649 uncertainty->dt_sol[ip*nwave+iw]*= (
geom->airmass_sph[ig] /
geom->airmass_plp[ig] *pow(
tmp,
geom->airmass_sph[ig] /
geom->airmass_plp[ig]-1) );
3658 for (iw = 0; iw < nwave; iw++) {
3659 ig = iw *
geom->gmult;
3660 taua[iw] = tauamin[iw]*(1.0 - modrat) + tauamax[iw] * modrat;
3661 tsen[iw] = tsenmin[iw]*(1.0 - modrat) + tsenmax[iw] * modrat;
3664 for(ib=0;ib<nbands_ac;ib++){
3665 inir=acbands_index[ib]-iwnir_s;
3666 uncertainty->derv_tsen_rhorc[iw][inir]=(tsenmax[iw]-tsenmin[iw])*uncertainty->derv_modrat_rhorc[inir];
3667 if(acbands_index[ib]==iwnir_l){
3668 uncertainty->derv_tsen_rhorc[iw][inir]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhorc_l[iwnir_l];
3669 uncertainty->derv_tsen_rhorc[iw][inir]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_rhorc_l[iwnir_l];
3672 uncertainty->derv_tsen_taua_l[iw]=(tsenmax[iw]-tsenmin[iw])*uncertainty->derv_modrat_taua_l;
3673 uncertainty->derv_tsen_taua_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_taua_l[iw];
3674 uncertainty->derv_tsen_taua_l[iw]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_taua_l[iw];
3676 uncertainty->derv_tsen_rhow_l[iw]=(tsenmax[iw]-tsenmin[iw])*uncertainty->derv_modrat_rhow_l;
3677 uncertainty->derv_tsen_rhow_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhow_l[iw];
3678 uncertainty->derv_tsen_rhow_l[iw]+=modrat *dtmax[iw]*uncertainty->derv_taua_max_rhow_l[iw];
3681 double tmp_pressure_diff = exp(-0.5 * taur[iw] /
geom->csenz[ig] *(
pr / p0 - 1));
3682 tsen[iw] = tsen[iw] * tmp_pressure_diff;
3685 for(ib=0;ib<nbands_ac;ib++){
3686 inir=acbands_index[ib]-iwnir_s;
3687 uncertainty->derv_tsen_rhorc[iw][inir]*=tmp_pressure_diff;
3689 uncertainty->derv_tsen_taua_l[iw]*=tmp_pressure_diff;
3690 uncertainty->derv_tsen_rhow_l[iw]*=tmp_pressure_diff;
3696 tsen[iw] = pow(tsen[iw],
geom->airmass_sph[ig] /
geom->airmass_plp[ig]);
3698 uncertainty->dt_sen[ip*nwave+iw]*= (
geom->airmass_sph[ig] /
geom->airmass_plp[ig] *pow(
tmp,
geom->airmass_sph[ig] /
geom->airmass_plp[ig]-1) );
3718 geom_str *
geom,
float wv,
float pr,
float taur[], int32_t modmin,
3719 int32_t modmax,
float modrat,
float rhoa[],
float taua[],
float tsol[],
3720 float tsen[],
float tauamin[],
float tauamax[],
int taua_opt,
int ip,uncertainty_t *uncertainty) {
3734 diff_tran_pca(
sensorID,wave,nwave,iwnir_l,
geom,wv,
pr,taur,modmin,modmax,modrat,rhoa,taua,tsol,tsen,\
3735 tauamin,tauamax,taua_opt,ip,uncertainty);
3739 if ((tsolmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3740 printf(
"Unable to allocate space for tsolmin.\n");
3743 if ((tsolmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3744 printf(
"Unable to allocate space for tsolmax.\n");
3747 if ((tsenmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3748 printf(
"Unable to allocate space for tsenmin.\n");
3751 if ((tsenmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
3752 printf(
"Unable to allocate space for tsenmax.\n");
3759 if ((dtmin = (
float *) calloc(nwave,
sizeof(
float))) ==
NULL) {
3760 printf(
"Unable to allocate space for dtmin.\n");
3763 if ((dtmax = (
float *) calloc(nwave,
sizeof(
float))) ==
NULL) {
3764 printf(
"Unable to allocate space for dtmax.\n");
3770 gmult =
geom->gmult;
3771 if (interpol == 1) gmult = 0;
3774 if (taua_opt == 0) {
3779 }
else if (taua_opt == 2) {
3790 for (iw = 0; iw < nwave; iw++) {
3791 ig = iw *
geom->gmult;
3792 tsol[iw] = tsolmin[iw]*(1.0 - modrat) + tsolmax[iw] * modrat;
3796 for(inir=iwnir_s;inir<=iwnir_l;inir++){
3797 uncertainty->derv_tsol_rhorc[iw][inir-iwnir_s]=(tsolmax[iw]-tsolmin[iw])*uncertainty->derv_modrat_rhorc[inir-iwnir_s];
3799 uncertainty->derv_tsol_rhorc[iw][inir-iwnir_s]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhorc_l[iw];
3800 uncertainty->derv_tsol_rhorc[iw][inir-iwnir_s]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_rhorc_l[iw];
3803 uncertainty->derv_tsol_taua_l[iw]=(tsolmax[iw]-tsolmin[iw])*uncertainty->derv_modrat_taua_l;
3804 uncertainty->derv_tsol_taua_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_taua_l[iw];
3805 uncertainty->derv_tsol_taua_l[iw]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_taua_l[iw];
3807 uncertainty->derv_tsol_rhow_l[iw]=(tsolmax[iw]-tsolmin[iw])*uncertainty->derv_modrat_rhow_l;
3808 uncertainty->derv_tsol_rhow_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhow_l[iw];
3809 uncertainty->derv_tsol_rhow_l[iw]+=modrat *dtmax[iw]*uncertainty->derv_taua_max_rhow_l[iw];
3813 double tmp_pressure_diff = exp(-0.5 * taur[iw] /
geom->csolz[ig]*(
pr / p0 - 1));
3814 tsol[iw] = tsol[iw] * tmp_pressure_diff;
3819 for(inir=iwnir_s;inir<=iwnir_l;inir++)
3820 uncertainty->derv_tsol_rhorc[iw][inir-iwnir_s]*=tmp_pressure_diff;
3821 uncertainty->derv_tsol_taua_l[iw]*=tmp_pressure_diff;
3822 uncertainty->derv_tsol_rhow_l[iw]*=tmp_pressure_diff;
3828 tsol[iw] = pow(tsol[iw],
geom->airmass_sph[ig] /
geom->airmass_plp[ig]);
3831 uncertainty->dt_sol[ip*nwave+iw]*= (
geom->airmass_sph[ig] /
geom->airmass_plp[ig] *pow(
tmp,
geom->airmass_sph[ig] /
geom->airmass_plp[ig]-1) );
3840 for (iw = 0; iw < nwave; iw++) {
3841 ig = iw *
geom->gmult;
3842 taua[iw] = tauamin[iw]*(1.0 - modrat) + tauamax[iw] * modrat;
3843 tsen[iw] = tsenmin[iw]*(1.0 - modrat) + tsenmax[iw] * modrat;
3846 for(inir=iwnir_s;inir<=iwnir_l;inir++){
3847 uncertainty->derv_tsen_rhorc[iw][inir-iwnir_s]=(tsenmax[iw]-tsenmin[iw])*uncertainty->derv_modrat_rhorc[inir-iwnir_s];
3849 uncertainty->derv_tsen_rhorc[iw][inir-iwnir_s]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhorc_l[iw];
3850 uncertainty->derv_tsen_rhorc[iw][inir-iwnir_s]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_rhorc_l[iw];
3853 uncertainty->derv_tsen_taua_l[iw]=(tsenmax[iw]-tsenmin[iw])*uncertainty->derv_modrat_taua_l;
3854 uncertainty->derv_tsen_taua_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_taua_l[iw];
3855 uncertainty->derv_tsen_taua_l[iw]+= modrat *dtmax[iw]*uncertainty->derv_taua_max_taua_l[iw];
3857 uncertainty->derv_tsen_rhow_l[iw]=(tsenmax[iw]-tsenmin[iw])*uncertainty->derv_modrat_rhow_l;
3858 uncertainty->derv_tsen_rhow_l[iw]+=(1.0 - modrat)*dtmin[iw]*uncertainty->derv_taua_min_rhow_l[iw];
3859 uncertainty->derv_tsen_rhow_l[iw]+=modrat *dtmax[iw]*uncertainty->derv_taua_max_rhow_l[iw];
3862 double tmp_pressure_diff = exp(-0.5 * taur[iw] /
geom->csenz[ig] *(
pr / p0 - 1));
3863 tsen[iw] = tsen[iw] * tmp_pressure_diff;
3866 for(inir=iwnir_s;inir<=iwnir_l;inir++)
3867 uncertainty->derv_tsen_rhorc[iw][inir-iwnir_s]*=tmp_pressure_diff;
3868 uncertainty->derv_tsen_taua_l[iw]*=tmp_pressure_diff;
3869 uncertainty->derv_tsen_rhow_l[iw]*=tmp_pressure_diff;
3875 tsen[iw] = pow(tsen[iw],
geom->airmass_sph[ig] /
geom->airmass_plp[ig]);
3877 uncertainty->dt_sen[ip*nwave+iw]*= (
geom->airmass_sph[ig] /
geom->airmass_plp[ig] *pow(
tmp,
geom->airmass_sph[ig] /
geom->airmass_plp[ig]-1) );
3901 int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *
geom,
float wv,
float rhoa[], int32_t *modmin, int32_t *modmax,
3902 float32 *modrat,
float *tau_pred_min,
float *tau_pred_max, int32_t ip,
float chi[],
float *mbac_w,uncertainty_t *uncertainty) {
3903 static float *tau_iwnir_l;
3904 static int firstcall=1, nir_l;
3905 static float **pc_rhoa;
3906 static int npc, ntau_870;
3908 static int wave_base_index;
3910 int iw, ib, itau,ipc,itemp;
3914 static float *rho_all_wav_pred,*rho_pred_min,*rho_pred_max;
3918 float diff_2 = 0.0,diff_1=0.0;
3919 float mbac_wsum = 0.0;
3920 static float *chi_old;
3924 static float **derv_chi_rhorc;
3925 static float *derv_chi_taua_l;
3926 static float **derv_rhoa_rhorc_l;
3927 float *derv_modrat_rhorc;
3928 float derv_modrat_taua_l=0.;
3929 float derv_modrat_rhow_l=0.;
3930 float *derv_taua_min_rhorc_l;
3931 float *derv_taua_min_taua_l;
3932 float *derv_taua_min_rhow_l;
3933 float *derv_taua_max_rhorc_l;
3934 float *derv_taua_max_taua_l;
3935 float *derv_taua_max_rhow_l;
3936 static float *derv_taua_rhorc_l;
3937 static float *derv_chi_obs_rhorc;
3938 float derv_chi_obs_taua_l;
3940 static float *derv_rhoa_min=
NULL,*derv_rhoa_max=
NULL;
3941 static float *derv_taua_min=
NULL,*derv_taua_max=
NULL;
3944 float derv_temp_rhorc;
3945 float derv_modrat_chi0,derv_modrat_chi1;
3947 static float **rhoa_modl;
3948 struct str chi_struct[nmodels];
3950 static int32_t nbands_ac, *acbands_index;
3957 printf(
"aer_wave_base is different from the aer_wave_long, which may break the uncertainty for mbac\n");
3961 ntau_870=aertab->ntau_870;
3963 rhoa_modl=(
float **)malloc(aertab->ntau_870*
sizeof(
float*));
3964 pc_rhoa=(
float **)malloc(aertab->ntau_870*
sizeof(
float *));
3966 for(itau=0;itau<aertab->ntau_870;itau++){
3967 pc_rhoa[itau]=(
float *)malloc(aertab->npc*
sizeof(
float));
3968 rhoa_modl[itau]=(
float *)malloc(nwave*
sizeof(
float));
3970 rho_all_wav_pred = (
float*) malloc(nwave *
sizeof (
float));
3971 rho_pred_min = (
float*) malloc(nwave *
sizeof (
float));
3972 rho_pred_max = (
float*) malloc(nwave *
sizeof (
float));
3973 chi_old=(
float*) malloc((nmodels) *
sizeof (
float));
3974 tau_iwnir_l=(
float*) malloc((nmodels) *
sizeof (
float));
3976 nbands_ac=
input->nbands_ac;
3977 acbands_index=
input->acbands_index;
3979 float *tempwave=(
float *)malloc(nbands_ac*
sizeof(
float));
3980 for(iw=0;iw<nbands_ac;iw++)
3981 tempwave[iw]=
input->mbac_wave[iw];
3983 wave_base_index=
windex(
input->aer_wave_base*1.,tempwave,nbands_ac);
3987 derv_chi_taua_l =(
float *)malloc(nmodels*
sizeof(
float));
3988 derv_chi_rhorc =(
float **)malloc(nmodels*
sizeof(
float *));
3989 derv_rhoa_rhorc_l=(
float **)malloc(nmodels*
sizeof(
float*));
3990 derv_taua_rhorc_l=(
float *)malloc(nmodels*
sizeof(
float));
3992 for(
im=0;
im<nmodels;
im++){
3993 derv_chi_rhorc[
im] =(
float *)malloc(nbands_ac*
sizeof(
float));
3994 derv_rhoa_rhorc_l[
im]=(
float *)malloc(nwave*
sizeof(
float));
3996 derv_chi_obs_rhorc=(
float *)malloc(nbands_ac*
sizeof(
float));
3997 derv_rhoa_min=(
float *)malloc(nwave*
sizeof(
float));
3998 derv_rhoa_max=(
float *)malloc(nwave*
sizeof(
float));
3999 derv_taua_min=(
float *)malloc(nwave*
sizeof(
float));
4000 derv_taua_max=(
float *)malloc(nwave*
sizeof(
float));
4005 derv_modrat_rhorc=uncertainty->derv_modrat_rhorc;
4006 derv_taua_min_rhorc_l=uncertainty->derv_taua_min_rhorc_l;
4007 derv_taua_min_taua_l=uncertainty->derv_taua_min_taua_l;
4008 derv_taua_min_rhow_l=uncertainty->derv_taua_min_rhow_l;
4009 derv_taua_max_rhorc_l=uncertainty->derv_taua_max_rhorc_l;
4010 derv_taua_max_taua_l=uncertainty->derv_taua_max_taua_l;
4011 derv_taua_max_rhow_l=uncertainty->derv_taua_max_rhow_l;
4016 for(iw = 0; iw <nbands_ac; iw++)
4017 derv_chi_obs_rhorc[iw]=0.;
4018 derv_chi_obs_taua_l=0.;
4021 for (ib = 0; ib <nbands_ac; ib++) {
4023 iw=acbands_index[ib];
4026 chi_obs += rhoa[iw]/rhoa[nir_l];
4027 mbac_wsum += mbac_w[iw];
4029 derv_chi_obs_rhorc[ib]=1/rhoa[nir_l];
4030 derv_chi_obs_rhorc[wave_base_index]+=(-rhoa[iw]/rhoa[nir_l]/rhoa[nir_l]);
4031 derv_chi_obs_taua_l+=(-1/rhoa[nir_l]*uncertainty->derv_Lg_taua[iw]+rhoa[iw]/rhoa[nir_l]/rhoa[nir_l]*uncertainty->derv_Lg_taua[nir_l]);
4036 derv_chi_obs_taua_l/=mbac_wsum;
4037 for (ib = 0; ib <nbands_ac; ib++)
4038 derv_chi_obs_rhorc[ib]/=mbac_wsum;
4041 for (
im = 0;
im < nmodels;
im++) {
4044 aermod=aertab->model[modl];
4049 for(itau=0;itau<ntau_870;itau++)
4050 for (ib = 0; ib <nbands_ac; ib++){
4052 iw=acbands_index[ib];
4053 rhoa_modl[itau][iw]=0;
4055 for(ipc=0;ipc<npc;ipc++)
4056 rhoa_modl[itau][iw]+=pc_rhoa[itau][ipc]*aermod->pc_components_rhoa[ipc][iw];
4058 rhoa_modl[itau][iw]+=aermod->pc_mean_rhoa[iw];
4064 if(rhoa[nir_l]<=rhoa_modl[0][nir_l])
4066 else if (rhoa[nir_l]>= rhoa_modl[ntau_870-1][nir_l])
4069 for(itau=1;itau<ntau_870;itau++){
4070 if(rhoa[nir_l]<rhoa_modl[itau][nir_l] && rhoa[nir_l]>=rhoa_modl[itau-1][nir_l])
4074 tau_iwnir_l[
im]=aermod->tau_870[itau]+(aermod->tau_870[itau]-aermod->tau_870[itau-1])*(rhoa[nir_l]-rhoa_modl[itau][nir_l])/(rhoa_modl[itau][nir_l]-rhoa_modl[itau-1][nir_l]);
4077 if(
comp_rhoa_pca(nwave, wave,
geom, tau_iwnir_l[
im], modl, tau_pred_min, rhoa,derv_rhoa_min,derv_taua_min) !=0)
4079 for(ib=0;ib<nwave;ib++)
4080 tau_pred_max[ib]=tau_pred_min[ib];
4081 *modmin=modl;*modmax=modl;*modrat=1.0;
4086 derv_taua_rhorc_l[
im]=(aermod->tau_870[itau]-aermod->tau_870[itau-1])/(rhoa_modl[itau][nir_l]-rhoa_modl[itau-1][nir_l]);
4089 for (ib = 0; ib <nbands_ac; ib++){
4090 iw=acbands_index[ib];
4091 rho_all_wav_pred[iw]=rhoa_modl[itau][iw]+(rhoa_modl[itau][iw]-rhoa_modl[itau-1][iw])*(tau_iwnir_l[
im]-aermod->tau_870[itau])/(aermod->tau_870[itau]-aermod->tau_870[itau-1]);
4095 derv_rhoa_rhorc_l[
im][iw]=derv_taua_rhorc_l[
im]*(rhoa_modl[itau][iw]-rhoa_modl[itau-1][iw])/(aermod->tau_870[itau]-aermod->tau_870[itau-1]);
4100 derv_chi_rhorc[
im][wave_base_index]=0.;
4101 derv_chi_taua_l[
im]=0.;
4104 for (ib = 0; ib <nbands_ac; ib++){
4105 iw=acbands_index[ib];
4109 diff_2 += rho_all_wav_pred[iw]/rhoa[nir_l];
4110 diff_1 +=pow((rhoa[iw] - rho_all_wav_pred[iw])/noise[iw], 2)*mbac_w[iw];
4111 mbac_wsum += mbac_w[iw];
4114 derv_temp_rhorc=(derv_rhoa_rhorc_l[
im][iw]/rhoa[nir_l]-rho_all_wav_pred[iw]/rhoa[nir_l]/rhoa[nir_l]);
4115 derv_chi_rhorc [
im][wave_base_index]+=derv_temp_rhorc;
4116 derv_chi_taua_l[
im]+=(-derv_temp_rhorc*uncertainty->derv_Lg_taua[nir_l]);
4120 chi[
im] = diff_2 / mbac_wsum;
4121 chi_old[
im]=diff_1/mbac_wsum;
4123 derv_chi_rhorc[
im][wave_base_index]/=mbac_wsum;
4124 derv_chi_taua_l [
im]/=mbac_wsum;
4134 qsort(chi_struct,nmodels,
sizeof(chi_struct[0]),cmp);
4142 *modrat=chi_struct[0].
value/(chi_struct[0].
value+chi_struct[1].
value);
4143 chi_old[0]=chi_struct[0].
value*(1-*modrat)+chi_struct[1].
value*(*modrat);
4146 for(
im=0;
im<nmodels;
im++)
4156 if(chi_obs>chi_struct[1].
value){
4158 diff_1=chi_struct[0].
value;
4163 chi_struct[1].
value=diff_1;
4166 *modmin=chi_struct[0].
index;
4167 *modmax=chi_struct[1].
index;
4169 *modrat=(chi_obs-chi_struct[0].
value)/(chi_struct[1].
value-chi_struct[0].
value);
4171 derv_modrat_chi0= (chi_obs-chi_struct[1].
value)/pow(chi_struct[1].
value-chi_struct[0].
value,2);
4172 derv_modrat_chi1=-(chi_obs-chi_struct[0].
value)/pow(chi_struct[1].
value-chi_struct[0].
value,2);
4174 for (ib = 0; ib <nbands_ac; ib++){
4175 iw=acbands_index[ib];
4177 derv_modrat_rhorc[ib]=derv_chi_obs_rhorc[ib]/(chi_struct[1].
value-chi_struct[0].
value);
4178 derv_modrat_rhow_l+=-derv_modrat_rhorc[ib]*uncertainty->ratio_rhow[ib];
4182 derv_modrat_rhorc[wave_base_index]=derv_modrat_chi0*derv_chi_rhorc[*modmin][wave_base_index]+derv_modrat_chi1*derv_chi_rhorc[*modmax][wave_base_index];
4183 derv_modrat_rhorc[wave_base_index]+=derv_chi_obs_rhorc[wave_base_index]/(chi_struct[1].
value-chi_struct[0].
value);
4184 derv_modrat_rhow_l+=-derv_modrat_rhorc[wave_base_index];
4186 derv_modrat_taua_l =derv_modrat_chi0*derv_chi_taua_l [*modmin]+derv_modrat_chi1*derv_chi_taua_l [*modmax];
4187 derv_modrat_taua_l +=derv_chi_obs_taua_l/(chi_struct[1].
value-chi_struct[0].
value);
4189 uncertainty->derv_modrat_rhow_l=derv_modrat_rhow_l;
4190 uncertainty->derv_modrat_taua_l=derv_modrat_taua_l;
4194 if(
comp_rhoa_pca(nwave, wave,
geom, tau_iwnir_l[*modmin], mindx[*modmin], tau_pred_min, rho_pred_min,derv_rhoa_min,derv_taua_min) !=0)
4196 if(
comp_rhoa_pca(nwave, wave,
geom, tau_iwnir_l[*modmax], mindx[*modmax], tau_pred_max, rho_pred_max,derv_rhoa_max,derv_taua_max) !=0)
4199 for (iw = 0; iw <nwave; iw++) {
4201 rhoa[iw]=rho_pred_min[iw]*(1-*modrat)+rho_pred_max[iw]*(*modrat);
4204 for(itemp=0;itemp<nbands_ac;itemp++){
4205 ib=acbands_index[itemp];
4207 uncertainty->derv_La_rhorc[iw][itemp]=(rho_pred_max[iw]-rho_pred_min[iw])*derv_modrat_rhorc[itemp];
4208 uncertainty->derv_La_rhow_l[iw]+=-uncertainty->derv_La_rhorc[iw][itemp]*uncertainty->ratio_rhow[itemp];
4211 uncertainty->derv_La_rhorc[iw][wave_base_index]=(rho_pred_max[iw]-rho_pred_min[iw])*derv_modrat_rhorc[wave_base_index];
4212 uncertainty->derv_La_rhorc[iw][wave_base_index]+= (1-*modrat)*derv_rhoa_min[iw]*derv_taua_rhorc_l[*modmin];
4213 uncertainty->derv_La_rhorc[iw][wave_base_index]+= (*modrat) *derv_rhoa_max[iw]*derv_taua_rhorc_l[*modmax];
4214 uncertainty->derv_La_rhow_l[iw]+=-uncertainty->derv_La_rhorc[iw][wave_base_index];
4216 uncertainty->derv_La_taua_l[iw]=(rho_pred_max[iw]-rho_pred_min[iw])*derv_modrat_taua_l;
4217 uncertainty->derv_La_taua_l[iw]+=(1-*modrat)*(-derv_rhoa_min[iw]*derv_taua_rhorc_l[*modmin]*uncertainty->derv_Lg_taua[iw]);
4218 uncertainty->derv_La_taua_l[iw]+=(*modrat) *(-derv_rhoa_max[iw]*derv_taua_rhorc_l[*modmax]*uncertainty->derv_Lg_taua[iw]);
4220 derv_taua_min_rhorc_l[iw]=derv_taua_min[iw]*derv_taua_rhorc_l[*modmin];
4221 derv_taua_max_rhorc_l[iw]=derv_taua_max[iw]*derv_taua_rhorc_l[*modmax];
4222 derv_taua_min_taua_l[iw]=-derv_taua_min_rhorc_l[iw]*uncertainty->derv_Lg_taua[iw];
4223 derv_taua_max_taua_l[iw]=-derv_taua_max_rhorc_l[iw]*uncertainty->derv_Lg_taua[iw];
4224 derv_taua_min_rhow_l[iw]=-derv_taua_min_rhorc_l[iw];
4225 derv_taua_max_rhow_l[iw]=-derv_taua_max_rhorc_l[iw];
4230 *modmin=mindx[*modmin];
4231 *modmax=mindx[*modmax];
4251 int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *
geom,
float wv,
float rhoa[], int32_t *modmin, int32_t *modmax,
4252 float32 *modrat,
float *tau_pred_min,
float *tau_pred_max, int32_t ip,
float chi[],
float *mbac_w,uncertainty_t *uncertainty) {
4254 if (!have_ms && !use_pca_lut) {
4255 printf(
"\nThe multi-scattering spectral matching atmospheric correction method requires\n");
4256 printf(
"ams_all, bms_all, cms_all, dms_all, ems in the aerosol model tables.\n");
4263 status=
smaer_pca(
sensorID,wave,nwave,iwnir_s,iwnir_l,nmodels,mindx,
geom,wv,rhoa,
4264 modmin,modmax,modrat,tau_pred_min,tau_pred_max,ip,chi,mbac_w,uncertainty);
4268 float *
ac, *
bc, *cc, *dc, *ec;
4269 float ext_iwnir_l,tau_iwnir_l;
4270 double ax, bx, cx, fx;
4271 static int firstcall=1, nir_l;
4275 static float **tau_all_wav, **rho_all_wav_pred;
4277 float ext_coef, lg_taua, lg_rhoa;
4279 int iwtab,
im, modl;
4281 float diff_2 = 0.0, diff_1 = 0.0;
4282 float mbac_wsum = 0.0;
4283 static float *chi_old;
4287 static int32_t nbands_ac, *acbands_index,wave_base_index;
4289 float **derv_chi_rhorc;
4290 float *derv_chi_taua_l;
4291 float **derv_rhoa_rhorc_l;
4292 float *derv_modrat_rhorc;
4293 float derv_modrat_taua_l=0.;
4294 float derv_modrat_rhow_l=0.;
4295 float *derv_taua_min_rhorc_l;
4296 float *derv_taua_min_taua_l;
4297 float *derv_taua_min_rhow_l;
4298 float *derv_taua_max_rhorc_l;
4299 float *derv_taua_max_taua_l;
4300 float *derv_taua_max_rhow_l;
4301 float **derv_taua_rhorc_l;
4302 float *derv_chi_obs_rhorc;
4303 float derv_chi_obs_taua_l;
4305 float derv_temp_rhorc;
4306 float derv_modrat_chi0,derv_modrat_chi1;
4309 derv_modrat_rhorc=uncertainty->derv_modrat_rhorc;
4310 derv_taua_min_rhorc_l=uncertainty->derv_taua_min_rhorc_l;
4311 derv_taua_min_taua_l=uncertainty->derv_taua_min_taua_l;
4312 derv_taua_min_rhow_l=uncertainty->derv_taua_min_rhow_l;
4313 derv_taua_max_rhorc_l=uncertainty->derv_taua_max_rhorc_l;
4314 derv_taua_max_taua_l=uncertainty->derv_taua_max_taua_l;
4315 derv_taua_max_rhow_l=uncertainty->derv_taua_max_rhow_l;
4317 derv_chi_taua_l =(
float *)malloc(nmodels*
sizeof(
float));
4318 derv_chi_rhorc =(
float **)malloc(nmodels*
sizeof(
float *));
4319 derv_rhoa_rhorc_l=(
float **)malloc(nmodels*
sizeof(
float*));
4320 derv_taua_rhorc_l=(
float **)malloc(nmodels*
sizeof(
float*));
4322 for(
im=0;
im<nmodels;
im++){
4323 derv_chi_rhorc[
im] =(
float *)malloc(nbands_ac*
sizeof(
float));
4324 derv_rhoa_rhorc_l[
im]=(
float *)malloc(nwave*
sizeof(
float));
4325 derv_taua_rhorc_l[
im]=(
float *)malloc(nwave*
sizeof(
float));
4327 for(
im=0;
im<nmodels;
im++)
4328 for(iw=0;iw<nbands_ac;iw++)
4329 derv_chi_rhorc[
im][iw]=0.;
4331 for (iw = 0; iw < nbands_ac; iw++)
4332 derv_modrat_rhorc[iw] = 0.;
4334 derv_chi_obs_rhorc = (
float *)malloc(nbands_ac *
sizeof(
float));
4341 tau_all_wav = (
float**) malloc(nwave *
sizeof (
float*));
4342 rho_all_wav_pred = (
float**) malloc(nwave *
sizeof (
float*));
4344 chi_old = (
float *)malloc((nmodels) *
sizeof(
float));
4346 for (
i = 0;
i < nwave;
i++) {
4347 tau_all_wav[
i] = (
float*) malloc((nmodels) *
sizeof (
float));
4348 rho_all_wav_pred[
i] = (
float*) malloc((nmodels) *
sizeof (
float));
4351 nbands_ac=
input->nbands_ac;
4352 acbands_index=
input->acbands_index;
4354 float *tempwave=(
float *)malloc(nbands_ac*
sizeof(
float));
4355 for(iw=0;iw<nbands_ac;iw++)
4356 tempwave[iw]=
input->mbac_wave[iw];
4358 wave_base_index=
windex(
input->aer_wave_base*1.,tempwave,nbands_ac);
4366 struct str chi_struct[nmodels];
4368 iwtab_l = iwatab[nir_l];
4370 float chi_obs = 0.0;
4372 for (iw = 0; iw <nbands_ac; iw++)
4373 derv_chi_obs_rhorc[iw ] = 0.;
4374 derv_chi_obs_taua_l = 0.;
4377 for (ib = 0; ib <nbands_ac; ib++) {
4379 iw=acbands_index[ib];
4380 if (mbac_w[iw] == 0. || iw == nir_l)
4382 chi_obs +=rhoa[iw] / rhoa[nir_l];
4383 mbac_wsum += mbac_w[iw];
4385 derv_chi_obs_rhorc[wave_base_index] = 1 / rhoa[nir_l];
4386 derv_chi_obs_rhorc[wave_base_index] += (-rhoa[iw] / rhoa[nir_l] / rhoa[nir_l]);
4387 derv_chi_obs_taua_l += (-1 / rhoa[nir_l] * uncertainty->derv_Lg_taua[iw] +rhoa[iw] / rhoa[nir_l] / rhoa[nir_l] * uncertainty->derv_Lg_taua[nir_l]);
4390 chi_obs /= mbac_wsum;
4392 derv_chi_obs_taua_l /= mbac_wsum;
4393 for (ib = 0; ib <nbands_ac; ib++)
4394 derv_chi_obs_rhorc[ib] /= mbac_wsum;
4397 for (
im = 0;
im < nmodels;
im++) {
4403 ax = (
double)
ac[iwtab_l] - log((
double) rhoa[nir_l]);
4405 cx = (
double) cc[iwtab_l];
4407 fx = bx * bx - 4.0 *
ax*cx;
4408 if (fx > 0.0 && cx != 0.0) {
4409 tau_iwnir_l = 0.5 * (-bx + sqrt(fx)) / cx;
4411 derv_temp_rhorc= 1 / rhoa[nir_l] / sqrt(fx);
4413 tau_iwnir_l = exp(tau_iwnir_l);
4415 derv_temp_rhorc *= tau_iwnir_l;
4424 ext_iwnir_l = aertab->model[modl]->extc[iwtab_l];
4425 for (iw = 0; iw <nwave ; iw++) {
4427 ext_coef = aertab->model[modl]->extc[iwtab];
4428 tau_all_wav[iw][
im] = (ext_coef / ext_iwnir_l) * tau_iwnir_l;
4429 lg_taua = log(tau_all_wav[iw][
im]);
4431 derv_taua_rhorc_l[
im][iw]=(ext_coef / ext_iwnir_l)*derv_temp_rhorc;
4432 derv_rhoa_rhorc_l[
im][iw]=1/tau_all_wav[iw][
im]*(ext_coef / ext_iwnir_l)*derv_temp_rhorc;
4435 lg_rhoa =
ac[iwtab] +
bc[iwtab] * lg_taua + cc[iwtab] * pow(lg_taua, 2) +
4436 dc[iwtab] * pow(lg_taua, 3) + ec[iwtab] * pow(lg_taua, 4);
4439 derv_rhoa_rhorc_l[
im][iw] *=
4440 (
bc[iwtab] + 2 * cc[iwtab] * lg_taua + 3 * dc[iwtab] * pow(lg_taua, 2) +
4441 4 * ec[iwtab] * pow(lg_taua, 3));
4443 rho_all_wav_pred[iw][
im] = exp(lg_rhoa);
4446 derv_rhoa_rhorc_l[
im][iw]*= rho_all_wav_pred[iw][
im];
4452 derv_chi_rhorc[
im][wave_base_index]=0.;
4453 derv_chi_taua_l[
im]=0.;
4456 for (ib = 0; ib <nbands_ac; ib++) {
4457 iw=acbands_index[ib];
4458 if(mbac_w[iw]==0.|| iw==nir_l)
4461 diff_2 += rho_all_wav_pred[iw][
im] / rhoa[nir_l];
4462 diff_1 += pow((rhoa[iw] - rho_all_wav_pred[iw][
im]) / noise[iw], 2) * mbac_w[iw];
4463 mbac_wsum += mbac_w[iw];
4466 derv_temp_rhorc = (derv_rhoa_rhorc_l[
im][iw] / rhoa[nir_l] -
4467 rho_all_wav_pred[iw][
im] / rhoa[nir_l] / rhoa[nir_l]);
4468 derv_chi_rhorc[
im][wave_base_index] += derv_temp_rhorc;
4469 derv_chi_taua_l[
im] += (-derv_temp_rhorc * uncertainty->derv_Lg_taua[nir_l]);
4473 chi[
im] = diff_2 / mbac_wsum;
4474 chi_old[
im] = diff_1 / mbac_wsum;
4476 derv_chi_rhorc[
im][wave_base_index] /= mbac_wsum;
4477 derv_chi_taua_l[
im] /= mbac_wsum;
4487 qsort(chi_struct,nmodels,
sizeof(chi_struct[0]),cmp);
4489 *modrat = chi_struct[0].
value / (chi_struct[0].
value + chi_struct[1].
value);
4491 chi_struct[0].
value * (1 - *modrat) +
4492 chi_struct[1].
value * (*modrat);
4496 for (
im = 0;
im < nmodels;
im++)
4499 if ((chi_obs - chi_struct[0].
value) * (chi_obs - chi_struct[1].
value) > 0) {
4501 *modmin = chi_struct[0].
index;
4505 if (chi_obs > chi_struct[1].
value) {
4507 diff_1 = chi_struct[0].
value;
4512 chi_struct[1].
value = diff_1;
4515 *modmin = chi_struct[0].
index;
4516 *modmax = chi_struct[1].
index;
4518 *modrat = (chi_obs - chi_struct[0].
value) / (chi_struct[1].
value - chi_struct[0].
value);
4521 (chi_obs - chi_struct[1].
value) / pow(chi_struct[1].
value - chi_struct[0].
value, 2);
4523 -(chi_obs - chi_struct[0].
value) / pow(chi_struct[1].
value - chi_struct[0].
value, 2);
4525 for (ib = 0; ib < nbands_ac; iw++) {
4526 iw=acbands_index[ib];
4528 derv_modrat_rhorc[ib] =
4529 derv_chi_obs_rhorc[ib] / (chi_struct[1].
value - chi_struct[0].
value);
4530 derv_modrat_rhow_l +=
4531 -derv_modrat_rhorc[ib] * uncertainty->ratio_rhow[ib];
4535 derv_modrat_rhorc[wave_base_index] = derv_modrat_chi0 * derv_chi_rhorc[*modmin][wave_base_index] +
4536 derv_modrat_chi1 * derv_chi_rhorc[*modmax][wave_base_index];
4537 derv_modrat_rhorc[wave_base_index] +=
4538 derv_chi_obs_rhorc[wave_base_index] / (chi_struct[1].
value - chi_struct[0].
value);
4539 derv_modrat_rhow_l += -derv_modrat_rhorc[wave_base_index];
4541 derv_modrat_taua_l =
4542 derv_modrat_chi0 * derv_chi_taua_l[*modmin] + derv_modrat_chi1 * derv_chi_taua_l[*modmax];
4543 derv_modrat_taua_l += derv_chi_obs_taua_l / (chi_struct[1].
value - chi_struct[0].
value);
4545 uncertainty->derv_modrat_rhow_l=derv_modrat_rhow_l;
4546 uncertainty->derv_modrat_taua_l=derv_modrat_taua_l;
4550 for (iw = 0; iw <nwave; iw++) {
4552 tau_pred_min[iw] = tau_all_wav[iw][*modmin];
4553 tau_pred_max[iw] = tau_all_wav[iw][*modmax];
4554 rhoa[iw]=rho_all_wav_pred[iw][*modmin]*(1-*modrat)+rho_all_wav_pred[iw][*modmax]*(*modrat);
4557 for (ib = 0; ib < nbands_ac; ib++) {
4558 i=acbands_index[ib];
4560 uncertainty->derv_La_rhorc[iw][ib] =
4561 (rho_all_wav_pred[iw][*modmax] - rho_all_wav_pred[iw][*modmin]) *
4562 derv_modrat_rhorc[ib];
4563 uncertainty->derv_La_rhow_l[iw] +=
4564 -uncertainty->derv_La_rhorc[iw][ib] * uncertainty->ratio_rhow[ib];
4567 uncertainty->derv_La_rhorc[iw][wave_base_index] =
4568 (rho_all_wav_pred[iw][*modmax] - rho_all_wav_pred[iw][*modmin]) *
4569 derv_modrat_rhorc[nir_l - iwnir_s];
4570 uncertainty->derv_La_rhorc[iw][wave_base_index] +=
4571 (1 - *modrat) * derv_rhoa_rhorc_l[*modmin][iw];
4572 uncertainty->derv_La_rhorc[iw][wave_base_index] += (*modrat) * derv_rhoa_rhorc_l[*modmax][iw];
4573 uncertainty->derv_La_rhow_l[iw] += -uncertainty->derv_La_rhorc[iw][wave_base_index];
4575 uncertainty->derv_La_taua_l[iw]=(rho_all_wav_pred[iw][*modmax]-rho_all_wav_pred[iw][*modmin])*derv_modrat_taua_l;
4576 uncertainty->derv_La_taua_l[iw]+=(1-*modrat)*(-derv_rhoa_rhorc_l[*modmin][iw]*uncertainty->derv_Lg_taua[iw]);
4577 uncertainty->derv_La_taua_l[iw]+=(*modrat) *(-derv_rhoa_rhorc_l[*modmax][iw]*uncertainty->derv_Lg_taua[iw]);
4579 derv_taua_min_rhorc_l[iw]=derv_taua_rhorc_l[*modmin][iw];
4580 derv_taua_max_rhorc_l[iw]=derv_taua_rhorc_l[*modmax][iw];
4581 derv_taua_min_taua_l[iw]=-derv_taua_rhorc_l[*modmin][iw]*uncertainty->derv_Lg_taua[iw];
4582 derv_taua_max_taua_l[iw]=-derv_taua_rhorc_l[*modmax][iw]*uncertainty->derv_Lg_taua[iw];
4583 derv_taua_min_rhow_l[iw]=-derv_taua_rhorc_l[*modmin][iw];
4584 derv_taua_max_rhow_l[iw]=-derv_taua_rhorc_l[*modmax][iw];
4589 *modmin = mindx[*modmin];
4590 *modmax = mindx[*modmax];
4591 chi[0] = chi_old[0];
4601 for(
im=0;
im<nmodels;
im++){
4602 free(derv_chi_rhorc[
im]);
4603 free(derv_rhoa_rhorc_l[
im]);
4604 free(derv_taua_rhorc_l[
im]);
4606 free(derv_chi_rhorc);
4607 free(derv_chi_taua_l);
4608 free(derv_rhoa_rhorc_l);
4609 free(derv_taua_rhorc_l);
4610 free(derv_chi_obs_rhorc);
4624 int32_t nmodels, int32_t mindx[],
4625 geom_str *
geom,
float wv,
float rhoa[],
4626 int32_t *modmin, int32_t *modmax,
float *modrat,
float *epsnir,
4627 float tau_pred_min[],
float tau_pred_max[],
int ip,uncertainty_t *uncertainty) {
4630 float rho_pred_min[nwave], rho_pred_max[nwave];
4631 float rho_aer[nwave], tau_aer[nwave];
4633 if (!have_ms && !use_pca_lut) {
4634 printf(
"\nThe multi-scattering epsilon atmospheric correction method requires\n");
4635 printf(
"ams_all, bms_all, cms_all, dms_all, ems in the aerosol model tables.\n");
4640 rhoa, modmin, modmax, modrat, epsnir,
4641 tau_pred_max, tau_pred_min, rho_pred_max, rho_pred_min, tau_aer, rho_aer) != 0)
4646 rhoa, modmin, modmax, modrat, epsnir,
4647 tau_pred_max, tau_pred_min, rho_pred_max, rho_pred_min, tau_aer, rho_aer,ip,uncertainty) != 0)
4651 for (iw = 0; iw < nwave; iw++) {
4652 rhoa[iw] = rho_aer[iw];
4667 int32_t iwnir_l, int32_t nmodels, int32_t mindx[], geom_str *
geom,
4668 float wv,
float rhoa[], int32_t *modmin, int32_t *modmax,
4669 float *modrat,
float *epsnir,
float tauamin[],
float tauamax[]) {
4685 if ((rhoasmin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4686 printf(
"Unable to allocate space for rhoasmin.\n");
4689 if ((rhoasmax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4690 printf(
"Unable to allocate space for rhoasmax.\n");
4693 if ((rhoamin = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4694 printf(
"Unable to allocate space for rhoamin.\n");
4697 if ((rhoamax = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4698 printf(
"Unable to allocate space for rhoasmax.\n");
4704 rhoa, iwnir_s, iwnir_l, modmin, modmax, modrat, epsnir);
4708 cc = log(*epsnir) / (wave[iwnir_l] - wave[iwnir_s]);
4730 for (iw = 0; iw < nwave; iw++) {
4732 epsmin1 = epsmin[iw];
4733 epsmax1 = epsmax[iw];
4736 epsmin1 = exp(cc * (wave[iwnir_l] - wave[iw]));
4740 rhoasmin[iw] = rhoasmin[iwnir_l] * epsmin1;
4741 rhoasmax[iw] = rhoasmax[iwnir_l] * epsmax1;
4749 for (iw = 0; iw < nwave; iw++) {
4750 rhoa[iw] = rhoamin[iw]*(1.0 - (*modrat)) + rhoamax[iw]*(*modrat);
4778 return (
x->rhoa <
y->rhoa ? -1 : 1);
4782 int32_t nmodel, int32_t mindx[], geom_str *
geom,
float wv,
4783 float rhoa[], int32_t iwnir_s, int32_t iwnir_l, int32_t *modmin,
4784 int32_t *modmax,
float *modrat,
float *epsnir) {
4800 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4801 printf(
"Unable to allocate space for rhoas.\n");
4804 if ((rhoa_tmp = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
4805 printf(
"Unable to allocate space for rhoas_tmp.\n");
4811 for (jm = 0; jm < nmodel; jm++) {
4825 rhoas[iwnir_s] = rhoas[iwnir_l] *
eps[iwnir_s];
4830 rhoa_tab[jm].modnum =
im;
4831 rhoa_tab[jm].rhoa = rhoa_tmp[iwnir_s];
4832 rhoa_tab[jm].eps =
eps[iwnir_s];
4836 qsort(rhoa_tab, nmodel,
sizeof (rhoaTstr), (
int (*)(
const void *,
const void *))
comp_rhoaT);
4839 for (jm = 0; jm < nmodel; jm++) {
4840 if (rhoa_tab[jm].rhoa > rhoa[iwnir_s])
4846 }
else if (jm == nmodel) {
4853 wt = (rhoa[iwnir_s] - rhoa_tab[jm1].rhoa) / (rhoa_tab[jm2].rhoa - rhoa_tab[jm1].rhoa);
4855 *modmin = rhoa_tab[jm1].modnum;
4856 *modmax = rhoa_tab[jm2].modnum;
4858 *epsnir = rhoa_tab[jm1].eps * (1.0 - wt) + rhoa_tab[jm2].
eps*wt;
4871 static int order_models(
const void *
p1,
const void *p2) {
4872 aermodstr *
x = *(aermodstr **)
p1;
4873 aermodstr *
y = *(aermodstr **) p2;
4875 if (
x->rh ==
y->rh) {
4894 int rhaer(int32_t
sensorID,
float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
4895 geom_str *
geom,
float wv,
float rh,
float pr,
float taur[],
float rhoa[],
4896 int32_t *modmin1, int32_t *modmax1,
float *modrat1, int32_t *modmin2, int32_t *modmax2,
float *modrat2,
4897 float *
eps,
float taua[],
float tsol[],
float tsen[], int32_t ip,
float *mbac_w,
float *chi, uncertainty_t *uncertainty) {
4898 static int firstCall = 1;
4919 float *tau_pred_min1;
4920 float *tau_pred_max1;
4921 float *tau_pred_min2;
4922 float *tau_pred_max2;
4926 int irh1, irh2, irh;
4930 static uncertainty_t *uncertainty1=
NULL, *uncertainty2=
NULL;
4937 float lastrh = -1.0;
4940 printf(
"-E- %s line %d: This aerosol selection method requires models with a Relative Humidity attribute and size distribution.\n",
4941 __FILE__, __LINE__);
4946 uncertainty1 = (uncertainty_t *) malloc(
sizeof(uncertainty_t));
4947 uncertainty2 = (uncertainty_t *) malloc(
sizeof(uncertainty_t));
4949 if (
alloc_uncertainty(uncertainty->nbands, uncertainty->nbands_ac, uncertainty->npix, uncertainty1) != 0) {
4950 printf(
"unable to allocate error record in rhaer()\n");
4953 if (
alloc_uncertainty(uncertainty->nbands, uncertainty->nbands_ac, uncertainty->npix, uncertainty2) != 0) {
4954 printf(
"unable to allocate error record in rhaer()\n");
4959 qsort(aertab->model, aertab->nmodel, sizeof (aermodstr*), (
int (*)(
const void *,
const void *)) order_models);
4967 for (
im = 0;
im < aertab->nmodel;
im++) {
4968 if (aertab->model[
im]->rh != lastrh) {
4969 rhtab[nrh] = aertab->model[
im]->rh;
4970 lastrh = rhtab[nrh];
4973 if (nrh == 1 && aertab->model[
im]->sd != lastsd) {
4974 sdtab[nsd] = aertab->model[
im]->sd;
4975 lastsd = sdtab[nsd];
4979 if (nrh * nsd != aertab->nmodel) {
4980 printf(
"-E- %s line %d: number of humidities (%d) x number of size distributions (%d) must equal number of models (%d).\n",
4981 __FILE__, __LINE__, nrh, nsd, aertab->nmodel);
4984 printf(
"%d aerosol models: %d humidities x %d size fractions\n", aertab->nmodel, nrh, nsd);
4985 for (irh = 0; irh < nrh; irh++) {
4986 for (isd = 0; isd < nsd; isd++) {
4987 im = irh * nsd + isd;
4988 printf(
"model %d, rh=%f, sd=%d, alpha=%f, name=%s\n",
4989 im, aertab->model[
im]->rh, aertab->model[
im]->sd, aertab->model[
im]->angstrom[0], aertab->model[
im]->name);
5001 printf(
"unable to copy the error record in rhaer()\n");
5005 printf(
"unable to copy the error record in rhaer()\n");
5011 if ((taua1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5012 printf(
"Unable to allocate space for taua1.\n");
5015 if ((taua2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5016 printf(
"Unable to allocate space for taua2.\n");
5019 if ((tsol1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5020 printf(
"Unable to allocate space for tsol1.\n");
5023 if ((tsol2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5024 printf(
"Unable to allocate space for tsol2.\n");
5027 if ((rhoa1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5028 printf(
"Unable to allocate space for rhoa1.\n");
5031 if ((rhoa2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5032 printf(
"Unable to allocate space for rhoa2.\n");
5035 if ((tsen1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5036 printf(
"Unable to allocate space for tsen1.\n");
5039 if ((tsen2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5040 printf(
"Unable to allocate space for tsen2.\n");
5043 if ((tau_pred_min1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5044 printf(
"Unable to allocate space for tau_pred_min1.\n");
5047 if ((tau_pred_min2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5048 printf(
"Unable to allocate space for tau_pred_min2.\n");
5051 if ((tau_pred_max1 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5052 printf(
"Unable to allocate space for tau_pred_max1.\n");
5055 if ((tau_pred_max2 = (
float *) malloc(nwave *
sizeof (
float))) ==
NULL) {
5056 printf(
"Unable to allocate space for tau_pred_max2.\n");
5060 if ((chi1 = (
float *) calloc(nmodels,
sizeof (
float))) ==
NULL) {
5061 printf(
"Unable to allocate space for chi1.\n");
5064 if ((chi2 = (
float *) calloc(nmodels,
sizeof (
float))) ==
NULL) {
5065 printf(
"Unable to allocate space for chi2.\n");
5070 for (iw = 0; iw < nwave; iw++) {
5074 rhoa1[iw] = rhoa[iw];
5075 rhoa2[iw] = rhoa[iw];
5089 if (nrh == 1 || rhtab[0] > rh) {
5094 }
else if (rhtab[nrh - 1] < rh) {
5100 for (irh = 0; irh < nrh; irh++) {
5101 if (rhtab[irh] > rh)
5104 irh1 =
MIN(
MAX(0, irh - 1), nrh - 2);
5106 wt = (rh - rhtab[irh1]) / (rhtab[irh2] - rhtab[irh1]);
5108 derv_wt_rh = 1.0 / (rhtab[irh2] - rhtab[irh1]);
5113 for (
im = 0;
im < nsd;
im++) {
5114 mindx1[
im] = irh1 * nsd +
im;
5115 mindx2[
im] = irh2 * nsd +
im;
5127 if (
smaer(
sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx1,
geom, wv, rhoa1, modmin1, modmax1, modrat1, tau_pred_min1, tau_pred_max1,ip, chi1, mbac_w,uncertainty1) != 0) {
5137 free(tau_pred_min1);
5138 free(tau_pred_max1);
5139 free(tau_pred_min2);
5140 free(tau_pred_max2);
5145 geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1,ip,uncertainty1) != 0) {
5155 free(tau_pred_min1);
5156 free(tau_pred_max1);
5157 free(tau_pred_min2);
5158 free(tau_pred_max2);
5164 geom, wv, rhoa1, modmin1, modmax1, modrat1, &eps1, tau_pred_min1, tau_pred_max1) != 0) {
5173 free(tau_pred_min1);
5174 free(tau_pred_max1);
5175 free(tau_pred_min2);
5176 free(tau_pred_max2);
5182 *modmin1, *modmax1, *modrat1, rhoa1, taua1, tsol1, tsen1, tau_pred_min1, tau_pred_max1, 1, ip,uncertainty1);
5190 if (
smaer(
sensorID, wave, nwave, iwnir_s, iwnir_l, nmodels, mindx2,
5191 geom, wv, rhoa2, modmin2, modmax2, modrat2, tau_pred_min2, tau_pred_max2, ip, chi2, mbac_w,uncertainty2) != 0) {
5201 free(tau_pred_min1);
5202 free(tau_pred_max1);
5203 free(tau_pred_min2);
5204 free(tau_pred_max2);
5213 geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2,ip,uncertainty2) != 0) {
5222 free(tau_pred_min1);
5223 free(tau_pred_max1);
5224 free(tau_pred_min2);
5225 free(tau_pred_max2);
5230 geom, wv, rhoa2, modmin2, modmax2, modrat2, &eps2, tau_pred_min2, tau_pred_max2) != 0) {
5239 free(tau_pred_min1);
5240 free(tau_pred_max1);
5241 free(tau_pred_min2);
5242 free(tau_pred_max2);
5248 *modmin2, *modmax2, *modrat2, rhoa2, taua2, tsol2, tsen2, tau_pred_min2, tau_pred_max2, 1,ip,uncertainty2);
5250 for (iw = 0; iw < nwave; iw++) {
5251 rhoa[iw] = rhoa1[iw]*(1 - wt) + rhoa2[iw] * wt;
5252 taua[iw] = taua1[iw]*(1 - wt) + taua2[iw] * wt;
5253 tsol[iw] = tsol1[iw]*(1 - wt) + tsol2[iw] * wt;
5254 tsen[iw] = tsen1[iw]*(1 - wt) + tsen2[iw] * wt;
5257 for(inir=0;inir<=iwnir_l-iwnir_s;inir++){
5258 uncertainty->derv_La_rhorc[iw][inir]=(1-wt)*uncertainty1->derv_La_rhorc[iw][inir]+ wt*uncertainty2->derv_La_rhorc[iw][inir];
5259 uncertainty->derv_taua_rhorc[iw][inir]=(1-wt)*uncertainty1->derv_taua_rhorc[iw][inir]+ wt*uncertainty2->derv_taua_rhorc[iw][inir];
5260 uncertainty->derv_tsen_rhorc[iw][inir]=(1-wt)*uncertainty1->derv_tsen_rhorc[iw][inir]+ wt*uncertainty2->derv_tsen_rhorc[iw][inir];
5261 uncertainty->derv_tsol_rhorc[iw][inir]=(1-wt)*uncertainty1->derv_tsol_rhorc[iw][inir]+ wt*uncertainty2->derv_tsol_rhorc[iw][inir];
5263 uncertainty->derv_La_taua_l[iw]=(1-wt)*uncertainty1->derv_La_taua_l[iw]+ wt*uncertainty2->derv_La_taua_l[iw];
5264 uncertainty->derv_La_rhow_l[iw]=(1-wt)*uncertainty1->derv_La_rhow_l[iw]+ wt*uncertainty2->derv_La_rhow_l[iw];
5265 uncertainty->derv_La_rh[iw]=(rhoa2[iw]-rhoa1[iw])*derv_wt_rh;
5268 uncertainty->derv_taua_taua_l[iw]=(1-wt)*uncertainty1->derv_taua_taua_l[iw]+ wt*uncertainty2->derv_taua_taua_l[iw];
5270 uncertainty->derv_taua_rhow_l[iw]=(1-wt)*uncertainty1->derv_taua_rhow_l[iw]+ wt*uncertainty2->derv_taua_rhow_l[iw];
5271 uncertainty->derv_taua_rh[iw]=(taua2[iw]-taua1[iw])*derv_wt_rh;
5276 uncertainty->derv_tsol_taua_l[iw]=(1-wt)*uncertainty1->derv_tsol_taua_l[iw]+ wt*uncertainty2->derv_tsol_taua_l[iw];
5277 uncertainty->derv_tsol_rhow_l[iw]=(1-wt)*uncertainty1->derv_tsol_rhow_l[iw]+ wt*uncertainty2->derv_tsol_rhow_l[iw];
5278 uncertainty->derv_tsol_rh[iw]=(tsol2[iw]-tsol1[iw])*derv_wt_rh;
5283 uncertainty->derv_tsen_taua_l[iw]=(1-wt)*uncertainty1->derv_tsen_taua_l[iw]+ wt*uncertainty2->derv_tsen_taua_l[iw];
5284 uncertainty->derv_tsen_rhow_l[iw]=(1-wt)*uncertainty1->derv_tsen_rhow_l[iw]+ wt*uncertainty2->derv_tsen_rhow_l[iw];
5285 uncertainty->derv_tsen_rh[iw]=(tsen2[iw]-tsen1[iw])*derv_wt_rh;
5289 *
eps = eps1 * (1 - wt) + eps2*wt;
5290 *chi =chi1[0]*(1-wt) +chi2[0]*wt;
5294 for (iw = 0; iw < nwave; iw++) {
5295 rhoa[iw] = rhoa1[iw];
5296 taua[iw] = taua1[iw];
5297 tsol[iw] = tsol1[iw];
5298 tsen[iw] = tsen1[iw];
5302 for(inir=0;inir<=iwnir_l-iwnir_s;inir++){
5303 uncertainty->derv_La_rhorc[iw][inir]=uncertainty1->derv_La_rhorc[iw][inir];
5304 uncertainty->derv_taua_rhorc[iw][inir]=uncertainty1->derv_taua_rhorc[iw][inir];
5305 uncertainty->derv_tsen_rhorc[iw][inir]=uncertainty1->derv_tsen_rhorc[iw][inir];
5306 uncertainty->derv_tsol_rhorc[iw][inir]=uncertainty1->derv_tsol_rhorc[iw][inir];
5308 uncertainty->derv_La_taua_l[iw]=uncertainty1->derv_La_taua_l[iw];
5310 uncertainty->derv_La_rhow_l[iw]=uncertainty1->derv_La_rhow_l[iw];
5315 uncertainty->derv_taua_taua_l[iw]=uncertainty1->derv_taua_taua_l[iw];
5317 uncertainty->derv_taua_rhow_l[iw]=uncertainty1->derv_taua_rhow_l[iw];
5322 uncertainty->derv_tsol_taua_l[iw]=uncertainty1->derv_tsol_taua_l[iw];
5323 uncertainty->derv_tsol_rhow_l[iw]=uncertainty1->derv_tsol_rhow_l[iw];
5328 uncertainty->derv_tsen_taua_l[iw]=uncertainty1->derv_tsen_taua_l[iw];
5329 uncertainty->derv_tsen_rhow_l[iw]=uncertainty1->derv_tsen_rhow_l[iw];
5344 free(tau_pred_min1);
5345 free(tau_pred_max1);
5346 free(tau_pred_min2);
5347 free(tau_pred_max2);
5361 int fixedaer(int32_t
sensorID, int32_t modnum,
float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
5362 char models[
MAXAERMOD][32], int32_t nmodels,
5363 geom_str *
geom,
float wv,
float rhoa[],
float *epsnir) {
5368 if ((rhoas = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5369 printf(
"Unable to allocate space for rhoas.\n");
5373 if (rhoa[iwnir_l] < 0.0) {
5386 printf(
"Error getting rhoas\n");
5392 for (iw = 0; iw < nwave; iw++) {
5393 rhoas[iw] = rhoas[iwnir_l] *
eps[iw];
5399 if (iwnir_s == iwnir_l)
5400 *epsnir =
eps[iwnir_l - 1];
5402 *epsnir =
eps[iwnir_s];
5417 int32_t iwnir_s, int32_t iwnir_l, geom_str *
geom,
float wv,
5418 int32_t modmin, int32_t modmax,
float modrat,
float rhoa[],
float *
eps) {
5426 if ((rhoa1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5427 printf(
"Unable to allocate space for rhoa1.\n");
5430 if ((rhoa2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5431 printf(
"Unable to allocate space for rhoa2.\n");
5435 if (modmin < 0 || modmin >=
input->naermodels ||
5436 modmax < 0 || modmax >=
input->naermodels ||
5437 modrat < 0.0 || modrat > 1.0) {
5438 printf(
"Bogus input for fixed model pair: %d %d %f\n", modmin + 1, modmax + 1, modrat);
5443 if (rhoa[iwnir_l] >
input->rhoamin) {
5445 rhoa2[iwnir_l] = rhoa1[iwnir_l] = rhoa[iwnir_l];
5455 for (iw = 0; iw < nwave; iw++) {
5457 rhoa[iw] =
MAX((1.0 - modrat) * rhoa1[iw] + modrat * rhoa2[iw], 0.0);
5459 *
eps = (1.0 - modrat) * eps1 + modrat*eps2;
5462 }
else if (rhoa[iwnir_l] > -(
input->rhoamin)) {
5466 for (iw = 0; iw < nwave; iw++) {
5467 rhoa[iw] =
MAX(rhoa[iwnir_l], 1e-6);
5493 int32_t iwnir_s, int32_t iwnir_l, geom_str *
geom,
float wv,
5494 int32_t *modmin, int32_t *modmax,
float *modrat,
float rhoa[],
5496 static int firstCall = 1;
5497 static int angst_band1 = -1;
5498 static int angst_band2 = -1;
5513 int ig, gmult, iw, iwtab;
5516 maxwave =
MAX(aertab->nwave, nwave);
5518 if ((rhoa1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5519 printf(
"Unable to allocate space for rhoa1.\n");
5522 if ((rhoa2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5523 printf(
"Unable to allocate space for rhoa2.\n");
5526 if ((rhoas1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5527 printf(
"Unable to allocate space for rhoas1.\n");
5530 if ((rhoas2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5531 printf(
"Unable to allocate space for rhoas2.\n");
5534 if ((
f1 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
5535 printf(
"Unable to allocate space for f1.\n");
5538 if ((
f2 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
5539 printf(
"Unable to allocate space for f2.\n");
5542 if ((lnf1 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
5543 printf(
"Unable to allocate space for lnf1.\n");
5546 if ((lnf2 = (
float *) calloc(maxwave,
sizeof (
float))) ==
NULL) {
5547 printf(
"Unable to allocate space for lnf2.\n");
5552 angst_band1 =
windex(520, wave, nwave);
5553 angst_band2 =
windex(865, wave, nwave);
5558 for (iw = 0; iw < nwave; iw++) {
5559 if (aot[iw] < 0.0) {
5573 if (aot[iwnir_l] > 0.0)
5574 angstrom = -log(aot[angst_band1] / aot[angst_band2]) /
5575 log(wave[angst_band1] / wave[angst_band2]);
5586 gmult = (interpol == 1) ? 0 :
geom->gmult;
5589 for (iw = 0; iw < aertab->nwave; iw++) {
5591 f1[iw] = aertab->model[*modmin]->albedo[iw] *
5592 phase1[iw] / 4.0 /
geom->csolz[ig] /
geom->csenz[ig];
5593 f2[iw] = aertab->model[*modmax]->albedo[iw] *
5594 phase2[iw] / 4.0 /
geom->csolz[ig] /
geom->csenz[ig];
5596 lnf1[iw] = log(
f1[iw]);
5597 lnf2[iw] = log(
f2[iw]);
5603 for (iw = 0; iw < nwave; iw++) {
5605 if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0) {
5606 rhoas1[iw] = aot[iw] * exp(
linterp(aertab->wave, lnf1, aertab->nwave, wave[iw]));
5607 rhoas2[iw] = aot[iw] * exp(
linterp(aertab->wave, lnf2, aertab->nwave, wave[iw]));
5609 rhoas1[iw] = aot[iw] *
f1[iwtab];
5610 rhoas2[iw] = aot[iw] *
f2[iwtab];
5614 for (iw = 0; iw < nwave; iw++) {
5616 rhoas1[iw] = aot[iw] *
f1[iwtab];
5617 rhoas2[iw] = aot[iw] *
f2[iwtab];
5620 eps1 = rhoas1[iwnir_s] / rhoas1[iwnir_l];
5621 eps2 = rhoas2[iwnir_s] / rhoas2[iwnir_l];
5627 for (iw = 0; iw < nwave; iw++) {
5628 rhoa[iw] = (1.0 - *modrat) * rhoa1[iw] + *modrat * rhoa2[iw];
5630 *epsnir = (1.0 - *modrat) * eps1 + *modrat * eps2;
5651 int aerosol(l2str *l2rec, int32_t aer_opt_in, aestr *aerec, int32_t ip,
5652 float wave[], int32_t nwave, int32_t iwnir_s_in, int32_t iwnir_l_in,
5653 float F0_in[],
float La1_in[],
float La2_out[],
5654 float t_sol_out[],
float t_sen_out[],
float *
eps,
float taua_out[],
5655 int32_t *modmin, int32_t *modmax,
float *modrat,
5656 int32_t *modmin2, int32_t *modmax2,
float *modrat2,
float *mbac_w) {
5657 static int firstCall = 1;
5676 float *taua_pred_min;
5677 float *taua_pred_max;
5679 static int32_t last_scan=-99;
5680 static int32_t aer_wave_base;
5684 l1str *
l1rec = l2rec->l1rec;
5685 uncertainty_t *uncertainty=
l1rec->uncertainty;
5687 static float *noise_temp;
5695 last_scan =
l1rec->iscan;
5701 float wv =
l1rec->wv [ip];
5702 float rh =
l1rec->rh [ip];
5708 Maxband = nwave + 1;
5711 if ((rhoa = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5712 printf(
"Unable to allocate space for rhoa.\n");
5715 if ((radref = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5716 printf(
"Unable to allocate space for raderef.\n");
5719 if ((
F0 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5720 printf(
"Unable to allocate space for F0.\n");
5723 if ((taur = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5724 printf(
"Unable to allocate space for taur.\n");
5727 if ((La1 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5728 printf(
"Unable to allocate space for La1.\n");
5731 if ((La2 = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5732 printf(
"Unable to allocate space for La2.\n");
5735 if ((t_sol = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5736 printf(
"Unable to allocate space for t_sol.\n");
5739 if ((t_sen = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5740 printf(
"Unable to allocate space for t_sen.\n");
5743 if ((taua = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5744 printf(
"Unable to allocate space for taua.\n");
5747 if ((taua_pred_min = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5748 printf(
"Unable to allocate space for taua_pred_min.\n");
5751 if ((taua_pred_max = (
float *) calloc(nwave,
sizeof (
float))) ==
NULL) {
5752 printf(
"Unable to allocate space for taua_pred_max.\n");
5766 geom.airmass = (
float *) malloc(
sizeof (
float));
5767 *
geom.airmass = 1. /
geom.csolz[0] + 1. /
geom.csenz[0];
5769 geom.airmass_sph = (
float *) malloc(
sizeof (
float));
5770 geom.airmass_plp = (
float *) malloc(
sizeof (
float));
5778 geom.senz = &
l1rec->geom_per_band->senz[ipb];
5779 geom.solz = &
l1rec->geom_per_band->solz[ipb];
5780 geom.phi = &
l1rec->geom_per_band->delphi[ipb];
5781 geom.csolz = &
l1rec->geom_per_band->csolz[ipb];
5782 geom.csenz = &
l1rec->geom_per_band->csenz[ipb];
5785 geom.airmass = (
float *) malloc(
Nbands *
sizeof (
float));
5786 for (iw = 0; iw <
Nbands; iw++) {
5787 geom.airmass[iw] = 1. /
geom.csolz[iw] + 1. /
geom.csenz[iw];
5790 geom.airmass_plp = (
float *) malloc(
Nbands *
sizeof (
float));
5791 geom.airmass_sph = (
float *) malloc(
Nbands *
sizeof (
float));
5792 for (iw = 0; iw <
Nbands; iw++) {
5803 aer_opt = aer_opt_in;
5806 for (iw = 0; iw < nwave; iw++) {
5807 F0 [iw] = F0_in [iw];
5808 taur[iw] =
l1rec->l1file->Tau_r[iw];
5809 La1 [iw] = La1_in[iw];
5810 if (iw == iwnir_s_in) iwnir_s = iw;
5811 if (iw == iwnir_l_in) iwnir_l = iw;
5815 mu0 = cos(
solz / radeg);
5816 mu = cos(
senz / radeg);
5817 airmass = 1.0 / mu0 + 1.0 / mu;
5831 for (iw = 0; iw < nwave; iw++) {
5836 radref[iw] =
pi /
F0[iw] /
geom.csolz[iw *
geom.gmult];
5840 ipb = ip * nwave + iw;
5851 for (
im = 0;
im < aertab->nmodel;
im++) mindx[
im] =
im;
5852 if (have_rh && aertab->nmodel >= 30) {
5853 printf(
"\nLimiting aerosol models based on RH.\n");
5860 if ((interpol == 1) && (
geom.gmult == 1)) {
5861 fprintf(
stderr,
"-E- %s line %d: Interpolated aerosol tables are\n",
5862 __FILE__, __LINE__);
5863 fprintf(
stderr,
" not permitted for use with band-dependent geometry, set geom_per_band=0\n");
5883 if (aer_opt > 0 && aer_opt <=
MAXAERMOD) {
5884 printf(
"\nUsing fixed aerosol model #%d (%s)\n", aer_opt,
input->aermodels[aer_opt - 1]);
5885 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5889 printf(
"\nUsing Spectral Matching of aerosols reflectance for\n");
5890 printf(
"wavelength from %4.1f nm to %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
5893 printf(
"\nUsing white-aerosol approximation\n");
5894 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5898 printf(
"\nUsing Gordon & Wang aerosol model selection\n");
5899 printf(
"Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
5900 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5904 printf(
"\nUsing Gordon & Wang aerosol model selection\n");
5905 printf(
" and NIR correction with up to %d iterations\n",
input->aer_iter_max);
5906 printf(
"Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
5907 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5911 printf(
"\nUsing Gordon & Wang aerosol model selection with NIR/SWIR switching.\n");
5912 printf(
"NIR correction with up to %d iterations\n",
input->aer_iter_max);
5913 printf(
"NIR bands at %d and %d nm\n",
input->aer_wave_short,
input->aer_wave_long);
5914 printf(
"SWIR bands at %d and %d nm\n\n",
input->aer_swir_short,
input->aer_swir_long);
5918 printf(
"\nUsing Gordon & Wang aerosol model selection\n");
5919 printf(
" and MUMM correction\n");
5920 printf(
"Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
5921 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5924 printf(
"\nUsing multi-scattering aerosol model selection.\n");
5925 printf(
" and NIR correction with up to %d iterations\n",
input->aer_iter_max);
5926 printf(
"Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
5927 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5930 printf(
"\nUsing multi-scattering aerosol model selection in linear space.\n");
5931 printf(
" and NIR correction with up to %d iterations\n",
input->aer_iter_max);
5932 printf(
"Using bands at %4.1f and %4.1f nm for model selection\n", wave[iwnir_s], wave[iwnir_l]);
5933 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5936 printf(
"\nUsing fixed, input aerosol optical thicknesses for aerosol selection.\n");
5939 printf(
"\nUsing multiple scattering aerosols with fixed model pair\n");
5942 printf(
"\nUsing multiple scattering aerosols with fixed model pair\n");
5943 printf(
" and NIR iteration with up to %d iterations\n",
input->aer_iter_max);
5944 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5947 printf(
"\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
5948 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5951 printf(
"\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
5952 printf(
" and NIR iteration with up to %d iterations\n",
input->aer_iter_max);
5953 printf(
"Extrapolating from %4.1f nm\n", wave[iwnir_l]);
5956 printf(
"\nErroneous aerosol option, %d\n", aer_opt);
5966 for (iw = iwnir_s; iw <= iwnir_l; iw +=
MAX(iwnir_l - iwnir_s, 1))
5967 rhoa[iw] = La1[iw] * radref[iw];
5977 if (iwnir_l <= iwnir_s || wave[iwnir_s] < 600) {
5978 printf(
"Aerosol selection bands must be greater than 600nm with short wave less than long wave (%d,%d)\n", iwnir_l, iwnir_s);
5986 for (iw = iwnir_s; iw <= iwnir_l; iw++) {
5987 rhoa[iw] = La1[iw] * radref[iw];
5991 if (rhoa[iwnir_s] >
input->rhoamin && rhoa[iwnir_l] >
input->rhoamin) {
5994 if (La1[iwnir_s] / La1[iwnir_l] > 0.1) {
5997 aertab->nmodel, mindx,
5998 &
geom, wv, rhoa, modmin, modmax, modrat,
eps, taua, taua);
6001 for (iw = 0; iw < nwave; iw++)
6002 La2[iw] = rhoa[iw] / radref[iw];
6005 }
else if (rhoa[iwnir_s] > -(
input->rhoamin) && rhoa[iwnir_l] > -(
input->rhoamin)) {
6009 *modmin = aertab->nmodel;
6010 *modmax = aertab->nmodel;
6012 temp =
MAX(rhoa[iwnir_l], 1e-6);
6013 for (iw = 0; iw < nwave; iw++) {
6015 La2 [iw] = rhoa[iw] / radref[iw];
6036 if (iwnir_l <= iwnir_s || wave[iwnir_s] < 600) {
6037 printf(
"Aerosol selection bands must be greater than 600nm with short wave less than long wave");
6045 for (iw = iwnir_s; iw <= iwnir_l; iw++) {
6046 rhoa[iw] = La1[iw] * radref[iw];
6050 if (rhoa[iwnir_s] >
input->rhoamin && rhoa[iwnir_l] >
input->rhoamin) {
6053 eps_tmp=rhoa[iwnir_s] / rhoa[iwnir_l];
6055 eps_tmp=rhoa[iwnir_s] / rhoa[aer_wave_base];
6056 if (eps_tmp > 0.1 && eps_tmp < 10.0) {
6059 &
geom, wv, rh,
pr, taur, rhoa,
6060 modmin, modmax, modrat, modmin2, modmax2, modrat2,
eps, taua, t_sol, t_sen, ip,mbac_w,&chi2,uncertainty);
6063 l2rec->chi2[ip]=chi2;
6064 for (iw = 0; iw < nwave; iw++){
6065 La2[iw] = rhoa[iw] / radref[iw];
6068 uncertainty->derv_La_taua_l[iw] /= radref[iw];
6069 uncertainty->derv_La_rhow_l[iw] /= radref[iw];
6070 uncertainty->derv_La_rh[iw] /= radref[iw];
6071 for(inir=0;inir<=iwnir_l-iwnir_s;inir++){
6072 uncertainty->derv_La_rhorc[iw][inir] /= radref[iw];
6079 }
else if (rhoa[iwnir_s] > -(
input->rhoamin) && rhoa[iwnir_l] > -(
input->rhoamin)) {
6083 *modmin = aertab->nmodel;
6084 *modmax = aertab->nmodel;
6086 *modmin2 = aertab->nmodel;
6087 *modmax2 = aertab->nmodel;
6089 temp =
MAX(rhoa[iwnir_l], 1e-6);
6090 for (iw = 0; iw < nwave; iw++) {
6092 La2 [iw] = rhoa[iw] / radref[iw];
6106 *modmin, *modmax, *modrat, rhoa, taua, t_sol, t_sen, taua_pred_min, taua_pred_max, taua_opt, ip, uncertainty);
6124 if (La1[iwnir_l] > 0.0) {
6131 for (iw = 0; iw < nwave; iw++) {
6132 La2 [iw] = *
eps *
F0[iw] /
F0[iwnir_l] * La1[iwnir_l];
6133 rhoa[iw] = La2[iw] * radref[iw];
6146 if (aerec !=
NULL && aerec->mode ==
ON) {
6147 *modmin = aerec->mod_min[ip] - 1;
6148 *modmax = aerec->mod_max[ip] - 1;
6149 *modrat = aerec->mod_rat[ip];
6151 *modmin =
input->aermodmin - 1;
6152 *modmax =
input->aermodmax - 1;
6153 *modrat =
input->aermodrat;
6157 *modmin, *modmax, *modrat, rhoa,
eps);
6159 for (iw = 0; iw < nwave; iw++)
6160 La2[iw] = rhoa[iw] / radref[iw];
6167 if (
input->aer_angstrom > -2) {
6168 angstrom =
input->aer_angstrom;
6172 unix2yds(l2rec->l1rec->scantime, &year, &
day, &sec);
6176 if (angstrom > -2) {
6181 *modmin, *modmax, *modrat, rhoa,
eps);
6186 for (iw = 0; iw < nwave; iw++)
6187 La2[iw] = rhoa[iw] / radref[iw];
6196 if (aerec !=
NULL && aerec->mode ==
ON)
6197 aot = &aerec->taua[ip *
Nbands];
6202 modmin, modmax, modrat, rhoa,
eps);
6204 for (iw = 0; iw < nwave; iw++)
6205 La2[iw] = rhoa[iw] / radref[iw];
6215 *modmin = aer_opt - 1;
6216 *modmax = aer_opt - 1;
6219 if (*modmin < 0 || *modmin >
input->naermodels - 1) {
6220 printf(
"Invalid aerosol option: %d\n", *modmin);
6225 rhoa[iwnir_l] = La1[iwnir_l] * radref[iwnir_l];
6234 for (iw = 0; iw < nwave; iw++)
6235 La2[iw] = rhoa[iw] / radref[iw];
6243 *modmin, *modmax, *modrat, rhoa, taua, t_sol, t_sen, taua_pred_min, taua_pred_max, 0,ip, uncertainty);
6247 for (iw = 0; iw < nwave; iw++) {
6248 La2_out [iw] = La2 [iw];
6249 taua_out [iw] = taua [iw];
6250 t_sol_out[iw] = t_sol[iw];
6251 t_sen_out[iw] = t_sen[iw];
6255 *modmin = *modmin + 1;
6256 *modmax = *modmax + 1;
6257 *modmin2 = *modmin2 + 1;
6258 *modmax2 = *modmax2 + 1;
6269 free(taua_pred_min);
6270 free(taua_pred_max);
6273 free(
geom.airmass_plp);
6274 free(
geom.airmass_sph);
6297 static int firstCall = 1;
6307 l1str *
l1rec = l2rec->l1rec;
6312 wave2 =
l1file->fwave[ib2];
6320 wave1 =
l1file->fwave[ib1];
6322 for (ip = 0; ip <
l1rec->npix; ip++) {
6323 aot1 = l2rec->taua[ip *
l1file->nbands + ib1];
6324 aot2 = l2rec->taua[ip *
l1file->nbands + ib2];
6325 if (aot1 > 0.0 && aot2 > 0.0)
6326 angst[ip] = -log(aot1 / aot2) / log(wave1 / wave2);
6327 else if (aot1 == 0.0 && aot2 == 0.0)
6345 int32_t ip, ipb1, ipb2;
6347 l1str *
l1rec = l2rec->l1rec;
6350 for (ip = 0; ip <
l1rec->npix; ip++) {
6351 ipb1 = ip *
l1file->nbands + iwnir_s;
6352 ipb2 = ip *
l1file->nbands + iwnir_l;
6353 if (l2rec->La[ipb2] > 0.0) {
6354 eps[ip] = l2rec->La[ipb1] / l2rec->La[ipb2]