sunscope/assets/js/physics.js

// ------------------------------------------------------------------------
//  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
//    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 (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.28-0.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 ~15-18% 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
//   20-26 -C - comfortable
//   26-32 -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. VERTICAL GLAZING GAIN (sun-angle dependent)
//     When solar elevation is between ~10- and ~60-, the sun's rays
//     cut through the side glass and windscreen at an angle that allows
//     significant transmission into the cabin. 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.
//     Outside that window a diffuse floor still applies: scattered sky
//     light enters the glass at any sun angle, so the tent function is
//     clamped rather than switched off (an unclamped tent produced a
//     discontinuity where a LOWER sun gave MORE glazing gain).
//
//  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.
  // Tent function peaks at 35- elevation, where the sun cuts squarely through
  // side glass and windscreen, and tapers either side. It is clamped at 0.30
  // (= the 0.15 diffuse floor once the 0.5 not-in-beam factor is applied) so
  // scattered sky light always gets in. Without that clamp the gain fell to
  // zero at 10- and 60- and then jumped back up outside the range.
  const tau = 0.70 * preset.glazingArea;
  const tent = (solElev != null && solElev > 0)
    ? Math.max(0, 1 - Math.abs(solElev - 35) / 25)
    : 0;
  const glazingGain = tau * globalRad * Math.max(0.30, tent) * 0.5; // occupant not in beam

  // -- 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);
  const adj = sunTerm - windTerm;
  // Local sun-trap nudge only - the env air-temp shift is already in Ta.
  return Ta + Math.max(-4, Math.min(5, adj));
}

// -- 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) {
  if (globalRad == null || Ta == null) return null;
  const elevDeg = (solElev != null) ? Math.max(5, solElev) : 45;
  const angleCorrection = Math.sin(elevDeg * Math.PI / 180);
  const absorbed = globalRad * (1 - furAlbedo) * angleCorrection;
  // 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) {
  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) ?? TaArr[0];
  for (let i = 0; i < n; i++) {
    const target = calcFurSurfaceTemp(TaArr[i], radArr[i], vaArr[i], elevArr[i], furAlbedo) ?? 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;
}