fix(weather): physically correct EPW radiation and sky fields (F1)

SAM/SAM.Tests/Resources/Weather/README.md

@@ -0,0 +1,36 @@
+# Reference EPW fixtures (EPW writer tests)
+
+Five TMYx EPW files, gzip-compressed (`gzip -9 -n`) without modification. `SHA256SUMS.txt` holds the
+SHA-256 of each **uncompressed** original. `WeatherEPWExportTests.ReferenceFixtures_MatchRecordedHashes`
+fails if a fixture changes.
+
+| File | Location | Lat / lon / time zone | Source | Why |
+|---|---|---|---|---|
+| `HUN_BU_Budapest.Met.Center.128400_TMYx.epw` | Budapest, HU | 47.517 / 19.033 / +1 | SRC-TMYx (ERA5-based) | The F1 incident file: the embedded weather of the SAM_UI PR #239 box model came from it |
+| `USA_AK_Fairbanks.Intl.AP.702610_TMYx.2004-2018.epw` | Fairbanks, AK | 64.804 / -147.876 / -9 | ISD-TMYx | High latitude and low sun; 13 degrees west of the zone meridian |
+| `AUS_NSW_Sydney.Intl.AP.947670_TMYx.2004-2018.epw` | Sydney, AU | -33.946 / 151.177 / +10 | ISD-TMYx | Southern hemisphere |
+| `SGP_SG_Changi.Intl.AP.486980_TMYx.2004-2018.epw` | Singapore | 1.350 / 103.994 / +8 | ISD-TMYx | Tropics; 14 degrees east of the zone meridian |
+| `USA_AZ_Phoenix-Sky.Harbor.Intl.AP.722780_TMYx.2004-2018.epw` | Phoenix, AZ | 33.433 / -112.017 / -7 | ISD-TMYx | Clear, high direct normal |
+
+## Provenance
+
+- Budapest: copied 2026-10-19 from the SAM_UI PR #239 validation inputs
+  (`C:\TasOut\licensed-e2e-20261019\inputs-box\Budapest.epw`, machine-local, byte-identical to the
+  `-EpwPath` file of that validation). Its `COMMENTS 2` header names climate.onebuilding.org as the source.
+- The other four: copied 2026-10-19 from the local Ladybug Tools bundle
+  (`C:\Program Files\ladybug_tools\resources\weather\<name>\<name>.epw`), which redistributes the TMYx
+  files published at https://climate.onebuilding.org/.
+
+## Licence
+
+TMYx files from climate.onebuilding.org are free to use with attribution: "Climatic data provided by
+climate.onebuilding.org".
+
+## Known properties used by the tests
+
+On the EPW time base (hour ending, local standard time), the Fairbanks, Sydney and Singapore direct
+normal closes with its own global - diffuse horizontal within about 1.3 %. The Budapest and Phoenix
+files do not (about 9 %). Their direct normal closes with global - diffuse when the solar geometry is
+taken about 30 minutes earlier, so within each file the direct normal and the horizontal components
+are on time bases about half an hour apart. Budapest's own extraterrestrial column is a further
+30 minutes later than standard-time geometry. The tests measure this rather than assume it.

SAM/SAM.Tests/Resources/Weather/SHA256SUMS.txt

@@ -0,0 +1,5 @@
+7641aaa476462f76c8598d0836205ccf3df8f8b79aa6343910349303a04d180d *AUS_NSW_Sydney.Intl.AP.947670_TMYx.2004-2018.epw
+472ffb4e2f10867600949a392ce46891fc1bd05f2fc8d8917a99bca7d26209c3 *HUN_BU_Budapest.Met.Center.128400_TMYx.epw
+d47089f83c4c35f74407050c45e54661f6148c6a35522c9a32af26e4b51144a9 *SGP_SG_Changi.Intl.AP.486980_TMYx.2004-2018.epw
+d1591f5ca40e07c229aa781d0514e84722a501678ff276f020698700a5da16b0 *USA_AK_Fairbanks.Intl.AP.702610_TMYx.2004-2018.epw
+dec64d7722ed677992f9023160c9c98a1043bb74cbcce696736cbf421a513a80 *USA_AZ_Phoenix-Sky.Harbor.Intl.AP.722780_TMYx.2004-2018.epw

