Files
sunscope/assets/js/physics.js
T
fraxle 862f631beb 5.5
New big and small scope sunset and sunrise simulation
New calculations for the day tabs
Added missing icons to simple view
New record popups
Week ahead description popups
Camping profile
Tweak to some of the indoor tems
2026-10-08 09:53:41 +01:00

1042 lines
53 KiB
JavaScript

// ------------------------------------------------------------------------
// physics.js - Physical constants and meteorological calculations.
//
// All numbers are peer-reviewed constants or coefficients. Nothing in
// here should need editing unless the underlying science changes.
//
// Exports (in order of appearance):
// SIGMA, EPSILON_P, A_K, ALBEDO_GRASS, ALBEDO_CONCRETE (radiation constants)
// LAG_HOURS_CONCRETE slab thermal time constant
// calcConcreteTemp(Ta, globalRad, windSpeed, ...) urban surface temp (instantaneous target)
// calcConcreteTempPass(arrays...) thermal-lag pass over full hourly arrays
// calcIndoorTempPass(TaArr, globArr, elevArr, buildingType) passive indoor temp
// calcManagedIndoorTempPass(TaArr, globArr, elevArr, buildingType) managed indoor temp
// calcVehicleInteriorTemp(Ta, globalRad, solElev, vehicleType, ventilated, speedMph, windMps) equilibrium cabin temp
// calcVehicleInteriorTempPass(TaArr, globArr, elevArr, vaArr, vehicleType, ventilated, speedMph) lagged cabin temp
// calcShadeAirTemp(TaEnv, effRadEnv, va, elev, shade) per-env shade air temp
// calcShadeFeltTemp(TaShade, difEnv, va, eh, elev, shade) felt temp in that shade
// vaporPressureHpa(Ta, RH) Magnus formula - hPa
// solarElevationDeg(lat, lon, dateUTC) NOAA simplified solar position (-)
// calcTmrt(Ta, dirRad, diffRad, globalRad, solElev) Mean radiant temperature
// utciApprox(Ta, Tmrt, va10, ehPa) Br-de et al. 2012 UTCI polynomial
// ------------------------------------------------------------------------
import { VEHICLE_TYPES, BUILDING_TYPES } from './utils.js';
// -------------------------------------------------------------------
// 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 - Stefan-Boltzmann 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
// multi-factor energy-balance approach. Six physical effects are
// modelled beyond the basic solar-gain / convective-loss pair:
//
// 1. SUN ANGLE CORRECTION
// At low solar elevation the sun hits the surface obliquely,
// spreading energy over a larger area. A sin(elev) factor
// reduces absorbed solar proportionally. Clamped at a 5-deg
// minimum so the result stays finite near the horizon.
//
// 2. CLOUD TYPE TRANSMITTANCE
// Cloud cover category drives a transmittance multiplier.
// Low cloud (stratus) is far more opaque than high cirrus:
// clear 1.00 - full beam reaches the surface
// wispy 0.92 - cirrus barely attenuates
// scattered 0.72 - broken cumulus, significant blocking
// overcast 0.28 - thick stratus, mostly diffuse remains
// The raw effectiveRad already dims with cloud cover from the
// API, but this adds the qualitative distinction between cloud
// types that the single radiation number does not capture.
//
// 3. UV CLARITY FACTOR
// UV index is a proxy for atmospheric clarity beyond cloud cover -
// aerosols, haze, and humidity all reduce it. A UV of 8+ indicates
// a very clean, dry atmosphere with maximum direct-beam intensity.
// Normalised to a 0.85-1.00 range so it modulates rather than
// dominates. No UV data defaults to neutral (1.0).
//
// 4. EVAPORATIVE COOLING FROM SOIL MOISTURE
// Wet concrete loses heat via evaporation. Soil moisture at 0-1cm
// is used as a proxy for surface wetness (0 = bone dry, 1 = fully
// saturated). A saturated surface loses up to ~8-C relative to
// the dry case - consistent with published wet-pavement studies.
//
// 5. RAIN-WET SURFACE
// Active precipitation forces surface wetness regardless of soil
// moisture data. Above 0.5 mm/h the surface is considered fully
// wet and the maximum evaporative penalty applies.
//
// 6. SNOW COVER
// Snow on concrete insulates the slab from solar gain AND strongly
// reflects incoming radiation (albedo ~0.80 for fresh snow vs 0.30
// for bare concrete). When snowfall is active or lying snow is
// implied (snow > 0), absorbed radiation is cut by 85% and a small
// insulating offset is applied instead.
//
// Colour thresholds (same as before - surface contact risk):
// < Ta - should not occur (clamped)
// Ta - 45 -C - warm but bearable contact
// 45 - 60 -C - pain threshold for bare skin contact
// > 60 -C - burns on contact (relevant for paws / bare feet)
// -------------------------------------------------------------------
export function calcConcreteTemp(Ta, globalRad, windSpeed, solElev, uv, cloudCat, soilM, precip, snow) {
if (globalRad == null || Ta == null) return null;
// -- 6. Snow short-circuit --------------------------------------------
// Snow-covered concrete behaves like a white reflective insulator.
// Absorbed solar collapses; slab temp stays close to air temp.
if (snow != null && snow > 0) {
const snowAbsorbed = globalRad * (1 - 0.80); // fresh snow albedo ~0.80
const hcSnow = 10 + 4.5 * Math.sqrt(Math.max(windSpeed || 0, 0));
const Ts = Ta + snowAbsorbed / hcSnow;
return Math.max(Ta - 1, Math.min(Ts, 40)); // snow-covered slab rarely exceeds 40-C
}
// -- 1. Sun angle correction ------------------------------------------
// Low-angle sun spreads energy across a larger surface area.
// sin(elev) = 1.0 at 90-deg (overhead), ~0.17 at 10-deg (grazing).
// Default to sin(45-deg) ~0.71 when elevation is unknown.
const elevDeg = (solElev != null) ? Math.max(5, solElev) : 45;
const angleCorrection = Math.sin(elevDeg * Math.PI / 180);
// -- 2. Cloud type transmittance --------------------------------------
// Modulates beam quality beyond what raw radiation already captures.
const cloudTransmit = cloudCat === 'clear' ? 1.00
: cloudCat === 'wispy' ? 0.92
: cloudCat === 'scattered' ? 0.72
: /* overcast */ 0.28;
// -- 3. UV clarity factor --------------------------------------------
// UV index as atmospheric-clarity proxy. Scaled to 0.85-1.00 range.
// uv=0 (night or heavy cloud) - neutral 1.0 (no adjustment needed,
// radiation already near-zero). uv=8+ - max clarity bonus of 1.0.
const uvClarity = (uv && uv > 0)
? 0.85 + 0.15 * Math.min(1, uv / 8)
: 1.0;
// -- Absorbed solar with all modifiers --------------------------------
const absorbed = globalRad * (1 - ALBEDO_CONCRETE)
* angleCorrection
* cloudTransmit
* uvClarity;
// -- Convective loss --------------------------------------------------
const hc = 10 + 4.5 * Math.sqrt(Math.max(windSpeed || 0, 0));
// -- Dry surface temperature ------------------------------------------
const Ts = Ta + absorbed / hc;
// -- 4 + 5. Evaporative cooling ---------------------------------------
// Rain-wet surface overrides soil moisture - surface is fully saturated.
const surfaceWet = (precip != null && precip >= 0.5)
? 1.0
: Math.max(0, Math.min(1, soilM ?? 0));
// Max evaporative delta ~8-C at full saturation (published wet-pavement data).
const evapCooling = surfaceWet * 8;
const TsCooled = Ts - evapCooling;
// -- Final clamp: no cooler than air, no hotter than 85-C -------------
return Math.max(Ta, Math.min(TsCooled, 85));
}
// -------------------------------------------------------------------
// CONCRETE THERMAL LAG TIME CONSTANT
// -------------------------------------------------------------------
// A standard urban pavement slab (~100 mm thick, exposed top surface,
// air below via sub-base) has moderate thermal mass. Real-world
// measurement studies put the e-folding time constant at 1.5-2 h for
// this geometry - we use 1.5 h as representative of the thin end of
// typical footway construction (block paving, tarmac-over-hardcore).
//
// If this ever needs to vary by surface type, pull it into a
// SURFACE_TYPES preset object mirroring BUILDING_TYPES in utils.js.
// For now a single named constant keeps the intent obvious.
// -------------------------------------------------------------------
export const LAG_HOURS_CONCRETE = 1.5;
// -------------------------------------------------------------------
// CONCRETE SURFACE TEMPERATURE - THERMAL LAG PASS
// -------------------------------------------------------------------
// Wraps calcConcreteTemp in the same exponential-blending pattern used
// by calcIndoorTempPass. Instead of snapping to the instantaneous
// target each hour, the slab temperature blends toward it at a rate
// controlled by LAG_HOURS_CONCRETE.
//
// Effect in practice:
// - A slab baking at 55-C when cloud rolls in will still read ~48-C
// an hour later and ~44-C two hours later - not instantly 22-C.
// - A cold slab at dawn takes 2-3 hours of strong sun to fully heat.
// - Post-rain cool-down persists into the next hour even if it stops.
//
// Call AFTER building hourlyRows (same pattern as calcIndoorTempPass).
// Returns a Float64-like plain Array of concrete temps, one per hour.
//
// Arguments are parallel arrays (one value per forecast hour):
// TaArr - air temperature (-C)
// radArr - effectiveRad (dir + dif*0.2) (W/m-)
// vaArr - wind speed (m/s)
// elevArr - solar elevation (degrees)
// uvArr - UV index
// cloudCatArr - cloud category string
// soilMArr - soil moisture 0-1cm (0-1 fraction)
// precipArr - precipitation (mm/h)
// snowArr - snowfall (cm/h)
// -------------------------------------------------------------------
export function calcConcreteTempPass(TaArr, radArr, vaArr, elevArr, uvArr, cloudCatArr, soilMArr, precipArr, snowArr) {
const n = TaArr.length;
const result = new Array(n);
const alpha = 1 - Math.exp(-1 / LAG_HOURS_CONCRETE);
// Seed from the first hours instantaneous value so we start somewhere
// physically reasonable rather than zero.
let Tc = calcConcreteTemp(
TaArr[0], radArr[0], vaArr[0],
elevArr[0], uvArr[0], cloudCatArr[0],
soilMArr[0], precipArr[0], snowArr[0]
) ?? TaArr[0];
for (let i = 0; i < n; i++) {
const target = calcConcreteTemp(
TaArr[i], radArr[i], vaArr[i],
elevArr[i], uvArr[i], cloudCatArr[i],
soilMArr[i], precipArr[i], snowArr[i]
) ?? TaArr[i];
// Blend slab temp toward this hours target.
Tc = Tc + alpha * (target - Tc);
// Slab cannot be cooler than air (no active cooling mechanism).
result[i] = Math.max(TaArr[i], Tc);
}
return result;
}
// -------------------------------------------------------------------
// UK HOUSE INDOOR TEMPERATURE - TWO-NODE THERMAL MODEL
// -------------------------------------------------------------------
// Estimates the free-running indoor air temperature of the selected
// BUILDING_TYPES preset when EMPTY: no heating, no air conditioning, no
// occupants - only standby appliances (fridge, router) as internal gain.
//
// Each hour builds an outdoor-driven TARGET - where the room would settle
// if this hour's conditions held forever:
//
// target = outdoor ref + solar lift + standby gain
//
// outdoor ref blend of outside air, the ground and heated neighbours:
// Ta x (1 - groundCoupling - neighbourCoupling)
// + ground temp x groundCoupling
// + neighbour temp x neighbourCoupling
// Heat lost through the floor goes to the ground, not the
// air, and UK ground under a building sits near 6.5-14.5 C
// through the year (ukGroundTempC). This is what keeps an
// empty house above freezing in a cold snap and a little
// cooler in a heatwave - with a FIXED gap per degree, not a
// thermostat-like one.
// Neighbours (party walls/floors of terraces and flats) are
// occupied homes: never below NEIGHBOUR_MIN_C (they heat), and
// otherwise assumed the same as this building, so no flow.
// solar lift glob x glazingRatio x gValue x orientFactor x solarScale
// (tent: glazingRatio stands in for the sunlit fabric)
// standby gain fridge/router/standby loads in an empty building (C)
//
// The room is then modelled as TWO coupled nodes rather than one lag:
//
// 1. ROOM AIR (Ti) - fast node, time constant lagHours. Its equilibrium is
// a blend of the outdoor-driven target and the thermal mass:
// airEq = (1 - massCoupling) x target + massCoupling x Tm
// massCoupling is how strongly the walls/floors hold the air compared
// with the outdoor pull through glass, fabric and leakage.
//
// 2. THERMAL MASS (Tm) - slow node, time constant massHours. Walls, floors
// and furniture exchange heat with the room air.
//
// Why two nodes: a single exponential lag lets ~70% of the outdoor daily
// swing through for a 4 h brick constant, so a brick house swung 21-34 C on
// an 18-32 C day. Monitored UK homes swing roughly a third of the outdoor
// range and instead RATCHET UP over a multi-day heatwave as the mass warms -
// the slow node reproduces both. The effective outdoor time constant of the
// mass is massHours / (1 - massCoupling) - about 1.5 days for brick.
//
// The forecast starts at midnight with no history, so the model is spun up
// by cycling the first day SPIN_UP_DAYS times (assumes the days before the
// forecast resembled the first one) instead of starting the mass cold.
//
// Indoor CAN be cooler than outdoor (a stone cottage on a hot afternoon);
// the only clamp is a 55 C physical cap.
//
// Colour thresholds:
// < 20 -C - cool, may need heating
// 20-26 -C - comfortable
// 26-32 -C - warm; WHO heatwave advisory threshold for sleeping
// > 32 -C - hot; risk for elderly and vulnerable occupants
// -------------------------------------------------------------------
// Days of the first forecast day replayed to warm up the thermal mass.
const SPIN_UP_DAYS = 5;
// Managed mode switches to the heatwave routine on any day forecast to reach
// this outdoor temperature, or whenever the room itself passes HOT_INDOOR_C.
const HOT_DAY_OUTDOOR_C = 24;
const HOT_INDOOR_C = 23;
// Windows are only opened to cool while the room is above this, so night
// ventilation never chills the house below a comfortable sleeping temp.
const VENT_FLOOR_C = 20;
// UK ground temperature at ~1 m under a building for a day of year.
// Annual mean ~10.5 C, swing +/-4 C, coldest early March, warmest early
// September (the ground lags the air by ~6 weeks). UK/northern-hemisphere
// climate assumed, matching the rest of the building model.
const GROUND_MEAN_C = 10.5;
const GROUND_SWING_C = 4;
const GROUND_COLDEST_DOY = 65;
// Occupied neighbours keep their homes at least this warm in cold weather.
const NEIGHBOUR_MIN_C = 18;
function ukGroundTempC(dayOfYear) {
return GROUND_MEAN_C - GROUND_SWING_C * Math.cos(2 * Math.PI * (dayOfYear - GROUND_COLDEST_DOY) / 365);
}
// One hour of the indoor model. Mutates and returns state { Ti, Tm }.
// managed - false: windows closed, curtains open all the time (passive)
// true : occupant follows the routine for the conditions -
// HEAT (hotDay): curtains closed while the sun is up,
// windows opened whenever outdoor air is cooler
// COOL: windows shut, curtains open by day to let the sun
// in, closed after dark to cut window heat loss
function indoorStep(state, Ta, glob, preset, managed, hotDay, groundTemp) {
const { lagHours, massHours, massCoupling, glazingRatio, gValue, orientFactor,
curtainBlock, curtainInsul, ventAlpha, solarScale, groundCoupling,
neighbourCoupling, internalGain } = preset;
const heatMode = managed && (hotDay || state.Ti >= HOT_INDOOR_C);
const coolMode = managed && !heatMode;
// Solar gain - only blocked when managing for heat and the sun is up.
const block = heatMode ? curtainBlock : 0;
const solarLift = glob * glazingRatio * gValue * orientFactor * solarScale * (1 - block);
// Outdoor pull is split between outside air, the ground under the floor
// and any heated neighbours through party walls/floors.
const neighbourTemp = Math.max(state.Ti, NEIGHBOUR_MIN_C);
const outdoorRef = (1 - groundCoupling - neighbourCoupling) * Ta
+ groundCoupling * groundTemp
+ neighbourCoupling * neighbourTemp;
let target = outdoorRef + solarLift + internalGain;
// Cool routine after dark: closed curtains cut window heat loss, so the
// room is pulled less hard toward outdoor - curtainInsul is the share of
// the outdoor heat loss they save.
if (coolMode && glob <= 0 && target < state.Ti) {
target += (state.Ti - target) * curtainInsul;
}
// Room air relaxes toward a blend of the outdoor pull and the mass.
const airEq = (1 - massCoupling) * target + massCoupling * state.Tm;
const aAir = 1 - Math.exp(-1 / lagHours);
const aMass = 1 - Math.exp(-1 / massHours);
state.Ti += aAir * (airEq - state.Ti);
// Heat routine: open windows only when outside is cooler and the room is
// still above the comfort floor; never ventilate below that floor.
if (heatMode && Ta < state.Ti && state.Ti > VENT_FLOOR_C) {
const vented = state.Ti + ventAlpha * (Ta - state.Ti);
state.Ti = Math.max(vented, VENT_FLOOR_C, Ta);
}
// Mass exchanges with the (possibly ventilated) room air - this is how
// night purging cools the fabric for the next day.
state.Tm += aMass * (state.Ti - state.Tm);
return state;
}
// Shared driver for the passive and managed passes.
// Hourly arrays are assumed to start at local midnight (forecast_days fetch),
// so hour i belongs to forecast day floor(i / 24).
// startIso is the first hour's ISO time (for the seasonal ground temp); if
// missing the annual mean ground temperature is used.
function runIndoorModel(TaArr, globArr, buildingType, managed, startIso) {
const preset = BUILDING_TYPES[buildingType] || BUILDING_TYPES.brick;
const n = TaArr.length;
const result = new Array(n);
// Ground changes only ~0.5 C over a 14-day forecast - one value is enough.
const start = startIso ? new Date(startIso.slice(0, 10) + 'T00:00:00Z') : null;
const doy = start && !isNaN(start)
? Math.floor((start - Date.UTC(start.getUTCFullYear(), 0, 1)) / 86400000) + 1
: null;
const groundTemp = doy != null ? ukGroundTempC(doy) : GROUND_MEAN_C;
// Daily outdoor maxima decide the managed routine for each forecast day.
const dayMax = [];
for (let i = 0; i < n; i++) {
const d = Math.floor(i / 24);
const Ta = TaArr[i];
if (Ta != null) dayMax[d] = Math.max(dayMax[d] ?? -Infinity, Ta);
}
const isHotDay = i => (dayMax[Math.floor(i / 24)] ?? -Infinity) >= HOT_DAY_OUTDOOR_C;
// Seed both nodes at the first hour's air temp, then spin up on day one.
const seed = TaArr[0] ?? groundTemp;
const state = { Ti: seed, Tm: seed };
const spinHours = Math.min(24, n);
for (let rep = 0; rep < SPIN_UP_DAYS; rep++) {
for (let i = 0; i < spinHours; i++) {
indoorStep(state, TaArr[i] ?? state.Ti, globArr[i] ?? 0, preset, managed, isHotDay(i), groundTemp);
}
}
for (let i = 0; i < n; i++) {
const Ta = TaArr[i] ?? state.Ti;
indoorStep(state, Ta, globArr[i] ?? 0, preset, managed, isHotDay(i), groundTemp);
result[i] = Math.min(state.Ti, 55);
}
return result;
}
// PASSIVE: empty building, windows closed, curtains open, no heating or cooling.
// Returns an array of indoor temperatures, one per hour.
// buildingType must be a key of BUILDING_TYPES; defaults to 'brick'.
// elevArr is accepted for signature compatibility but not needed - the
// orientFactor average already covers sun angle on the glazing.
// startIso (optional) - first hour's ISO time, sets the seasonal ground temp.
export function calcIndoorTempPass(TaArr, globArr, elevArr, buildingType = 'brick', startIso) {
return runIndoorModel(TaArr, globArr, buildingType, false, startIso);
}
// -------------------------------------------------------------------
// UK HOUSE INDOOR TEMPERATURE - MANAGED
// -------------------------------------------------------------------
// The same empty building, but with its curtains and windows set sensibly
// for the conditions (by a visiting owner or smart blinds - see indoorStep):
//
// HEAT routine (day forecast >= HOT_DAY_OUTDOOR_C, or room >= HOT_INDOOR_C)
// Standard UK heatwave advice. Curtains/blinds closed while the sun is up
// (curtainBlock - internal curtains only stop ~40-55% because the heat
// is already through the glass; external shutters would be ~80%+).
// Windows opened whenever outdoor air is cooler than inside, down to a
// VENT_FLOOR_C comfort floor - so it never makes things worse.
//
// COOL routine (all other days)
// Windows shut. Curtains open in daylight to bank free solar gain,
// closed after dark to cut window heat loss (curtainInsul). Managed is
// therefore slightly WARMER than passive on cool days, not colder.
// -------------------------------------------------------------------
export function calcManagedIndoorTempPass(TaArr, globArr, elevArr, buildingType = 'brick', startIso) {
return runIndoorModel(TaArr, globArr, buildingType, true, startIso);
}
// -------------------------------------------------------------------
// 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. VERTICAL GLAZING GAIN (sun-angle dependent)
// Sun cuts through the side glass and windscreen and 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. The beam on
// vertical glass goes as cos(elevation) averaged over azimuth, over a
// diffuse floor that applies at any sun angle. A parked vehicle is not
// aimed at the sun, so only orientFactor of the glazing is catching it.
//
// 3. HORIZONTAL GLAZING GAIN (rooflights / panoramic roofs)
// Rooflights collect most strongly when the sun is high - the exact
// condition under which the vertical-glass term is tailing off. Scales
// with sin(elevation) and uses a lower transmission (0.40) because
// rooflights are typically smoked/tinted acrylic rather than clear glass.
//
// Ambient wind scrubs the bodywork whether or not the windows are open, so
// it feeds the shell convection coefficient in every case; road speed does
// the same job and the larger of the two wins. No evaporative cooling.
// Cars warm quickly; motorhomes and caravans are treated as insulated living
// spaces with 25-35 mm sandwich panels, so panel heat gain is much smaller.
// Their interiors are SLOWER, not COOLER - the insulation that keeps heat out
// also keeps it in, so a closed-up van ends up about as hot as a closed-up car
// once it has had a few hours to get there. Because they are occupied living
// spaces, they also retain a little warmth from previous hours, people,
// appliances, and background heating.
//
// Colour thresholds in the table:
// < 35 -C - warm but tolerable for short periods
// 35-45 -C - dangerous for children/pets (hyperthermia risk)
// > 45 -C - potentially fatal within minutes
// -------------------------------------------------------------------
export function calcVehicleInteriorTemp(Ta, globalRad, solElev, vehicleType = 'car', ventilated = false, speedMph = 0, windMps = 0) {
if (globalRad == null || Ta == null) return null;
// Look up vehicle preset; fall back to a standard car if key unknown.
const preset = VEHICLE_TYPES[vehicleType] || VEHICLE_TYPES.car;
// Road speed in m/s. speedMph = 0 is "Static" (parked); higher speeds scrub
// the shell with forced airflow and (windows down) flush the cabin toward
// ambient.
const vMs = Math.max(0, speedMph) * 0.447; // mph -> m/s
// -- 1. Panel conduction ------------------------------------------
const panelAbsorbed = globalRad * (1 - preset.albedo); // W/m- absorbed by bodywork
// Panel surface temp: absorbed solar / convective loss to outside air.
// Convection over the shell is driven by whichever airflow is stronger -
// the ambient wind or the car's own road speed. A parked vehicle in a
// 15 mph wind has bodywork far closer to ambient than one in still air,
// so ambient wind must not be ignored just because the windows are shut.
// The 1.5 m/s floor is not fudge: hOut is steep near calm, and a reported
// wind of 0 m/s does not mean still air at the bodywork - thermal plumes off
// hot panels and ordinary gustiness keep it moving. Without the floor a dead
// calm hour sent a sealed car to ~67 -C. It only bites below ~3 mph, so it
// leaves the calibration anchors untouched.
const vShell = Math.max(1.5, vMs, Math.max(0, windMps || 0));
const hOut = 10 + 5.5 * Math.sqrt(vShell);
const panelSurfaceTemp = Ta + panelAbsorbed / hOut;
// Conductive gain into cabin. Cars are thin metal + trim; motorhomes and
// caravans use insulated sandwich panels, so their bodyU is much lower.
const hCabin = preset.bodyU ?? 4;
const conductionGain = hCabin * (panelSurfaceTemp - Ta); // W/m-
// -- 2. Vertical glazing gain (angle-dependent) --------------------
// Glazing transmission for auto glass ~0.70; scaled by vehicle glazing area.
//
// The beam landing on VERTICAL glass goes as cos(elevation) once averaged
// over azimuth, plus a diffuse floor of scattered sky light that gets in at
// any sun angle. orientFactor then accounts for the glazing that is NOT
// pointing at the sun - a parked vehicle is not aimed, and a motorhome's
// windscreen faces wherever it happened to park.
//
// This replaced a tent function that peaked at 35- elevation. That shape
// assumed the glass was always squarely aimed at the sun, so its factor
// nearly TRIPLED between 3pm and 5pm as the sun dropped toward the peak -
// sending a ventilated van climbing to +9 over ambient in the late
// afternoon when the measured excess stays flat around +5..6. It also had
// a discontinuity at its 10-/60- cut-offs where a LOWER sun gave MORE gain.
const tau = 0.70 * preset.glazingArea;
const angleFactor = (solElev != null && solElev > 0)
? 0.30 + 0.70 * Math.cos(solElev * Math.PI / 180)
: 0.30;
const glazingGain = tau * globalRad * angleFactor * (preset.orientFactor ?? 0.65) * 0.5;
// -- 3. Horizontal glazing gain (rooflights, panoramic roof) -------
// Rooflights collect in proportion to sin(elevation), so they peak at midday
// just as the vertical-glass term is falling away. Smoked acrylic transmits
// roughly 0.40 rather than the 0.70 of clear auto glass.
const roofGain = (solElev != null && solElev > 0)
? 0.40 * (preset.roofGlazing ?? 0) * globalRad * Math.sin(solElev * Math.PI / 180)
: 0;
// -- Combine into cabin air temperature ---------------------------
// Total heat input per m- of cabin surface
const totalGain = conductionGain + glazingGain + roofGain;
// Cabin heat rejection: an effective blend of leakage, internal air volume,
// and surfaces exchanging heat with the outside. Opening the windows when
// parked multiplies it by the preset's ventMult - a car with every window
// down flushes far harder per unit volume than a van with two windows and a
// rooflight open, so that multiplier is per-vehicle, not a shared constant.
// On the move the through-draught multiplies this further (a 70 mph open
// window flushes the cabin almost to ambient). Sealed but moving rejects a
// little faster too, because the cooler shell pulls cabin heat out.
const speedFactor = 1 + vMs / 12; // grows with road speed (windows-open draught)
const lossMult = ventilated ? (preset.ventMult ?? 4) * speedFactor : 1 + vMs / 40;
const effectiveHLoss = preset.hCabinLoss * lossMult;
const solarRise = totalGain / effectiveHLoss;
// Motorhomes/caravans behave more like small insulated rooms than parked
// cars. This term captures retained living-space warmth: strongest on cool
// days, tapering away as outdoor air warms, and reduced when ventilated.
// It used to be much larger (8 - 0.25*Ta) because it was silently standing
// in for solar gain the model was throwing away; now that the glazing terms
// are right it only has to cover occupancy and residual warmth.
const retainedWarmth = preset.retainedWarmth
? Math.max(0, 5 - 0.20 * Ta) * (ventilated ? 0.35 : 1)
: 0;
const internalGain = (preset.internalGain ?? 0) * (ventilated ? 0.35 : 1);
const Ti = Ta + solarRise + retainedWarmth + internalGain;
// Clamp: can't be cooler than outside air; physical cap at 90 -C
return Math.max(Ta, Math.min(Ti, 90));
}
// -------------------------------------------------------------------
// VEHICLE INTERIOR TEMPERATURE - TWO-PASS (thermal lag)
// -------------------------------------------------------------------
// calcVehicleInteriorTemp above returns the EQUILIBRIUM cabin temperature
// for one hour's conditions - where the interior would settle if those
// conditions held. Real cabins take time to get there, so this pass relaxes
// toward that target with a per-vehicle time constant, exactly as
// calcConcreteTempPass and calcIndoorTempPass do.
//
// This replaces an earlier 'thermalMass' multiplier that scaled the
// equilibrium rise down (motorhomes ran at 0.35x). That conflated two
// different things: thermal mass delays how fast you reach equilibrium, it
// does not lower the equilibrium itself. A van parked in the sun since
// breakfast is close to equilibrium by mid-afternoon, so the multiplier
// under-predicted every long parked spell - badly enough that a ventilated
// motorhome came out barely a degree above ambient at peak sun.
//
// Cars use a short constant (~0.5 h - a car is hot within the half hour);
// motorhomes and caravans a longer one (~1.5-1.7 h) reflecting their larger
// air volume and heavier interior fit-out.
//
// vaArr is ambient wind in m/s (Open-Meteo wind_speed_10m with
// wind_speed_unit=ms) and may be null/omitted, in which case still air is
// assumed. speedMph is the selected road-speed class, not per-hour data.
// -------------------------------------------------------------------
export function calcVehicleInteriorTempPass(TaArr, globArr, elevArr, vaArr, vehicleType = 'car', ventilated = false, speedMph = 0) {
const preset = VEHICLE_TYPES[vehicleType] || VEHICLE_TYPES.car;
const n = TaArr.length;
const result = new Array(n);
const alpha = 1 - Math.exp(-1 / (preset.lagHours ?? 0.5));
// Seed at the first hour's air temp - a vehicle left overnight has
// equalised with the outside air.
let Ti = TaArr[0] ?? 15;
for (let i = 0; i < n; i++) {
const Ta = TaArr[i];
if (Ta == null) { result[i] = null; continue; }
const target = calcVehicleInteriorTemp(
Ta, globArr[i] ?? 0, elevArr[i], vehicleType, ventilated, speedMph,
vaArr ? (vaArr[i] ?? 0) : 0
);
if (target == null) { result[i] = null; continue; }
Ti = Ti + alpha * (target - Ti);
// Can't be cooler than outside air; physical cap at 90 -C.
result[i] = Math.max(Ta, Math.min(Ti, 90));
}
return result;
}
// -------------------------------------------------------------------
// SHADE AIR TEMPERATURE (per-environment microclimate)
// -------------------------------------------------------------------
// Estimates the AIR temperature you'd actually experience sitting in
// the typical shade of the selected Solar Model - NOT the felt temp.
//
// Baseline: this builds on the SAME environment air temp that SunSoak
// uses - TaEnv = Ta + env.taOffset - so every Solar Model's air-temp
// shift (urban heat island, forest/river evapotranspiration cooling,
// desert baking) flows through here automatically and the two columns
// can never drift out of sync. On top of that env baseline we add only
// the LOCAL sun-trap effects: surrounding surfaces (walls, sand, rock)
// re-radiate heat into the still air pocket, while wind mixing pulls
// the pocket back toward the ambient reading. Crucially the sun-trap
// warming is driven by the environment-REDUCED radiation, so a forest
// canopy shades the sun-trap just as it shades the person.
//
// IMPORTANT: this is an air-temperature nudge, deliberately gentle.
// The radiant load of direct sun is already handled by Tmrt/SunSoak,
// so we do NOT re-add it here - that would double-count the sun.
//
// The environment air-temp shift is NOT re-specified here - it lives in
// env.taOffset and reaches us via the TaEnv passed in as `Ta`. The shade
// block only carries the two local micro-effects:
//
// shade.shelter - 0-1 wind shelter. High = enclosed (canyon, canopy),
// so wind mixing has little effect. Low = breezy
// (open water, single beach umbrella).
// shade.sun - 0-~1.2 surface solar gain. How sun-baked the
// surroundings are: hot concrete/sand high, shaded
// forest floor / cool water low. Drives how much the
// still air warms on a sunny hour.
//
// Behaviour: a calm sunny hour in a sheltered, sun-baked environment
// reads a degree or two above the environment air temp; a windy or
// overcast hour collapses back toward it. Cooling models (forest,
// river) sit below official Air because their taOffset already has.
//
// TaEnv - environment air temp (Ta + env.taOffset), the shared
// SunSoak baseline. Passed in as the `Ta` argument.
// effectiveRad - environment-REDUCED radiation reaching the sun-trap.
// -------------------------------------------------------------------
export const SHADE_K_SUN = 1.5; // surface-warming strength (per full-sun, fully sun-baked)
export const SHADE_K_WIND = 0.15; // wind-mixing strength (per m/s, fully exposed)
export function calcShadeAirTemp(Ta, effectiveRad, va, elev, shade) {
if (Ta == null || !shade) return Ta ?? null;
const rad = Math.max(0, Math.min(1, (effectiveRad || 0) / 800)); // 0-1, saturates ~800 W/m-
const daytime = (elev != null && elev > 0) ? 1 : 0;
const sunTerm = SHADE_K_SUN * (shade.sun ?? 0) * rad * daytime;
const windTerm = SHADE_K_WIND * (1 - (shade.shelter ?? 0)) * Math.max(0, va || 0);
// Wind can only cancel the sun-trap warming, never push the pocket BELOW
// ambient: mixing drags trapped air back toward the open-air reading, it
// does not manufacture cold. Without the floor a breezy night read a
// degree under Air for no physical reason (and Pet Shade with it).
// Genuinely cool environments still sit below Air via env.taOffset.
const adj = Math.max(0, sunTerm - windTerm);
// Local sun-trap nudge only - the env air-temp shift is already in Ta.
return Ta + Math.min(5, adj);
}
// -------------------------------------------------------------------
// SHADE FELT TEMPERATURE - "SunSoak, but out of the sun"
// -------------------------------------------------------------------
// calcShadeAirTemp() above answers "what would a thermometer in the shade
// read". That is NOT what someone standing in the shade feels: they are
// still outdoors, so wind still strips heat off them and humidity still
// blunts their sweating, exactly as in the open. Reading a shade AIR temp
// beside SunSoak's FELT temp made Shade look hotter than full sun on a
// breezy night - both numbers were right, but they were different
// quantities sitting in adjacent columns.
//
// So Shade runs the same UTCI machinery as SunSoak, changing one input:
// the radiant load. Out of the beam a person sees no direct sun, only the
// part of the sky their shade does not block, so:
//
// direct -> 0 (that is what shade IS)
// diffuse -> difEnv * skyView
// global -> the same diffuse (no beam reaching the ground to reflect)
//
// skyView derives from shade.shelter, reusing it as an enclosure proxy: a
// forest canopy or building canyon that stops the wind also hides most of
// the sky dome, while an umbrella on open sand blocks little besides the
// sun itself. It is a proxy, not a measurement - like the rest of the
// shade block, this is an estimate of a typical spot in that environment.
//
// Consequences worth knowing: at night there is no beam to remove, so
// Shade converges on SunSoak (they differ only by the sun-trap nudge in
// the air temp). In daytime sun Shade sits several degrees under SunSoak,
// which is the whole point of standing in it. On a windy hour BOTH fall
// together, because both now carry the same wind term.
// -------------------------------------------------------------------
export const SHADE_SKY_BLOCK = 0.6; // fraction of the sky an enclosed shade hides, per unit shelter
export function calcShadeFeltTemp(TaShade, difRad, va, eh, elev, shade) {
if (TaShade == null) return null;
const skyView = 1 - SHADE_SKY_BLOCK * (shade?.shelter ?? 0);
const difShade = Math.max(0, (difRad || 0) * skyView);
const Tmrt = calcTmrt(TaShade, 0, difShade, difShade, elev);
return utciApprox(TaShade, Tmrt, va, eh);
}
// -- Pets Profile - Fur Surface Temperature -----------------------------------
// Cats and small dogs experience solar heat differently from bare pavement:
// a coat absorbs radiation like calcConcreteTemp()'s slab, but a thin/short
// coat over a living, blood-perfused body insulates less than a heavy double
// coat and never bakes as hot as dead concrete. furAlbedo varies by coat
// colour (FUR_COLORS in utils.js): black ~0.05, brown ~0.15, golden ~0.25,
// white ~0.40. Paw burn risk reuses calcConcreteTemp() output directly (the
// ground/pavement surface temp is already the right number for paw contact):
// < 40 -C - safe
// 40-52 -C - discomfort / possible burn (hold-your-hand-for-7-seconds test)
// > 52 -C - burns within 60 seconds
// A breed-specific multiplier (brachycephalic/senior) is deliberately not
// modelled - this profile targets cats and small dogs as the reference animal.
// -------------------------------------------------------------------------
export function calcFurSurfaceTemp(Ta, globalRad, windSpeed, solElev, furAlbedo = 0.15, uv, cloudCat) {
if (globalRad == null || Ta == null) return null;
const elevDeg = (solElev != null) ? Math.max(5, solElev) : 45;
const angleCorrection = Math.sin(elevDeg * Math.PI / 180);
// Same cloud/UV-clarity derating calcConcreteTemp() applies on top of the
// raw radiation figure - both models share the same input array
// (hourlyRows.effectiveRad), so without these fur was absorbing the full
// beam while concrete was further discounted for haze/cloud, making fur
// read hotter than pavement on anything but a clear sky (backwards from
// the "never bakes as hot as dead concrete" intent above).
const cloudTransmit = cloudCat === 'clear' ? 1.00
: cloudCat === 'wispy' ? 0.92
: cloudCat === 'scattered' ? 0.72
: cloudCat === 'overcast' ? 0.28
: 1.00;
const uvClarity = (uv && uv > 0)
? 0.85 + 0.15 * Math.min(1, uv / 8)
: 1.0;
const absorbed = globalRad * (1 - furAlbedo) * angleCorrection * cloudTransmit * uvClarity;
// Lower convective coefficient than bare concrete (hc=10+4.5*sqrt(v)) -
// trapped air in the coat insulates, running the surface hotter for a
// given wind speed than pavement, but not as insulated as a thick double
// coat, so it sits between "bare skin" and "heavy fur".
const hc = 7 + 3.5 * Math.sqrt(Math.max(windSpeed || 0, 0));
const Ts = Ta + absorbed / hc;
// Floor at air temp; cap below concrete's 85-C ceiling - a thin coat over
// a living body doesn't bake as hot as dead pavement. No wet-fur/rain
// evaporative term yet (future refinement, mirrors calcConcreteTemp's).
return Math.max(Ta, Math.min(Ts, 70));
}
export const LAG_HOURS_FUR = 0.4; // ~24 min - low thermal mass, animal moves in/out of sun
export function calcFurSurfaceTempPass(TaArr, radArr, vaArr, elevArr, furAlbedo = 0.15, uvArr, cloudCatArr) {
const n = TaArr.length;
const result = new Array(n);
const alpha = 1 - Math.exp(-1 / LAG_HOURS_FUR);
let Tf = calcFurSurfaceTemp(TaArr[0], radArr[0], vaArr[0], elevArr[0], furAlbedo, uvArr?.[0], cloudCatArr?.[0]) ?? TaArr[0];
for (let i = 0; i < n; i++) {
const target = calcFurSurfaceTemp(TaArr[i], radArr[i], vaArr[i], elevArr[i], furAlbedo, uvArr?.[i], cloudCatArr?.[i]) ?? TaArr[i];
Tf = Tf + alpha * (target - Tf);
result[i] = Math.max(TaArr[i], Tf);
}
return result;
}
// -- Vehicle-at-speed thermal model (IMPLEMENTED) ------------------------------
// calcVehicleInteriorTemp now takes a speedMph argument (Static / 20 / 50 / 70
// mph in the UI). Forced convection over the shell scales with road speed, and
// windows-open through-flow scales further with speed, so a moving cabin runs
// cooler than the same parked car. speedMph = 0 is the parked case - which is
// no longer "no airflow", since ambient wind now drives shell convection too.
//
// -- CYCLIST AT SPEED (FUTURE) ------------------------------------------------
// A cyclist generates their own headwind, so the felt temperature (UTCI) is
// very different from a stationary person - even without ambient wind.
// Could share the same speed-class approach as the vehicle model:
// Slow (10 mph) / Moderate (15 mph) / Fast (25 mph)
//
// Key differences from the vehicle model:
// - The cyclist IS the exposed person - use the UTCI polynomial directly
// with va = max(ambientWind, cyclingSpeed * conversionFactor)
// - No cabin heating effect - the cyclist gets wind chill, not solar entrapment
// - High metabolic heat generation raises core temp (links to the activity
// modifier planned in compute.js - cycling at speed is a combined effect)
// - UVA/UVB exposure is unchanged - still fully exposed to the sun
//
// Could be a sub-option within the existing Cycling activity variant rather
// than a separate column - e.g. a speed picker in the variant controls.
// -----------------------------------------------------------------------------
// -------------------------------------------------------------------
// VAPOUR PRESSURE - Magnus formula.
// -------------------------------------------------------------------
// Converts air temperature (-C) and relative humidity (%) to vapour
// pressure in hPa. Used as the humidity input to utciApprox().
// -------------------------------------------------------------------
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 algorithm (degrees above horizon).
// -------------------------------------------------------------------
// Accurate to within ~0.01- for most practical purposes. Returns a
// negative value when the sun is below the horizon (civil twilight
// starts at -6-, nautical at -12-, astronomical at -18-).
// -------------------------------------------------------------------
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 (Tmrt)
// -------------------------------------------------------------------
// Tmrt is the uniform temperature of an imaginary enclosure that would
// cause the same net radiation exchange as the actual environment.
// It accounts for:
// - Direct solar beam (DNI), scaled by the projected-area factor fp
// - Diffuse sky radiation (scattered and cloud-reflected)
// - Ground-reflected shortwave (albedo - global radiation)
// - Longwave thermal emission from surrounding surfaces (- blackbody at Ta)
//
// The fabric index 0.308 (fp formula from ISO 7933) projects the sun
// onto a standing person's silhouette as a function of solar elevation.
// Output feeds directly into utciApprox() as the Tmrt argument.
// -------------------------------------------------------------------
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));
// Clamp Tmrt-Ta to the polynomial's validated domain. The UTCI fit is only
// valid for a mean-radiant offset of -30..+70 -C; on hot urban/concrete
// profiles the modelled Tmrt can exceed Ta by more than +70-, where the
// polynomial extrapolates and loses accuracy. Clamping keeps the felt
// temperature inside the range the coefficients were derived for.
const D_Tmrt = Math.max(-30, Math.min(70, 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;
}