src/gen/gensdaymtx.c

#ifndef lint
static const char RCSid[] = "$Id: gensdaymtx.c,v 1.4 2024/08/08 02:00:20 greg Exp $";
#endif

#include <stdlib.h>
#include <ctype.h>
#ifdef _WIN32
#include <windows.h>
#else
#include <errno.h>
#include <sys/stat.h>
#include <sys/types.h>
#endif

#include "atmos.h"
#include "copyright.h"
#include "data.h"
#include "platform.h"
#include "rtio.h"
#include "rtmath.h"
#include "sun.h"
#include "loadEPW.h"


const   double  SUN_ANG_DEG             = 0.533;                 /* sun full-angle in degrees */
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;

/* 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};

enum {
        NSUNPATCH = 4      /* max. # patches to spread sun into */
};

char    *progname;
double  altitude;          /* Solar altitude (radians) */
double  azimuth;                /* Solar azimuth (radians) */
int             julian_date;    /* Julian date */
double  sun_zenith;      /* Sun zenith angle (radians) */
int             nskypatch;        /* number of Reinhart patches */
float   *rh_palt;          /* sky patch altitudes (radians) */
float   *rh_pazi;          /* sky patch azimuths (radians) */
float   *rh_dom;                /* sky patch solid angle (sr) */
FVECT   sundir;
double  sun_ct;          /* cos(theta) of sun altitude angle */

int             input                   = 0;                                    /* Input type */
int             output                  = 0;                                    /* Output type */
int             nsuns                   = NSUNPATCH;                    /* number of sun patches to use */
double  fixed_sun_sa    = -1;                              /* fixed solid angle per sun? */
int             verbose                 = 0;                                    /* progress reports to stderr? */
int             outfmt                  = 'a';                            /* output format */
int             rhsubdiv                = 1;                                    /* Reinhart sky subdivisions */
COLOR   skycolor                = {.96, 1.004, 1.118};  /* sky coloration */
COLOR   suncolor                = {1., 1., 1.};          /* sun color */
double  grefl                   = .2;                              /* ground reflectance */


static inline double deg_to_rad(double deg)
{
        return deg * (PI / 180.);
}

static inline double rad_to_deg(double rad)
{
        return rad * (180. / PI);
}

static inline void vectorize(double altitude, double azimuth, FVECT v)
{
        v[1] = cos(altitude);
        v[0] = (v)[1] * sin(azimuth);
        v[1] *= cos(azimuth);
        v[2] = sin(altitude);
}

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 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
}

static const char *getfmtname(int fmt)
{
        switch (fmt) {
        case 'a':
                return ("ascii");
        case 'f':
                return ("float");
        case 'd':
                return ("double");
        }
        return ("unknown");
}


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 / D65EFFICACY;
        }
        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);
}

double solar_sunset(int month, int day)
{
        float W;
        W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
        return 12 + (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
}


double solar_sunrise(int month, int day)
{
        float W;
        W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day))));
        return 12 - (PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (PI * 15);
}

int rh_init(void) 
{
#define NROW 7
        static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6};
        const double alpha = (PI / 2.) / (NROW * rhsubdiv + .5);
        int p, i, j;
        /* allocate patch angle arrays */
        nskypatch = 0;
        for (p = 0; p < NROW; p++)
                nskypatch += tnaz[p];
        nskypatch *= rhsubdiv * rhsubdiv;
        nskypatch += 2;
        rh_palt = (float *)malloc(sizeof(float) * nskypatch);
        rh_pazi = (float *)malloc(sizeof(float) * nskypatch);
        rh_dom = (float *)malloc(sizeof(float) * nskypatch);
        if ((rh_palt == NULL) | (rh_pazi == NULL) | (rh_dom == NULL)) {
                fprintf(stderr, "%s: out of memory in rh_init()\n", progname);
                exit(1);
        }
        rh_palt[0] = -PI / 2.; /* ground & zenith patches */
        rh_pazi[0] = 0.;
        rh_dom[0] = 2. * PI;
        rh_palt[nskypatch - 1] = PI / 2.;
        rh_pazi[nskypatch - 1] = 0.;
        rh_dom[nskypatch - 1] = 2. * PI * (1. - cos(alpha * .5));
        p = 1; /* "normal" patches */
        for (i = 0; i < NROW * rhsubdiv; i++) {
                const float ralt = alpha * (i + .5);
                const int ninrow = tnaz[i / rhsubdiv] * rhsubdiv;
                const float dom =
                2. * PI * (sin(alpha * (i + 1)) - sin(alpha * i)) / (double)ninrow;
                for (j = 0; j < ninrow; j++) {
                        rh_palt[p] = ralt;
                        rh_pazi[p] = 2. * PI * j / (double)ninrow;
                        rh_dom[p++] = dom;
                }
        }
        return nskypatch;
