diff --git a/scripts/bake-sky.mjs b/scripts/bake-sky.mjs index 4e9d015a..827110c2 100644 --- a/scripts/bake-sky.mjs +++ b/scripts/bake-sky.mjs @@ -1,10 +1,16 @@ -// Bakes Midnight Sky and Noonday Sky from src/sky/sun.js into src/sky/baked.ts. -// Run `npm run bake:sky` after changing sun.js, skyTokens.js or a moment. +// Bakes Midnight Sky and Noonday Sky from src/sky/sun.js into src/sky/baked.ts, with the +// atmospheric-scattering tables (src/sky/scatter.js) as `scatter`. +// Run `npm run bake:sky` after changing sun.js, skyTokens.js, scatter.js or a moment. import { writeFileSync } from 'node:fs'; import { bakeSky } from '../src/sky/skyTokens.js'; +import { bakeScatter } from '../src/sky/scatter.js'; const out = new URL('../src/sky/baked.ts', import.meta.url); -const baked = { midnight: bakeSky('midnight'), noonday: bakeSky('noonday') }; +const baked = { + midnight: bakeSky('midnight'), + noonday: bakeSky('noonday'), + scatter: bakeScatter(), +}; writeFileSync( out, `// Generated by scripts/bake-sky.mjs from src/sky/sun.js. Do not edit.\n` + diff --git a/src/sky/README.md b/src/sky/README.md new file mode 100644 index 00000000..7afc80a7 --- /dev/null +++ b/src/sky/README.md @@ -0,0 +1,83 @@ +# Baked sky + +`npm run bake:sky` (`scripts/bake-sky.mjs`) lights two fixed moments with `sun.js`, the single source of +the sun and the air, and writes `baked.ts`. `baked.midnight`, `baked.noonday` and `baked.scatter` are its +three keys. Everything here is linear RGB unless a name says `hex`. + +## Per sky: `baked.midnight`, `baked.noonday` + +| field | meaning | +| ------------------------ | -------------------------------------------------------------------------------------------------------------------------------------------- | +| `moment` | `{ time, lat, lon }` fed to `sun.js` (equinox, 45N, prime meridian; 00:00 and 12:00 UTC) | +| `elevation` | sun elevation in degrees at that moment: -44.6 midnight, +45.2 noonday | +| `sun.elevation_deg` | the same, for scene code; fade stars by this (and see below for twilight) | +| `sun.azimuth_deg` | degrees clockwise from north, 0..360 | +| `sun.color_linear` | direct sun colour at normal incidence, clear zenith ~ 1; multiply by `max(0, dot(n, toSun))` | +| `ambient_linear` | sky ambient on an upward surface, with a night floor of `AMBIENT_FLOOR` (a moonlit skyglow, not black); daylight is `sun.js`'s own `diffuse` | +| `horizon_hex` | colour of the last gradient stop, the horizon | +| `gradient` | 7 sRGB hex stops, zenith to horizon | +| `tokens`, `bubbles`, ... | theme colours (see `skyTokens.js`) | + +The two skies are single moments. To fade stars through twilight, take the elevation from `sun.js` at the +viewer's real time (`sunPosition(time, lat, lon).elevation`) rather than from `baked..sun`, which only +ever says -44.6 or +45.2. The scattering table below goes dark at -12 degrees and below (the planet's shadow +covers every sunlit air path), so a star fade between about -6 and -12 degrees follows it. + +## `baked.scatter`: atmospheric scattering for a WebGL2 shader + +A Hillaire-style (2020) single-scattering atmosphere, reduced to two small tables. The air is `sun.js`'s: +`TAU` (its zenith optical depth per channel) is the vertical optical depth here, split into a grey aerosol +part (0.02) and a Rayleigh part, so the transmittance table equals `exp(-TAU)` at the zenith. Planet radius +6360 km, atmosphere 100 km, scale heights 8 km (Rayleigh) and 1.2 km (Mie, g = 0.8, albedo 0.9), observer +100 m above the ground. Single scattering only: no ground bounce, no multiple scattering, no sun disc +(the consumer draws it from `sun.color_linear`). That makes dusk a little darker than the real thing. + +Constants: `planet_radius_km`, `atmosphere_height_km`, `observer_altitude_km`, `extinction_per_km` (RGB at the +observer, for aerial perspective: `exp(-k * distance_km)`). + +### Encoding (both tables) + +IEEE 754 half floats (16 bit), little-endian, base64 of the raw bytes, no header. The stored value is +`value * store_scale` (`store_scale` = 64; it lifts the dark end of twilight clear of the half-float +subnormals). Exact decode: for a 16-bit word `h`, `e = (h >> 10) & 31`, `m = h & 1023`, +`half = e == 0 ? m * 2^-24 : (1 + m / 1024) * 2^(e - 15)` (no infinities or NaNs are stored), and +`value = half / store_scale`. In WebGL2 upload as `RGB16F` / `HALF_FLOAT` from a `Uint16Array` over the +decoded bytes, and the GPU does the decode: sample, then divide by `store_scale`. Precision is 11 bits +(about 0.05%): the horizon glow does not step. + +### `transmittance`: `width` 32, `height` 1, 3 channels + +Transmittance from the observer to space along a view elevation, per channel. Texel `i` is the elevation +`-4 + 94 * (i / 31)^2` degrees (`viewElevation(i)` in `scatter.js`; the inverse is `viewIndex(deg)`): +finest at the horizon. Elevations that look at the ground (-4 to about -0.3) are 0. Sample it for the sun's +direct colour: at the sun's elevation, `sun.color_linear`-like attenuation, redder as it sets. + +### `sky_view`: `width` 7, `height` 32, `depth` 12, 3 channels + +Sky radiance per steradian for a sun of unit irradiance (multiply by your own sun colour and exposure). +One 3D texture, `x` fastest: + +- `x`: azimuth of the view away from the sun, 0, 30, ... 180 degrees (`azimuth_step_deg`). The sky is + symmetric about the sun's vertical plane, so use `abs` of the relative azimuth. +- `y`: view elevation, same axis as the transmittance table (`viewElevation(i)`). +- `z`: sun elevation, `sun_elevations_deg` = -18, -12, -8, -5, -2, 0, 3, 8, 15, 25, 45, 90. Between them + interpolate linearly; `sunIndex(deg)` in `scatter.js` gives the fractional index. + +Element `((z * 32 + y) * 7 + x) * 3 + c` of the decoded array. Texel centres, for linear filtering: + +```glsl +uniform highp sampler3D uSky; // RGB16F, 7 x 32 x 12, LINEAR, CLAMP_TO_EDGE +uniform float uSunIndex; // sunIndex(sunElevationDeg), 0..11, from the CPU +vec3 skyRadiance(float viewElDeg, float azFromSunDeg) { + float y = sqrt(clamp((viewElDeg + 4.0) / 94.0, 0.0, 1.0)) * 31.0; + float x = clamp(abs(azFromSunDeg) / 180.0, 0.0, 1.0) * 6.0; + return texture(uSky, vec3((x + 0.5) / 7.0, (y + 0.5) / 32.0, (uSunIndex + 0.5) / 12.0)).rgb / 64.0; +} +``` + +Row width is 21 bytes, so set `gl.pixelStorei(gl.UNPACK_ALIGNMENT, 1)` before `texImage3D`. + +Size: 16,128 bytes for `sky_view` (21.5 KB as base64) and 192 bytes for `transmittance`; about +twice the "few KB" of an 8-bit table, the price of not quantising the glow. `scatter.test.ts` pins the +layout, the encoding, the checksums and a set of physical anchors (blue noon sky, red sunset toward the sun, +black below the planet's shadow limit). diff --git a/src/sky/baked.ts b/src/sky/baked.ts index ecaffbe5..b732fbc1 100644 --- a/src/sky/baked.ts +++ b/src/sky/baked.ts @@ -233,4 +233,32 @@ export const baked = { ], link: '#2a4f27', }, + scatter: { + planet_radius_km: 6360, + atmosphere_height_km: 100, + observer_altitude_km: 0.1, + extinction_per_km: [0.020271980981493083, 0.02767884208799963, 0.044961518003181586], + format: 'half-float', + store_scale: 64, + view_elevation_deg: { + min: -4, + max: 90, + samples: 32, + }, + transmittance: { + width: 32, + height: 1, + channels: 3, + data: 'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAArkjBP0wqR02AR145qk9vS4BBtlCITWBGW1EcT3dJ1FE1UN5LL1K9UBdNdVIrUSdOq1KFURlP1lLOUe5P+VILUlRQFlM+UqVQLVNoUupQQFOLUidRUFOpUltRXlPCUohRaVPYUq5Rc1PqUs9Re1P4UupRgVMFUwFShlMPUxRSilMWUyJSjVMcUyxSj1MfUzNSj1MgUzVS', + }, + sky_view: { + sun_elevations_deg: [-18, -12, -8, -5, -2, 0, 3, 8, 15, 25, 45, 90], + azimuth_step_deg: 30, + width: 7, + height: 32, + depth: 12, + channels: 3, + data: 'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAHwAHAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAyhRWDusADhBFCV8AUAFcAAMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABRKEGYE5g40C8wBGQGDAA0AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAqQ++DQsGtQrbCHkCiABRAA4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAVgq5CZkEmAXWBNoBMgAbAAQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADwX3BIkCOAIMAvUACwAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA/AH2ARkB1wC/AFcAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAwACvAFcASgA1ABAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAQQAvAA4AFAAIAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAEQAHAAAAAwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAApCUCH6sLDiMuHHoIqRhIEVABTwFNAAIAAAAAAAAAAAAAAAAAAAAAAAAAKiaHIl4WviPsH4wT1BmGFdsI7wIZARcAAAAAAAAAAAAAAAAAAAAAAAAASSQ2IkMahCGpH4sXmhjnFVgN0wP+AVYAAAAAAAAAAAAAAAAAAAAAAAAAOCGPIBgb3R61HVgYTBbZFMgOPQSpArIAAAAAAAAAAAAAAAAAAAAAAAAAXh4LHisaDxywG68XJBQsE5oOMgQCAw0BAAAAAAAAAAAAAAAAAAAAAAAAoBvAG7oYARkJGQEWcxEiEagNCwQcA1IBAAAAAAAAAAAAAAAAAAAAAAAAoRj/GNwWZxabFnAUdA9cD5MMswMMA30BCQABAAAAAAAAAAAAAAAAAAAAERaMFuEUExRpFHESVw1aDT0LVQPmApIBIQAMAAEAAAAAAAAAAAAAAAAA3xNjFO0SiRELEqsQzwvwC7QJAAO2ApYBQwAiAAUAAAAAAAAAAAAAAAAAWhEIEvQQxg87ENAOygkDCoQItwKDApABaABAABAAAgAAAAAAAAAAAAAAnw88ECgPdA0RDggNZQioCDcHdwJTAoQBjABhACEAEAAEAAAAAgAAAAAAYQ0SDjwN9AtyDIUL6AZaB9QFPwImAnUBqQCAADUAKAARAAIAEAAEAAAA0AtwDMALCAqjCrAJlQXnBcIEDwL+AWUBwgCcAEsARQAnAAcAKAARAAIA6AmWCswJlAgICVsIkwTRBO4D5QHaAVQB1wC0AGAAYwBCABMARAAmAAgAeAjzCGEI+wa9B7wGyQP8A0gDwQG6AUQB4wDHAHMAgQBeACQAYgBBABQAwAaJB6wGcAUKBkMFKQNXA8QCogGfATUB7QDXAIQAmQB6ADgAgABeACUAMwXOBSEFWgTFBCYErALXAlwCiAGHASgB9wDjAJMAsQCTAE4AmAB6ADoAHQSFBPkDhwPPA04DSwJxAggCcQFzARwB/wDuAKAAwwCqAGMAsQCVAFAASgOPAxwD3gIUA6gCAAIhAsQBXwFhAREBAwH2AKwA0QC+AHcAwwCsAGYAogLWAnUCWAKFAikCxgHgAY0BTwFTAQgBBgH+ALcA4ADQAIoA0gDBAHwAIAJLAvcB8wEZAscBmQGsAWABQwFHAQEBCQEFAcEA7QDhAJ0A4wDVAJEAwwHjAZgBqwHGAXwBdQGDATwBOQE+AfsADQENAcsA+ADxAK8A8QDoAKUAgwGWAVEBdwGHAUMBWAFhAR8BMgE4AfcAEwEWAdYAAgECAcIA/QD7ALsAUwFdARwBTgFXARYBQAFHAQcBLgE0AfQAHgEiAeMAEwEXAdgAEAETAdQALQEyAfMALQEyAfMALQEyAfMALQEyAfMALQEyAfMALQEyAfMALQEyAfMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAVy9pKJIUVi03JksSdCcFIIELkRq8ERgBCwABAAAAAAAAAAAAAAAAAAAAxDDfLH4g/C4LK1kedykLJTAY9R5OGS4LwwFcAAMAAAAAAAAAAAAAAAAAHzCmLV4lKS4pLL4jMCl9JpUdhiDFHMASHwxjBX0AAAAAAAAAAAAAAAAAQS4eLX8nyyyqK4AlTyhvJkog0CAxHuoWzhPcDqYEAAAAAAAAAAAAAAAAnSwyLAgoLithKgMmxya1JQYhpSDLHh0ZEBhtFJAMugRDARIADAABAAAA2CqjKp8nYykdKcQlUyXYJCAhQyDUHksa9xmgF0MRkhDfC+ACCAkeA0QAKik1KcAmGigRKC8lPSQOJN8glx+KHvgaRBsxGWsUuhVaEnILVBL6DUoF8CcbKMwlWyZ/JoMk2yLJIncgxR4dHj8b5RstGh4WkhjzFZkQoRbcEzsNNCaKJukkAiU8Jb4joiG2IQUgDx6lHUAbFBzDGoQX2xlCGO0T5RisFqUR6yRDJSIkAiRDJJsitCDZIC0fcB0wHRcbFBwMG0YYsRpEGZYVExqLGIQU5iNHJPQigiIEI6Eh9R8nIGQe5hzEHNYa5xsdG54YHRv1GfkWyBp3GRkWayIGI9shWyHUIc0g0R4sH7EdbxxiHIkaixsHG9EYVBtfGgYYJhsVGmgXRiHQIfEgdSDiIBwg5B1AHhQdChwLHDkaJxvaGukYUBuRGmcYURtxGjIYYCDYIC4ggR8gIAwfKB19HYwcZxt9G+oZxhqeGuwYHBuZGqQYQxuYGokYUx8QIBgfXx4DHxEejhzcHBcc0xr4Gp8ZahpcGuAY0BqBGsMYBhuVGrsYLh7bHggedR0AHj0dDxxWHGMbVhqFGloZExoYGssYehpUGsoYtBp0Gs8YQR3QHSYdtxwsHYscTBvNG7Ya7BkhGhwZwxnVGa8YIxoZGr8YWBo/Gs0YgBz2HGkcGBx8HOwbnRoSGyEakxnNGeQYeRmWGZEYzRnXGagY/Bn8GbkYwxtEHJkbMxvZG/QaDRp4GqQZSBmFGbQYNxlcGXMYfBmTGYkYohm0GZoYxxpmG5gaZBrvGicalxn5GToZChlJGYoY/RgpGVcYMBlSGWcYThlrGXUY/xmAGsoZvhkzGoAZOBmSGeQY2RgZGWgYzBj9GD8Y7hgXGUYYAhkoGU8YYxnOGSkZPhmhGf0Y7hhBGZ4YshjzGE0YphjcGCwYtxjmGCsYwxjvGC8Y8RhLGbAY3xg1GZoYuBgEGWgYlxjYGDoYixjFGCAYkBjFGBoYlBjHGBoYphjzGFwYoRjrGFQYkxjaGEEYhhjIGC0Yfhi8GB8YfBi3GBcYfBi2GBUYgBjCGCkYgBjCGCkYgBjCGCkYgBjCGCkYgBjCGCkYgBjCGCkYgBjCGCkYAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAoTzFMZAc7zXcLqga4jGFKqwV7ixVJAsOHCRXGLsBAAAAAAAAAAAAAAAAzTuBNX0oDzgONPUmmzRTML8i0zCDK8sc0Cz0JNET5h9hEUYAAAAAAAAAKTqbNsEtwTc0NZUsrTTzMdcomzEmLjwkLTDjKqceHy31JP0TUyrKH0cIyjiHNm4wrzZMNS4vLDRYMgYsZTFVLxAo2zBqLawkiTClK0sgKDDlKcgcbzfsNU4xkzXkNFswJjMQMhot2jCOL7MpzzCCLhkoQDEDLuklWTGYLcAk3jUvNZMxnzRWNJ8wEDKFMZwtRDBCL8AqdzDdLpEpQzE6L/8olDFCL58oujR6NHYxsjOPM5EwJjHsMLUtdy+/LkQrCzDFLo8q9DCpL6MqWzH+L5EqwjO0MyQxdzKOMlYwZjBcMI8tii4tLmcrQS9yLhorizCeL7kr+DAPMO8rdzKjMrwwezGwMQYwlS+1L0gtwS2fLU0rfS4HLlArHjBPLygsiDDkL14sdTG8MU0wrzDzMFwvli7VLvEsGi0dLRArzi2ULUkrai/fLkQsGDB6L4ospzD5MMEvCDBRMK4uxS0ULpYsjiyoLMIqNS0kLR0rqS5iLjwsXS/5Lo0s/S9TMPQu/S6QLwkuGC1uLT0sGSxCLGwqsSy7LNsq+i3iLR4smy5wLnIs5y6LLzYuFy6lLm8thyzgLNIrbSvRKxYqPixaLIwqXi1lLeMr6y3mLUMs/i2ZLoktVi3bLeIsDixmLDYryCo3K8MptCsCLDgq0izuLHgrTC1iLQksOS3ILewssSwtLWMsTiv5K6gqPCqxKnYpBytnK+MpViyALAQrvSzkLJErkiwTLWAsJSyXLOIroCpFKygqxik+Ki8pcCraKpEpzysaLI4qPixvLAorAyx3LMcrWysWLBYrDSqqKrYpYyncKfAo7yleKkQpCyt6KxkqmCsELIMqESvjK+sqkCpOK2IqkCkmKlIpDymIKbgofynwKf0oXyrTKqopzipEKwAqQyr+Kiwq5CmRKsYpKCm2KfsoyihCKYcoHymSKb0oyik/KkMpHiqVKoYplyk/KogpVSnyKUEp0ShZKbEokSgHKV4ozyhBKYQoSSm+KeYohin9KRcpCymiKQEp4ChwKdIojCgNKXMoYyjYKDwojCj+KFQo3ihRKZUoBil7KbYonignKZUohigLKXooVijSKEMoPyizKCEoWCjJKC4oiCj5KFMonygSKWYoTyjOKEUoRCjBKDkoLyinKB4oJSiYKA0oMSiiKBIoSCi4KCIoUyjDKCsoICiXKBIoHSiTKA8oGCiMKAcoFiiIKAEoGSiKKAEoHyiQKAQoIiiTKAYoECiDKPonECiDKPonECiDKPonECiDKPonECiDKPonECiDKPonECiDKPonsDGjJYIJViNRGGsBMRuhEZ4AyRcpD3AAgRfSD4IAphgnEa8AJxnNEcYA3jHQJa4JhyNuGHUBXxvFEaIA+hdWD3IAsBcBEIUAwxhHEbMARxnxEcsAdzJnJkAKFCTMGJQB+Rs9Eq8ATxjsD3sAJhhSEI8AJBmxEcEAshloEtoAsTOdJ20LvSSRGdQBnBw1E8oA+RiREI4AyRj5EKQA7BmMEt0AjxpdE/oAGTULKekMGSYpG1kC5x2cFAIBWBrREbQAFhpQEs8AhRslFBYBKhyoFDoBJTgZLPcP0SilHbADnyAyF5AB7hx/FBMBshzWFDgBxh1OFp4BYR4RF9IBvz+iM10XXjAUJYUKCCgoHh8FFyQ0Gy8DLCMMGzUDHyRFHOoDWySUHDcEkkq6PwYn1z0KNq4hnTjFMUYdfDUjLiQYvDT3KxgSEzRdKNEI0TLsJSoFFEk1QNYuNz5+OYcsADpBNgkpCDjVMzAlKjjXMl4izTg7Moge/DirMeEbMUeYPz4zeT2DOoIxvjn1N6AuDzhlNeErfThnNVwqpjkBNlcpKjonNqoo+ERfPkQ1kTx8Ojk0AjkdOFAxRzfWNd8uKThENmcubTmZN6AuCDoUOJMuxEIrPTk2hDv6OSE1QTjFN6UySjauNYgwUTdXNpUw2zgDOD4xcTlmOH0xl0AmPJA2MTpVOYE1MTcPN0ozZDVFNSgxVTYENnoxPTjMN50yxDhSOB0zTD61On82IjmwOIM1HTZONm8zojTJNGgxeDWLNewxWDdGN2EzIzgOOAc0bDx+OTM2SjgbOE81PTWbNUYzAzRQNGkxvjQJNQsyXjalNqszLTdvNzg0ZzqROMo1ODcvN/40hzT6NPIyATPEM0UxIzSNNPQxijUDNqEzPTa7Njc0zTi3N1U1HzZMNqE05DNvNIoyLDIBMwwxRTMaNL4x1jRpNWIzcDUONhc0cDeYNt80OTWJNUE07jLsMxsyfTFYMsswbTJpM3cxPjTcNAUzwjRuNcgz8DWxNWw0fTTgNMYzIjIcM6wx7DDFMYgwuTG1MikxejNcNJgyLTTcNE0z2jTyNP4zwTNPNBMzdjFoMkIxczBHMUYwIDEXMtkwnDLTMyQyWzNYNMgyCzRTNDMzvDKjM2wy5jDMMd8wDjDbMAgwnjCMMYow3TEHM68xfTLEM0Ay2DKXM3ky4THJMtExazBGMYQwcS9/MJ8vMDATMUAwNzFRMjwxvDHvMrox3DGvMs4xKjEMMkUxAjDSMDIw4i4vMDgvoy+pMPYvqDCvMdAwFDEyMjoxEDHnMTQxjjBnMcYwUi9uMM8vai7YL90uAi9OMHcvLDAgMWowgjCJMcIwaTA6MaowCzDZMFQwui4YMEsvBS5lL4wueS79LwYvgS+iMA0wAzD1MFMwvy+nMDEwOi9fMN8vOi6eL9gusS0FL0cuBi51L6Euyi42MHMvLS90MN0v3i4qMI8vgy7vLzIvzy0jL3YubC21Lgwupi0EL0suLy6yL+EudS4FMC0vKC6HL94u7S1EL6AueC2+LiMuNS11LtstWC2nLgMusS0YL2Uu3i1RL5cumS3kLk8ueC29LisuMy1vLuEtCi1DLrUtHC1eLsotTi2eLgEuZy2/Lh4uMS1rLuQtIS1ZLtMtAC01LrAt7CwfLpkt8iwpLqEtBi1GLrktEi1ULsUt7iwfLp8t6iwaLpot4CwQLpAt2SwJLogt2SwKLogt3SwQLo0t3ywULpAt0iwBLoIt0iwBLoIt0iwBLoIt0iwBLoIt0iwBLoIt0iwBLoIt0iwBLoItuz3qOHAtkTCoLOsioyhFJjMeFSUMJHQc8yR4JEEdKybuJSsf2iaxJhIg7j0VOZ4tsDDHLBgjwChtJloeNiUlJJAcEiWUJGEdUiYTJlcfBSfbJisglj6jOTguFjEtLa0jIynxJtseoSV6JOwceiXxJM0d1CaPJuoflSdmJ34g8z/MOngv6TECLnAk8CkCKOYffyYpJakdUiayJase4SePJ40gYChEKCohYkGXPAMxnTO3L6slmCsdKQQhJiiTJi0fCChAJzkgBinPKMIhlCltKYkifESbPxo0+zUHMmAo6i3tK7EjcioVKXshQiqYKW4iyytqK2AkUywuLPckCkylRrI6ED3xOLMu8jRqMsopWzESMA0oKzFwMKsoaTLXMUYqGzOTMhYr8FB4SpY7d0N/PhIxMz2tOfIssTqkN3UqOzsgOKAqpTwtOfkrKD2rOU4sSk/ESSA9hUL4P6I1fT0zPIcyojsMOoUwWzzkOtIwzT2APOgxhT4CPWQyJ02oSPs9QkEKQK843DyiPOI1+zriOkE0EjwHPL80gD1jPQ02Nz4NPqk2J0suR0I+MEB8Px86Fjx8PP83+DnGOvE1CTsFPNQ2yDxwPXg4aj0jPvo43EhaRQo+tj6nPuQ6zzoaPJs4BjlCOvc28jl6OxE4DDwVPWg5lzzAPQ86eEbPQ4o9cT3QPRw7szlRO9Q4PTihOWM3BTm9Ols41DqWPNU5vjsyPY86QUSsQfQ8dzwNPfg60Th6Oss4NjcDOWk3Qzj/OWI4yTkUPOQ5jjqfPKM6n0EsQGQ8aDtjPKI6HTi5OZ44MzZ1ODI3TTdROUI47zg1O7g5ljkUPHI6lD9MPss7MTqiOzI6GzcSOV04ZDXzN9w2UDa4OAw4PThcOms5yjgwOxs6Rj3oPPA6NzmoOrk5MTaCOBU4vDQfN3g2gTUzOJo3VTeeOQ85IjhVOrA5oDvhOzI6azjTOUA5cTUHOJc3NDRqNhM21zSANxU3Yjb4OKw4KzeVOT45uzl/Oow5iDcdOco40TQ7Nwk3hjPQNbA1SjS6NpM2lTVnOEo4PzbuOMs4djh3Ofo4dTaAOFw4TDSFNoQ2yzJNNVQ1qDMQNhc26DTRN9Y3djVbOFw4KTerOHc4kTX0N+w3tTPnNQg2MDLeNAA14TJ+NaQ1VTT1NiQ3yzS1N+Y34jUJOAI41DQNNzA39DJfNZc1rTF/NLQ0OTIANTw1rjM1Nn42NzTTNiI37jQJNzE3NjRENoU2TzLqNDA1QDEuNHA0qzGTNN001jKNNec1cTMONm42MTQrNnQ2YzOXNek1wjGDNNQ05DDTMzU0MjE2NIk0HDL8NF41ljJhNcs1NzNwNcs1hTICNV41STErNII0lzBfMwE0zTDOMz80fjF+NOQ02zHMNDk1RzLSNDU1yzGBNOM04jC/Mzo0VzD+MqkzeTBJM/8z+DASNHk0PDFMNLg0hDFNNLI0MTEUNHc0izA/M/gzIzCuMl8zNDDbMpMzijBxMx40uDDBM0o06DC8M0M0sjBxMxs0QzDUMo4z8S9tMiEz+y+DMjwzMTDgMqUzTjASM90zbDAKM88zTDDdMp8zCTB+MjgzsC87MvEyqS9AMvky3C9wMjAz+i+LMk8zDzCDMkIz/i9tMiozui88MvYygi8XMs8ycy8SMssygi8kMt8yjS8uMusyni8nMuIylC8gMtoyfC8QMskyZS8BMrkyWS/7MbMyVi/7MbQyVy/8MbYyWy/5MbIyWy/5MbIyWy/5MbIyWy/5MbIyWy/5MbIyWy/5MbIyWy/5MbIyej2nO5k2xTKgMW0ucisPLBEqJyhRKWooECjlKTspDynQKyErnSlnLAwsrD3tO9M29DLGMZcuoispLDcqQihyKYYoKigKKlspLykALE0rwCmDLCUsVT5rPJQ3kDNGMiQvISx/LLYqmijhKeEofyiEKsYpmSlRLN8rNireLHgsuT9fPZQ4azROMyQwyCwyLbwrUCnIKp0pMCmDK6MqdSr6LIYsKiucLSMtUkFePzg6ujW3NE0xHy6jLuosySpSLB4rnirHLDQsHSxULrotkSwiL4AujkRCQjE9iThoNxo0xzAmMYcvRy2uLnAtIi1fL2guYi7hMF0wFS+AMfMwQkygSUlEtz8QPi06/TcmOJk1YzRbNQQ0QTTiNbQ0SDXHN2M23DVhOD4341DNTHZD5kQYQho64T5YPTU2fjyFO4M09DxTPEE1iz7LPQw3XT+JPvI3u09xTJhEI0SwQrw83z7kPng5zDwYPRs4hD0HPtw4bj8YQJI6NECkQGw7FE7ASxJFmEJpQkE+8j0jP8E7QzxiPew5+TxuPhQ7vz5lQNE8pj/9QHQ91UyeSi9FPUHFQTU/8jyyPqM8MDsYPSU7NjwaPkk8uT0tQNo9fj6/QKM+mUt+SfNENkAQQZI/GjwFPvk8ADqZPLM7BzuAPZ48xjyHP08+aj1HQCg/wUleSHBE4T5mQIA/2zpTPQA9CTkUPMM74TncPKQ89zuiPlQ+hDyIPy4/IEiYRodDsT2iPy4/zTmyPNc8Rjg0O4g79jhEPHo8sTrLPRY+kzuSPuY+nUXGRCpCwjyiPro++DglPJY8VjdfOiY7Ojh7Ozg8qzkOPbQ9aDq5PXU+aUPOQgFBAjzGPTk+TThYO0s8YDaoObU6RjeVOtg71zhpPEM9dTn8PPI94kDpQB1AyToLPbU9hTeKOvs7mTUOOUI6UjbROT07KTi5O888rjhaPGw9lT5PP+Y+wzlpPDY9ojbaOWQ7+DSMONI5ijUrOac6NTfHOl88CjibO+g8nTynPek96ji+O7485TVEOdY6czQeOGs55TSdOBs6Rzb3Oek7ADeqOmw8wDqJPCU9NTjPOk48RjXCOFI6BjSCNw05XDQkOJs5fzVEOSQ7GDbaOfE7Jzl+O4g8PDcCOs87vzRSONo5VTPlNrk40jN7Nyg51TSqOHE6UzUmOR47FDhUOgU8QDZSORU7TDTiN2w5vTJfNm84EjPMNsI4RjQlONA5qzSKOF86pjZvOSo7bjW6OGs60DM5Nwk5PTLtNS04cTI4Nmg4mDNkN0A5HDQDOLM5izW4OGc6vzQ2ONI5JDOlNq840jGNNeg36jG6NRk4yTKfNsE4RTMdNxw5tTQjOLw5LDSIN0k5kDIkNl84eTE8NYU3ejFQNaw3GzL5NVI4dzJWNpc4DTRMNyU5YjPENs44DzK0NRg4LzH5NDE3HzH6NDs3jDFtNek3zDGwNSY4EDN7NqE4lDIbNmE4oTFTNbQ38jDBNOs21zC0NN82GDH8NE43QjEoNY83OTLMNS846jGMNQM4RDEBNUg3wjCUNLE2oDB+NJg2vzCkNNM21zC9NPk2jzE9NZ43XzEVNWc39jC8NO02nDBxNIU2ejBYNGU2gTBlNHs2izBwNIw2CjHLNAQ38TC2NOc2uDCGNKM2gTBYNGQ2ZDBCNEY2XTA/NEU2XjBCNEk2qTB4NJE2oDBwNIY2iTBdNGw2cTBJNFA2XzA7ND02VTA0NDQ2UzAyNDI2azBENEo2azBENEo2azBENEo2azBENEo2azBENEo2azBENEo2azBENEo27DkmOUI3djKfMvIxJixwLSsuyihTK5gstCgQLG8tzilWLVcvbCr/LScwITpTOX43pTLOMhoyQSyTLVIu6SiBK7Qs0SgpLJAt8yl3LYQvlCokLkAw0zrpOSM4PzNnM58ynCwILtQuTykNLBMtMyl9LP4tayrlLQwwGSugLpQwJDwkO/Q4QDRTNLEzWC37LuAvIiqtLNct/yksLeMuZivJLqYwFyygL0MxrT3cPKc6jDWcNfE03C54MAMx1yv2LWcvpSuSLlswtixQMOExNi3YMKcy10AaQHw9bjhwOK03XzHyMrEzGi6cMKgx7i0RMaIyTC+kMnk0CTB2Mw81hUhWR3hEnj9QP8s9gTibObc5FTVlNy446jQLOOA4CjZLOZA6rDbyOWw7QU1ISoRD2kQjQ049ez+5PsE59zzTPE44cT2MPRA5Ij9sP9Q6AkAwQLg7VkwXSr5EKkTAQ6M/Yj8qQI88ND07Pv067D1TPyc85T/zQKE9dUCaQWA+PkulSXNFvUJoQ8FAXj47QP09ljxrPqM8Rj2VP4Y9ET8jQYE//D/SQUBAV0pBSexFcEGuQkhBTD3cP9Y+tDv+PU89bjwQP08+8D3IQEhAtz5rQdpAzUn4SDZGc0DlQW9BZzwLPyU/azpfPYs9WTtLPow+6zxAQG5Ajj3QQAVBcUm+SFdGaj8rQVRBYzs7PhM/Yjm+PHs9HjqDPWo+FDxpP1NAmjwxQOVACklzSD9GRT6JQBVBRzqAPcw+kTgtPEA9JDnOPBU+zzprPhNArjtBP5tAZEjkR89FXT36P8RAZznfPGo+2TdhO/E8XTgxPKc9vTmOPYU/dDpEPj5A70ZoRgNFnjwIP2tAtThWPP891DaOOpw8fTdXOzQ93zjRPNo+eDlrPbY/70TFRA5E+js1PhJAJDjCO5Q9AjbcOUk8fTZ2OsI8LDgxPDI+qjiyPPM+jULGQmJC5Tp8PXU/Wzf6Oi09VzVHOfY7rTW4OVk8MTdQO5M9ATgVPDo+REDPQBRB8znZPM8+lTZQOs08zDTJOGc7ATUYOfA7PDZoOgI96DYcO5A9rj0QPyhAITlJPDQ+7jW9OXU8WTRgOOc6dDSROEM7cTWhOX48/TU2OvY86DtvPQU/aziWO6Q9YDU9OSQ89TMHOHY6/zMgOKs6xjT4OAg8NzVxOWs8yzlcPBE+oDe6OiE95TTOOLQ7VzN3NxI6PTOANyU6NzRpOEA7kDTJOOA7dDg+O1Q9mjb6Oak8ejRsOC470zL4Nrw5nTLgNrI5fzPdN4o6AzQ5OAc7ITcpOro8vjVUOT08HDQWOLQ6ZTKNNnE5GjJcNlE5uDIQN+05HTOAN0w64jVVOTg8BjXEOLU7lDOUN0Y6CTIzNjE5sDHwNQA5FTJnNmk5WzK2Nq059zSrOJA7bDRHOAQ7AzMNN+I5vjHpNfo4XTGaNb84kzHeNf04wDESNio5RTQgOMw61zO2N2Y6hDKXNok5gDGsNc04HjFYNYw4LzF1Nag4SDGSNcM4cTNZNyE6ADP7Nto5FjIwNjo5TzF7Nak48jApNWg46TApNWs48jA1NXc4kjKYNo45TjJdNmA5tzHXNfU4KTFVNYw42DANNVM4wDD6NEU4vjD7NEc44jH6NRM5vjHZNfk4aDGMNbo4DjE6NXg40TACNUs4tDDpNDc4rTDjNDM4WjF9NbA4TDFxNaU4KDFQNYs4/jApNWs42zAKNVI4xjD3NEI4vzDwND04+DAkNWc4+DAkNWc4+DAkNWc4+DAkNWc4+DAkNWc4+DAkNWc4+DAkNWc4QzV3NXI1/TDtMfIyHSz2LSkwFSk8LHgu+yipLJAvBCr6LfwwmCqoLpcxbjWjNZs1ITEWMiEzOCweLkQwNSlXLKAuGinGLL4vKSofLhoxwSrRLrkx/zUzNiI2mjGhMrwzkyyfLpwwoimwLCUvgikmLSswpSqZLn4xSStZLykyLjdiNzs3lTLBM380US2tL1MxgyppLR0wWirsLckwpyuWL00yMSw5MBEz1TjsOL84TzQKNcg12S7sMMkyKyzmLjcxCyyGLwwy3izQMPQzVS1bMXY0BjwKPJo76jb9N4A4YTGtMzg1fy5XMfgzRS7LMZk0hi9lMwo2HzAeNMU2TkP+QvZB+z2cPs0+iTg3OsY7azVJOOQ5MTWeOL86NzbiOWw8zTaMOvU8b0hURiBBBkSuQjY+iT9QP5I7Mz1tPdo5mD0gPsQ6Ij/6P3Q88T93QP88hUeCRtBCMEN3Q3JAdT96QMM9aT3ePog8Aj7hP0Q9yD8tQfU+WUDOQdA/RkYpRgpE+EFBQ3xBcz6GQFk/wTwAP849Tz0EQLo+5j5IQXFAuT/tQf5AaUW1RWxE7EChQg9CZT0zQCBA+jt/Pn8+bzxoP3o/wj3cQOxAcT50QYdB5URXRZ5EHUDtQTtCgjyLP0VApTrOPbE+VDuNPp0/vDxHQP1ATD3LQJlBr0QhRbpEAj9FQSNCnTuzPjdAlTkgPY8+FDq0PV0/zjtmP8tAXDwkQF9BvUQcRddEFz60QOVBhTryPQxAvziEPEE+FznyPOk+fDpdPnZAOjscP/1ACkVKRQJFYT04QJVBqTlNPak/Fzj9O+A9TzhNPGE+bjl5PRNACjoZPotAhEWgRT5F0DyePz9B+jjAPDI/IzcdO3s9YTeCO9Y9lzi3PFw/Fjk9PRdA8EXtRXFFVzzoPudAbThJPLs+TjZhOho9YzaZOlE90zcVPJs+UTiDPFA/10XMRVVF2DtFPpBA8zfJO0o+ozXFOcA8ljXWOdg8uDYXO+k9XzfNO4I+6UTzRLZEETuxPTxAMzccO+E9GDVDOW887zQ0OWw80DUxOkg9VjbEOsc9AkNjQ4xDUTolPdM/kTaHOoA9pjTWOCg8aDStOA08ETVvObo8ezXlOSI9i0AoQdZBljmhPDM/BjYGOig9STR7ONI78jM9OHY7dDTNOD48xTQpOZA80j1IP5xA5DgmPJs+jTWUOdc8+DMvOGY7PzPDN+o65TNGOKU7LzSLOBI8xjtqPZg/PzhpOws+IzUwOYw8eTPgNwk7rzItN3Y6EzOtN/A6aTMJOE07hzlDPHw+WjedOoY9xTTWOEg8ETN2N7k6PDK2Nhc6ajL5Nls6ozI8N5o6LjgFO6w9XzbpOQw9cDSGOAo8vDIfN3Y65TFZNs055zFqNuM5BzKSNgo6pTb5OQo9jDVNOZ48IzQ9OKA7djLYNj86pTEUNpY5hTEANog5kTERNps5fTUxOYU83jTFODk8uzP2NzY7PzKeNhE6ezHlNXA5RDG2NUo5QTG3NUw5qTSUOBM8TTRQOL47PDN/N9Q6EzJwNu05ZDHLNVw5IDGNNSc5EzGCNR45CjQVOGA7rDPVNx07yzIUN3w68jFONtE5YTHENVg5GzGCNR85CTFxNRA5HjNXN7U65zIoN446ZzK2Ni062zE1Nr05cTHQNWU5NDGXNTI5IjGFNSM5XTKqNiM6RzKXNhM6EDJlNuc5zTEmNrE5kzHvNYE5bTHMNWE5YDG/NVc5yTEhNq05yTEhNq05yTEhNq05yTEhNq05yTEhNq05yTEhNq05yTEhNq05Ri/uMBIzpi0OMB0yKSu3LagwPSmjLAcwECnpLHowqSnLLWwxBCpBLuUxey8RMUIzzy0qMEYyWSvdLccwXynALCAwMCkHLZUwzCnuLYwxKSpnLggyFzCFMeIzVC6HMM0y+itbLiox0CkiLXMwmSlrLe4wQCpjLvcxoyrkLnwy0TB3MpY0ay9IMeUzpCxfL/kxuirtLR8xdCo7LqcxMCtVL9QyoivpL2wzUzI1NOk11jDXMhE1/i29MKAzUCyQL34yHCznLyIzkSykMEw02CwBMas0GTWyNpw4uTNmNd03ujBoM+A1vC7cMfY0Xi4UMmo1DS8dM4I2eC+qMxA3cDyQPf0+pDpxPOw9ATgFOmM8ojW1OFY7RjXXOO870DWkOcA8JjYSOic9cULrQXg+HUHvQJM95D4lPzk8XD3cPSM7ez07Pr87YD5nP6c83D4CQA09CULbQqxA7EDbQRhA7z5kQFQ+jz1TP2A9zT3aP9Q90z6wQAE/YD8WQZw/KUHOQslBPkDeQRxBED5zQPc/4jxrP74+Fj3pP0E//T22QFhAeD4cQbdAXEBUQmxCJj93QbBBHj0nQHFAGTzgPnQ/Ozw+P+w/9TxKQLZAWj2lQBtBdz/BQadCED70QOVBUTyBP5hA2jomPqA/+TphPvo/DTx+P7RAXTwMQBZBiz45QaFCOz12QNtBXzu5PoxAxzlwPXU/xzmLPak/pTpwPndAIzvwPtFA5j3GQHhCmzwIQK9BZDoHPmZA8TjQPB8/2DjPPCo/fzmDPSBA4jnqPW5AfD1uQEJCIjxXP3RBoTlwPTNASzhJPLg+HTgwPJw+mTi6PH4/5jgNPQJASD0wQAxCkjvAPjRBCTnxPPk/kTezO08+FjdZOxE+xzcTPMA+HzhUPC8/TD0OQONBEztGPvZAkDiHPI8/xDb7Ouw9LzaBOo89pzYROxE+ADd2O2o+lT0QQM5BvTrnPb1AMDgwPCw/ITZkOpM9eTXROR09wTUuOnU9AjZ6Ork9Pz5CQNxBizqePYpAyDfOO9I+oTXpOUU96TRCObo8CjV1Oe08ODWrOSA9fT+5QCFCeTpqPV5ATTdUO4E+OzWFOQM9dzTQOGg8eTTfOHo8ljQDOZ08xkCSQbRCfjpEPThA6jbtOjk+6zQ1Ocs8HzR2OCQ8BzRnOBw8FjR7ODA8KULCQo9DjDonPRZAmjaWOvg9rDT0OJ08uDMvON47XTMKOJ87ZDMROLA7IUOTQxFEiToGPeQ/WTZMOr09fDTCOHg8VTP2N4472jKFNys7zjJ/Nyg7YELVQnVDVjrRPI4/IDYLOog9WDSaOFo8ETOsN1U7fTIgN9g6YzIJN8Q6YUD1QNtB3jl+PCA/7TXQOVY9PzR8OEM86DJ9NzM7RDLfNqM6HzK8NoQ6YD22PndALDkSPJ4+uTWYOSc9LjRmODI82DJnNyU7KzLANow6/jGVNmQ6yzrKPC0/Zjg+OxU+gjVgOfg8JTRXOCY83jJmNyk7MTLANpE6ADKQNmU6rjhsOxo+YTdoOpM9RDUnOco8ITRNOB08+zJ6Nz47VTLeNq86IzKuNoM6BzclOmM9OjazOSA9/zTtOJs8ITRHOBg8LjOiN2M7mDIYN+c6aTLsNr86pzVQOds8VzUeObo8tjSzOGw8JDREOBU8dzPcN5Y7/TJxNzg71jJNNxg7xDS4OG48qzSmOGE8bTR5OD88JjRDOBM81zMUONg7izPpN6I7cjPTN487KDRDOBM8KDRDOBM8KDRDOBM8KDRDOBM8KDRDOBM8KDRDOBM8KDRDOBM8SSnQLG4wSSnQLG4wSSnQLG4wSSnQLG4wSSnQLG4wSSnQLG4wSSnQLG4waynvLIowaynvLIowaynvLIowaynvLIowaynvLIowaynvLIowaynvLIow3SlULeUw3SlULeUw3SlULeUw3SlULeUw3SlULeUw3SlULeUw3SlULeUwzCooLqIxzCooLqIxzCooLqIxzCooLqIxzCooLqIxzCooLqIxzCooLqIxXCzbLyQzXCzbLyQzXCzbLyQzXCzbLyQzXCzbLyQzXCzbLyQzXCzbLyQz0y4XMnY10y4XMnY10y4XMnY10y4XMnY10y4XMnY10y4XMnY10y4XMnY1ujXnOAo8ujXnOAo8ujXnOAo8ujXnOAo8ujXnOAo8ujXnOAo8ujXnOAo8cD0PPsk7cD0PPsk7cD0PPsk7cD0PPsk7cD0PPsk7cD0PPsk7cD0PPsk7oz2LP8w9oz2LP8w9oz2LP8w9oz2LP8w9oz2LP8w9oz2LP8w9oz2LP8w99TyiPzg/9TyiPzg/9TyiPzg/9TyiPzg/9TyiPzg/9TyiPzg/9TyiPzg/LDwUP/Q/LDwUP/Q/LDwUP/Q/LDwUP/Q/LDwUP/Q/LDwUP/Q/LDwUP/Q/AjtaPhJAAjtaPhJAAjtaPhJAAjtaPhJAAjtaPhJAAjtaPhJAAjtaPhJA9DmnPf0/9DmnPf0/9DmnPf0/9DmnPf0/9DmnPf0/9DmnPf0/9DmnPf0/IzkLPa0/IzkLPa0/IzkLPa0/IzkLPa0/IzkLPa0/IzkLPa0/IzkLPa0/gziKPE4/gziKPE4/gziKPE4/gziKPE4/gziKPE4/gziKPE4/gziKPE4/CTghPPA+CTghPPA+CTghPPA+CTghPPA+CTghPPA+CTghPPA+CTghPPA+VzeZO5s+VzeZO5s+VzeZO5s+VzeZO5s+VzeZO5s+VzeZO5s+VzeZO5s+yjYVO1I+yjYVO1I+yjYVO1I+yjYVO1I+yjYVO1I+yjYVO1I+yjYVO1I+YTavOhc+YTavOhc+YTavOhc+YTavOhc+YTavOhc+YTavOhc+YTavOhc+FjZjOuk9FjZjOuk9FjZjOuk9FjZjOuk9FjZjOuk9FjZjOuk9FjZjOuk95jUsOsg95jUsOsg95jUsOsg95jUsOsg95jUsOsg95jUsOsg95jUsOsg9zzUKOrI9zzUKOrI9zzUKOrI9zzUKOrI9zzUKOrI9zzUKOrI9zzUKOrI90DX5Oac90DX5Oac90DX5Oac90DX5Oac90DX5Oac90DX5Oac90DX5Oac97jX8Oag97jX8Oag97jX8Oag97jX8Oag97jX8Oag97jX8Oag97jX8Oag9MzYVOrQ9MzYVOrQ9MzYVOrQ9MzYVOrQ9MzYVOrQ9MzYVOrQ9MzYVOrQ9sjZMOs89sjZMOs89sjZMOs89sjZMOs89sjZMOs89sjZMOs89sjZMOs89lje2OgE+lje2OgE+lje2OgE+lje2OgE+lje2OgE+lje2OgE+lje2OgE+oTh+O1s+oTh+O1s+oTh+O1s+oTh+O1s+oTh+O1s+oTh+O1s+oTh+O1s+SDqFPAo/SDqFPAo/SDqFPAo/SDqFPAo/SDqFPAo/SDqFPAo/SDqFPAo/2TwfPjZA2TwfPjZA2TwfPjZA2TwfPjZA2TwfPjZA2TwfPjZA2TwfPjZA6z9/QHJB6z9/QHJB6z9/QHJB6z9/QHJB6z9/QHJB6z9/QHJB6z9/QHJBGUGQQVtCGUGQQVtCGUGQQVtCGUGQQVtCGUGQQVtCGUGQQVtCGUGQQVtC', + }, + }, } as const; diff --git a/src/sky/scatter.js b/src/sky/scatter.js new file mode 100644 index 00000000..5781c089 --- /dev/null +++ b/src/sky/scatter.js @@ -0,0 +1,249 @@ +// Atmospheric scattering tables, baked for a WebGL2 shader (layout: src/sky/README.md). +// Dependency-free apart from sun.js, so scripts/bake-sky.mjs and the tests run it under plain node. +// +// bakeScatter() -> { transmittance, sky_view, ... } with base64 RGB16F tables (half floats). +// +// A Hillaire-style single-scattering model (Rayleigh + Mie, spherical planet, observer 100 m up), +// reduced to two small tables: the transmittance of the whole atmosphere along a view direction, and +// the sky's radiance by sun elevation, view elevation and azimuth from the sun. sun.js stays the +// single source of the air: its zenith optical depths (TAU, per channel) are the vertical optical +// depth of this atmosphere, split into a grey aerosol part and a Rayleigh part. So the direct sun +// sun.js returns and the transmittance here agree at the zenith; only the air-mass law differs +// (Kasten-Young there, the spherical integral here). + +import { TAU } from './sun.js'; + +const RAD = Math.PI / 180; +const clamp = (x, lo, hi) => Math.min(hi, Math.max(lo, x)); + +// Geometry (km) and scale heights. +export const PLANET_RADIUS_KM = 6360; +export const ATMOSPHERE_HEIGHT_KM = 100; +export const OBSERVER_ALTITUDE_KM = 0.1; +const H_RAYLEIGH = 8; +const H_MIE = 1.2; +const MIE_ALBEDO = 0.9; // single-scattering albedo of the aerosol (EYEBALLED) +const MIE_G = 0.8; // Henyey-Greenstein asymmetry (Hillaire 2020) +const MIE_TAU = 0.02; // grey aerosol vertical optical depth (EYEBALLED); the rest of TAU is Rayleigh + +// Vertical optical depth of an exponential layer of scale height H over the atmosphere. +const column = (H) => H * (1 - Math.exp(-ATMOSPHERE_HEIGHT_KM / H)); +const BETA_RAYLEIGH = TAU.map((tau) => (tau - MIE_TAU) / column(H_RAYLEIGH)); // per km at sea level +const BETA_MIE_EXT = MIE_TAU / column(H_MIE); +const BETA_MIE_SCA = BETA_MIE_EXT * MIE_ALBEDO; + +// Table axes. +export const SUN_ELEVATIONS_DEG = [-18, -12, -8, -5, -2, 0, 3, 8, 15, 25, 45, 90]; +export const VIEW_SAMPLES = 32; +export const AZIMUTH_SAMPLES = 7; // 0, 30, ... 180 degrees from the sun +export const AZIMUTH_STEP_DEG = 30; +export const VIEW_MIN_DEG = -4; +export const VIEW_MAX_DEG = 90; +// Tables hold IEEE 754 half floats, little-endian, holding value * STORE_SCALE. The scale lifts the +// dark end (twilight sky radiance is ~1e-4) clear of the half-float subnormals; the brightest +// stored value, a Mie peak near 0.3 x 64, is far below the format's 65504. +export const STORE_SCALE = 64; + +/** View elevation in degrees of sample index i of VIEW_SAMPLES: finest at the horizon. */ +export const viewElevation = (i) => + VIEW_MIN_DEG + (VIEW_MAX_DEG - VIEW_MIN_DEG) * (i / (VIEW_SAMPLES - 1)) ** 2; + +/** Fractional index into SUN_ELEVATIONS_DEG of a sun elevation (clamped to the table): the z texel coordinate. */ +export const sunIndex = (elevationDeg) => { + const e = clamp( + elevationDeg, + SUN_ELEVATIONS_DEG[0], + SUN_ELEVATIONS_DEG[SUN_ELEVATIONS_DEG.length - 1] + ); + const k = Math.min(SUN_ELEVATIONS_DEG.findIndex((x) => x > e) - 1, SUN_ELEVATIONS_DEG.length - 2); + const i = k < 0 ? SUN_ELEVATIONS_DEG.length - 2 : k; + return i + (e - SUN_ELEVATIONS_DEG[i]) / (SUN_ELEVATIONS_DEG[i + 1] - SUN_ELEVATIONS_DEG[i]); +}; + +/** Inverse of viewElevation: the fractional sample index of an elevation. */ +export const viewIndex = (elevationDeg) => + Math.sqrt(clamp((elevationDeg - VIEW_MIN_DEG) / (VIEW_MAX_DEG - VIEW_MIN_DEG), 0, 1)) * + (VIEW_SAMPLES - 1); + +/** Round a non-negative number to the nearest IEEE 754 half float (ties to even); returns its 16 bits. */ +export function toHalf(x) { + if (!(x > 0)) return 0; + const u = new Uint32Array(new Float32Array([x]).buffer)[0]; + const exp = ((u >> 23) & 0xff) - 127 + 15; + const mant = u & 0x7fffff; + if (exp >= 31) return 0x7bff; + let half; + let rem; + let halfway; + if (exp <= 0) { + if (exp < -10) return 0; + const m = mant | 0x800000; + const shift = 14 - exp; + half = m >> shift; + rem = m & ((1 << shift) - 1); + halfway = 1 << (shift - 1); + } else { + half = (exp << 10) | (mant >> 13); + rem = mant & 0x1fff; + halfway = 0x1000; + } + return rem > halfway || (rem === halfway && half & 1) ? half + 1 : half; +} + +/** The exact value of 16 half-float bits (finite values only). */ +export const fromHalf = (h) => { + const e = (h >> 10) & 31; + const m = h & 1023; + return e === 0 ? m * 2 ** -24 : (1 + m / 1024) * 2 ** (e - 15); +}; + +/** Stored 16 bits of a value, and the value back out of them. */ +export const encode = (v) => toHalf(v * STORE_SCALE); +export const decode = (h) => fromHalf(h) / STORE_SCALE; + +const extinction = (h) => { + const r = Math.exp(-h / H_RAYLEIGH); + const m = Math.exp(-h / H_MIE); + return BETA_RAYLEIGH.map((b) => b * r + BETA_MIE_EXT * m); +}; + +// Distance along ray (origin o, unit direction d) to the sphere of radius r centred on the origin +// of the planet: the far root, or the near root if `near`. Infinity when it misses or is behind. +function sphereHit(o, d, r, near) { + const b = o[0] * d[0] + o[1] * d[1] + o[2] * d[2]; + const c = o[0] * o[0] + o[1] * o[1] + o[2] * o[2] - r * r; + const disc = b * b - c; + if (disc < 0) return Infinity; + const sq = Math.sqrt(disc); + const t = near ? -b - sq : -b + sq; + return t > 0 ? t : Infinity; +} + +// Optical depth per channel from p along d to the top of the atmosphere; Infinity if the planet is in the way. +function opticalDepthToSpace(p, d, steps) { + if (sphereHit(p, d, PLANET_RADIUS_KM, true) !== Infinity) return null; + const top = sphereHit(p, d, PLANET_RADIUS_KM + ATMOSPHERE_HEIGHT_KM, false); + const dt = top / steps; + const depth = [0, 0, 0]; + for (let k = 0; k < steps; k += 1) { + const t = (k + 0.5) * dt; + const q = [p[0] + d[0] * t, p[1] + d[1] * t, p[2] + d[2] * t]; + const e = extinction(Math.hypot(...q) - PLANET_RADIUS_KM); + for (let c = 0; c < 3; c += 1) depth[c] += e[c] * dt; + } + return depth; +} + +const observer = [0, PLANET_RADIUS_KM + OBSERVER_ALTITUDE_KM, 0]; +const viewDirection = (elevationDeg, azimuthDeg) => { + const e = elevationDeg * RAD; + const a = azimuthDeg * RAD; + return [Math.cos(e) * Math.cos(a), Math.sin(e), Math.cos(e) * Math.sin(a)]; +}; +const sunDirection = (elevationDeg) => [ + Math.cos(elevationDeg * RAD), + Math.sin(elevationDeg * RAD), + 0, +]; + +/** Transmittance per channel from the observer along an elevation to space; zeros where the ground is in the way. */ +export function transmittance(elevationDeg) { + const depth = opticalDepthToSpace(observer, viewDirection(elevationDeg, 0), 64); + return depth ? depth.map((x) => Math.exp(-x)) : [0, 0, 0]; +} + +const rayleighPhase = (c) => (3 / (16 * Math.PI)) * (1 + c * c); +const miePhase = (c) => + (1 - MIE_G * MIE_G) / (4 * Math.PI * (1 + MIE_G * MIE_G - 2 * MIE_G * c) ** 1.5); + +/** + * Sky radiance per channel, per steradian, for a sun of unit irradiance (so the consumer multiplies by + * its own sun colour), in a view direction at an elevation and an azimuth measured from the sun. + * Single scattering, no ground bounce, no sun disc. + */ +export function skyRadiance(sunElevationDeg, viewElevationDeg, azimuthDeg) { + const v = viewDirection(viewElevationDeg, azimuthDeg); + const s = sunDirection(sunElevationDeg); + const cosTheta = v[0] * s[0] + v[1] * s[1] + v[2] * s[2]; + const pr = rayleighPhase(cosTheta); + const pm = miePhase(cosTheta); + const ground = sphereHit(observer, v, PLANET_RADIUS_KM, true); + const top = sphereHit(observer, v, PLANET_RADIUS_KM + ATMOSPHERE_HEIGHT_KM, false); + const end = Math.min(ground, top); + const steps = 24; + const dt = end / steps; + const along = [0, 0, 0]; + const radiance = [0, 0, 0]; + for (let k = 0; k < steps; k += 1) { + const t = (k + 0.5) * dt; + const q = [observer[0] + v[0] * t, observer[1] + v[1] * t, observer[2] + v[2] * t]; + const h = Math.hypot(...q) - PLANET_RADIUS_KM; + const e = extinction(h); + for (let c = 0; c < 3; c += 1) along[c] += e[c] * dt; + const toSun = opticalDepthToSpace(q, s, 12); + if (toSun) { + const r = Math.exp(-h / H_RAYLEIGH); + const m = Math.exp(-h / H_MIE); + for (let c = 0; c < 3; c += 1) { + const scatter = BETA_RAYLEIGH[c] * r * pr + BETA_MIE_SCA * m * pm; + radiance[c] += Math.exp(-along[c] - toSun[c]) * scatter * dt; + } + } + } + return radiance; +} + +const base64 = (halves) => { + const bytes = new Uint8Array(halves.length * 2); + halves.forEach((h, i) => { + bytes[i * 2] = h & 0xff; // little-endian + bytes[i * 2 + 1] = h >> 8; + }); + let s = ''; + for (let i = 0; i < bytes.length; i += 0x8000) { + s += String.fromCharCode(...bytes.subarray(i, i + 0x8000)); + } + return btoa(s); +}; + +export function bakeScatter() { + const view = Array.from({ length: VIEW_SAMPLES }, (_, i) => viewElevation(i)); + + const t = new Uint16Array(VIEW_SAMPLES * 3); + view.forEach((e, i) => transmittance(e).forEach((v, c) => (t[i * 3 + c] = encode(v)))); + + // x = azimuth, y = view elevation, z = sun elevation, three half floats per texel: a 3D texture. + const sky = new Uint16Array(SUN_ELEVATIONS_DEG.length * VIEW_SAMPLES * AZIMUTH_SAMPLES * 3); + let at = 0; + SUN_ELEVATIONS_DEG.forEach((sun) => + view.forEach((e) => { + for (let a = 0; a < AZIMUTH_SAMPLES; a += 1) { + skyRadiance(sun, e, a * AZIMUTH_STEP_DEG).forEach((v) => { + sky[at] = encode(v); + at += 1; + }); + } + }) + ); + + return { + planet_radius_km: PLANET_RADIUS_KM, + atmosphere_height_km: ATMOSPHERE_HEIGHT_KM, + observer_altitude_km: OBSERVER_ALTITUDE_KM, + // Extinction per km at the observer, per channel, for aerial perspective: exp(-k * distance_km). + extinction_per_km: extinction(OBSERVER_ALTITUDE_KM), + format: 'half-float', + store_scale: STORE_SCALE, + view_elevation_deg: { min: VIEW_MIN_DEG, max: VIEW_MAX_DEG, samples: VIEW_SAMPLES }, + transmittance: { width: VIEW_SAMPLES, height: 1, channels: 3, data: base64(t) }, + sky_view: { + sun_elevations_deg: SUN_ELEVATIONS_DEG, + azimuth_step_deg: AZIMUTH_STEP_DEG, + width: AZIMUTH_SAMPLES, + height: VIEW_SAMPLES, + depth: SUN_ELEVATIONS_DEG.length, + channels: 3, + data: base64(sky), + }, + }; +} diff --git a/src/sky/scatter.test.ts b/src/sky/scatter.test.ts new file mode 100644 index 00000000..9ebcb014 --- /dev/null +++ b/src/sky/scatter.test.ts @@ -0,0 +1,160 @@ +import { describe, expect, it } from 'vitest'; +import { + AZIMUTH_SAMPLES, + AZIMUTH_STEP_DEG, + STORE_SCALE, + SUN_ELEVATIONS_DEG, + VIEW_SAMPLES, + bakeScatter, + decode, + fromHalf, + skyRadiance, + sunIndex, + toHalf, + transmittance, + viewElevation, + viewIndex, +} from './scatter.js'; +import { TAU } from './sun.js'; +import { baked } from './baked'; + +const table = (b64: string): Uint16Array => { + const bin = atob(b64); + const bytes = Uint8Array.from(bin, (c) => c.charCodeAt(0)); + return new Uint16Array(bytes.buffer); +}; +// Checksums of the baked halves, so a change to the model or the encoding is a deliberate re-pin. +const SKY_CHECKSUM = 622927107; +const TRANS_CHECKSUM = 1498461058; +// eslint-disable-next-line @typescript-eslint/no-explicit-any +const scatter = baked.scatter as any; +const sky = table(scatter.sky_view.data); +const trans = table(scatter.transmittance.data); +const sunIndex_ = (deg: number) => SUN_ELEVATIONS_DEG.indexOf(deg); +// Texel of the sky table: x = azimuth, y = view elevation, z = sun elevation; value = three halves. +const texel = (s: number, v: number, a: number) => + [0, 1, 2].map((c) => decode(sky[((s * VIEW_SAMPLES + v) * AZIMUTH_SAMPLES + a) * 3 + c])); + +describe('half floats', () => { + it('round-trip within half precision across the table range, ties to even', () => { + [6.1e-5, 1e-4, 0.0123, 0.3, 1, 20, 1000].forEach((x) => + expect(Math.abs(fromHalf(toHalf(x)) - x) / x).toBeLessThan(5.1e-4) + ); + expect(toHalf(0)).toBe(0); + expect(toHalf(1)).toBe(0x3c00); + expect(toHalf(2 ** -24)).toBe(1); // smallest subnormal + expect(toHalf(1 + 2 ** -11)).toBe(0x3c00); // exact tie goes to the even neighbour + expect(toHalf(1 + 3 * 2 ** -11)).toBe(0x3c02); + expect(fromHalf(0x7bff)).toBe(65504); + }); + + it('keeps twilight radiance, ~1e-4, within 0.1% thanks to the store scale', () => { + const x = 1.234e-4; + expect(Math.abs(decode(toHalf(x * STORE_SCALE)) - x) / x).toBeLessThan(1e-3); + }); +}); + +describe('sunIndex', () => { + it('is the fractional position of an elevation on the sun axis, clamped to the table', () => { + expect(sunIndex(-18)).toBe(0); + expect(sunIndex(90)).toBe(SUN_ELEVATIONS_DEG.length - 1); + expect(sunIndex(-100)).toBe(0); + expect(sunIndex(100)).toBe(SUN_ELEVATIONS_DEG.length - 1); + expect(sunIndex(0)).toBe(sunIndex_(0)); + expect(sunIndex(1.5)).toBeCloseTo(sunIndex_(0) + 0.5, 12); // halfway between 0 and 3 + expect(sunIndex(35)).toBeCloseTo(sunIndex_(25) + 0.5, 12); // halfway between 25 and 45 + }); +}); + +describe('scattering tables', () => { + it('are what scripts/bake-sky.mjs would write now (run npm run bake:sky if not)', () => { + expect(JSON.parse(JSON.stringify(bakeScatter()))).toEqual(JSON.parse(JSON.stringify(scatter))); + }); + + it('have the documented layout and size', () => { + expect(scatter.format).toBe('half-float'); + expect(scatter.sky_view).toMatchObject({ + width: AZIMUTH_SAMPLES, + height: VIEW_SAMPLES, + depth: SUN_ELEVATIONS_DEG.length, + channels: 3, + azimuth_step_deg: AZIMUTH_STEP_DEG, + }); + expect(sky.length).toBe(AZIMUTH_SAMPLES * VIEW_SAMPLES * SUN_ELEVATIONS_DEG.length * 3); + expect(trans.length).toBe(VIEW_SAMPLES * 3); + expect(sky.length * 2).toBe(16128); // bytes: the budget the README states + }); + + it('pin the view-elevation axis and its inverse', () => { + expect(viewElevation(0)).toBe(-4); + expect(viewElevation(VIEW_SAMPLES - 1)).toBe(90); + expect(viewElevation(6)).toBeCloseTo(-4 + 94 * (6 / 31) ** 2, 10); + [0, 5, 6, 17, 31].forEach((i) => expect(viewIndex(viewElevation(i))).toBeCloseTo(i, 9)); + }); + + it('transmittance to space agrees with sun.js at the zenith and falls toward the horizon', () => { + const top = VIEW_SAMPLES - 1; + [0, 1, 2].forEach((c) => + expect(Math.abs(decode(trans[top * 3 + c]) / Math.exp(-TAU[c]) - 1)).toBeLessThan(0.01) + ); + expect(transmittance(90).map((v, c) => v / Math.exp(-TAU[c]))).toSatisfy((r: number[]) => + r.every((x) => Math.abs(x - 1) < 0.01) + ); + // blue goes first: lower view, lower blue, and redder light at the horizon + let prev = Infinity; + for (let i = top; i >= 8; i -= 1) { + const blue = decode(trans[i * 3 + 2]); + expect(blue).toBeLessThanOrEqual(prev); + prev = blue; + } + expect(decode(trans[8 * 3]) / decode(trans[8 * 3 + 2])).toBeGreaterThan(10); // ~2 degrees up: sunset red + expect(decode(trans[0])).toBe(0); // -4 degrees looks at the ground + }); + + it('colour the daytime sky blue, and a low sun redder toward it than away from it', () => { + const noon = texel(sunIndex_(45), VIEW_SAMPLES - 1, 0); + expect(noon[2]).toBeGreaterThan(noon[1]); + expect(noon[1]).toBeGreaterThan(noon[0]); + const horizon = Math.round(viewIndex(5)); + const toward = texel(sunIndex_(3), horizon, 0); + const away = texel(sunIndex_(3), horizon, AZIMUTH_SAMPLES - 1); + expect(toward[0] / toward[2]).toBeGreaterThan(5); + expect(toward[0] / toward[2]).toBeGreaterThan(2 * (away[0] / away[2])); + expect(toward[0]).toBeGreaterThan(5 * away[0]); // the glow is around the sun + }); + + it('fade through twilight to black below the planet-shadow limit', () => { + const mid = Math.round(viewIndex(30)); + const lum = (s: number) => texel(sunIndex_(s), mid, 3).reduce((x, y) => x + y, 0); + expect(lum(8)).toBeGreaterThan(lum(3)); + expect(lum(3)).toBeGreaterThan(lum(-2)); + expect(lum(-2)).toBeGreaterThan(lum(-5)); + expect(lum(-5)).toBeGreaterThan(0); + expect(lum(-12)).toBe(0); + expect(lum(-18)).toBe(0); + }); + + it('hold the model to half-float precision, with no quantised steps along the horizon glow', () => { + // Every table texel off a handful of sampled ones is the model's own value to < 0.2%. + [ + [sunIndex_(3), 6, 0], + [sunIndex_(3), 9, 2], + [sunIndex_(0), 8, 0], + [sunIndex_(45), 20, 3], + [sunIndex_(-2), 12, 1], + ].forEach(([s, v, a]) => { + const direct = skyRadiance(SUN_ELEVATIONS_DEG[s], viewElevation(v), a * AZIMUTH_STEP_DEG); + texel(s, v, a).forEach((got, c) => expect(Math.abs(got / direct[c] - 1)).toBeLessThan(2e-3)); + }); + // Moving up from the horizon toward the sun, red drops smoothly: distinct, strictly decreasing. + const s = sunIndex_(3); + const reds = Array.from({ length: 8 }, (_, k) => texel(s, 7 + k, 0)[0]); + reds.slice(1).forEach((r, k) => expect(r).toBeLessThan(reds[k])); + expect(new Set(reds).size).toBe(reds.length); + }); + + it('are the same bytes as when pinned', () => { + const sum = (a: Uint16Array) => a.reduce((x, y) => (x * 31 + y) % 2147483647, 7); + expect([sum(sky), sum(trans)]).toEqual([SKY_CHECKSUM, TRANS_CHECKSUM]); + }); +}); diff --git a/src/sky/sun.js b/src/sky/sun.js index 55cad1ac..3173eb4d 100644 --- a/src/sky/sun.js +++ b/src/sky/sun.js @@ -57,7 +57,7 @@ export const sunElevation = (time, lat, lon) => sunPosition(time, lat, lon).elev // Optical depth of the clear atmosphere at the zenith, per channel (R, G, B). // Rayleigh at 650/550/450 nm plus a little ozone and haze (EYEBALLED to keep // low sun orange rather than brown). -const TAU = [0.06, 0.12, 0.26]; +export const TAU = [0.06, 0.12, 0.26]; // Relative air mass, Kasten & Young 1989. Elevation in degrees, >= 0. const airMass = (el) => 1 / (Math.sin(el * RAD) + 0.50572 * Math.pow(el + 6.07995, -1.6364));