src/gen/genssky.c

#ifndef lint
static const char RCSid[] = "$Id: genssky.c,v 2.8 2025/06/07 05:09:45 greg Exp $";
#endif
/* Main function for generating spectral sky */
/* Cloudy sky computed as weight average of clear and cie overcast sky */

#include "atmos.h"
#include "copyright.h"
#include "color.h"
#include "paths.h"
#include "resolu.h"
#include "rtio.h"
#include <ctype.h>
#ifdef _WIN32
#include <windows.h>
#else
#include <errno.h>
#include <sys/stat.h>
#include <sys/types.h>
#endif

const double ARCTIC_LAT = 67.;
const double TROPIC_LAT = 23.;
const int SUMMER_START = 4;
const int SUMMER_END = 9;
const double GNORM = 0.777778;

const double D65EFF = 203.; /* standard illuminant D65 */

/* Mean normalized relative daylight spectra where CCT = 6415K for overcast; */
const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
                                                          1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
                                                          1.0434,  1.05547, 0.98212, 0.94445, 0.9722,
                                                          0.82387, 0.87853, 0.82559, 0.75111, 0.78925};

/* European and North American zones */
struct
{
        char zname[8];                 /* time zone name (all caps) */
        float zmer;                 /* standard meridian */
} tzone[] = {{"YST", 135},       {"YDT", 120},   {"PST", 120},  {"PDT", 105},
                         {"MST", 105},   {"MDT", 90},    {"CST", 90},   {"CDT", 75},
                         {"EST", 75},    {"EDT", 60},    {"AST", 60},   {"ADT", 45},
                         {"NST", 52.5},  {"NDT", 37.5},  {"GMT", 0},    {"BST", -15},
                         {"CET", -15},   {"CEST", -30},  {"EET", -30},  {"EEST", -45},
                         {"AST", -45},   {"ADT", -60},   {"GST", -60},  {"GDT", -75},
                         {"IST", -82.5}, {"IDT", -97.5}, {"JST", -135}, {"NDT", -150},
                         {"NZST", -180}, {"NZDT", -195}, {"", 0}};

static int
make_directory
(
        const char *path
)
{
#ifdef _WIN32
        if (CreateDirectory(path, NULL) || GetLastError() == ERROR_ALREADY_EXISTS) {
                return 1;
        }
        return 0;

#else
        if (mkdir(path, 0777) == 0 || errno == EEXIST) {
                return 1;
        }
        return 0;

#endif
}

inline static float
deg2rad
(
        float deg
)
{
        return deg * (PI / 180.);
}

static int
cvthour
(
        char *hs,
        int *tsolar,
        double *hour
)
{
        char *cp = hs;
        int i, j;

        if ((*tsolar = *cp == '+')) {
                cp++;                                 /* solar time? */
        }
        while (isdigit(*cp)) {
                cp++;
        }
        if (*cp == ':') {
                *hour = atoi(hs) + atoi(++cp) / 60.0;
        }else{
                *hour = atof(hs);
                if (*cp == '.') {
                        cp++;
                }
        }
        while (isdigit(*cp)) {
                cp++;
        }
        if (!*cp) {
                return (0);
        }
        if (*tsolar || !isalpha(*cp)) {
                fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
                exit(1);
        }
        i = 0;
        do {
                for (j = 0; cp[j]; j++) {
                        if (toupper(cp[j]) != tzone[i].zname[j]) {
                                break;
                        }
                }
                if (!cp[j] && !tzone[i].zname[j]) {
                        s_meridian = tzone[i].zmer * (PI / 180);
                        return (1);
                }
        } while (tzone[i++].zname[0]);

        fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
        fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
        for (i = 1; tzone[i].zname[0]; i++) {
                fprintf(stderr, " %s", tzone[i].zname);
        }
        putc('\n', stderr);
        exit(1);
}

