Files
sunscope/assets/js/physics.js
T

581 lines
33 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.
//
// 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)
// calcConcreteTemp(Ta, globalRad, windSpeed) urban surface temp
// calcIndoorTempPass(TaArr, globArr, elevArr, buildingType) passive indoor temp
// calcManagedIndoorTempPass(TaArr, globArr, elevArr, buildingType) managed indoor temp
// calcVehicleInteriorTemp(Ta, globalRad, solElev, vehicleType, ventilated)
// 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 — 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: ~10 W/m²K still air, rises with wind
// (10 reflects realistic natural convection; 5 was too low and ran too hot)
const hc = 10 + 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. Each hour builds a realistic
// target temperature from outdoor air, window solar gain, retained
// warmth, and internal gains, then the room temperature lags toward
// that target. This avoids runaway accumulation while still giving
// the characteristic late-day indoor 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
// ═══════════════════════════════════════════════════════════════════
// Two-pass function: call with the full arrays of hourly Ta and globalRad.
// Returns an array of indoor temperatures, one per hour.
// buildingType must be a key of BUILDING_TYPES; defaults to 'brick'.
export function calcIndoorTempPass(TaArr, globArr, elevArr, buildingType = 'brick') {
const preset = BUILDING_TYPES[buildingType] || BUILDING_TYPES.brick;
const { lagHours, glazingRatio, gValue, orientFactor, solarScale, baseTemp, internalGain, retainedScale } = preset;
const n = TaArr.length;
const result = new Array(n);
const alpha = 1 - Math.exp(-1 / lagHours); // per-hour blending weight
// Seed to a plausible occupied indoor baseline rather than outdoor air.
let Ti = Math.max(TaArr[0] ?? 15, baseTemp ?? 16);
for (let i = 0; i < n; i++) {
const Ta = TaArr[i] ?? Ti;
const glob = globArr[i] ?? 0;
// Solar gain through windows (W/m² effective)
const solarGain = glob * glazingRatio * gValue * orientFactor;
const retainedWarmth = Math.max(0, (baseTemp ?? 16) - Ta) * (retainedScale ?? 0.35);
const target = Ta + solarGain * (solarScale ?? 0.1) + retainedWarmth + (internalGain ?? 0.7);
// Apply thermal lag: blend toward the hour's target instead of adding
// solar gain repeatedly onto the previous indoor temperature.
Ti = Ti + alpha * (target - Ti);
// Can't be colder than outdoor (house doesn't actively cool)
result[i] = Math.max(Math.min(Ti, 55), Math.min(Ta, Ti));
}
return result;
}
// ═══════════════════════════════════════════════════════════════════
// UK HOUSE INDOOR TEMPERATURE — MANAGED (curtains closed, windows open)
// ───────────────────────────────────────────────────────────────────
// Models the same UK brick house as calcIndoorTempPass but with two
// behavioural interventions that reflect standard heatwave advice:
//
// 1. CURTAINS CLOSED
// Thick curtains block ~80% of window solar gain before it enters
// the room. The small remaining fraction is diffuse light through
// the curtain fabric. Wall conduction is unchanged.
//
// 2. WINDOWS OPEN (smart ventilation)
// When outdoor air is cooler than the indoor air, windows are open
// and a ventilation heat exchange pulls the indoor temp toward
// outdoor. When outdoor is hotter than indoor, windows are kept
// shut — so this strategy never makes things worse, only better.
// Ventilation rate ~2 air changes/hour for a well-opened house.
//
// Same thermal lag model as calcIndoorTempPass (4 h brick time constant).
// ═══════════════════════════════════════════════════════════════════
// buildingType must be a key of BUILDING_TYPES; defaults to 'brick'.
export function calcManagedIndoorTempPass(TaArr, globArr, elevArr, buildingType = 'brick') {
const preset = BUILDING_TYPES[buildingType] || BUILDING_TYPES.brick;
const { lagHours, glazingRatio, gValue, orientFactor, curtainBlock, ventAlpha, solarScale, baseTemp, internalGain, retainedScale } = preset;
const n = TaArr.length;
const result = new Array(n);
const alpha = 1 - Math.exp(-1 / lagHours);
let Ti = Math.max(TaArr[0] ?? 15, baseTemp ?? 16);
for (let i = 0; i < n; i++) {
const Ta = TaArr[i] ?? Ti;
const glob = globArr[i] ?? 0;
// Solar gain — curtains block curtainBlock fraction
const solarGain = glob * glazingRatio * gValue * orientFactor * (1 - curtainBlock);
const retainedWarmth = Math.max(0, (baseTemp ?? 16) - Ta) * (retainedScale ?? 0.35);
const target = Ta + solarGain * (solarScale ?? 0.1) + retainedWarmth + (internalGain ?? 0.7);
// Apply thermal lag toward the managed target.
Ti = Ti + alpha * (target - Ti);
// Smart ventilation: only open windows when outside is cooler
if (Ta < Ti) {
Ti = Ti + ventAlpha * (Ta - Ti);
}
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 unless ventilation is enabled. No evaporative cooling.
// Cars warm quickly; motorhomes and caravans are treated as insulated living
// spaces with 2535 mm sandwich panels, so panel heat gain is much smaller
// and the interior response is slower than a car cabin. Because they are
// occupied living spaces, they also retain warmth from previous hours, people,
// appliances, and background heating; without that, cool sunny days are
// under-estimated badly.
//
// 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, vehicleType = 'car', ventilated = false) {
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;
// ── 1. Panel conduction ──────────────────────────────────────────
const panelAbsorbed = globalRad * (1 - preset.albedo); // 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. 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. Side-window glazing gain (angle-dependent) ─────────────────
// Glazing transmission for auto glass ~0.70; scaled by vehicle glazing area.
const tau = 0.70 * preset.glazingArea;
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 rejection: an effective blend of leakage, internal air volume,
// and surfaces exchanging heat with the outside. With windows open it is
// roughly 5× higher — air moves freely
// through the cabin, flushing heat out and capping interior temperature much
// closer to ambient. Cabin temp still rises a little due to panel/roof solar gain.
const effectiveHLoss = ventilated ? preset.hCabinLoss * 5 : preset.hCabinLoss;
const thermalMass = preset.thermalMass ?? 1;
const solarRise = (totalGain / effectiveHLoss) * thermalMass;
// 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.
const retainedWarmth = preset.retainedWarmth
? Math.max(0, 8 - 0.25 * 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));
}
// ── FUTURE FEATURE v2: Pets Profile — Fur Temperature & Heat Stress ──────────
// Dogs and cats experience heat very differently from humans:
//
// FUR SURFACE TEMPERATURE
// Dark/thick fur absorbs solar radiation and can run significantly hotter
// than air temperature — similar in principle to calcConcreteTemp() but
// with fur-specific albedo values:
// Black fur: albedo ~0.05 (almost all radiation absorbed)
// Brown fur: albedo ~0.15
// Golden fur: albedo ~0.25
// White fur: albedo ~0.40
// Add a calcFurSurfaceTemp(Ta, globalRad, windSpeed, furAlbedo) function
// mirroring calcConcreteTemp but with appropriate convective coefficients
// for fur (lower hc than concrete — fur insulates).
//
// PAW BURN RISK
// Paw pads are highly sensitive to surface temperature. Use calcConcreteTemp()
// output directly — the existing concrete surface temp is already the right
// number. Threshold guidance:
// < 40 °C — safe
// 4052 °C — discomfort / possible burn (hold-your-hand-for-7-seconds test)
// > 52 °C — burns within 60 seconds
// Display as a simple traffic-light column in the Pets profile.
//
// BREED-SPECIFIC HEAT STRESS
// Brachycephalic breeds (Bulldogs, Pugs, French Bulldogs, Persians) have
// severely impaired thermoregulation — their safe UTCI ceiling is much lower.
// A breed selector (Normal / Brachycephalic / Senior) applies a risk multiplier
// to the UTCI thresholds, similar to the planned activity modifier.
//
// SUGGESTED COLUMNS FOR PETS PROFILE
// Hour | Air Temp | UTCI+P | Fur Surface Temp | Paw Burn Risk | Shade Advised
//
// ─────────────────────────────────────────────────────────────────────────────
// ── FUTURE FEATURE: Vehicle-at-speed thermal model ────────────────────────────
// Idea: extend calcVehicleInteriorTemp (or add a companion function) to model
// cabin temperature for a vehicle travelling at speed, not just parked.
//
// Key physics differences from the static model:
// • Forced convection over the shell scales with vehicle speed (v²), so
// hOut rises significantly — shell cools much faster than when parked.
// • Above ~30 mph the vehicle's own forward motion dominates airflow, so
// ambient wind direction becomes largely irrelevant (simplifies the model).
// • Speed classes to model: urban (~20 mph), dual carriageway (~50 mph),
// motorway (~70 mph) — each with a derived hOut multiplier.
// • Windows-open behaviour changes completely at speed: at 70 mph open
// windows create high-velocity through-flow, dramatically cutting cabin
// temp vs. the sealed-car case (but much less pleasant than AC!).
// • Roof and bonnet solar gain stays the same; side-glass gain is unchanged.
// • AC-off vs AC-on would be the primary user toggle alongside speed class.
//
// Suggested signature:
// calcVehicleInteriorTempAtSpeed(Ta, globalRad, solElev, vehicleType, speedMph, windowsOpen)
//
// This would be a useful companion output column (e.g. "Vehicle (moving)")
// for road-trip planning, dog-in-car safety at a rest stop vs. motorway, etc.
//
// ── CYCLIST AT SPEED (companion to the above) ────────────────────────────────
// 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));
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;
}