#undef NROW
}

/* Resize daylight matrix (GW) */
float *resize_dmatrix(float *mtx_data, int nsteps, int npatch)
{
        if (mtx_data == NULL)
                mtx_data = (float * ) malloc(sizeof(float) * NSSAMP * nsteps * npatch);
        else
                mtx_data = (float * ) realloc(mtx_data,
                                sizeof(float) * NSSAMP * nsteps * npatch);
        if (mtx_data == NULL) {
                fprintf(stderr,
                                "%s: out of memory in resize_dmatrix(%d,%d)\n",
                                progname, nsteps, npatch);
                exit(1);
  }
  return mtx_data;
}

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;
}

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)
{
        /* Set rayleigh density profile */
        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';
}


/* Add in solar direct to nearest sky patches (GW) */
void add_direct(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m, DATARRAY *irrad,
                double ccover, double dirnorm, float *parr)
{
        FVECT svec;
        double near_dprod[NSUNPATCH];
        int near_patch[NSUNPATCH];
        double wta[NSUNPATCH], wtot;
        int i, j, p;

        /* identify nsuns closest patches */
        for (i = nsuns; i--;)
                near_dprod[i] = -1.;
        vectorize(altitude, azimuth, svec);
        for (p = 1; p < nskypatch; p++) {
                FVECT pvec;
                double dprod;
                vectorize(rh_palt[p], rh_pazi[p], pvec);
                dprod = DOT(pvec, svec);
                for (i = 0; i < nsuns; i++)
                        if (dprod > near_dprod[i]) {
                                for (j = nsuns; --j > i;) {
                                        near_dprod[j] = near_dprod[j - 1];
                                        near_patch[j] = near_patch[j - 1];
                                }
                                near_dprod[i] = dprod;
                                near_patch[i] = p;
                                break;
                        }
        }
        /* Get solar radiance */
        double sun_radiance[NSSAMP] = {0};
        get_solar_radiance(tau, scat, scat1m, sundir, ER, sun_ct, sun_radiance);
        if (ccover > 0) {
                double zenithbr = get_overcast_zenith_brightness(sundir);
                double skybr = get_overcast_brightness(sundir[2], zenithbr);
                int l;
                for (l = 0; l < NSSAMP; ++l) {
                        sun_radiance[l] = wmean2(sun_radiance[l], D6415[l] * skybr / WVLSPAN, ccover);
                }
        }
        /* Normalize */
        double sum = 0.0;
        for (i = 0; i < NSSAMP; ++i) {
                sum += sun_radiance[i];
        }
        double mean = sum / NSSAMP;

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

        /* weight by proximity */
        wtot = 0;
        for (i = nsuns; i--;) 
                wtot += wta[i] = 1. / (1.002 - near_dprod[i]);
        /* add to nearest patch radiances */
        for (i = nsuns; i--;) {
                float *pdest = parr + NSSAMP * near_patch[i];
                int k;
                for (k = 0; k < NSSAMP; k++) {
                        *pdest++ = sun_radiance[k] * wta[i] / wtot;
                }
        }
}


void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m,
                DATARRAY *irrad_clear, double ccover, double dif_ratio,
                double overcast_zenithbr, FVECT view_point, float *parr)
{
        double  mu_sky; /* Sun-sky point azimuthal angle */
        double  sspa;   /* Sun-sky point angle */
        int      i;
        for (i = 1; i < nskypatch; i++) {
                FVECT   rdir_sky;
                vectorize(rh_palt[i], rh_pazi[i], rdir_sky);
                mu_sky = fdot(view_point, rdir_sky) / ER;
                sspa = fdot(rdir_sky, sundir);

                SCOLOR  sky_radiance = {0};
                get_sky_radiance(scat, scat1m, ER, mu_sky, sun_ct, sspa, sky_radiance);
                int     k;
                for (k = 0; k < NSSAMP; ++k) {
                        sky_radiance[k] *= WVLSPAN;
                }

                if (ccover > 0) {
                        double skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr);
                        for (k = 0; k < NSSAMP; ++k) {
                                sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover);
                        }
                }

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

                for (k = 0; k < NSSAMP; ++k) {
                        parr[NSSAMP * i + k] = sky_radiance[k];
                }
        }
}