static void
basename
(
        const char *path,
        char *output,
        size_t outsize
)
{
        const char *last_slash = strrchr(path, '/');
        const char *last_backslash = strrchr(path, '\\');
        const char *filename = path;
        const char *last_dot;

        if (last_slash && last_backslash) {
                filename =
                        (last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
        } else if (last_slash) {
                filename = last_slash + 1;
        } else if (last_backslash) {
                filename = last_backslash + 1;
        }

        last_dot = strrchr(filename, '.');
        if (last_dot) {
                size_t length = last_dot - filename;
                if (length < outsize) {
                        strncpy(output, filename, length);
                        output[length] = '\0';
                } else {
                        strncpy(output, filename, outsize - 1);
                        output[outsize - 1] = '\0';
                }
        }
}

static char *
join_paths
(
        const char *path1,
        const char *path2
)
{
        size_t len1 = strlen(path1);
        size_t len2 = strlen(path2);
        int need_separator = (path1[len1 - 1] != DIRSEP);

        char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
        if (!result) {
                return NULL;
        }

        strcpy(result, path1);
        if (need_separator) {
                result[len1] = DIRSEP;
                len1++;
        }
        strcpy(result + len1, path2);

        return result;
}

static inline double
wmean2
(
        const double a,
        const double b,
        const double x
)
{
        return a * (1 - x) + b * x;
}

static inline double
wmean
(
        const double a,
        const double x,
        const double b,
        const double y
)
{
        return (a * x + b * y) / (a + b);
}

static double
get_overcast_zenith_brightness
(
        const double sundir[3]
)
{
        double zenithbr;
        if (sundir[2] < 0) {
                zenithbr = 0;
        } else {
                zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
        }
        return zenithbr;
}

/* from gensky.c */
static double
get_overcast_brightness
(
        const double dz,
        const double zenithbr
)
{
        double groundbr = zenithbr * GNORM;
        return wmean(
                pow(dz + 1.01, 10),
                zenithbr * (1 + 2 * dz) / 3,
                pow(dz + 1.01, -10),
                groundbr);
}

static void
write_header
(
        const int argc,
        char **argv,
        const double cloud_cover,
        const double grefl,
        const int res
)
{
        int i;
        printf("# ");
        for (i = 0; i < argc; i++) {
                printf("%s ", argv[i]);
        }
        printf("\n");
        printf(
                "#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: "
                "%d\n\n",
                cloud_cover,
                grefl,
                res);
}

static void
write_rad
(
        const double *sun_radiance,
        const double intensity,
        const FVECT sundir,
        const char *ddir,
        const char *skyfile
)
{
        if (sundir[2] > 0) {
                printf("void spectrum sunrad\n0\n0\n22 380 780 ");
                int i;
                for (i = 0; i < NSSAMP; ++i) {
                        printf("%.3f ", sun_radiance[i]);
                }
                printf(
                        "\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n",
                        intensity,
                        intensity,
                        intensity);
                printf(
                        "solar source sun\n0\n0\n4 %f %f %f 0.533\n\n",
                        sundir[0],
                        sundir[1],
                        sundir[2]);
        }
        printf(
                "void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
                "'1-Acos(Dz)/PI'\n0\n0\n\n",
                skyfile);
}

static void
write_hsr_header
(
        FILE *fp,
        RESOLU *res
)
{
        newheader("RADIANCE", fp);
        fputncomp(NSSAMP, fp);
        fputwlsplit(WLPART, fp);
        fputformat(SPECFMT, fp);
        fputc('\n', fp);
        fputsresolu(res, fp);
}

static void
reverse_array_float
(
        float arr[],
        int size
)
{
        int start = 0;
        int end = size - 1;

        while (start < end) {
                float temp = arr[start];
                arr[start] = arr[end];
                arr[end] = temp;
                start++;
                end--;
        }
}

int
gen_spect_sky
(
        DATARRAY *tau_clear,
        DATARRAY *scat_clear,
        DATARRAY *scat1m_clear,
        DATARRAY *irrad_clear,
        const double cloud_cover,
        const FVECT sundir,
        const double grefl,
        const int res,
        const char *outname,
        const char *ddir,
        const double dirnorm,
        const double difhor
)
{
        char skyfile[PATH_MAX];
        if (!snprintf(
                        skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP, outname)) {
                fprintf(stderr, "Error setting sky file name\n");
                return 0;
        }
        ;
        int xres = res;
        int yres = xres / 2;
        RESOLU rs = {PIXSTANDARD, xres, yres};
        FILE *skyfp = fopen(skyfile, "w");
        write_hsr_header(skyfp, &rs);

        CNDX[3] = NSSAMP;

        FVECT view_point = {0, 0, ER + 10};
        const double radius = VLEN(view_point);
        const double sun_ct = fdot(view_point, sundir) / radius;

        double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
        double overcast_grndbr = overcast_zenithbr * GNORM;

        double dif_ratio = 1;
        if (difhor > 0) {
                DATARRAY *indirect_irradiance_clear =
                        get_indirect_irradiance(irrad_clear, radius, sun_ct);
                double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
                double diffuse_irradiance = 0;
                int l;
                for (l = 0; l < NSSAMP; ++l) {
                        diffuse_irradiance +=
                                indirect_irradiance_clear->arr.d[l] * 20;                                                                 /* 20nm interval */
                }
                free(indirect_irradiance_clear);
                diffuse_irradiance =
                        wmean2(diffuse_irradiance, overcast_ghi, cloud_cover);
                if (diffuse_irradiance > 0) {
                        dif_ratio =
                                difhor / WHTEFFICACY / diffuse_irradiance / 1.15;                                                                 /* fudge */
                }
        }
        int i, j, k;
        for (j = 0; j < yres; ++j) {
                for (i = 0; i < xres; ++i) {
                        SCOLOR radiance = {0};
                        SCOLR sky_sclr = {0};

                        float px = i / (xres - 1.0);
                        float py = j / (yres - 1.0);
                        float lambda = ((1 - py) * PI) - (PI / 2.0);
                        float phi = (px * 2.0 * PI) - PI;

                        FVECT rdir = {
                                cos(lambda) * cos(phi), cos(lambda) * sin(phi), sin(lambda)
                        };

                        const double mu = fdot(view_point, rdir) / radius;
                        const double nu = fdot(rdir, sundir);

                        /* hit ground */
                        if (rdir[2] < 0) {
                                get_ground_radiance(
                                        tau_clear,
                                        scat_clear,
                                        scat1m_clear,
                                        irrad_clear,
                                        view_point,
                                        rdir,
                                        radius,
                                        mu,
                                        sun_ct,
                                        nu,
                                        grefl,
                                        sundir,
                                        radiance);
                        } else {
                                get_sky_radiance(
                                        scat_clear, scat1m_clear, radius, mu, sun_ct, nu, radiance);
                        }

                        for (k = 0; k < NSSAMP; ++k) {
                                radiance[k] *= WVLSPAN;
                        }

                        if (cloud_cover > 0) {
                                double skybr =
                                        get_overcast_brightness(rdir[2], overcast_zenithbr);
                                if (rdir[2] < 0) {
                                        for (k = 0; k < NSSAMP; ++k) {
                                                radiance[k] = wmean2(
                                                        radiance[k],
                                                        overcast_grndbr * D6415[k],
                                                        cloud_cover);
                                        }
                                } else {
                                        for (k = 0; k < NSSAMP; ++k) {
                                                radiance[k] =
                                                        wmean2(radiance[k], skybr * D6415[k], cloud_cover);
                                        }
                                }
                        }

                        for (k = 0; k < NSSAMP; ++k) {
                                radiance[k] *= dif_ratio;
                        }

                        reverse_array_float(radiance, NSSAMP);

                        scolor2scolr(sky_sclr, radiance, NSSAMP);
                        putbinary(sky_sclr, LSCOLR, 1, skyfp);
                }
        }
        fclose(skyfp);

        /* Get solar radiance */
        double sun_radiance[NSSAMP] = {0};
        get_solar_radiance(
                tau_clear,
                scat_clear,
                scat1m_clear,
                sundir,
                radius,
                sun_ct,
                sun_radiance);
        if (cloud_cover > 0) {
                double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
                int i;
                for (i = 0; i < NSSAMP; ++i) {
                        sun_radiance[i] = wmean2(
                                sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
                }
        }

        /* Normalize */
        double sum = 0.0;
        for (i = 0; i < NSSAMP; ++i) {
                sum += sun_radiance[i];
        }
        double mean = sum / NSSAMP;
        for (i = 0; i < NSSAMP; ++i) {
                sun_radiance[i] /= mean;
        }
        double intensity = mean * WVLSPAN;
        if (dirnorm > 0) {
                intensity = dirnorm / SOLOMG / WHTEFFICACY;
        }

        write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
        return 1;
}

static DpPaths
get_dppaths
(
        const char *dir,
        const double aod,
        const char *mname,
        const char *tag
)
{
        DpPaths paths;

        snprintf(
                paths.tau,
                PATH_MAX,
                "%s%ctau_%s_%s_%.2f.dat",
                dir,
                DIRSEP,
                tag,
                mname,
                aod);
        snprintf(
                paths.scat,
                PATH_MAX,
                "%s%cscat_%s_%s_%.2f.dat",
                dir,
                DIRSEP,
                tag,
                mname,
                aod);
        snprintf(
                paths.scat1m,
                PATH_MAX,
                "%s%cscat1m_%s_%s_%.2f.dat",
                dir,
                DIRSEP,
                tag,
                mname,
                aod);
        snprintf(
                paths.irrad,
                PATH_MAX,
                "%s%cirrad_%s_%s_%.2f.dat",
                dir,
                DIRSEP,
                tag,
                mname,
                aod);

        return paths;
}

static void
set_rayleigh_density_profile
(
        Atmosphere *atmos,
        char *tag,
        const int is_summer,
        const double s_latitude
)
{
        if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
                tag[0] = 's';
                if (is_summer) {
                        tag[1] = 's';
                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
                        atmos->beta_r0 = BR0_SS;
                } else {
                        tag[1] = 'w';
                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
                        atmos->beta_r0 = BR0_SW;
                }
        } else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
                tag[0] = 'm';
                if (is_summer) {
                        tag[1] = 's';
                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
                        atmos->beta_r0 = BR0_MS;
                } else {
                        tag[1] = 'w';
                        atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
                        atmos->beta_r0 = BR0_MW;
                }
        } else {
                tag[0] = 't';
                tag[1] = 'r';
                atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
                atmos->beta_r0 = BR0_T;
        }
        tag[2] = '\0';
}

