--- ray/src/gen/gensdaymtx.c 2024/08/08 02:00:20 1.4 +++ ray/src/gen/gensdaymtx.c 2025/04/10 23:30:58 1.5 @@ -1,13 +1,9 @@ #ifndef lint -static const char RCSid[] = "$Id: gensdaymtx.c,v 1.4 2024/08/08 02:00:20 greg Exp $"; +static const char RCSid[] = "$Id: gensdaymtx.c,v 1.5 2025/04/10 23:30:58 greg Exp $"; #endif -#include "atmos.h" -#include "copyright.h" -#include "data.h" -#include "platform.h" -#include "rtio.h" -#include + #include +#include #ifdef _WIN32 #include #else @@ -16,6 +12,30 @@ static const char RCSid[] = "$Id: gensdaymtx.c,v 1.4 2 #include #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" + +#ifndef M_PI + #define M_PI 3.14159265358579 +#endif + +#define vector(v, alt, azi) \ + ((v)[1] = cos(alt), (v)[0] = (v)[1] * sin(azi), (v)[1] *= cos(azi), \ + (v)[2] = sin(alt)) + +#define rh_vector(v, i) vector(v, rh_palt[i], rh_pazi[i]) + +#define rh_cos(i) tsin(rh_palt[i]) + +#define solar_minute(jd, hr) ((24 * 60) * ((jd) - 1) + (int)((hr) * 60. + .5)) + + char *progname; double altitude; /* Solar altitude (radians) */ @@ -71,16 +91,6 @@ float *rh_dom; /* sky patch solid angle (sr) */ double sun_ct; -#define vector(v, alt, azi) \ - ((v)[1] = cos(alt), (v)[0] = (v)[1] * sin(azi), (v)[1] *= cos(azi), \ - (v)[2] = sin(alt)) - -#define rh_vector(v, i) vector(v, rh_palt[i], rh_pazi[i]) - -#define rh_cos(i) tsin(rh_palt[i]) - -#define solar_minute(jd, hr) ((24 * 60) * ((jd) - 1) + (int)((hr) * 60. + .5)) - inline void vectorize(double altitude, double azimuth, FVECT v) { v[1] = cos(altitude); v[0] = (v)[1] * sin(azimuth); @@ -123,7 +133,8 @@ static inline double wmean(const double a, const doubl return (a * x + b * y) / (a + b); } -static double get_zenith_brightness(const double sundir[3]) { + +static double get_overcast_zenith_brightness(const double sundir[3]) { double zenithbr; if (sundir[2] < 0) { zenithbr = 0; @@ -140,6 +151,23 @@ static double get_overcast_brightness(const double dz, 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 + (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15)); +} + + +double +solar_sunrise(int month, int day) +{ + float W; + W = -1 * (tan(s_latitude) * tan(sdec(jdate(month, day)))); + return(12 - (M_PI / 2 - atan2(W, sqrt(1 - W * W))) * 180 / (M_PI * 15)); +} + int rh_init(void) { #define NROW 7 static const int tnaz[NROW] = {30, 30, 24, 24, 18, 12, 6}; @@ -196,7 +224,7 @@ float *resize_dmatrix(float *mtx_data, int nsteps, int } static Atmosphere init_atmos(const double aod, const double grefl) { - Atmosphere atmos = {.ozone_density = {.layers = + Atmosphere atmos = {.ozone_density = {.layers = { {.width = 25000.0, .exp_term = 0.0, @@ -221,146 +249,168 @@ static Atmosphere init_atmos(const double aod, const d .beta_scale = aod / AOD0_CA, .beta_m = NULL, .grefl = grefl}; - return atmos; + return atmos; } static DpPaths get_dppaths(const char *dir, const double aod, const char *mname, const char *tag) { - DpPaths paths; + DpPaths paths; - snprintf(paths.tau, PATH_MAX, "%s%ctau_%s_%s_%.2f.dat", dir, DIRSEP, tag, + 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, + 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, + 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, + snprintf(paths.irrad, PATH_MAX, "%s%cirrad_%s_%s_%.2f.dat", dir, DIRSEP, tag, mname, aod); - return paths; + 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; + /* 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[1] = 'w'; - atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW; - atmos->beta_r0 = BR0_SW; + tag[0] = 't'; + tag[1] = 'r'; + atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T; + atmos->beta_r0 = BR0_T; } - } 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'; + 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, float *parr) { - FVECT svec; - double near_dprod[NSUNPATCH]; - int near_patch[NSUNPATCH]; - double wta[NSUNPATCH], wtot; - int i, j, p; + 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]; + /* 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); } - 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_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); } - } - /* 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; + /* Normalize */ + double sum = 0.0; + for (i = 0; i < NSSAMP; ++i) { + sum += sun_radiance[i]; } - } -} + double mean = sum / NSSAMP; -void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m, double ccover, - float *parr) { - int i; - double mu_sky; /* Sun-sky point azimuthal angle */ - double sspa; /* Sun-sky point angle */ - FVECT view_point = {0, 0, ER}; - for (i = 1; i < nskypatch; i++) { - FVECT rdir_sky; - int k; - 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); - for (k = 0; k < NSSAMP; ++k) { - sky_radiance[k] *= WVLSPAN; + 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; + } - if (ccover > 0) { - double zenithbr = get_zenith_brightness(sundir); - double grndbr = zenithbr * GNORM; - double skybr = get_overcast_brightness(rdir_sky[2], zenithbr); - int k; - for (k = 0; k < NSSAMP; ++k) { - sky_radiance[k] = wmean2(sky_radiance[k], skybr * D6415[k], ccover); - } + /* 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; + } } +} - for (k = 0; k < NSSAMP; ++k) { - parr[NSSAMP * i + k] = sky_radiance[k]; +void calc_sky_patch_radiance(DATARRAY *scat, DATARRAY *scat1m, DATARRAY *irrad_clear, double ccover, double dif_ratio, + double overcast_zenithbr, float *parr) { + int i; + double mu_sky; /* Sun-sky point azimuthal angle */ + double sspa; /* Sun-sky point angle */ + FVECT view_point = {0, 0, ER}; + for (i = 1; i < nskypatch; i++) { + FVECT rdir_sky; + int k; + 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); + for (k = 0; k < NSSAMP; ++k) { + sky_radiance[k] *= WVLSPAN; + } + + if (ccover > 0) { + double skybr = get_overcast_brightness(rdir_sky[2], overcast_zenithbr); + int k; + 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]; + } } - } } /* Return maximum of two doubles */ @@ -368,354 +418,400 @@ static inline double dmax(double a, double b) { return /* Compute sky patch radiance values (modified by GW) */ void compute_sky(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m, - DATARRAY *irrad, double ccover, float *parr) { - int index; /* Category index */ - int i; - float sun_zenith; - SCOLOR sky_radiance = {0}; - SCOLOR ground_radiance = {0}; - SCOLR sky_sclr = {0}; - SCOLR ground_sclr = {0}; - FVECT view_point = {0, 0, ER}; - 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); - int j; + DATARRAY *irrad, double ccover, double difhor, float *parr) { + float sun_zenith; + SCOLOR sky_radiance = {0}; + SCOLOR ground_radiance = {0}; + SCOLR sky_sclr = {0}; + SCOLR ground_sclr = {0}; + FVECT view_point = {0, 0, ER}; + 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 = DegToRad(90.0); - else if (altitude >= DegToRad(87.0)) - sun_zenith = DegToRad(3.0); - else - sun_zenith = DegToRad(90.0) - altitude; + /* 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 = DegToRad(90.0); + else if (altitude >= DegToRad(87.0)) + sun_zenith = DegToRad(3.0); + else + sun_zenith = DegToRad(90.0) - altitude; - /* 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); - for (j = 0; j < NSSAMP; j++) { - parr[j] *= WVLSPAN; - } - calc_sky_patch_radiance(scat, scat1m, ccover, parr); + 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, parr); } int main(int argc, char *argv[]) { - char buf[256]; - int doheader = 1; /* output header? */ - double rotation = 0.0; - double elevation = 0; - int leap_day = 0; /* add leap day? */ - int sun_hours_only = 0; /* only output sun hours? */ - float *mtx_data = NULL; - 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 */ - int mo, da; - double hr, aod, cc; - double dni, dhi; - int mtx_offset = 0; - int i, j; - 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; - FVECT view_point = {0, 0, ER}; + EPWheader *epw = NULL; /* EPW/WEA input file */ + EPWrecord erec; /* current EPW/WEA input record */ + int doheader = 1; /* output header? */ + double rotation = 0.0; + double elevation = 0; + 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? */ + float *mtx_data = NULL; + 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, dhi; + 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; + FVECT view_point = {0, 0, ER}; + int i, j; - progname = argv[0]; + 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') + 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; - 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); + + epw = EPWopen(argv[i]); + if (epw == NULL) exit(1); - } - fixed_sun_sa *= fixed_sun_sa * PI; - break; - case 'o': /* output format */ - switch (argv[i][2]) { - case 'f': - case 'd': - case 'a': - outfmt = argv[i][2]; + 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 \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; - default: - goto userr; - } - 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 userr; + goto fmterr; } - } - if (i < argc - 1) - goto userr; - 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 \n", progname); - } - /* read weather tape header */ - if (scanf("place %[^\r\n] ", buf) != 1) - goto fmterr; - if (scanf("latitude %lf\n", &s_latitude) != 1) - goto fmterr; - if (scanf("longitude %lf\n", &s_longitude) != 1) - goto fmterr; - if (scanf("time_zone %lf\n", &s_meridian) != 1) - goto fmterr; - if (scanf("site_elevation %lf\n", &elevation) != 1) - goto fmterr; - if (scanf("weather_data_file_units %d\n", &input) != 1) - goto fmterr; - rh_init(); - if (verbose) { - fprintf(stderr, "%s: location '%s'\n", progname, buf); - fprintf(stderr, "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n", + 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); - } + if (rotation != 0) + fprintf(stderr, "%s: rotating output %.0f degrees\n", progname, rotation); + } - s_latitude = DegToRad(s_latitude); - s_longitude = DegToRad(s_longitude); - s_meridian = DegToRad(s_meridian); - /* initial allocation */ - mtx_data = resize_dmatrix(mtx_data, tstorage = 2, nskypatch); + s_latitude = DegToRad(s_latitude); + s_longitude = DegToRad(s_longitude); + s_meridian = DegToRad(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; - } + /* Load mie density data */ + DATARRAY *mie_dp = getdata(mie_path); + if (mie_dp == NULL) { + fprintf(stderr, "Error reading mie data\n"); + return 0; + } - while (scanf("%d %d %lf %lf %lf %lf %lf\n", &mo, &da, &hr, &dni, &dhi, &aod, - &cc) == 7) { - if (aod == 0.0) { - aod = AOD0_CA; - fprintf(stderr, "aod is zero, using default value %.3f\n", AOD0_CA); + if (epw->isWEA == WEAnot) { + fprintf(stderr, "EPW input\n"); + } else if (epw->isWEA != WEAphotnorm) { + fprintf(stderr, "need WEA in photopic unit\n"); + exit(1); } - double sda, sta; - 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); - altitude = salt(sda, hr + sta); - sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG / 2.)); - azimuth = sazi(sda, hr + sta) + PI - DegToRad(rotation); + 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; - vectorize(altitude, azimuth, sundir); - if (sun_hours_only && sundir[2] <= 0.) { - continue; /* skipping nighttime points */ - } - sun_ct = fdot(view_point, sundir) / ER; + 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 > -DegToRad(SUN_ANG_DEG / 2.)); - mtx_offset = NSSAMP * nskypatch * nstored; - nstored += 1; - printf("mtx_offset = %d nstored = %d nskypatch = %d\n", mtx_offset, nstored, - nskypatch); - /* make space for next row */ - if (nstored > tstorage) { - printf("make space for next row nstored = %d tstorage = %d\n", nstored, - tstorage); - tstorage += (tstorage >> 1) + nstored + 7; - mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch); - } - ntsteps++; /* keep count of time steps */ + azimuth = sazi(sda, st) + PI - DegToRad(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; + printf("%d %d %f %f %f %f %f\n", mo, da, hr, dni, dhi, aod, cc); + + 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); + 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; + 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); + 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; - } - } + 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); + 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, mtx_data + mtx_offset); - if (!sky_only) - add_direct(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp, - cc, mtx_data + mtx_offset); - /* monthly reporting */ - if (verbose && mo != last_monthly) - fprintf(stderr, "%s: stepping through month %d...\n", progname, - last_monthly = mo); - } - freedata(mie_dp); - if (!ntsteps) { - fprintf(stderr, "%s: no valid time steps on input\n", progname); - exit(1); - } - /* check for junk at end */ - while ((i = fgetc(stdin)) != EOF) - if (!isspace(i)) { - fprintf(stderr, "%s: warning - unexpected data past EOT: ", progname); - buf[0] = i; - buf[1] = '\0'; - fgets(buf + 1, sizeof(buf) - 1, stdin); - fputs(buf, stderr); - fputc('\n', stderr); - break; + if (!solar_only) + compute_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp, irrad_clear_dp, + cc, dhi, 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); } - /* write out matrix */ - if (outfmt != 'a') - SET_FILE_BINARY(stdout); + 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); + flockfile(stdout); #endif - if (verbose) - fprintf(stderr, "%s: writing %smatrix with %d time steps...\n", progname, + 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", RadToDeg(s_latitude), - -RadToDeg(s_longitude)); - printf("NROWS=%d\n", nskypatch); - printf("NCOLS=%d\n", nstored); - printf("NCOMP=%d\n", NSSAMP); - 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]); + if (doheader) { + newheader("RADIANCE", stdout); + printargs(argc, argv, stdout); + printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude), + -RadToDeg(s_longitude)); + printf("NROWS=%d\n", nskypatch); + printf("NCOLS=%d\n", nstored); + printf("NCOMP=%d\n", NSSAMP); + 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; } - 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; } - if (ferror(stdout)) - goto writerr; - } -alldone: - if (fflush(NULL) == EOF) - goto writerr; - if (verbose) - fprintf(stderr, "%s: done.\n", progname); - exit(0); userr: - fprintf(stderr, + 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); + exit(1); fmterr: - fprintf(stderr, "%s: weather tape format error in header\n", progname); - exit(1); + 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); + fprintf(stderr, "%s: write error on output\n", progname); + exit(1); }