/* Compute sky patch radiance values (modified by GW) */
void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m,
                DATARRAY *irrad, double ccover, double difhor, FVECT view_point, float *parr)
{
        float sun_zenith;
        SCOLOR sky_radiance = {0};
        SCOLOR ground_radiance = {0};
        SCOLR sky_sclr = {0};
        SCOLR ground_sclr = {0};
        const double radius = VLEN(view_point);
        const double sun_ct = fdot(view_point, sundir) / radius;
        const FVECT rdir_grnd = {0, 0, -1};
        const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
        const double nu_grnd = fdot(rdir_grnd, sundir);

        /* Calculate sun zenith angle (don't let it dip below horizon) */
        /* Also limit minimum angle to keep circumsolar off zenith */
        if (altitude <= 0.0)
                sun_zenith = deg_to_rad(90.0);
        else if (altitude >= deg_to_rad(87.0))
                sun_zenith = deg_to_rad(3.0);
        else
                sun_zenith = deg_to_rad(90.0) - altitude;

        double overcast_zenithbr = get_overcast_zenith_brightness(sundir);

        /* diffuse calibration factor */
        double dif_ratio = 1;
        if (difhor > 0) {
                DATARRAY *indirect_irradiance_clear = get_indirect_irradiance(irrad, 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, ccover);
                if (diffuse_irradiance > 0) {
                        dif_ratio = difhor / WHTEFFICACY / diffuse_irradiance / 1.15;      /* fudge */
                }
        }

        /* Compute ground radiance (include solar contribution if any) */
        get_ground_radiance(tau, scat, scat1m, irrad, view_point, rdir_grnd, radius,
                                                mu_grnd, sun_ct, nu_grnd, grefl, sundir, parr);
        int j;
        for (j = 0; j < NSSAMP; j++) {
                parr[j] *= WVLSPAN;
        }
        calc_sky_patch_radiance(scat, scat1m, irrad, ccover, dif_ratio, overcast_zenithbr, view_point, parr);
}