static Atmosphere
init_atmos
(
        const double aod,
        const double grefl
)
{
        Atmosphere atmos = {
                .ozone_density =
                {.layers =
                 {
                         {.width = 25000.0,
                                          .exp_term = 0.0,
                                          .exp_scale = 0.0,
                                          .linear_term = 1.0 / 15000.0,
                                          .constant_term = -2.0 / 3.0},
                         {.width = AH,
                                          .exp_term = 0.0,
                                          .exp_scale = 0.0,
                                          .linear_term = -1.0 / 15000.0,
                                          .constant_term = 8.0 / 3.0},
                 }},
                .rayleigh_density =
                {.layers =
                 {
                         {.width = AH,
                                          .exp_term = 1.0,
                                          .exp_scale = -1.0 / HR_MS,
                                          .linear_term = 0.0,
                                          .constant_term = 0.0},
                 }},
                .beta_r0 = BR0_MS,
                .beta_scale = aod / AOD0_CA,
                .beta_m = NULL,
                .grefl = grefl
        };
        return atmos;
}

int
main
(
        int argc,
        char *argv[]
)
{
        int month, day;
        double hour;
        FVECT sundir;
        int num_threads = 1;
        int sorder = 4;
        int year = 0;
        int tsolar = 0;
        int got_meridian = 0;
        double grefl = 0.2;
        double ccover = 0.0;
        int res = 64;
        double aod = AOD0_CA;
        char *outname = "out";
        char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
        char mie_name[20] = "mie_ca";
        char lstag[3];
        char *ddir = ".";
        int i;
        double dirnorm = 0;                 /* direct normal illuminance */
        double difhor = 0;                 /* diffuse horizontal illuminance */

        fixargv0(argv[0]);
        if (argc == 2 && !strcmp(argv[1], "-defaults")) {
                printf("-i %d\t\t\t\t#scattering order\n", sorder);
                printf("-g %f\t\t\t#ground reflectance\n", grefl);
                printf("-c %f\t\t\t#cloud cover\n", ccover);
                printf("-r %d\t\t\t\t#image resolution\n", res);
                printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
                printf("-f %s\t\t\t\t#output name (-f)\n", outname);
                printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
                exit(0);
        }

        if (argc < 4) {
                fprintf(
                        stderr,
                        "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod "
                        "-r "
                        "res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
                        "-g grefl -f outpath\n",
                        argv[0]);
                return 0;
        }

        month = atoi(argv[1]);
        if (month < 1 || month > 12) {
                fprintf(stderr, "bad month");
                exit(1);
        }
        day = atoi(argv[2]);
        if (day < 1 || day > 31) {
                fprintf(stderr, "bad month");
                exit(1);
        }
        got_meridian = cvthour(argv[3], &tsolar, &hour);

        if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
                fprintf(stderr, "Cannot compute solar angle\n");
                exit(1);
        }

        for (i = 4; i < argc; i++) {
                if (argv[i][0] == '-') {
                        switch (argv[i][1]) {
                        case 'a':
                                s_latitude = atof(argv[++i]) * (PI / 180.0);
                                break;
                        case 'c':
                                ccover = atof(argv[++i]);
                                break;
                        case 'd':
                                aod = atof(argv[++i]);
                                break;
                        case 'f':
                                outname = argv[++i];
                                break;
                        case 'g':
                                grefl = atof(argv[++i]);
                                break;
                        case 'i':
                                sorder = atoi(argv[++i]);
                                break;
                        case 'l':
                                mie_path = argv[++i];
                                basename(mie_path, mie_name, sizeof(mie_name));
                                break;
                        case 'm':
                                if (got_meridian) {
                                        ++i;
                                        break;
                                }
                                s_meridian = atof(argv[++i]) * (PI / 180.0);
                                break;
                        case 'n':
                                num_threads = atoi(argv[++i]);
                                break;
                        case 'o':
                                s_longitude = atof(argv[++i]) * (PI / 180.0);
                                break;
                        case 'L':
                                dirnorm = atof(argv[++i]);
                                difhor = atof(argv[++i]);
                                break;
                        case 'p':
                                ddir = argv[++i];
                                break;
                        case 'r':
                                res = atoi(argv[++i]);
                                break;
                        case 'y':
                                year = atoi(argv[++i]);
                                break;
                        default:
                                fprintf(stderr, "Unknown option %s\n", argv[i]);
                                exit(1);
                        }
                }
        }
        if (year && (year < 1950) | (year > 2050)) {
                fprintf(
                        stderr,
                        "%s: warning - year should be in range 1950-2050\n",
                        progname);
        }
        if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180) {
                fprintf(
                        stderr,
                        "%s: warning - %.1f hours btwn. standard meridian and "
                        "longitude\n",
                        progname,
                        (s_longitude - s_meridian) * 12 / PI);
        }

        Atmosphere clear_atmos = init_atmos(aod, grefl);

        int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
        if (s_latitude < 0) {
                is_summer = !is_summer;
        }
        set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);

        /* Load mie density data */
        DATARRAY *mie_dp = getdata(mie_path);
        if (mie_dp == NULL) {
                fprintf(stderr, "Error reading mie data\n");
                return 0;
        }
        clear_atmos.beta_m = mie_dp;

        char gsdir[PATH_MAX];
        size_t siz = strlen(ddir);
        if (ISDIRSEP(ddir[siz - 1])) {
                ddir[siz - 1] = '\0';
        }
        snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
        if (!make_directory(gsdir)) {
                fprintf(stderr, "Failed creating atmos_data directory");
                exit(1);
        }
        DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);

        if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
                getpath(clear_paths.scat, ".", R_OK) == NULL ||
                getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
                getpath(clear_paths.irrad, ".", R_OK) == NULL) {
                printf("# Pre-computing...\n");
                if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
                        fprintf(stderr, "Pre-compute failed\n");
                        return 0;
                }
        }

        DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
        DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
        DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
        DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);

        write_header(argc, argv, ccover, grefl, res);

        if (!gen_spect_sky(
                        tau_clear_dp,
                        scat_clear_dp,
                        scat1m_clear_dp,
                        irrad_clear_dp,
                        ccover,
                        sundir,
                        grefl,
                        res,
                        outname,
                        ddir,
                        dirnorm,
                        difhor)) {
                fprintf(stderr, "gen_spect_sky failed\n");
                exit(1);
        }

        freedata(mie_dp);
        freedata(tau_clear_dp);
        freedata(scat_clear_dp);
        freedata(irrad_clear_dp);
        freedata(scat1m_clear_dp);

        return 1;
}
Revision:	2.9
Committed:	Fri Jul 11 18:12:25 2025 UTC (3 weeks, 2 days ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad6R0, HEAD
Changes since 2.8:	+811 -574 lines
Log Message:	fix(genssky,gensdaymtx): Spectral data ordering was reversed in output HSR's (TW)
#	User	Rev	Content
1	greg	2.2	#ifndef lint
2	greg	2.9	static const char RCSid[] = "$Id: genssky.c,v 2.8 2025/06/07 05:09:45 greg Exp $";
3	greg	2.2	#endif
4			/* Main function for generating spectral sky */
5			/* Cloudy sky computed as weight average of clear and cie overcast sky */
6	greg	2.1
7	greg	2.9	#include "atmos.h"
8	greg	2.8	#include "copyright.h"
9	greg	2.9	#include "color.h"
10	greg	2.8	#include "paths.h"
11	greg	2.1	#include "resolu.h"
12	greg	2.2	#include "rtio.h"
13			#include <ctype.h>
14			#ifdef _WIN32
15			#include <windows.h>
16			#else
17			#include <errno.h>
18			#include <sys/stat.h>
19			#include <sys/types.h>
20			#endif
21	greg	2.1
22			const double ARCTIC_LAT = 67.;
23			const double TROPIC_LAT = 23.;
24			const int SUMMER_START = 4;
25			const int SUMMER_END = 9;
26			const double GNORM = 0.777778;
27
28			const double D65EFF = 203.; /* standard illuminant D65 */
29
30	greg	2.2	/* Mean normalized relative daylight spectra where CCT = 6415K for overcast; */
31	greg	2.1	const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
32	greg	2.9	1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
33			1.0434, 1.05547, 0.98212, 0.94445, 0.9722,
34			0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
35	greg	2.1
36	greg	2.2	/* European and North American zones */
37	greg	2.9	struct
38			{
39			char zname[8]; /* time zone name (all caps) */
40			float zmer; /* standard meridian */
41			} tzone[] = {{"YST", 135}, {"YDT", 120}, {"PST", 120}, {"PDT", 105},
42			{"MST", 105}, {"MDT", 90}, {"CST", 90}, {"CDT", 75},
43			{"EST", 75}, {"EDT", 60}, {"AST", 60}, {"ADT", 45},
44			{"NST", 52.5}, {"NDT", 37.5}, {"GMT", 0}, {"BST", -15},
45			{"CET", -15}, {"CEST", -30}, {"EET", -30}, {"EEST", -45},
46			{"AST", -45}, {"ADT", -60}, {"GST", -60}, {"GDT", -75},
47			{"IST", -82.5}, {"IDT", -97.5}, {"JST", -135}, {"NDT", -150},
48			{"NZST", -180}, {"NZDT", -195}, {"", 0}};
49
50			static int
51			make_directory
52			(
53			const char *path
54			)
55			{
56			#ifdef _WIN32
57			if (CreateDirectory(path, NULL) \|\| GetLastError() == ERROR_ALREADY_EXISTS) {
58			return 1;
59			}
60			return 0;
61	greg	2.2
62			#else
63	greg	2.9	if (mkdir(path, 0777) == 0 \|\| errno == EEXIST) {
64			return 1;
65			}
66			return 0;
67
68	greg	2.2	#endif
69			}
70
71	greg	2.9	inline static float
72			deg2rad
73			(
74			float deg
75			)
76			{
77			return deg * (PI / 180.);
78			}
79
80			static int
81			cvthour
82			(
83			char *hs,
84			int *tsolar,
85			double *hour
86			)
87			{
88			char *cp = hs;
89			int i, j;
90
91			if ((tsolar = cp == '+')) {
92			cp++; /* solar time? */
93			}
94			while (isdigit(*cp)) {
95			cp++;
96			}
97			if (*cp == ':') {
98			*hour = atoi(hs) + atoi(++cp) / 60.0;
99			}else{
100			*hour = atof(hs);
101			if (*cp == '.') {
102			cp++;
103			}
104			}
105			while (isdigit(*cp)) {
106			cp++;
107			}
108			if (!*cp) {
109			return (0);
110			}
111			if (tsolar \|\| !isalpha(cp)) {
112			fprintf(stderr, "%s: bad time format: %s\n", progname, hs);
113			exit(1);
114			}
115			i = 0;
116			do {
117			for (j = 0; cp[j]; j++) {
118			if (toupper(cp[j]) != tzone[i].zname[j]) {
119			break;
120			}
121			}
122			if (!cp[j] && !tzone[i].zname[j]) {
123			s_meridian = tzone[i].zmer * (PI / 180);
124			return (1);
125			}
126			} while (tzone[i++].zname[0]);
127
128			fprintf(stderr, "%s: unknown time zone: %s\n", progname, cp);
129			fprintf(stderr, "Known time zones:\n\t%s", tzone[0].zname);
130			for (i = 1; tzone[i].zname[0]; i++) {
131			fprintf(stderr, " %s", tzone[i].zname);
132			}
133			putc('\n', stderr);
134			exit(1);
135			}
136
137			static void
138			basename
139			(
140			const char *path,
141			char *output,
142			size_t outsize
143			)
144			{
145			const char *last_slash = strrchr(path, '/');
146			const char *last_backslash = strrchr(path, '\\');
147			const char *filename = path;
148			const char *last_dot;
149
150			if (last_slash && last_backslash) {
151			filename =
152			(last_slash > last_backslash) ? last_slash + 1 : last_backslash + 1;
153			} else if (last_slash) {
154			filename = last_slash + 1;
155			} else if (last_backslash) {
156			filename = last_backslash + 1;
157			}
158
159			last_dot = strrchr(filename, '.');
160			if (last_dot) {
161			size_t length = last_dot - filename;
162			if (length < outsize) {
163			strncpy(output, filename, length);
164			output[length] = '\0';
165			} else {
166			strncpy(output, filename, outsize - 1);
167			output[outsize - 1] = '\0';
168			}
169			}
170			}
171
172			static char *
173			join_paths
174			(
175			const char *path1,
176			const char *path2
177			)
178			{
179			size_t len1 = strlen(path1);
180			size_t len2 = strlen(path2);
181			int need_separator = (path1[len1 - 1] != DIRSEP);
182
183			char *result = malloc(len1 + len2 + (need_separator ? 2 : 1));
184			if (!result) {
185			return NULL;
186			}
187
188			strcpy(result, path1);
189			if (need_separator) {
190			result[len1] = DIRSEP;
191			len1++;
192			}
193			strcpy(result + len1, path2);
194
195			return result;
196			}
197
198			static inline double
199			wmean2
200			(
201			const double a,
202			const double b,
203			const double x
204			)
205			{
206			return a * (1 - x) + b * x;
207			}
208
209			static inline double
210			wmean
211			(
212			const double a,
213			const double x,
214			const double b,
215			const double y
216			)
217			{
218			return (a * x + b * y) / (a + b);
219			}
220
221			static double
222			get_overcast_zenith_brightness
223			(
224			const double sundir[3]
225			)
226			{
227			double zenithbr;
228			if (sundir[2] < 0) {
229			zenithbr = 0;
230			} else {
231			zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
232			}
233			return zenithbr;
234	greg	2.1	}
235
236	greg	2.2	/* from gensky.c */
237	greg	2.9	static double
238			get_overcast_brightness
239			(
240			const double dz,
241			const double zenithbr
242			)
243			{
244			double groundbr = zenithbr * GNORM;
245			return wmean(
246			pow(dz + 1.01, 10),
247			zenithbr * (1 + 2 * dz) / 3,
248			pow(dz + 1.01, -10),
249			groundbr);
250			}
251
252			static void
253			write_header
254			(
255			const int argc,
256			char **argv,
257			const double cloud_cover,
258			const double grefl,
259			const int res
260			)
261			{
262			int i;
263			printf("# ");
264			for (i = 0; i < argc; i++) {
265			printf("%s ", argv[i]);
266			}
267			printf("\n");
268			printf(
269			"#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: "
270			"%d\n\n",
271			cloud_cover,
272			grefl,
273			res);
274			}
275
276			static void
277			write_rad
278			(
279			const double *sun_radiance,
280			const double intensity,
281			const FVECT sundir,
282			const char *ddir,
283			const char *skyfile
284			)
285			{
286			if (sundir[2] > 0) {
287			printf("void spectrum sunrad\n0\n0\n22 380 780 ");
288			int i;
289			for (i = 0; i < NSSAMP; ++i) {
290			printf("%.3f ", sun_radiance[i]);
291			}
292			printf(
293			"\n\nsunrad light solar\n0\n0\n3 %.1f %.1f %.1f\n\n",
294			intensity,
295			intensity,
296			intensity);
297			printf(
298			"solar source sun\n0\n0\n4 %f %f %f 0.533\n\n",
299			sundir[0],
300			sundir[1],
301			sundir[2]);
302			}
303			printf(
304			"void specpict skyfunc\n5 noop %s . 'Atan2(Dy,Dx)/PI+1' "
305			"'1-Acos(Dz)/PI'\n0\n0\n\n",
306			skyfile);
307			}
308
309			static void
310			write_hsr_header
311			(
312			FILE *fp,
313			RESOLU *res
314			)
315			{
316			newheader("RADIANCE", fp);
317			fputncomp(NSSAMP, fp);
318			fputwlsplit(WLPART, fp);
319			fputformat(SPECFMT, fp);
320			fputc('\n', fp);
321			fputsresolu(res, fp);
322			}
323
324			static void
325			reverse_array_float
326			(
327			float arr[],
328			int size
329			)
330			{
331			int start = 0;
332			int end = size - 1;
333
334			while (start < end) {
335			float temp = arr[start];
336			arr[start] = arr[end];
337			arr[end] = temp;
338			start++;
339			end--;
340			}
341			}
342
343			int
344			gen_spect_sky
345			(
346			DATARRAY *tau_clear,
347			DATARRAY *scat_clear,
348			DATARRAY *scat1m_clear,
349			DATARRAY *irrad_clear,
350			const double cloud_cover,
351			const FVECT sundir,
352			const double grefl,
353			const int res,
354			const char *outname,
355			const char *ddir,
356			const double dirnorm,
357			const double difhor
358			)
359			{
360			char skyfile[PATH_MAX];
361			if (!snprintf(
362			skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP, outname)) {
363			fprintf(stderr, "Error setting sky file name\n");
364			return 0;
365			}
366			;
367			int xres = res;
368			int yres = xres / 2;
369			RESOLU rs = {PIXSTANDARD, xres, yres};
370			FILE *skyfp = fopen(skyfile, "w");
371			write_hsr_header(skyfp, &rs);
372
373			CNDX[3] = NSSAMP;
374
375			FVECT view_point = {0, 0, ER + 10};
376			const double radius = VLEN(view_point);
377			const double sun_ct = fdot(view_point, sundir) / radius;
378
379			double overcast_zenithbr = get_overcast_zenith_brightness(sundir);
380			double overcast_grndbr = overcast_zenithbr * GNORM;
381
382			double dif_ratio = 1;
383			if (difhor > 0) {
384			DATARRAY *indirect_irradiance_clear =
385			get_indirect_irradiance(irrad_clear, radius, sun_ct);
386			double overcast_ghi = overcast_zenithbr * 7.0 * PI / 9.0;
387			double diffuse_irradiance = 0;
388			int l;
389			for (l = 0; l < NSSAMP; ++l) {
390			diffuse_irradiance +=
391			indirect_irradiance_clear->arr.d[l] * 20; /* 20nm interval */
392			}
393			free(indirect_irradiance_clear);
394			diffuse_irradiance =
395			wmean2(diffuse_irradiance, overcast_ghi, cloud_cover);
396			if (diffuse_irradiance > 0) {
397			dif_ratio =
398			difhor / WHTEFFICACY / diffuse_irradiance / 1.15; /* fudge */
399			}
400			}
401			int i, j, k;
402			for (j = 0; j < yres; ++j) {
403			for (i = 0; i < xres; ++i) {
404			SCOLOR radiance = {0};
405			SCOLR sky_sclr = {0};
406
407			float px = i / (xres - 1.0);
408			float py = j / (yres - 1.0);
409			float lambda = ((1 - py) * PI) - (PI / 2.0);
410			float phi = (px * 2.0 * PI) - PI;
411
412			FVECT rdir = {
413			cos(lambda) * cos(phi), cos(lambda) * sin(phi), sin(lambda)
414			};
415
416			const double mu = fdot(view_point, rdir) / radius;
417			const double nu = fdot(rdir, sundir);
418
419			/* hit ground */
420			if (rdir[2] < 0) {
421			get_ground_radiance(
422			tau_clear,
423			scat_clear,
424			scat1m_clear,
425			irrad_clear,
426			view_point,
427			rdir,
428			radius,
429			mu,
430			sun_ct,
431			nu,
432			grefl,
433			sundir,
434			radiance);
435			} else {
436			get_sky_radiance(
437			scat_clear, scat1m_clear, radius, mu, sun_ct, nu, radiance);
438			}
439
440			for (k = 0; k < NSSAMP; ++k) {
441			radiance[k] *= WVLSPAN;
442			}
443
444			if (cloud_cover > 0) {
445			double skybr =
446			get_overcast_brightness(rdir[2], overcast_zenithbr);
447			if (rdir[2] < 0) {
448			for (k = 0; k < NSSAMP; ++k) {
449			radiance[k] = wmean2(
450			radiance[k],
451			overcast_grndbr * D6415[k],
452			cloud_cover);
453			}
454			} else {
455			for (k = 0; k < NSSAMP; ++k) {
456			radiance[k] =
457			wmean2(radiance[k], skybr * D6415[k], cloud_cover);
458			}
459			}
460			}
461
462			for (k = 0; k < NSSAMP; ++k) {
463			radiance[k] *= dif_ratio;
464			}
465
466			reverse_array_float(radiance, NSSAMP);
467
468			scolor2scolr(sky_sclr, radiance, NSSAMP);
469			putbinary(sky_sclr, LSCOLR, 1, skyfp);
470			}
471			}
472			fclose(skyfp);
473
474			/* Get solar radiance */
475			double sun_radiance[NSSAMP] = {0};
476			get_solar_radiance(
477			tau_clear,
478			scat_clear,
479			scat1m_clear,
480			sundir,
481			radius,
482			sun_ct,
483			sun_radiance);
484			if (cloud_cover > 0) {
485			double skybr = get_overcast_brightness(sundir[2], overcast_zenithbr);
486			int i;
487			for (i = 0; i < NSSAMP; ++i) {
488			sun_radiance[i] = wmean2(
489			sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
490			}
491			}
492
493			/* Normalize */
494			double sum = 0.0;
495			for (i = 0; i < NSSAMP; ++i) {
496			sum += sun_radiance[i];
497			}
498			double mean = sum / NSSAMP;
499			for (i = 0; i < NSSAMP; ++i) {
500			sun_radiance[i] /= mean;
501			}
502			double intensity = mean * WVLSPAN;
503			if (dirnorm > 0) {
504			intensity = dirnorm / SOLOMG / WHTEFFICACY;
505			}
506
507			write_rad(sun_radiance, intensity, sundir, ddir, skyfile);
508			return 1;
509			}
510
511			static DpPaths
512			get_dppaths
513			(
514			const char *dir,
515			const double aod,
516			const char *mname,
517			const char *tag
518			)
519			{
520			DpPaths paths;
521
522			snprintf(
523			paths.tau,
524			PATH_MAX,
525			"%s%ctau_%s_%s_%.2f.dat",
526			dir,
527			DIRSEP,
528			tag,
529			mname,
530			aod);
531			snprintf(
532			paths.scat,
533			PATH_MAX,
534			"%s%cscat_%s_%s_%.2f.dat",
535			dir,
536			DIRSEP,
537			tag,
538			mname,
539			aod);
540			snprintf(
541			paths.scat1m,
542			PATH_MAX,
543			"%s%cscat1m_%s_%s_%.2f.dat",
544			dir,
545			DIRSEP,
546			tag,
547			mname,
548			aod);
549			snprintf(
550			paths.irrad,
551			PATH_MAX,
552			"%s%cirrad_%s_%s_%.2f.dat",
553			dir,
554			DIRSEP,
555			tag,
556			mname,
557			aod);
558
559			return paths;
560			}
561
562			static void
563			set_rayleigh_density_profile
564			(
565			Atmosphere *atmos,
566			char *tag,
567			const int is_summer,
568			const double s_latitude
569			)
570			{
571			if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
572			tag[0] = 's';
573			if (is_summer) {
574			tag[1] = 's';
575			atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
576			atmos->beta_r0 = BR0_SS;
577			} else {
578			tag[1] = 'w';
579			atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
580			atmos->beta_r0 = BR0_SW;
581			}
582			} else if (fabs(s_latitude * 180.0 / PI) > TROPIC_LAT) {
583			tag[0] = 'm';
584			if (is_summer) {
585			tag[1] = 's';
586			atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
587			atmos->beta_r0 = BR0_MS;
588			} else {
589			tag[1] = 'w';
590			atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
591			atmos->beta_r0 = BR0_MW;
592			}
593			} else {
594			tag[0] = 't';
595			tag[1] = 'r';
596			atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
597			atmos->beta_r0 = BR0_T;
598			}
599			tag[2] = '\0';
600			}
601
602			static Atmosphere
603			init_atmos
604			(
605			const double aod,
606			const double grefl
607			)
608			{
609			Atmosphere atmos = {
610			.ozone_density =
611			{.layers =
612			{
613			{.width = 25000.0,
614			.exp_term = 0.0,
615			.exp_scale = 0.0,
616			.linear_term = 1.0 / 15000.0,
617			.constant_term = -2.0 / 3.0},
618			{.width = AH,
619			.exp_term = 0.0,
620			.exp_scale = 0.0,
621			.linear_term = -1.0 / 15000.0,
622			.constant_term = 8.0 / 3.0},
623			}},
624			.rayleigh_density =
625			{.layers =
626			{
627			{.width = AH,
628			.exp_term = 1.0,
629			.exp_scale = -1.0 / HR_MS,
630			.linear_term = 0.0,
631			.constant_term = 0.0},
632			}},
633			.beta_r0 = BR0_MS,
634			.beta_scale = aod / AOD0_CA,
635			.beta_m = NULL,
636			.grefl = grefl
637			};
638			return atmos;
639			}
640
641			int
642			main
643			(
644			int argc,
645			char *argv[]
646			)
647			{
648			int month, day;
649			double hour;
650			FVECT sundir;
651			int num_threads = 1;
652			int sorder = 4;
653			int year = 0;
654			int tsolar = 0;
655			int got_meridian = 0;
656			double grefl = 0.2;
657			double ccover = 0.0;
658			int res = 64;
659			double aod = AOD0_CA;
660			char *outname = "out";
661			char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
662			char mie_name[20] = "mie_ca";
663			char lstag[3];
664			char *ddir = ".";
665			int i;
666			double dirnorm = 0; /* direct normal illuminance */
667			double difhor = 0; /* diffuse horizontal illuminance */
668
669			fixargv0(argv[0]);
670			if (argc == 2 && !strcmp(argv[1], "-defaults")) {
671			printf("-i %d\t\t\t\t#scattering order\n", sorder);
672			printf("-g %f\t\t\t#ground reflectance\n", grefl);
673			printf("-c %f\t\t\t#cloud cover\n", ccover);
674			printf("-r %d\t\t\t\t#image resolution\n", res);
675			printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
676			printf("-f %s\t\t\t\t#output name (-f)\n", outname);
677			printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
678			exit(0);
679			}
680
681			if (argc < 4) {
682			fprintf(
683			stderr,
684			"Usage: %s month day hour -y year -a lat -o lon -m tz -d aod "
685			"-r "
686			"res -n nproc -c ccover -l mie -L dirnorm_illum difhor_illum "
687			"-g grefl -f outpath\n",
688			argv[0]);
689			return 0;
690			}
691
692			month = atoi(argv[1]);
693			if (month < 1 \|\| month > 12) {
694			fprintf(stderr, "bad month");
695			exit(1);
696			}
697			day = atoi(argv[2]);
698			if (day < 1 \|\| day > 31) {
699			fprintf(stderr, "bad month");
700			exit(1);
701			}
702			got_meridian = cvthour(argv[3], &tsolar, &hour);
703
704			if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
705			fprintf(stderr, "Cannot compute solar angle\n");
706			exit(1);
707			}
708
709			for (i = 4; i < argc; i++) {
710			if (argv[i][0] == '-') {
711			switch (argv[i][1]) {
712			case 'a':
713			s_latitude = atof(argv[++i]) * (PI / 180.0);
714			break;
715			case 'c':
716			ccover = atof(argv[++i]);
717			break;
718			case 'd':
719			aod = atof(argv[++i]);
720			break;
721			case 'f':
722			outname = argv[++i];
723			break;
724			case 'g':
725			grefl = atof(argv[++i]);
726			break;
727			case 'i':
728			sorder = atoi(argv[++i]);
729			break;
730			case 'l':
731			mie_path = argv[++i];
732			basename(mie_path, mie_name, sizeof(mie_name));
733			break;
734			case 'm':
735			if (got_meridian) {
736			++i;
737			break;
738			}
739			s_meridian = atof(argv[++i]) * (PI / 180.0);
740			break;
741			case 'n':
742			num_threads = atoi(argv[++i]);
743			break;
744			case 'o':
745			s_longitude = atof(argv[++i]) * (PI / 180.0);
746			break;
747			case 'L':
748			dirnorm = atof(argv[++i]);
749			difhor = atof(argv[++i]);
750			break;
751			case 'p':
752			ddir = argv[++i];
753			break;
754			case 'r':
755			res = atoi(argv[++i]);
756			break;
757			case 'y':
758			year = atoi(argv[++i]);
759			break;
760			default:
761			fprintf(stderr, "Unknown option %s\n", argv[i]);
762			exit(1);
763			}
764			}
765			}
766			if (year && (year < 1950) \| (year > 2050)) {
767			fprintf(
768			stderr,
769			"%s: warning - year should be in range 1950-2050\n",
770			progname);
771			}
772			if (month && !tsolar && fabs(s_meridian - s_longitude) > 45 * PI / 180) {
773			fprintf(
774			stderr,
775			"%s: warning - %.1f hours btwn. standard meridian and "
776			"longitude\n",
777			progname,
778			(s_longitude - s_meridian) * 12 / PI);
779			}
780
781			Atmosphere clear_atmos = init_atmos(aod, grefl);
782
783			int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
784			if (s_latitude < 0) {
785			is_summer = !is_summer;
786			}
787			set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
788
789			/* Load mie density data */
790			DATARRAY *mie_dp = getdata(mie_path);
791			if (mie_dp == NULL) {
792			fprintf(stderr, "Error reading mie data\n");
793			return 0;
794			}
795			clear_atmos.beta_m = mie_dp;
796
797			char gsdir[PATH_MAX];
798			size_t siz = strlen(ddir);
799			if (ISDIRSEP(ddir[siz - 1])) {
800			ddir[siz - 1] = '\0';
801			}
802			snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
803			if (!make_directory(gsdir)) {
804			fprintf(stderr, "Failed creating atmos_data directory");
805			exit(1);
806			}
807			DpPaths clear_paths = get_dppaths(gsdir, aod, mie_name, lstag);
808
809			if (getpath(clear_paths.tau, ".", R_OK) == NULL \|\|
810			getpath(clear_paths.scat, ".", R_OK) == NULL \|\|
811			getpath(clear_paths.scat1m, ".", R_OK) == NULL \|\|
812			getpath(clear_paths.irrad, ".", R_OK) == NULL) {
813			printf("# Pre-computing...\n");
814			if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
815			fprintf(stderr, "Pre-compute failed\n");
816			return 0;
817			}
818			}
819
820			DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
821			DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
822			DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
823			DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
824
825			write_header(argc, argv, ccover, grefl, res);
826
827			if (!gen_spect_sky(
828			tau_clear_dp,
829			scat_clear_dp,
830			scat1m_clear_dp,
831			irrad_clear_dp,
832			ccover,
833			sundir,
834			grefl,
835			res,
836			outname,
837			ddir,
838			dirnorm,
839			difhor)) {
840			fprintf(stderr, "gen_spect_sky failed\n");
841			exit(1);
842			}
843
844			freedata(mie_dp);
845			freedata(tau_clear_dp);
846			freedata(scat_clear_dp);
847			freedata(irrad_clear_dp);
848			freedata(scat1m_clear_dp);
849	greg	2.1
850	greg	2.9	return 1;
851	greg	2.1	}