src/gen/gensdaymtx.c

#ifndef lint
static const char RCSid[] = "$Id: gensdaymtx.c,v 1.12 2025/08/04 23:45: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;
const double M_PER_KM        = 1e3;

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

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 = -1.0;                                                 /* direct normal illuminance */
        double dhi = -1.0;                                                 /* 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, k;
        FVECT view_point      = {0, 0, ER};

        fixargv0(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 1;
        }

        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 * M_PER_KM;
                if (aod >= 999.0) {
                        aod = AOD0_CA;
                        fprintf(stderr, "aod is not set, using default value %.3f\n", AOD0_CA);
                }
                double cc = erec.skycover;
                if (cc >= 99.0) {
                        cc = 0.0;
                        fprintf(stderr, "skycover is not set, using default value %.3f\n", 0.0);
                }
                double sda, sta, st;
                int sun_in_sky;

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

                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);
                fputwlsplit(WLPART, 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++) {
                                for (k = NSSAMP - 1; k >= 0; 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++) {
                                float ment[NSSAMP];
                                for (k = 0; k < NSSAMP; k++) {
                                        ment[NSSAMP-1 - k] = mtx_data[mtx_offset + k];
                                }
                                putbinary(ment, sizeof(float), NSSAMP, stdout);
                                mtx_offset += NSSAMP * nskypatch;
                        }
                        break;
                case 'd':
                        for (j = 0; j < nstored; j++) {
                                double ment[NSSAMP];
                                for (k = 0; k < NSSAMP; k++) {
                                        ment[NSSAMP-1 - k] = mtx_data[mtx_offset + k];
                                }
                                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.13
Committed:	Fri Aug 22 01:04:30 2025 UTC (2 months, 2 weeks ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.12:	+6 -6 lines
Log Message:	fix(gensdaymtx): Corrected handling of zero values