int main(int argc, char *argv[])
{
        EPWheader       *epw = NULL;                    /* EPW/WEA input file */
        EPWrecord       erec;                                   /* current EPW/WEA input record */
        int                     doheader = 1;                   /* output header? */
        double          rotation = 0.0;                 /* site rotation (degrees) */
        double          elevation = 0;                  /* site elevation (meters) */
        int                     leap_day = 0;                   /* add leap day? */
        int                     sun_hours_only = 0;             /* only output sun hours? */
        int                     dir_is_horiz;                   /* direct is meas. on horizontal? */
        int                     ntsteps = 0;                    /* number of time steps */
        int                     tstorage = 0;                   /* number of allocated time steps */
        int                     nstored = 0;                    /* number of time steps in matrix */
        int                     last_monthly = 0;               /* month of last report */
        double          dni;                                    /* direct normal illuminance */
        double          dhi;                                    /* diffuse horizontal illuminance */

        float           *mtx_data = NULL;
        int                     mtx_offset = 0;
        double          timeinterval = 0;
        char            lstag[3];
        char            *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
        char            *ddir = ".";
        char            mie_name[20] = "mie_ca";
        int                     num_threads = 1;
        int                     sorder = 4;
        int                     solar_only = 0;
        int                     sky_only = 0;
        int                     i, j;
        FVECT           view_point      = {0, 0, ER};

        progname = argv[0];

        for (i = 1; i < argc && argv[i][0] == '-'; i++) {
                switch (argv[i][1]) {
                case 'd': /* solar (direct) only */
                        solar_only = 1;
                        break;
                case 's': /* sky only (no direct) */
                        sky_only = 1;
                        break;
                case 'g':
                        grefl = atof(argv[++i]);
                        break;
                case 'm':
                        rhsubdiv = atoi(argv[++i]);
                        break;
                case 'n':
                        num_threads = atoi(argv[++i]);
                        break;
                case 'r': /* rotate distribution */
                        if (argv[i][2] && argv[i][2] != 'z')
                                goto userr;
                        rotation = atof(argv[++i]);
                        break;
                case 'u': /* solar hours only */
                        sun_hours_only = 1;
                        break;
                case 'p':
                        ddir = argv[++i];
                        break;
                case 'v': /* verbose progress reports */
                        verbose++;
                        break;
                case 'h': /* turn off header */
                        doheader = 0;
                        break;
                case '5': /* 5-phase calculation */
                        nsuns = 1;
                        fixed_sun_sa = PI / 360. * atof(argv[++i]);
                        if (fixed_sun_sa <= 0) {
                                fprintf(
                                        stderr, 
                                        "%s: missing solar disk size argument for '-5' option\n",
                                        progname);
                                exit(1);
                        }
                        fixed_sun_sa *= fixed_sun_sa * PI;
                        break;
                case 'i':
                        timeinterval = atof(argv[++i]);
                        break;
                case 'o': /* output format */
                        switch (argv[i][2]) {
                        case 'f':
                        case 'd':
                        case 'a':
                                outfmt = argv[i][2];
                                break;
                        default:
                                goto userr;
                        }
                        break;
                default:
                        goto userr;
                }
        }
        if (i < argc - 1)
                goto userr;

        epw = EPWopen(argv[i]);
        if (epw == NULL)
                exit(1);
        if (i == argc - 1 && freopen(argv[i], "r", stdin) == NULL) {
                fprintf(stderr, "%s: cannot open '%s' for input\n", progname, argv[i]);
                exit(1);
        }
        if (verbose) {
                if (i == argc - 1)
                        fprintf(stderr, "%s: reading weather tape '%s'\n", progname, argv[i]);
                else
                        fprintf(stderr, "%s: reading weather tape from <stdin>\n", progname);
        }
        s_latitude = epw->loc.latitude;
        s_longitude = -epw->loc.longitude;
        s_meridian = -15.*epw->loc.timezone;
        elevation = epw->loc.elevation;
        switch (epw->isWEA) {           /* translate units */
        case WEAnot:
        case WEAradnorm:
                input = 1;              /* radiometric quantities */
                dir_is_horiz = 0;       /* direct is perpendicular meas. */
                break;
        case WEAradhoriz:
                input = 1;              /* radiometric quantities */
                dir_is_horiz = 1;       /* solar measured horizontally */
                break;
        case WEAphotnorm:
                input = 2;              /* photometric quantities */
                dir_is_horiz = 0;       /* direct is perpendicular meas. */
                break;
        default:
                goto fmterr;
        }

        rh_init();

        if (verbose) {
                fprintf(stderr, "%s: location '%s %s'\n", progname, epw->loc.city, epw->loc.country);
                fprintf(
                        stderr, 
                        "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
                        progname, s_latitude, s_longitude);
                if (rotation != 0)
                        fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation);
        }

        s_latitude = deg_to_rad(s_latitude);
        s_longitude = deg_to_rad(s_longitude);
        s_meridian = deg_to_rad(s_meridian);
        /* initial allocation */
        mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch);

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

        if (epw->isWEA == WEAnot) {
                fprintf(stderr, "EPW input\n");
        } else if (epw->isWEA != WEAphotnorm) {
                fprintf(stderr, "need WEA in photopic unit\n");
                exit(1);
        }

        while ((j = EPWread(epw, &erec)) > 0) {
                const int          mo = erec.date.month+1;
                const int          da = erec.date.day;
                const double    hr = erec.date.hour;
                double aod = erec.optdepth;
                double cc = erec.skycover;
                double            sda, sta, st;
                int                      sun_in_sky;

                if (aod == 0.0) {
                        aod = AOD0_CA;
                        fprintf(stderr, "aod is zero, using default value %.3f\n", AOD0_CA);
                }
                /* compute solar position */
                if ((mo == 2) & (da == 29)) {
                        julian_date = 60;
                        leap_day = 1;
                } else
                        julian_date = jdate(mo, da) + leap_day;
                sda = sdec(julian_date);
                sta = stadj(julian_date);
                st = hr + sta;
                if (timeinterval > 0) {
                        if (fabs(solar_sunrise(mo, da) - st) <= timeinterval/120)
                                st = (st + timeinterval/120 + solar_sunrise(mo, da))/2;
                        else if (fabs(solar_sunset(mo, da) - st) < timeinterval/120)
                                st = (st - timeinterval/120 + solar_sunset(mo, da))/2;
                }
                altitude = salt(sda, st);
                sun_in_sky = (altitude > -deg_to_rad(SUN_ANG_DEG / 2.));

                azimuth = sazi(sda, st) + PI - deg_to_rad(rotation);

                vectorize(altitude, azimuth, sundir);
                if (sun_hours_only && !sun_in_sky) {
                        continue; /* skipping nighttime points */
                }
                sun_ct = fdot(view_point, sundir) / ER;

                dni = erec.dirillum;
                dhi = erec.diffillum;

                mtx_offset = NSSAMP * nskypatch * nstored;
                nstored += 1;

                /* make space for next row */
                if (nstored > tstorage) {
                        tstorage += (tstorage >> 1) + nstored + 7;
                        mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
                }
                ntsteps++; /* keep count of time steps */

                /* compute sky patch values */
                Atmosphere clear_atmos = init_atmos(aod, grefl);
                int is_summer = (mo >= SUMMER_START && mo <= SUMMER_END);
                if (s_latitude < 0) {
                        is_summer = !is_summer;
                }
                set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);

                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) {
                        fprintf(stderr, "# 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);

                if (!solar_only)
                        compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
                                  cc, dhi, view_point, mtx_data + mtx_offset);
                if (!sky_only)
                        add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp,
                                 cc, dni, mtx_data + mtx_offset);
                /* monthly reporting */
                if (verbose && mo != last_monthly)
                        fprintf(stderr, "%s: stepping through month %d...\n", progname,
                                last_monthly = mo);
        }
        if (j != EOF) {
                fprintf(stderr, "%s: error on input\n", progname);
                exit(1);
        }
        EPWclose(epw); epw = NULL;
        freedata(mie_dp);
        if (!ntsteps) {
                fprintf(stderr, "%s: no valid time steps on input\n", progname);
                exit(1);
        }
        /* write out matrix */
        if (outfmt != 'a')
                SET_FILE_BINARY(stdout);
