Files
sunscope/assets/js/physics.js
T
fraxle e0712300e5 Beta 6
Now with concrete, vehicle and home temps and profiles
2026-05-14 22:57:42 +01:00

426 lines
22 KiB
JavaScript
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
// ════════════════════════════════════════════════════════════════════════
// physics.js — Physical constants and meteorological calculations.
//
// Exports:
// SIGMA, EPSILON_P, A_K, ALBEDO_GRASS (radiation constants)
// vaporPressureHpa(Ta, RH) Magnus formula → hPa
// solarElevationDeg(lat, lon, dateUTC) NOAA simplified solar position
// calcTmrt(Ta, dirRad, diffRad, globalRad, solElev) Mean radiant temp
// utciApprox(Ta, Tmrt, va10, ehPa) Bröde et al. 2012 polynomial
//
// Nothing in here should need editing unless the underlying science
// changes. All numbers are peer-reviewed constants or coefficients.
// ════════════════════════════════════════════════════════════════════════
// ═══════════════════════════════════════════════════════════════════
// PHYSICAL CONSTANTS
// ───────────────────────────────────────────────────────────────────
// These are real-world physics values — don't change them unless you
// have a peer-reviewed reason. They're used by calcTmrt() below to
// work out how much heat your skin actually absorbs from the sun.
// SIGMA — StefanBoltzmann constant (radiates heat)
// EPSILON_P — emissivity of human skin (~0.97)
// A_K — short-wave absorption coefficient for clothing
// ALBEDO_GRASS — how much sun grass reflects back at you (23%)
// ALBEDO_CONCRETE — how much sun concrete reflects back (30%)
// Concrete absorbs more net solar than grass and has
// no evaporative cooling, so its surface runs hot.
// ═══════════════════════════════════════════════════════════════════
export const SIGMA = 5.670374419e-8;
export const EPSILON_P = 0.97;
export const A_K = 0.7;
export const ALBEDO_GRASS = 0.23;
export const ALBEDO_CONCRETE = 0.30;
// ═══════════════════════════════════════════════════════════════════
// CONCRETE SURFACE TEMPERATURE (Urban profile)
// ───────────────────────────────────────────────────────────────────
// Estimates the surface temperature of exposed concrete using a
// simplified energy-balance approach:
// • Absorbed solar = globalRad × (1 albedo)
// • No latent heat (no evaporation — concrete is dry)
// • Convective loss to air proportional to wind speed
// • Result is clamped to a physically plausible range
//
// This is what matters for contact heat stress in cities — the UTCI
// standard uses grass, which runs ~515 °C cooler than urban concrete
// on a sunny day because grass sweats (transpires).
// ═══════════════════════════════════════════════════════════════════
export function calcConcreteTemp(Ta, globalRad, windSpeed) {
if (globalRad == null || Ta == null) return null;
const absorbed = globalRad * (1 - ALBEDO_CONCRETE); // W/m²
// Convective heat transfer coefficient: ~5 W/m²K still air, rises with wind
const hc = 5 + 4.5 * Math.sqrt(Math.max(windSpeed || 0, 0));
// Surface temp: Ta + solar gain / convective loss
const Ts = Ta + absorbed / hc;
// Clamp: can't be cooler than air, cap at 85 °C (melting asphalt territory)
return Math.max(Ta, Math.min(Ts, 85));
}
// ═══════════════════════════════════════════════════════════════════
// UK HOUSE INDOOR TEMPERATURE (windows closed, no active cooling)
// ───────────────────────────────────────────────────────────────────
// Estimates the ambient indoor air temperature of a typical UK brick
// house with windows closed and no air conditioning.
//
// Two heat pathways are modelled:
//
// 1. WALL CONDUCTION
// Heat conducts through brick cavity walls and roof. UK Part L
// compliant walls have a U-value around 0.280.45 W/m²K; older
// solid-brick stock runs higher. A representative mid-stock value
// is used. This drives a slow, steady heat transfer proportional
// to the difference between outdoor and indoor air temperature.
//
// 2. WINDOW SOLAR GAIN
// A typical UK semi has ~1518% glazing ratio. Solar energy
// transmits through glass, is absorbed by floors and furniture,
// and heats the indoor air. Gain is averaged across orientations
// (not all windows face south). Diffuse radiation contributes
// regardless of sun angle.
//
// THERMAL LAG
// Brick and concrete have high thermal mass — the house responds
// slowly to outdoor temperature swings. This function takes a
// weighted average of the current hour's heat load and the
// previous few hours', giving the characteristic lag where indoor
// temperature peaks 24 hours after the outdoor peak.
// Call calcIndoorTempPass() on the full hourly arrays after
// building rows — it returns a per-hour indoor temp array.
//
// No mechanical cooling. Minimal infiltration (windows closed).
// Internal heat gains (people, appliances) are not modelled.
//
// Colour thresholds:
// < 20 °C — cool, may need heating
// 2026 °C — comfortable
// 2632 °C — warm; WHO heatwave advisory threshold for sleeping
// > 32 °C — hot; risk for elderly and vulnerable occupants
// ═══════════════════════════════════════════════════════════════════
// Single-hour instantaneous heat load (W/m² effective, indoor side).
// Used internally by calcIndoorTempPass — not exported.
function _houseHeatLoad(Ta, globalRad, solElev) {
// Wall + roof conduction: U-value ~0.35 W/m²K × effective envelope area ratio
// We express as a gain-per-degree-delta — applied against Ti later in the pass.
// (See calcIndoorTempPass for how this feeds the lag model.)
// Window solar gain: glazing ratio 0.16, g-value 0.63 (standard double glazing),
// averaged across orientations (0.5 factor — not all windows face the sun).
const glazingRatio = 0.16;
const gValue = 0.63;
const orientFactor = 0.50;
const solarGain = globalRad * glazingRatio * gValue * orientFactor;
return { conductionDelta: Ta, solarGain };
}
// Two-pass function: call with the full arrays of hourly Ta and globalRad.
// Returns an array of indoor temperatures, one per hour.
export function calcIndoorTempPass(TaArr, globArr, elevArr) {
const n = TaArr.length;
const result = new Array(n);
// Thermal resistance of the building envelope (°C per W/m² of heat load)
// Lower = faster response to outdoor changes. UK brick mid-stock ~0.45.
const uWall = 0.35; // W/m²K effective wall U-value
// Thermal mass time constant: heavier = longer lag.
// ~4 h lag for typical UK brick semi (expressed as exponential decay weight).
const lagHours = 4;
const alpha = 1 - Math.exp(-1 / lagHours); // per-hour blending weight
// Seed indoor temp to first outdoor temp
let Ti = TaArr[0] ?? 15;
for (let i = 0; i < n; i++) {
const Ta = TaArr[i] ?? Ti;
const glob = globArr[i] ?? 0;
const solElev = elevArr[i] ?? 0;
// Solar gain through windows (W/m² effective)
const glazingRatio = 0.16;
const gValue = 0.63;
const orientFactor = 0.50;
const solarGain = glob * glazingRatio * gValue * orientFactor;
// Conductive heat flow through walls: proportional to (Ta - Ti)
// uWall drives how quickly the indoor temp chases outdoor temp.
const conductionGain = uWall * (Ta - Ti);
// Target indoor temp this hour if there were no thermal mass:
// Ti_instant = Ti + conduction + solar load / heat capacity proxy
// heatCap proxy: how many degrees does 1 W/m² raise the indoor air?
// For a typical 90 m² house, ~0.15 °C per W/m² effective.
const heatCapProxy = 0.15;
const Ti_instant = Ti + (conductionGain + solarGain) * heatCapProxy;
// Apply thermal lag: blend toward Ti_instant slowly
Ti = Ti + alpha * (Ti_instant - Ti);
// Can't be colder than outdoor (house doesn't actively cool)
// Cap at 55 °C (physically implausible above this for a house interior)
result[i] = Math.max(Math.min(Ti, 55), Math.min(Ta, Ti));
}
return result;
}
// ═══════════════════════════════════════════════════════════════════
// VEHICLE INTERIOR CABIN TEMPERATURE (seated occupant, not in sunbeam)
// ───────────────────────────────────────────────────────────────────
// Models the ambient cabin air temperature experienced by an occupant
// seated out of direct sunlight inside a sealed, parked vehicle.
// Two heat sources are combined:
//
// 1. PANEL CONDUCTION
// Aluminium body panels absorb solar radiation and conduct heat
// into the cabin. Albedo ~0.25 (mid-point for typical mixed-colour
// fleet; dark paint ~0.10, silver/white ~0.40).
// Panel surface temp → conductive gain into cabin air.
//
// 2. SIDE-WINDOW GLAZING GAIN (sun-angle dependent)
// When solar elevation is between ~10° and ~60°, the sun's rays
// cut through the side glass at an angle that allows significant
// transmission into the cabin (rather than hitting the roof or
// reflecting off at a shallow angle). This warms the cabin air
// but the occupant is modelled as NOT sitting in the beam —
// so it adds to ambient cabin temp, not direct radiant load.
// Above 60° the sun mostly hits the roof; below 10° it reflects.
//
// Wind is ignored (vehicle is sealed). No evaporative cooling.
// Steady-state reached after ~4560 min of parking.
//
// Colour thresholds in the table:
// < 35 °C — warm but tolerable for short periods
// 3545 °C — dangerous for children/pets (hyperthermia risk)
// > 45 °C — potentially fatal within minutes
// ═══════════════════════════════════════════════════════════════════
export function calcVehicleInteriorTemp(Ta, globalRad, solElev) {
if (globalRad == null || Ta == null) return null;
// ── 1. Panel conduction ──────────────────────────────────────────
const albedoPanel = 0.25; // typical mixed fleet
const panelAbsorbed = globalRad * (1 - albedoPanel); // W/m² absorbed by bodywork
// Panel surface temp: absorbed solar / convective loss to outside air
// hOut ~10 W/m²K (light breeze over panel surface even when parked)
const hOut = 10;
const panelSurfaceTemp = Ta + panelAbsorbed / hOut;
// Conductive gain into cabin: panel-to-cabin air, ~4 W/m²K through metal + trim
const hCabin = 4;
const conductionGain = hCabin * (panelSurfaceTemp - Ta); // W/m²
// ── 2. Side-window glazing gain (angle-dependent) ─────────────────
// Glazing transmission for auto glass ~0.70
const tau = 0.70;
let glazingGain = 0;
if (solElev != null && solElev > 10 && solElev < 60) {
// Scale factor: peaks around 3040° elevation (sun cuts squarely
// through side glass), tapers off toward 10° (shallow/reflected)
// and 60° (sun increasingly hitting roof not side glass).
// Use a simple tent function peaking at 35°.
const peak = 35;
const halfWidth = 25; // degrees either side
const factor = Math.max(0, 1 - Math.abs(solElev - peak) / halfWidth);
// Diffuse radiation also enters through glass regardless of angle
glazingGain = tau * globalRad * factor * 0.5; // occupant not in beam → 50% ambient
} else {
// Outside the side-window zone: diffuse only (scattered sky light)
glazingGain = tau * (globalRad * 0.15); // ~15% diffuse fraction
}
// ── Combine into cabin air temperature ───────────────────────────
// Total heat input per m² of cabin surface
const totalGain = conductionGain + glazingGain;
// Cabin heat loss: poor natural ventilation in sealed car ~2.0 W/m²K
const hCabinLoss = 2.0;
const Ti = Ta + totalGain / hCabinLoss;
// Clamp: can't be cooler than outside air; physical cap at 90 °C
return Math.max(Ta, Math.min(Ti, 90));
}
// Vapour pressure (Magnus → hPa)
export function vaporPressureHpa(Ta, RH) {
const es = 6.105 * Math.exp((17.27 * Ta) / (237.7 + Ta));
return es * (RH / 100);
}
// Solar elevation (NOAA simplified, degrees)
export function solarElevationDeg(lat, lon, dateUTC) {
const start = Date.UTC(dateUTC.getUTCFullYear(), 0, 0);
const diff = dateUTC - start;
const DOY = Math.floor(diff / 86400000);
const hourUTC =
dateUTC.getUTCHours() +
dateUTC.getUTCMinutes() / 60 +
dateUTC.getUTCSeconds() / 3600;
const gamma = ((2 * Math.PI) / 365) * (DOY - 1 + (hourUTC - 12) / 24);
const eqtime =
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));
const decl =
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);
const timeOffset = eqtime + 4 * lon;
const tst = hourUTC * 60 + timeOffset;
const ha = (((tst / 4) - 180) * Math.PI) / 180;
const latRad = (lat * Math.PI) / 180;
const cosZenith =
Math.sin(latRad) * Math.sin(decl) +
Math.cos(latRad) * Math.cos(decl) * Math.cos(ha);
const zenith = Math.acos(Math.max(-1, Math.min(1, cosZenith)));
return (Math.PI / 2 - zenith) * (180 / Math.PI);
}
// Mean radiant temperature
export function calcTmrt(Ta, dirRad, diffRad, globalRad, solElev) {
const TaK = Ta + 273.15;
let fp = 0;
if (solElev > 0) {
const h2 = solElev;
fp = 0.308 * Math.cos((Math.PI / 180) * h2 * (0.998 - (h2 * h2) / 50000));
}
let DNI = 0;
if (solElev > 1) {
DNI = dirRad / Math.sin((solElev * Math.PI) / 180);
DNI = Math.min(DNI, 1100);
}
const Sshort = A_K * (fp * DNI + 0.5 * diffRad + 0.5 * ALBEDO_GRASS * globalRad);
const Slong = EPSILON_P * SIGMA * Math.pow(TaK, 4);
const TmrtK = Math.pow((Sshort + Slong) / (EPSILON_P * SIGMA), 0.25);
return TmrtK - 273.15;
}
// ═══════════════════════════════════════════════════════════════════
// UTCI POLYNOMIAL — DO NOT EDIT (or be VERY careful if you do)
// ───────────────────────────────────────────────────────────────────
// This is the official 210-term Bröde et al. (2012) approximation.
// It takes air temp, mean radiant temp, wind, and humidity and returns
// the "felt" temperature. Every number is a peer-reviewed coefficient.
// Scroll past it — there's nothing here you'll want to change.
// ═══════════════════════════════════════════════════════════════════
export function utciApprox(Ta, Tmrt, va10, ehPa) {
const va = Math.max(0.5, Math.min(17, va10));
const D_Tmrt = Tmrt - Ta;
const Pa = ehPa / 10;
const T = Ta, V = va, D = D_Tmrt, P = Pa;
const T2=T*T, T3=T2*T, T4=T3*T, T5=T4*T, T6=T5*T;
const V2=V*V, V3=V2*V, V4=V3*V, V5=V4*V, V6=V5*V;
const D2=D*D, D3=D2*D, D4=D3*D, D5=D4*D, D6=D5*D;
const P2=P*P, P3=P2*P, P4=P3*P, P5=P4*P, P6=P5*P;
return T +
6.07562052e-1 +
-2.27712343e-2 * T + 8.06470249e-4 * T2 + -1.54271372e-4 * T3 +
-3.24651735e-6 * T4 + 7.32602852e-8 * T5 + 1.35959073e-9 * T6 +
-2.25836520e0 * V + 8.80326035e-2 * T*V + 2.16844454e-3 * T2*V +
-1.53347087e-5 * T3*V + -5.72983704e-7 * T4*V + -2.55090145e-9 * T5*V +
-7.51269505e-1 * V2 + -4.08350271e-3 * T*V2 + -5.21670675e-5 * T2*V2 +
1.94544667e-6 * T3*V2 + 1.14099531e-8 * T4*V2 +
1.58137256e-1 * V3 + -6.57263143e-5 * T*V3 + 2.22697524e-7 * T2*V3 +
-4.16117031e-8 * T3*V3 +
-1.27762753e-2 * V4 + 9.66891875e-6 * T*V4 + 2.52785852e-9 * T2*V4 +
4.56306672e-4 * V5 + -1.74202546e-7 * T*V5 +
-5.91491269e-6 * V6 +
3.98374029e-1 * D + 1.83945314e-4 * T*D + -1.73754510e-4 * T2*D +
-7.60781159e-7 * T3*D + 3.77830287e-8 * T4*D + 5.43079673e-10 * T5*D +
-2.00518269e-2 * V*D + 8.92859837e-4 * T*V*D + 3.45433048e-6 * T2*V*D +
-3.77925774e-7 * T3*V*D + -1.69699377e-9 * T4*V*D +
1.69992415e-4 * V2*D + -4.99204314e-5 * T*V2*D + 2.47417178e-7 * T2*V2*D +
1.07596466e-8 * T3*V2*D +
8.49242932e-5 * V3*D + 1.35191328e-6 * T*V3*D + -6.21531254e-9 * T2*V3*D +
-4.99410301e-6 * V4*D + -1.89489258e-8 * T*V4*D +
8.15300114e-8 * V5*D +
7.55043090e-4 * D2 + -5.65095215e-5 * T*D2 + -4.52166564e-7 * T2*D2 +
2.46688878e-8 * T3*D2 + 2.42674348e-10 * T4*D2 +
1.54547250e-4 * V*D2 + 5.24110970e-6 * T*V*D2 + -8.75874982e-8 * T2*V*D2 +
-1.50743064e-9 * T3*V*D2 +
-1.56236307e-5 * V2*D2 + -1.33895614e-7 * T*V2*D2 + 2.49709824e-9 * T2*V2*D2 +
6.51711721e-7 * V3*D2 + 1.94960053e-9 * T*V3*D2 +
-1.00361113e-8 * V4*D2 +
-1.21206673e-5 * D3 + -2.18203660e-7 * T*D3 + 7.51269482e-9 * T2*D3 +
9.79063848e-11 * T3*D3 +
1.25006734e-6 * V*D3 + -1.81584736e-9 * T*V*D3 + -3.52197671e-10 * T2*V*D3 +
-3.36514630e-8 * V2*D3 + 1.35908359e-10 * T*V2*D3 +
4.17032620e-10 * V3*D3 +
-1.30369025e-9 * D4 + 4.13908461e-10 * T*D4 + 9.22652254e-12 * T2*D4 +
-5.08220384e-9 * V*D4 + -2.24730961e-11 * T*V*D4 +
1.17139133e-10 * V2*D4 +
6.62154879e-10 * D5 + 4.03863260e-13 * T*D5 + 1.95087203e-12 * V*D5 +
-4.73602469e-12 * D6 +
5.12733497e0 * P + -3.12788561e-1 * T*P + -1.96701861e-2 * T2*P +
9.99690870e-4 * T3*P + 9.51738512e-6 * T4*P + -4.66426341e-7 * T5*P +
5.48050612e-1 * V*P + -3.30552823e-3 * T*V*P + -1.64119440e-3 * T2*V*P +
-5.16670694e-6 * T3*V*P + 9.52692432e-7 * T4*V*P +
-4.29223622e-2 * V2*P + 5.00845667e-3 * T*V2*P + 1.00601257e-6 * T2*V2*P +
-1.81748644e-6 * T3*V2*P +
-1.25813502e-3 * V3*P + -1.79330391e-4 * T*V3*P + 2.34994441e-6 * T2*V3*P +
1.29735808e-4 * V4*P + 1.29064870e-6 * T*V4*P +
-2.28558686e-6 * V5*P +
-3.69476348e-2 * D*P + 1.62325322e-3 * T*D*P + -3.14279680e-5 * T2*D*P +
2.59835559e-6 * T3*D*P + -4.77136523e-8 * T4*D*P +
8.64203390e-3 * V*D*P + -6.87405181e-4 * T*V*D*P + -9.13863872e-6 * T2*V*D*P +
5.15916806e-7 * T3*V*D*P +
-3.59217476e-5 * V2*D*P + 3.28696511e-5 * T*V2*D*P + -7.10542454e-7 * T2*V2*D*P +
-1.24382300e-5 * V3*D*P + -7.38584400e-9 * T*V3*D*P +
2.20609296e-7 * V4*D*P +
-7.32469180e-4 * D2*P + -1.87381964e-5 * T*D2*P + 4.80925239e-6 * T2*D2*P +
-8.75492040e-8 * T3*D2*P +
2.77862930e-5 * V*D2*P + -5.06004592e-6 * T*V*D2*P + 1.14325367e-7 * T2*V*D2*P +
2.53016723e-6 * V2*D2*P + -1.72857035e-8 * T*V2*D2*P +
-3.95079398e-8 * V3*D2*P +
-3.59413173e-7 * D3*P + 7.04388046e-7 * T*D3*P + -1.89309167e-8 * T2*D3*P +
-4.79768731e-7 * V*D3*P + 7.96079978e-9 * T*V*D3*P +
1.62897058e-9 * V2*D3*P +
3.94367674e-8 * D4*P + -1.18566247e-9 * T*D4*P +
3.34678041e-10 * V*D4*P +
-1.15606447e-10 * D5*P +
-2.80626406e0 * P2 + 5.48712484e-1 * T*P2 + -3.99428410e-3 * T2*P2 +
-9.54009191e-4 * T3*P2 + 1.93090978e-5 * T4*P2 +
-3.08806365e-1 * V*P2 + 1.16952364e-2 * T*V*P2 + 4.95271903e-4 * T2*V*P2 +
-1.90710882e-5 * T3*V*P2 +
2.10787756e-3 * V2*P2 + -6.98445738e-4 * T*V2*P2 + 2.30109073e-5 * T2*V2*P2 +
4.17856590e-4 * V3*P2 + -1.27043871e-5 * T*V3*P2 +
-3.04620472e-6 * V4*P2 +
5.14507424e-2 * D*P2 + -4.32510997e-3 * T*D*P2 + 8.99281156e-5 * T2*D*P2 +
-7.14663943e-7 * T3*D*P2 +
-2.66016305e-4 * V*D*P2 + 2.63789586e-4 * T*V*D*P2 + -7.01199003e-6 * T2*V*D*P2 +
-1.06823306e-4 * V2*D*P2 + 3.61341136e-6 * T*V2*D*P2 +
2.29748967e-7 * V3*D*P2 +
3.04788893e-4 * D2*P2 + -6.42070836e-5 * T*D2*P2 + 1.16257971e-6 * T2*D2*P2 +
7.68023384e-6 * V*D2*P2 + -5.47446896e-7 * T*V*D2*P2 +
-3.59937910e-8 * V2*D2*P2 +
-4.36497725e-6 * D3*P2 + 1.68737969e-7 * T*D3*P2 +
2.67489271e-8 * V*D3*P2 +
3.23926897e-9 * D4*P2 +
-3.53874123e-2 * P3 + -2.21201190e-1 * T*P3 + 1.55126038e-2 * T2*P3 +
-2.63917279e-4 * T3*P3 +
4.53433455e-2 * V*P3 + -4.32943862e-3 * T*V*P3 + 1.45389826e-4 * T2*V*P3 +
2.17508610e-4 * V2*P3 + -6.66724702e-5 * T*V2*P3 +
3.33217140e-5 * V3*P3 +
-2.26921615e-3 * D*P3 + 3.80261982e-4 * T*D*P3 + -5.45314314e-9 * T2*D*P3 +
-7.96355448e-4 * V*D*P3 + 2.53458034e-5 * T*V*D*P3 +
-6.31223658e-6 * V2*D*P3 +
3.02122035e-4 * D2*P3 + -4.77403547e-6 * T*D2*P3 +
1.73825715e-6 * V*D2*P3 +
-4.09087898e-7 * D3*P3 +
6.14155345e-1 * P4 + -6.16755931e-2 * T*P4 + 1.33374846e-3 * T2*P4 +
3.55375387e-3 * V*P4 + -5.13027851e-4 * T*V*P4 +
1.02449757e-4 * V2*P4 +
-1.48526421e-3 * D*P4 + -4.11469183e-5 * T*D*P4 +
-6.80434415e-6 * V*D*P4 +
-9.77675906e-6 * D2*P4 +
8.82773108e-2 * P5 + -3.01859306e-3 * T*P5 +
1.04452989e-3 * V*P5 +
2.47090539e-4 * D*P5 +
1.48348065e-3 * P6;
}