SAM/SAM.Tests/SAM.Tests.csproj

@@ -82,6 +82,10 @@
     <None Include="Golden\GenOpt\*.tsv">
       <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
     </None>
+    <!-- TMYx reference EPWs, gzip-compressed (see Resources\Weather\README.md and WeatherEPWExportTests). -->
+    <None Include="Resources\Weather\*.epw.gz;Resources\Weather\SHA256SUMS.txt">
+      <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
+    </None>
   </ItemGroup>
 
 </Project>

SAM/SAM.Weather/Query/EPWSolarGeometry.cs

@@ -0,0 +1,188 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (c) 2020–2026 Michal Dengusiak & Jakub Ziolkowski and contributors
+
+using System;
+
+namespace SAM.Weather
+{
+    public static partial class Query
+    {
+        /// <summary>
+        /// Solar constant [W/m2] used for the extraterrestrial irradiance of the EPW record
+        /// (the value the EnergyPlus weather converter and the TMY2/TMY3 EPW sources use).
+        /// </summary>
+        private const double EPWSolarConstant = 1367.0;
+
+        /// <summary>
+        /// Number of equal sub-intervals used to average the solar geometry over one EPW hour.
+        /// </summary>
+        private const int EPWSolarSubIntervals = 60;
+
+        /// <summary>
+        /// Hour-averaged solar geometry of one EPW record, used to convert SAM horizontal beam
+        /// radiation into EPW direct-normal radiation.
+        /// <para>
+        /// An EPW hour <c>h</c> (1-24) holds the radiation received during the hour that ends at
+        /// <c>h:00</c> local standard time (EnergyPlus Auxiliary Programs, EPW data dictionary:
+        /// "during the number of minutes preceding the time indicated"); SAM stores it at hour index
+        /// <c>h - 1</c> (see <see cref="TryGetData"/>). The interval is therefore
+        /// <c>[hour, hour + 1)</c> local standard time for SAM <paramref name="hour"/>.
+        /// </para>
+        /// <para>
+        /// <paramref name="meanCosZenith"/> is the mean over that interval of max(0, cos zenith), so a
+        /// constant direct-normal irradiance <c>DNI</c> over the hour puts <c>DNI * meanCosZenith</c>
+        /// on the horizontal - the energy relation the conversion inverts. Sun position: NOAA
+        /// general solar position equations (Spencer 1971 Fourier series for declination and
+        /// equation of time), no refraction; accuracy is a few arc-minutes, well below the
+        /// sub-hour averaging.
+        /// </para>
+        /// <para>
+        /// <paramref name="limitDirectNormalIrradiance"/> is the largest direct-normal irradiance
+        /// that can reach the ground in the hour: the beam of a pure Rayleigh (aerosol-, water- and
+        /// ozone-free) atmosphere, <c>I0 exp(-dR(m') m')</c>, with the Kasten and Young (1989)
+        /// relative air mass <c>m</c>, <c>m' = m p / 101325</c> and the Kasten (1996) integral Rayleigh
+        /// optical thickness <c>dR</c>, averaged over the hour with cos-zenith weights (the constant
+        /// direct normal that would give that beam on the horizontal). A real atmosphere only
+        /// attenuates more, so the value is a physical bound, not a clear-sky model.
+        /// </para>
+        /// </summary>
+        /// <param name="dayOfYear">Day of year, 1-based.</param>
+        /// <param name="daysInYear">Days in the calendar year of the record (365 or 366).</param>
+        /// <param name="hour">SAM hour index 0-23 (EPW hour minus one).</param>
+        /// <param name="latitude">Latitude [degrees], north positive.</param>
+        /// <param name="longitude">Longitude [degrees], east positive.</param>
+        /// <param name="timeZone">Standard time zone [hours from UTC], east positive (EPW LOCATION convention).</param>
+        /// <param name="pressure">Station pressure [Pa] for the air-mass correction; NaN uses 101325 Pa.</param>
+        /// <param name="meanCosZenith">Hour mean of max(0, cos zenith) [-].</param>
+        /// <param name="extraterrestrialNormalIrradiance">Extraterrestrial direct-normal irradiance for the day [W/m2].</param>
+        /// <param name="limitDirectNormalIrradiance">Rayleigh-atmosphere bound on the hour's direct normal [W/m2]; 0 when the sun is down all hour.</param>
+        /// <returns>False when the location or time zone is not stated (NaN) or out of range.</returns>
+        private static bool TryGetEPWSolarGeometry(int dayOfYear, int daysInYear, int hour, double latitude, double longitude, double timeZone, double pressure, out double meanCosZenith, out double extraterrestrialNormalIrradiance, out double limitDirectNormalIrradiance)
+        {
+            meanCosZenith = double.NaN;
+            extraterrestrialNormalIrradiance = double.NaN;
+            limitDirectNormalIrradiance = double.NaN;
+
+            if (double.IsNaN(latitude) || double.IsNaN(longitude) || double.IsNaN(timeZone))
+                return false;
+
+            if (latitude < -90 || latitude > 90 || longitude < -180 || longitude > 180 || timeZone < -14 || timeZone > 14)
+                return false;
+
+            if (hour < 0 || hour > 23 || dayOfYear < 1 || dayOfYear > daysInYear)
+                return false;
+
+            double pressureRatio = double.IsNaN(pressure) || pressure <= 0 ? 1.0 : pressure / 101325.0;
+
+            double dayAngle = 2.0 * Math.PI * (dayOfYear - 1) / daysInYear;
+            extraterrestrialNormalIrradiance = EPWSolarConstant * (1.000110 + 0.034221 * Math.Cos(dayAngle) + 0.001280 * Math.Sin(dayAngle) + 0.000719 * Math.Cos(2 * dayAngle) + 0.000077 * Math.Sin(2 * dayAngle));
+
+            double latitudeRadians = latitude * Math.PI / 180.0;
+            double sumCosZenith = 0;
+            double sumLimitBeam = 0;
+            for (int i = 0; i < EPWSolarSubIntervals; i++)
+            {
+                // Mid-point of the sub-interval, in hours of local standard time.
+                double time = hour + (i + 0.5) / EPWSolarSubIntervals;
+
+                double gamma = 2.0 * Math.PI / daysInYear * (dayOfYear - 1 + (time - 12.0) / 24.0);
+
+                double equationOfTime = 229.18 * (0.000075 + 0.001868 * Math.Cos(gamma) - 0.032077 * Math.Sin(gamma) - 0.014615 * Math.Cos(2 * gamma) - 0.040849 * Math.Sin(2 * gamma));
+
+                double declination = 0.006918 - 0.399912 * Math.Cos(gamma) + 0.070257 * Math.Sin(gamma) - 0.006758 * Math.Cos(2 * gamma) + 0.000907 * Math.Sin(2 * gamma) - 0.002697 * Math.Cos(3 * gamma) + 0.00148 * Math.Sin(3 * gamma);
+
+                // True solar time [min]: standard meridian is 15 degrees per hour of time zone.
+                double trueSolarTime = time * 60.0 + equationOfTime + 4.0 * longitude - 60.0 * timeZone;
+                double hourAngle = (trueSolarTime / 4.0 - 180.0) * Math.PI / 180.0;
+
+                double cosZenith = Math.Sin(latitudeRadians) * Math.Sin(declination) + Math.Cos(latitudeRadians) * Math.Cos(declination) * Math.Cos(hourAngle);
+                if (cosZenith <= 0)
+                    continue;
+
+                sumCosZenith += cosZenith;
+
+                double airMass = RelativeAirMass(cosZenith) * pressureRatio;
+                sumLimitBeam += extraterrestrialNormalIrradiance * Math.Exp(-RayleighOpticalThickness(airMass) * airMass) * cosZenith;
+            }
+
+            meanCosZenith = sumCosZenith / EPWSolarSubIntervals;
+            limitDirectNormalIrradiance = sumCosZenith > 0 ? sumLimitBeam / sumCosZenith : 0;
+            return true;
+        }
+
+        /// <summary>
+        /// Relative optical air mass, Kasten and Young (1989), from the cosine of the solar zenith
+        /// (sun above the horizon).
+        /// </summary>
+        private static double RelativeAirMass(double cosZenith)
+        {
+            double zenith = Math.Acos(Math.Min(1.0, cosZenith)) * 180.0 / Math.PI;
+            return 1.0 / (cosZenith + 0.50572 * Math.Pow(96.07995 - zenith, -1.6364));
+        }
+
+        /// <summary>
+        /// Integral Rayleigh optical thickness, Kasten (1996), for the pressure-corrected air mass.
+        /// </summary>
+        private static double RayleighOpticalThickness(double airMass)
+        {
+            if (airMass <= 20)
+                return 1.0 / (6.6296 + 1.7513 * airMass - 0.1202 * Math.Pow(airMass, 2) + 0.0065 * Math.Pow(airMass, 3) - 0.00013 * Math.Pow(airMass, 4));
+
+            return 1.0 / (10.4 + 0.718 * airMass);
+        }
+
+        /// <summary>
+        /// Direct-normal radiation [Wh/m2] for one EPW hour from the horizontal beam (global minus
+        /// diffuse) and the hour-averaged solar geometry: <c>DNI = beam / meanCosZenith</c>, limited
+        /// to the Rayleigh-atmosphere bound of the hour. Horizontal beam in an hour where the sun is
+        /// below the horizon for the whole interval cannot be represented as a direct-normal value
+        /// and gives 0. NaN when the horizontal beam or the geometry is not known (written as the EPW
+        /// missing value).
+        /// </summary>
+        private static double EPWDirectNormalRadiation(double horizontalBeamRadiation, double meanCosZenith, double limitDirectNormalIrradiance)
+        {
+            if (double.IsNaN(horizontalBeamRadiation) || double.IsNaN(meanCosZenith))
+                return double.NaN;
+
+            if (horizontalBeamRadiation <= 0 || meanCosZenith <= 0)
+                return 0;
+
+            double result = horizontalBeamRadiation / meanCosZenith;
+            if (!double.IsNaN(limitDirectNormalIrradiance) && result > limitDirectNormalIrradiance)
+                result = limitDirectNormalIrradiance;
+
+            return result;
+        }
+
+        /// <summary>
+        /// Dew-point temperature [C] from dry-bulb temperature [C] and relative humidity [%]
+        /// (Magnus form over water, Alduchov and Eskridge 1996 coefficients; within about 0.35 K
+        /// from -40 C to 50 C). NaN when either input is missing or the humidity is not positive.
+        /// </summary>
+        private static double EPWDewPointTemperature(double dryBulbTemperature, double relativeHumidity)
+        {
+            if (double.IsNaN(dryBulbTemperature) || double.IsNaN(relativeHumidity) || relativeHumidity <= 0)
+                return double.NaN;
+
+            const double a = 17.625;
+            const double b = 243.04;
+
+            double value = Math.Log(Math.Min(relativeHumidity, 100.0) / 100.0) + a * dryBulbTemperature / (b + dryBulbTemperature);
+            return b * value / (a - value);
+        }
+
+        /// <summary>
+        /// Standard time zone [hours from UTC] stated on the weather data, or NaN when it is not stated.
+        /// </summary>
+        private static double EPWTimeZone(this WeatherData weatherData)
+        {
+            if (weatherData == null)
+                return double.NaN;
+
+            if (!weatherData.TryGetValue(WeatherDataParameter.TimeZone, out string timeZone) || string.IsNullOrWhiteSpace(timeZone))
+                return double.NaN;
+
+            return Core.Query.Double(Core.Query.UTC(timeZone));
+        }
+    }
+}

THIRD_PARTY.md

@@ -34,3 +34,10 @@
   listing comments and termination text. It has no GenOpt banner, copyright header, console transcript or code.
 - Distribution: no GenOpt code, `genopt.jar` or Java runtime is included in this repository or its builds. The oracle
   tool `SAM/SAM.Tests.GenOptOracle` runs a locally supplied `genopt.jar`, which is never committed.
+
+## TMYx weather files (test fixtures only)
+- Files: `SAM/SAM.Tests/Resources/Weather/*.epw.gz`, five TMYx EPW files from https://climate.onebuilding.org/,
+  gzip-compressed without modification. The SHA-256 of each original `.epw` is in `SHA256SUMS.txt` next to them, and
+  provenance per file is in that folder's `README.md`.
+- License: free to use with attribution: "Climatic data provided by climate.onebuilding.org".
+- Use: reference data for the EPW writer tests (`WeatherEPWExportTests`). They are not part of any build output.