#ifdef getc_unlocked
        flockfile(stdout);
#endif
        if (verbose)
                fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname,
                        outfmt == 'a' ? "" : "binary ", nstored);
        if (doheader) {
                newheader("RADIANCE", stdout);
                printargs(argc, argv, stdout);
                printf("LATLONG= %.8f %.8f\n", rad_to_deg(s_latitude),
                        -rad_to_deg(s_longitude));
                printf("NROWS=%d\n", nskypatch);
                printf("NCOLS=%d\n", nstored);
                printf("NCOMP=%d\n", NSSAMP);
                float wvsplit[4] = {380, 480, 588, 780};
                fputwlsplit(wvsplit, stdout);
                if ((outfmt == 'f') | (outfmt == 'd'))
                        fputendian(stdout);
                fputformat((char *)getfmtname(outfmt), stdout);
                putchar('\n');
        }
        /* patches are rows (outer sort) */
        for (i = 0; i < nskypatch; i++) {
                mtx_offset = NSSAMP * i;
                switch (outfmt) {
                case 'a':
                        for (j = 0; j < nstored; j++) {
                                int k;
                                for (k = 0; k < NSSAMP; k++) {
                                        printf("%.3g ", mtx_data[mtx_offset + k]);
                                }
                                printf("\n");
                                mtx_offset += NSSAMP * nskypatch;
                        }
                        if (nstored > 1)
                                fputc('\n', stdout);
                        break;
                case 'f':
                        for (j = 0; j < nstored; j++) {
                                putbinary(mtx_data + mtx_offset, sizeof(float), NSSAMP, stdout);
                                mtx_offset += NSSAMP * nskypatch;
                        }
                        break;
                case 'd':
                        for (j = 0; j < nstored; j++) {
                                double ment[NSSAMP];
                                for (j = 0; j < NSSAMP; j++)
                                        ment[j] = mtx_data[mtx_offset + j];
                                putbinary(ment, sizeof(double), NSSAMP, stdout);
                                mtx_offset += NSSAMP * nskypatch;
                        }
                        break;
                }
                if (ferror(stdout))
                        goto writerr;
        }
        return 0;
userr:
        fprintf(stderr,
                  "Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-r "
                  "deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
                  progname);
        exit(1);
fmterr:
        fprintf(stderr, "%s: weather tape format error in header\n", progname);
        exit(1);
writerr:
        fprintf(stderr, "%s: write error on output\n", progname);
        exit(1);
}
Revision:	1.6
Committed:	Tue Apr 15 20:15:50 2025 UTC (2 weeks, 3 days ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.5:	+739 -727 lines
Log Message:	fix(gensdaymtx): Added WAVELENGTH_SPLITS= to header output and code cleanup by TW