Sunday, September 27, 2026

Audit File for OpenSpace from Claude Code

I finally got some traction on making changes to OpenSpace.

I had Claude Code do a revision of my little lunar retroreflection hack and it created the following audit file.

Here are a few pictures from the real time simulation.

CT18.B1.Openspace — AUDIT

Project: CT18.B1.Openspace · Target: OpenSpace 0.21.4 (installed at D:\OpenSpace-21.4) · Date: 2026-09-26 Purpose of this file: let a reader who does not trust AI-written code check, claim by claim, where every number and every equation in the shader comes from, what was actually verified, and what was assumed.

0. How to read this audit

Every physical number carries a source tag (the same tags appear in tools/ct18b1_config.py, the only place numbers are defined):

tag	meaning
[V:xxx]	The value was read from the cited source during this build (web fetch of the paper, report, data page or search-result text).
[L:xxx]	Standard literature value or equation that was not re-fetched during this build (textbook constants, the Hapke equations, IAU rotation model). Treat with slightly less confidence; cross-checks are listed where they exist.
[D:xxx]	Derived in this project from [V] sources; the derivation is written out in section 3.4.
[A]	Assumption or design choice with no measured lunar source. Every [A] is listed in section 9.

"Verified" in this document means one of three different things, and the text always says which:

Code-transcription verified — the GLSL that OpenSpace will run was executed (software OpenGL 4.5) and compared number-for-number with an independent Python implementation of the published equations.
Physics-anchor verified — a model output that was not used to fit any parameter agrees with a published measurement (e.g. normal albedo).
Not verifiable here — OpenSpace itself could not be run in the build environment. The reference images come from an emulator that compiles the same shader files through unmodified copies of OpenSpace's own shader wrapper code (section 5.6). You must run the test keys in OpenSpace to close this gap (section 8).

The automated results are in verification/VERIFICATION_REPORT.md (33 of 33 checks passed on the delivered build).

1. What "retroreflective" means for the Moon, physically

The lunar surface is not a mirror-like retroreflector (like a bicycle reflector or the Apollo laser arrays). It is a porous, dark, fine-grained regolith that sends back noticeably more light toward the direction it came from than in any other direction. This opposition effect has two established physical causes, both modelled here:

Shadow-hiding opposition effect (SHOE). Every grain casts a shadow on the grains behind it. Seen from the direction of the light source, those shadows are hidden behind the grains that cast them, so the surface looks brighter. The peak width is set by porosity and grain-size distribution. Hapke (1986) gives the analytic form used here ([E3]).
Coherent-backscatter opposition effect (CBOE). Light scattered along two time-reversed paths through the grains interferes constructively exactly in the backward direction, giving a narrow peak. Hapke (2002) gives the form used here ([E4]).

Observations that constrain these on the real Moon, and how each is used:

Clementine found the Moon's disk-integrated brightness rises by more than 40 % between phase angles of 4° and 0°, and attributed most of it to shadow hiding [V:BUR96]. Used to derive the SHOE width h_S (section 3.4).
Helfenstein, Veverka & Hillier (1997) found the lunar opposition effect is best described by a narrow coherent- backscatter peak (most pronounced below ~2°) plus a broad shadow-hiding peak (up to ~20°) [V:HV97]. That is exactly the two-component structure of the model.
LROC WAC data give a CBOE maximum amplitude near 8 % and widths of 1–4°, wider for maria and in blue light [V:VEL16]. Used for B_C0 and h_C.

The rover lamp sits 0.1016 m (4 in) above the camera, so a ground point at distance d is seen at a phase angle g ≈ 0.1016/d rad (1° at 5.8 m, 0.25° at 23 m). The whole lamp-lit scene is therefore inside the opposition peak, and its brightness is enhanced up to (1+B_S0)(1+B_C0) ≈ 4.1× over what the same surface would give without the opposition effect. Test T02/T03 measures this factor directly (check C05: ×3.45 at 3 m rising to ×3.97 at 40 m).

2. References

Status column: V = content read during this build (what was read is stated); L = not re-fetched.

tag	reference	status	what was taken from it
HV87	Helfenstein, P. & Veverka, J. (1987) Photometric properties of lunar terrains derived from Hapke's equation, Icarus 72, 342–357. NASA NTRS copy: https://ntrs.nasa.gov/api/citations/19870013999/downloads/19870013999.pdf	V	Table 1, "dark terrains" column (maria): w = 0.12, h = 0.12, S(0) = 0.51, b = 0.41, c = 0.10, θ̄ = 8.1°. Relation B0 = S(0)/[w(1+b+c)]. Legendre form of P(g). Data span 0–150° phase, visual band.
SATO14	Sato, H. et al. (2014) Resolved Hapke parameter maps of the Moon, JGR Planets, doi:10.1002/2013JE004580. LROC data page: https://data.lroc.im-ldi.com/lroc/view_rdr/WAC_HAPKEPARAMMAP	V (data page; paper paywalled)	Fixed values in the WAC Hapke maps: θ̄ = 23.657°, filling factor 1.0, B_C0 = 0, h_C = 1.0. Used: θ̄ as the alternative roughness (checked in A03), K = 1. The per-pixel w, b, c maps could not be downloaded (server blocked).
BUR96	Buratti, B. J., Hillier, J. K. & Wang, M. (1996) The lunar opposition surge: observations by Clementine, Icarus 124, 490–499. https://www.sciencedirect.com/science/article/abs/pii/S0019103596902250	V (abstract)	Brightness rises more than 40 % from g = 4° to 0°; surge ~3–4 % larger at 0.41 µm than at 1.00 µm; ~10 % greater in highlands; shadow hiding is the principal cause.
HV97	Helfenstein, P., Veverka, J. & Hillier, J. (1997) The lunar opposition effect: a test of alternative models, Icarus 128, 2–14. https://www.sciencedirect.com/science/article/pii/S0019103597957262	V (abstract)	Two-component opposition effect: narrow CBOE (α < 2°) + broad SHOE (α < 20°).
VEL16	Velikodsky, Yu. I. et al. (2016) Opposition effect of the Moon from LROC WAC data, Icarus 275, 1 (ADS 2016Icar..275....1V). https://www.sciencedirect.com/science/article/abs/pii/S0019103516300392	V (abstract)	CBOE maximum amplitude near 8 %; widths 1–4° (highlands red ≈1.2°, maria blue ≈3.9°); coherent backscatter is probably the main visible-light mechanism near zero phase.
SP8023	NASA SP-8023 (1969) Lunar Surface Models (NASA Space Vehicle Design Criteria). https://ntrs.nasa.gov/api/citations/19700009596/downloads/19700009596.pdf	V	Equilibrium crater law N₀ = K D^n with n ≈ −2, K = 10⁻²; mare normal albedo 0.07–0.12, average 0.095.
ROS11	Rosenburg, M. A. et al. (2011) Global surface slopes and roughness of the Moon from the Lunar Orbiter Laser Altimeter, JGR 116, E02001. https://authors.library.caltech.edu/22811/	V (abstract)	Median Hurst exponent: maria 0.76, highlands 0.95, over ~17 m – 2.7 km; transition often near 1 km.
ROW71	Rowan, L. C., McCauley, J. F. & Holm, E. A. (1971) Lunar terrain mapping and relative-roughness analysis, USGS Professional Paper 599-G. https://ntrs.nasa.gov/api/citations/19720008109/downloads/19720008109.pdf	V	Mare median slope component ≈ 3.6° at 1 m (Ranger VII); typical maria below 1° over 0.75 km.
DI16	Di, K. et al. (2016) Rock size-frequency distribution analysis at the Chang'E-3 landing site, Planet. Space Sci. 120, 103–112. http://www.pmrslab.cn/publications/publications/Rock-size-frequency-distribution-analysis-at-the-Chang-E-3-landing-site_2016_Planetary-and-Space-Science.pdf	V	Cumulative fractional rock area F(D) = k exp(−qD), k = 0.0125, q = 1.743 m⁻¹; rock height H = 0.2347 D + 0.0039 m; N(D ≥ 0.1 m) = 0.17 m⁻² at CE-3, 0.20 Surveyor VI, 0.34 Surveyor I.
DAR05	Darula, S., Kittler, R. & Gueymard, C. A. (2005) Reference luminous solar constant and solar luminance for illuminance calculations, Solar Energy 79, 559. https://www.sciencedirect.com/science/article/abs/pii/S0038092X05000320	V (search-result text)	Luminous solar constant 133.8 klx.
ROB26	Inferring and interpreting the visual geometric albedo and phase function of Earth, PSJ 7, 12 (2026); arXiv:2507.22258. https://arxiv.org/html/2507.22258v1	V	Earth V-band geometric albedo 0.226 (B 0.277, R 0.221); visual geometric albedo 0.242; phase integral 1.22; 30–40 % below the long-quoted 0.367.
OS	OpenSpace v0.21.4 source (tag `releases/v0.21.4`) and Ghoul (commit b737947). https://github.com/OpenSpace/OpenSpace	V (source read)	RenderableModel custom-shader interface, uniform names, SceneGraphLightSource math, texture loader (vertical flip), HDR tone map, depth encoding, default window FOV (80° × 50.534°), default exposure 3.7 and gamma 0.95.
H81, H84, H86, H02, H12	Hapke, B. (1981) JGR 86, 3039; (1984) Icarus 59, 41 (roughness); (1986) Icarus 67, 264 (SHOE); (2002) Icarus 157, 523 (CBOE, H function); (2012) Theory of Reflectance and Emittance Spectroscopy, 2nd ed., Cambridge.	L	Equations [E1]–[E6]. Cross-checked in the build: the H-function approximation against an exact numerical solution (A04), the roughness correction against its analytic limits (θ̄ = 0 gives S = 1; i = e = g = 0 gives S = 1), and the complete GLSL against an independent Python implementation (A03).
IAU2009	Archinal, B. A. et al. (2011) Report of the IAU WGCCRE: 2009, Celest. Mech. Dyn. Astr. 109, 101.	L	Lunar rotation model (α₀, δ₀, W series). Cross-checked against DE421's integrated librations: agreement 0.003° (A09). Mean lunar radius 1737.4 km.
DE421	Folkner, W. M., Williams, J. G. & Boggs, D. H. (2009) The Planetary and Lunar Ephemeris DE421, IPN PR 42-178; data from the `de421` Python package.	L (data used directly)	Sun, Earth, Moon positions and lunar librations for the test epochs.
PIKE77	Pike, R. J. (1977) size-dependence of crater morphology; fresh simple craters depth/diameter ≈ 0.2, rim height ≈ 0.036 D.	L	Crater shape (fresh end of the degradation range).
A11	Apollo 11 lunar module landing coordinates, 0.67409° N, 23.47298° E (LROC).	L	Site location only. The terrain is statistical, not the real Apollo 11 terrain.

Sources that could not be reached from the build environment (NASA PDS, LROC data servers, USGS Astrogeology, JPL NAIF and the Wiley/AGU full texts were blocked by the network allowlist): real LOLA/LROC DEMs, the full Sato et al. parameter maps, Hapke et al. (2012) JGR, and Hapke, Nelson & Smythe (1998). Section 7.7 explains how to plug in a real DEM if you download one yourself.

3. The reflectance model and where it lives in the source

3.1 Equations

The shader computes the radiance factor RADF = I/F = π·r (r = Hapke bidirectional reflectance):

[E1] r = K·w/(4π) · μ₀ₑ/(μ₀ₑ+μₑ) · [ P(g)·(1 + B_S(g)) + H(μ₀ₑ/K)·H(μₑ/K) − 1 ] · (1 + B_C(g)) · S(i,e,ψ) — Hapke (2012) book form [H12]. Sato et al. (2014) and the MDPI/Chang'E-1 photometric-correction papers use the same structure.
[E2] P(g) = 1 + b·cos g + c·(3cos²g − 1)/2 — Legendre phase function, exactly as in HV87 so their b, c can be used unchanged.
[E3] B_S(g) = B_S0 / (1 + tan(g/2)/h_S) — SHOE [H86].
[E4] B_C(g) = B_C0 · [1 + (1 − e^(−y))/y] / [2(1+y)²], y = tan(g/2)/h_C — CBOE [H02].
[E5] H(x) = [1 − w·x·(r₀ + ½(1 − 2r₀x)·ln((1+x)/x))]⁻¹, r₀ = (1−γ)/(1+γ), γ = √(1−w) — [H02].
[E6] μ₀ₑ, μₑ, S(i,e,ψ): macroscopic roughness with mean slope angle θ̄, including the i ≤ e and i > e branches and the azimuth factor f(ψ) = exp(−2 tan(ψ/2)) — [H84].

3.2 Correspondence: equation → GLSL → independent Python

eq.	GLSL (shipped, `shaders/ct18b1_hapke.glsl`)	Python reference (`tools/hapke_reference.py`)	verification
E1	`hapkeRadf()`; vector front end `hapkeRadfVec()`	`radf()`, `radf_vectors()`	A03: 60,000 geometries × 5 parameter sets, worst relative difference 4.2e-6
E2	`hapkeP()`	`P_legendre()`	A03 (probe outputs P via RADF)
E3	`hapkeBS()`	`B_shoe()`	A03 (probe outputs B_S separately)
E4	`hapkeBC()`	`B_cboe()`	A03 (probe outputs B_C separately); A07 HWHM
E5	`hapkeH()`	`H_2002()`; exact `H_exact()`	A03; A04: approximation vs exact Chandrasekhar solution, max error 7.7e-5 at w = 0.12
E6	`hapkeRoughness()`, `hapkeE1()`, `hapkeE2()`	`roughness()`	A03 including θ̄ = 0, 8.1°, 23.657°

Independence caveat, stated plainly: both implementations were written by the same author (the AI) from the same understanding of the equations. A03 therefore proves the GLSL does what the Python does, including float32 behaviour on a GPU. It does not by itself prove the equations were remembered correctly. That is what the physics anchors in 3.5 are for: they compare model outputs with measurements that played no part in choosing the equations.

3.3 Parameter values used by the shader

symbol	GLSL constant	value	tag	note
w	`CT18B1_W`	0.12	V:HV87	dark terrains (maria)
b	`CT18B1_B`	0.41	V:HV87	b > 0: back-scattering particles
c	`CT18B1_C`	0.10	V:HV87
B_S0	`CT18B1_BS0`	2.8146	D:BS0	= 0.51 / (0.12 × 1.51), see 3.4
h_S	`CT18B1_HS`	0.066877	D:HS	solved from Clementine, see 3.4
B_C0	`CT18B1_BC0`	0.08	V:VEL16
h_C	`CT18B1_HC`	0.048915	D:HC	HWHM 2.0° [A] inside the 1–4° range of VEL16
θ̄	`CT18B1_THETA_BAR`	8.1° (0.14137 rad)	V:HV87	alternative 23.657° [V:SATO14], see 3.6
K	`CT18B1_K`	1.0	V:SATO14	no porosity correction (as in the WAC maps)

Check A01 re-reads the generated ct18b1_params.glsl and confirms every constant equals the configuration value (worst relative difference 2.9e-9, i.e. only decimal printing).

3.4 Derivations

[D:BS0] HV87 report the product S(0) = B₀·w·P(0) and give B₀ = S(0)/[w(1+b+c)]. With their dark-terrain values, B_S0 = 0.51/(0.12 × 1.51) = 2.815. In Hapke's 1986 theory B_S0 ≤ 1 for pure shadow hiding. Values above 1 are what an empirical fit produces when coherent backscatter and other effects are folded into the one SHOE term. This is kept deliberately because it is the published fit to real lunar data (check A08 reports it rather than hiding it). The consequence is that the model's total surge is realistic (A06) even though the split between SHOE and CBOE is not unique.
[D:HC] h_C is found by bisection so that B_C(g)/B_C0 = ½ at g = 2.0°. The 2.0° half width is an [A] choice inside VEL16's 1–4° range (maria, visual band). hapke_reference.solve_hc(); A02 re-solves it every run.
[D:HS] HV87's h = 0.12 was fitted to Earth-based data that cannot reach below about 2° phase. With it, and with CBOE added, the model's disk-integrated surge from 4° to 0° is only ×1.275 (×1.208 without CBOE; computed with hapke_reference.disk_integrated), contradicting Clementine's >40 % [V:BUR96]. h_S is therefore re-solved (bisection on a numerically integrated full disk, hapke_reference.solve_hs()) so that I(0°)/I(4°) = 1.42 with every other parameter held at the values above. Result: h_S = 0.0669, which is close to HV87's own disk-integrated h of 0.07. A02 re-solves it every run; A06 confirms the ratio 1.420. verification/plots/phase_curve.png shows the resulting phase curve next to the unmodified HV87 set.
[D:ROCKK] "More rocks" (the brief's B1 change) was implemented as the rockiest site in DI16's comparison, Surveyor I: k = 0.0125 × 0.34/0.17 = 0.025. The analytic N(D ≥ 0.1 m) of the resulting law is 0.337 m⁻² (A13), and the generated corridor has 0.327 m⁻² (A12).

3.5 Physics anchors (outputs that were not fitted)

check	model output	measurement	result
A05	normal albedo π·r(0,0,0) = 0.0948	SP-8023 mare normal albedo 0.07–0.12, average 0.095	agrees to 0.2 % of the average (not used in any fit)
A06	disk-integrated I(0°)/I(4°) = 1.420	Clementine: > 1.40	pass (h_S was fitted to this; listed for completeness)
A07	CBOE amplitude 8 %, HWHM 2.0°	VEL16: ~8 %, 1–4°	pass (input, listed for completeness)
C07b	full disk from infinity at g = 0.5°: brightness at 0.8–0.9 radius / centre = 1.00 (Lambert: 0.53)	The full Moon shows essentially no limb darkening (classical observation)	pass
C06	ground brightness around the anti-solar point in T08 falls 30 % from 0–0.5° to 4–6°	Heiligenschein seen around astronauts' shadows; surge shape per BUR96/HV97	pass (qualitative anchor)
C09	earthshine at full Earth = 6.2e-5 × sunlight	order 1e-4 is the textbook value	pass

3.6 Choices between published parameter sets

Roughness θ̄: HV87's 8.1° belongs with their w, b, c and is used by default because the set was fitted together. The LROC WAC maps fix θ̄ = 23.657° [V:SATO14]. Test T19 turns roughness off entirely to show its size; the GLSL is verified at both 8.1° and 23.657° (A03). If you prefer the LROC value, change HAPKE_THETA_DEG in tools/ct18b1_config.py and rerun the generator (section 10).
Wavelength: all anchors are visual-band (V, ~0.55 µm), so the shader outputs a single grey V-band value in R = G = B. The Moon is slightly red, and the surge is ~3–4 % larger in blue [V:BUR96]; neither is modelled (section 9).

4. Light sources and exposure

All light is expressed as irradiance relative to sunlight at 1 AU (dimensionless), so the output is the radiance factor the surface would have under 1 AU sunlight, scaled by a camera gain. There is no ambient term: check C03 shows shadowed ground in T01 is at 0.63 % of the sunlit median, and equals the predicted earthshine + lamp contribution to within 0.2 %.

source	model (GLSL function)	numbers	tags
Sun	parallel light, finite disk; `ct18b1SourceVisibility()` gives the visible fraction of a disk of radius 0.2666° above the terrain horizon (exact disk-over-chord area)	irradiance 1 at 1 AU; the real Sun–Moon distance varies by ±1.7 % (±3.4 % irradiance), not applied	L:IAU (959.63″ radius)
Earth	Lambert-sphere phase law, `ct18b1EarthIrradiance()`: E/E_sun = p_V · Φ(α) · (R_E/d)², α = Earth's phase angle computed in the shader from the Sun and Earth directions	p_V = 0.226; R_E/d = 6371/384400; Φ(α) = [sin α + (π−α)cos α]/π; Earth disk radius 0.95° for shadows	V:ROB26, L:IUGG, L:IAU
Rover lamp	point source 0.1016 m above the camera along the camera's up axis; Gaussian beam 40° FWHM on the view axis; `ct18b1LampIrradiance()`: E/E_sun = I·beam/(d² · 133 800 lx)	20 000 cd peak	[A] design, V:DAR05
Exposure	camera gain = 10^(6a), a = RenderableModel `AmbientIntensity` (re-purposed, 5.2)	daylight a = 0.10 (×3.98), night lamp a = 0.50 (×1000), night earthshine a = 0.72 (×2.1e4)	[A] like choosing an f-stop

End-to-end checks recompute every lit pixel in Python from the GPU's surface positions and normals: lamp (C01, 99.7 % of 568k pixels within 3 %, median ratio 0.9999), Sun (C02, 99.8 % of 451k fully sunlit pixels within 3 %, median 1.0001), earthshine (C08, 75th percentile 0.9999), whole body (C07, 100 % within 3 %).

Note on Earth's phase law: ROB26 give Earth's phase integral as 1.22, while a Lambert sphere has 1.5. The Lambert shape is therefore only an approximation away from full Earth. At the night-test epoch Earth is 3.8° from full, where the approximation matters least.

5. OpenSpace integration (read from the 0.21.4 source)

5.1 How a custom shader gets into OpenSpace without recompiling

RenderableModel accepts VertexShader and FragmentShader file paths in the asset (modules/base/rendering/renderablemodel.cpp). The fragment file is wrapped by OpenSpace's shaders/render.frag → shaders/framebuffer/renderframebuffer.frag, which calls our getFragment(). Both CT18.B1 renderables (CT18B1_Globe, CT18B1_Site) use this route, so no OpenSpace binaries are modified.

5.2 Re-purposed uniforms (the only runtime knobs OpenSpace offers a custom model shader)

RenderableModel property	uniform	CT18.B1 meaning
`AmbientIntensity`	`ambientIntensity`	camera exposure a → gain 10^(6a) (not an ambient light)
`DiffuseIntensity`	`diffuseIntensity`	lamp dimmer, 0 = off, 1 = 20 000 cd
`SpecularIntensity`	`specularIntensity`	mode = round(10 s): 10 physics, 0 OE off, 1 SHOE only, 2 CBOE only, 3 Lambert, 4 lamp phase-angle map, 5 Sun phase-angle map, 6 Sun visibility, 7 level map, 8 rock mask, 9 roughness off
`LightSources.Sun.Intensity`	`lightIntensities[0]`	Sun on (1) / off (0)
`LightSources.Earth.Intensity`	`lightIntensities[1]`	Earth on (1) / off (0)

Do not disable a light source with its Enabled flag: RenderableModel then compacts the arrays and the Earth would arrive as light 0. Set its Intensity to 0 instead (the test keys do this).

5.3 A bug in OpenSpace 0.21.4's SceneGraphLightSource, and the work-around

modules/base/lightsource/scenegraphlightsource.cpp, directionViewSpace(), transforms the direction (light position − node position) with combinedViewMatrix * vec4(direction, 1.0). Because w = 1, the camera translation is applied to a direction vector. The result is proportional to R(L − N − C), where C is the camera's world position measured from the solar-system barycentre (~1.5 × 10¹¹ m when near the Moon):

Sun: L ≈ barycentre, N ≈ C ≈ Moon, so R(L − 2·Moon) is nearly parallel to the truth. The error is about 0.2–0.4°, comparable to the Sun's radius.
Earth: L − N is about 3.8 × 10⁸ m while C is about 1.5 × 10¹¹ m, so the "Earth" light points roughly toward the Sun. Earthshine would be completely wrong.

Work-around (no OpenSpace code changed): each light source follows a helper node (CT18B1_SunHelper, CT18B1_EarthHelper) whose LuaTranslation places it at P = L + C. Then R(P − N − C) = R(L − N), exactly the true direction. C is read each frame from openspace.navigation.getNavigationState("Root") (anchor world position + camera offset). If that call fails, C = position of CT18B1_Site, which is correct to better than 0.001° for any camera near the site. Source: lua/ct18b1_light_helper_sun.lua. Check A17 executes both helper scripts in Lua 5.4 with a mocked API and confirms P = L + C exactly.

Residual risk (cannot be tested here): the helper uses the camera position of the previous frame, which is irrelevant at rover speeds but lags during fast camera flights. If OpenSpace's Lua table layout for vectors differs from both forms handled ({x=..} and {1,2,3}), the fallback path is used. Section 8, step 4 tells you how to see which.

5.4 Camera-mounted lamp

In OpenSpace the camera sits at the origin of view space, with +y up and −z the viewing direction. The shader places the lamp at view-space (0, 0.1016, 0) and aims it along −z, then transforms both into the model frame with inverse(modelViewTransform) (ct18b1Frame() in ct18b1_lighting.glsl). The lamp therefore stays "affixed to the top of the camera" however the camera moves or turns. C04 checks the resulting phase-angle geometry in the image: every band edge in T05 is bracketed by adjacent pixels whose GPU-derived phase angles straddle the limit (for example 1.004° → 0.997° at 5.8 m, predicted 5.82 m).

5.5 Textures, precision and depth (Ghoul details that matter)

Image rows are flipped on load (ghoul/src/io/texture/texturereaderstb.cpp). The shader undoes this (ct18b1Texel(): a.y = CT18B1_ATLAS_H - 1 - a.y). A16 runs ct18b1TerrainZ() on the GPU through a flipped upload of the real atlas and matches the Python decode to 0.19 mm.
Heights are 16-bit, split over R (high byte) and G (low byte) of an 8-bit RGB PNG, and read with texelFetch and manual bilinear interpolation, so OpenSpace's automatic mip-mapping and anisotropic filtering cannot corrupt them. The worst normal tilt the quantisation can cause is 0.26° (A15).
OpenSpace writes gl_FragDepth = distance/10³⁰ (floatoperations.glsl), which needs a 32-bit float depth buffer. The emulator allocates one; a 24-bit buffer would silently break occlusion.

5.6 The reference renderer ("emulator") and what it does not reproduce

tools/os_harness.py compiles the shipped shader files through unmodified copies of OpenSpace's render.frag, renderframebuffer.frag, fragment.glsl, floatoperations.glsl and PowerScaling headers (in tools/openspace_ref/, MIT licence). It sets the same uniforms RenderableModel sets, uses the default 1280×720 window with 50.534° vertical FOV, and applies the tone map of hdrAndFiltering.frag (1 − 2^(−3.7c), then c^(1/0.95)). Not reproduced: stars, the Sun and Earth disks and OpenSpace's own Earth, FXAA, the GUI, and OpenSpace's frame timing. Your screenshots will therefore contain extra sky objects and may differ by anti-aliasing, while the terrain itself should match.

6. The surface ("large, roughly spherical, bumpy")

6.1 Two objects

CT18B1_Globe: a sphere of radius 1737.4 km (IAU mean lunar radius) with synthetic relief of 1.5 km RMS (a sum of 600 random waves on the sphere with a k^−1.76 amplitude spectrum [A]), 163,842 vertices; vertex normals are generated by OpenSpace's model loader (Assimp smooth normals), which the emulator reproduces. It is a child of OpenSpace's stock Moon node, so it has the Moon's real position and IAU_MOON orientation. The stock Moon's renderable is switched off while the asset is loaded and switched back on when it is unloaded.
CT18B1_Site: a 6.4 km × 6.4 km patch at the Apollo 11 coordinates, lying on the same sphere (exact curvature drop r²/(R + √(R² − r²))). Its 40 m skirt and the globe's 30 m dimple under the patch stop gaps showing from orbit. From the rover the geometric horizon is at √(2Rh) = 2.3 km for h = 1.524 m, so the patch covers everything visible (C11: the GPU render and a Python line-of-sight ray cast agree on 3199 m vs 3208 m).

6.2 Where "real data" could and could not be used

No real lunar DEM was reachable from the build environment (section 2). The surface is therefore a statistical mare surface whose every ingredient comes from measurements:

ingredient	law	parameters	tags	checks
relief	Gaussian fractal, isotropic, power spectrum ∝ k^−(2H+2)	H = 0.76; amplitude: median slope component 3.6° over 1 m	V:ROS11, V:ROW71	A14 (spectrum reproduced band by band), A14b (3.55° at 1 m), A14c (0.74° at 0.75 km, maria < 1°)
craters	equilibrium N(>D) = 0.01·D⁻²	0.2 m – 1 km, depth/diameter log-uniform 0.02–0.20, rim height 0.18 × depth, rim decays as (r/R)⁻³ to 3R	V:SP8023, L:PIKE77, [A] degradation mix	A11 (counts within Poisson expectation)
rocks	Golombek–Rapp area law F(D) = k·e^(−qD), height 0.2347D + 0.0039 m	q = 1.743 m⁻¹, k = 0.025 (Surveyor I density), 0.1–8 m, lumpy elliptical caps [A]	V:DI16, D:ROCKK	A12, A13

The same three statistics agree with each other: ROW71's two slope values imply H ≈ 0.73, close to LOLA's 0.76. That agreement was not engineered.

Terrain seed: the random seed was chosen from 1800–1899 so that the rover path has an open view West, the way a real traverse would be planned to crest a rise. Every seed obeys the same statistics; with the first seed tried, the rover sat 12 m below a ridge 300 m ahead. This is a site-selection choice, not a parameter change (RANDOM_SEED in the config, with this note).

6.3 Resolution pyramid and why shading does not depend on the mesh

Heights live in three levels: near 2.5 cm (128 m × 25.6 m along the rover lane), mid 12.5 cm (256 m square) and far 4 m (6.4 km square), packed into one 5120 × 3072 atlas. Each finer level equals the coarser one plus extra detail that fades to zero at its border, so level blending is seamless. The shader picks the finest level a pixel can resolve (full weight at ≤ 1 texel per pixel, zero at ≥ 2), which is proper prefiltering. Unresolved roughness is exactly what Hapke's θ̄ describes.

The triangle mesh (580,301 vertices, spacing 5 cm at the lane centre growing with distance) is only a bounding proxy. The vertex shader displaces it from the same heightfield, and the fragment shader then intersects the view ray with the heightfield itself (ct18b1RefineHit()) before shading. This was added during this build after test T18 showed a black crescent in front of a rock. Where the mesh is coarser than a rock, the pixel was being shaded on the rock's back slope. The fix was verified visually (T18) and by all radiometric checks.

6.4 Shadows

Cast shadows are computed per pixel by marching from the surface point toward the Sun and the Earth over the heightfield: geometric steps from 4 cm out to 3.2 km, bilinear heights, sphere curvature included. The result is the horizon angle in the source's vertical plane, from which the visible fraction of the solar or terrestrial disk gives the penumbra. Lamp shadows march along the straight segment to the lamp. The first version sampled heights with nearest-texel lookups, which produced 2.5 cm "staircase" teeth on shadow edges seen close up. The march now uses bilinear sampling (fixed in this build, visible difference in T01's foreground).

6.5 Plugging in a real DEM

python3 tools/generate_terrain.py --dem your_far_level.npy replaces the statistical far level (1600 × 1600 float metres at 4 m spacing, centred on the site, x = East, y = North) with your data. Craters are then not added at that level; mid/near detail is still synthesised on top. An LROC NAC DTM of the Apollo 11 site (2 m/pixel), resampled to 4 m, would make the whole 6.4 km patch real at 4 m resolution.

7. Test cases (one keystroke each in OpenSpace)

Load the profile CT18.B1.Openspace, then press a key. The dashboard (top left) shows the test id, the UTC time and the settings. Reference images are in verification/images/<id>_f0.png, rendered by the emulator at the same time and camera. Press F12 to save an OpenSpace screenshot for comparison. Shift+0 lists all keys.

key	id	what it shows	what it verifies	what you should see	automated check
`1`	T01	Rover traverse - lunar morning, Sun low behind the rover (UTC 2026-10-17T12:00:16)	Main requested scene: camera 1.524 m (5 ft) above the ground, lamp 0.1016 m (4 in) above the camera, rover driving West at 1 m/s. Sun 12 deg high in the East, behind the camera, so the anti-solar point lies on the ground ahead. Checks: rocks and craters cast fully black shadows (no ambient light); the ground brightens toward the anti-solar point (Sun opposition surge, retroreflection of sunlight); rocks at Surveyor-I density; horizon at the distance set by lunar curvature (~2.2 km).	Black shadows; bright patch on the ground ~12 deg below the horizon straight ahead; rover advances 30 m every 30 s.	C02, C03, C10, C11
`2`	T02	Rover traverse - lunar night under a full Earth, headlamp on (UTC 2026-10-10T15:30:16)	The retroreflective headlamp: the lamp sits 10 cm above the eye, so every lit ground point is seen at a phase angle of only ~0.1 m / distance. The ground stays bright much farther out than a 1/r^2 fall-off suggests; the beam is a Gaussian 40 deg FWHM so brightness also falls off to the sides. Earthshine (full Earth, ~1e-4 of sunlight) is too faint to see at this exposure.	Wide bright pool that fades slowly with distance and looks almost flat (relief hidden at zero phase).	C05
`3`	T03	As T02 but opposition effect OFF (A/B control) (UTC 2026-10-10T15:30:16)	Same frame as T02 with B_S0 = B_C0 = 0. Distant ground (small lamp phase angle) must be much darker than in T02 while ground right in front (larger phase angle) changes less. The ratio T02/T03 per pixel is the retroreflection factor and is printed by the verification suite.	Same frame as 2 but far ground clearly darker.	C05
`4`	T04	As T02 but Lambert shading ('computer graphics playbook') (UTC 2026-10-10T15:30:16)	What a conventional renderer would show: Lambert (cosine) reflection with the same normal albedo. No retroreflection: the scene falls off as 1/r^2 and surfaces facing away from the lamp darken. Compare with T02.	Relief strongly shaded, far ground much darker than in 2.	C05
`5`	T05	Headlamp phase-angle map (false colour) (UTC 2026-10-10T15:30:00)	Geometry check. Colour bands show the lamp phase angle g at each ground point: white <0.25 deg, red <0.5, orange <1, yellow <2, green <4, blue <8, purple >=8. For a lamp 0.1016 m above the eye, g ~ 0.1016/d rad, so the 1 deg boundary must lie ~5.8 m away, 0.5 deg ~11.6 m, 0.25 deg ~23 m (checked numerically in the suite).	Colour bands; yellow/orange edge ~5.8 m, orange/red ~11.6 m, red/white ~23 m.	C04
`6`	T06	Headlamp only, physical shading (same view as T05) (UTC 2026-10-10T15:30:00)	Pairs with T05: brightness along the centre line vs distance is compared with the analytic prediction (lamp inverse-square x Gaussian beam x Hapke with the opposition surge) in the verification report.	Smooth lamp pool; compare with 5.	C01
`7`	T07	Sun phase-angle map around the anti-solar point (false colour) (UTC 2026-10-17T12:00:16)	Camera looks straight at the anti-solar point (azimuth = Sun azimuth + 180, depression = Sun elevation). The white/red bands must be centred on the image centre, i.e. where the shadow of the observer's head would fall.	Concentric colour rings centred in the image.	visual
`8`	T08	Heiligenschein: physical image of the T07 view (UTC 2026-10-17T12:00:16)	The retroreflective glow around the anti-solar point that Apollo astronauts photographed around their own shadows. Brightness must peak at the image centre and fall ~30 % within ~4 deg (shape set by h_S, h_C).	Brightest point at the image centre, fading outward.	C06
`9`	T09	Sun visibility (shadow / penumbra) map (UTC 2026-10-17T12:00:16)	White = Sun fully visible, dark blue = fully hidden. Penumbrae are narrow because the solar disk is only 0.53 deg across; their width must equal (distance from the shadow-casting edge) x 0.0093 rad.	White ground, dark-blue shadow shapes with thin soft edges.	used by C02, C03
`0`	T10	Rover view at lunar noon (Sun 88 deg high) (UTC 2026-10-23T22:00:16)	Near-vertical Sun: shadows shrink to small black patches under rocks; the scene looks flat and washed out (well documented by Apollo crews).	Flat-looking, low-contrast ground; small shadows under rocks.	visual
`Shift+1`	T11	Evening - looking West into the Sun (forward scatter) (UTC 2026-10-30T07:49:10)	Large phase angles: the regolith is dark toward the Sun (no forward-scattering glare, b > 0 back-scattering particles); rocks show bright rims only where their faces tilt toward the Sun.	Dark ground toward the Sun; rock rims lit.	visual
`Shift+2`	T12	Whole body at full phase (Hapke) (UTC 2026-10-17T12:00:16)	The famous flatness of the full Moon: a Hapke/Lommel-Seeliger surface is almost as bright at the limb as at the centre (no limb darkening).	Nearly uniformly bright disk with a slightly brighter centre.	C07, C07c
`Shift+3`	T13	Whole body at full phase (Lambert control) (UTC 2026-10-17T12:00:16)	Same view with Lambert shading: strong limb darkening, which the real full Moon does not show. Radial profiles of T12 and T13 are compared in the report.	Disk clearly darker toward the edge than 12.	C07c
`Shift+4`	T14	Whole body at quarter phase (UTC 2026-10-17T12:00:16)	Terminator and relief of the (synthetic) global shape; brightness rises toward the bright limb, as on the real quarter Moon.	Half-lit disk; terminator; bright limb side.	visual
`Shift+5`	T15	Heightfield level diagnostic (rover view) (UTC 2026-10-17T12:00:16)	Red = 2.5 cm near level, green = 12.5 cm mid level, blue = 4 m far level. Transitions must be smooth and happen where one pixel covers 1-2 texels (no aliasing, no popping).	Red near, green middle, blue far; smooth transitions.	visual
`Shift+6`	T16	Rock mask over the physical image (UTC 2026-10-17T12:00:16)	Rocks tinted red, so their density and sizes can be judged against the Surveyor I statistic (0.34 rocks >= 10 cm per m^2).	Red-tinted rocks scattered densely near the rover.	A12 (visual)
`Shift+7`	T17	Lunar night, headlamp OFF: earthshine only (UTC 2026-10-10T15:30:16)	Full Earth 71 deg high lights the surface at ~1e-4 of sunlight (Earth V-band geometric albedo 0.226). Needs 10^4 more exposure than daylight; shadows point away from the Earth.	Faint grey ground, shadows pointing away from the Earth.	C08
`Shift+8`	T18	Close-up of a rock at night, headlamp on (UTC 2026-10-10T15:30:00)	Lamp shadows are almost invisible because the lamp is only 10 cm above the eye (the same geometry that hides the observer's shadow in sunlight); only a thin dark fringe shows above the rock's upper edge.	Rock and ground lit; no black crescent; at most a thin dark fringe above the rock.	visual (rock lit; thin lamp-shadow fringe only)
`Shift+9`	T19	As T01 frame 0 with macroscopic roughness OFF (UTC 2026-10-17T12:00:16)	Sensitivity control: removing Hapke's sub-pixel roughness (theta = 8.1 deg) brightens grazing-incidence ground and weakens the anti-solar brightening less than removing the opposition effect does.	Like 1 but grazing ground brighter.	visual

At the night epoch the Earth is at azimuth 284 deg, elevation 71 deg, 3.8 deg from full; the Sun is 72 deg below the horizon.

Additional automated checks (no key; run python3 tools/run_verification.py): A01–A18 as listed in verification/VERIFICATION_REPORT.md, covering parameter bookkeeping, GPU-vs-Python equality, literature anchors, ephemeris frames, terrain statistics, Lua scripts and texture layout.

8. What only you can verify (OpenSpace on your machine)

Copy openspace_user/ into D:\OpenSpace-21.4\user\ (see README.md). Start OpenSpace and choose the profile CT18.B1.Openspace.
Check the log (D:\OpenSpace-21.4\logs\ or on screen) for the line "CT18.B1.Openspace loaded" and for any shader compile error mentioning ct18b1_. The emulator compiled the shaders with Mesa's GLSL compiler. AMD, NVIDIA and Intel compilers are stricter or looser in different places, and a message there would be the first thing to send back.
Press 1 and compare with verification/images/T01_f0.png: same horizon, same rocks, same shadow shapes, a glow on the ground at the lower centre. Stars, the Earth and the Sun disk will be extra in OpenSpace.
Press Shift+7 (earthshine only) and look at the shadows: they must fall away from the Earth (the Earth is 71° high at azimuth 284°, i.e. West-north-west). If they fall away from where the Sun would be (below the horizon in the East), the light-helper fix in 5.3 is not working on your installation.
Press 2, then 3, then 4 (headlamp with and without the opposition effect, and Lambert): the far ground must be clearly brighter in 2 than in 3.

9. Limitations, assumptions and known deviations (complete list)

Assumptions [A]: lamp 20 000 cd, 40° Gaussian beam; eye height 1.524 m and lamp offset 0.1016 m (user's "few feet"/"few inches"); rover speed 1 m/s; 2.3 m chassis smoothing; CBOE half width 2.0° (inside VEL16's 1–4°); crater degradation mix; rock plan shapes; globe relief spectrum and 1.5 km RMS; exposure gains; ray-march step counts.

Physics not modelled:

Terrain inter-reflection. Light scattered from sunlit slopes into shadows is ignored, so shadows next to sunlit rock faces are slightly too dark. On a dark (normal albedo 0.095), mostly flat mare this is a small effect, but it is not zero.
Rock photometry equals regolith photometry (no measured difference was found). Real rocks lack the fine regolith's porosity, so their opposition surge is probably weaker.
Colour. Grey V-band only: no lunar reddening and no wavelength dependence of the surge.
Distances. Fixed at 1 AU (Sun, ±3.4 %) and the mean Earth distance (±~11 % in earthshine).
Earth shadows are traced only where the Sun is not fully visible. In sunlight earthshine is < 2e-4 of the light, so this cannot change a pixel.
Eclipses of the site by the Earth are not modelled.
Rover motion. The rover drives a straight line and the camera stays level (it does not pitch or roll with the chassis).
No distant mountains. Nothing lies beyond the 3.2 km patch.

Engineering limits: the reference images are from an emulator, not OpenSpace (5.6); the SceneGraphLightSource work-around depends on OpenSpace's Lua API behaving as read from the source (5.3); the patch is static (the rover cannot leave the 128 m high-resolution lane without the near level ending; mid and far levels continue out to 3.2 km).

Previous version: CT18.A1.Openspace was not available to this build (it was on your computer and was not uploaded). B1 was therefore built from scratch, and nothing from A1 was reused or compared.

10. Reproducing everything

# Linux / WSL, Python 3.10+
pip install numpy scipy pillow matplotlib moderngl PyOpenGL lupa jplephem
pip download --no-deps --no-binary :all: de421 && tar xzf de421-2008.1.tar.gz -C tools   # ephemeris data
python3 tools/hapke_reference.py --solve   # re-derive h_C, h_S (values are frozen in the config)
python3 tools/generate_terrain.py          # atlas, meshes, parameter header, rover paths (~1 min)
python3 tools/make_assets.py               # OpenSpace assets, test keys, profile
python3 tools/run_verification.py          # 33 checks + reference images (~1-2 min, software OpenGL)
python3 tools/make_audit.py                # refresh the appendices of this file

Appendix A — Shader source (as shipped)

`shaders/ct18b1_hapke.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - Hapke (2012) reflectance model, GLSL implementation.
   3   *
   4   * Pure functions only (no uniforms), so the verification suite can compile this file
   5   * unchanged into a probe shader and compare it against tools/hapke_reference.py.
   6   *
   7   * Equation numbers [E1]..[E6] refer to AUDIT.md section 3 and to the docstring of
   8   * tools/hapke_reference.py.
   9   ****************************************************************************************/
  10  #ifndef CT18B1_HAPKE_GLSL
  11  #define CT18B1_HAPKE_GLSL
  12  
  13  struct HapkeParams {
  14    float w;       // single-scattering albedo
  15    float b;       // Legendre P1 coefficient
  16    float c;       // Legendre P2 coefficient
  17    float bs0;     // shadow-hiding opposition amplitude  (SHOE)
  18    float hs;      // shadow-hiding opposition width
  19    float bc0;     // coherent-backscatter opposition amplitude (CBOE)
  20    float hc;      // coherent-backscatter opposition width
  21    float theta;   // macroscopic roughness, radians
  22    float K;       // porosity factor
  23  };
  24  
  25  // [E5] Ambartsumian-Chandrasekhar H function, Hapke (2002) approximation
  26  float hapkeH(float x, float w) {
  27    x = max(x, 1e-6);
  28    float gamma = sqrt(1.0 - w);
  29    float r0 = (1.0 - gamma) / (1.0 + gamma);
  30    return 1.0 / (1.0 - w * x * (r0 + 0.5 * (1.0 - 2.0 * r0 * x) * log((1.0 + x) / x)));
  31  }
  32  
  33  // [E2] single-particle phase function, Legendre form used by Helfenstein & Veverka (1987)
  34  float hapkeP(float g, float b, float c) {
  35    float cg = cos(g);
  36    return 1.0 + b * cg + c * (1.5 * cg * cg - 0.5);
  37  }
  38  
  39  // [E3] shadow-hiding opposition effect
  40  float hapkeBS(float g, float B0, float h) {
  41    if (B0 == 0.0) return 0.0;
  42    return B0 / (1.0 + tan(0.5 * g) / h);
  43  }
  44  
  45  // [E4] coherent-backscatter opposition effect
  46  float hapkeBC(float g, float B0, float h) {
  47    if (B0 == 0.0) return 0.0;
  48    float y = tan(0.5 * g) / h;
  49    float term = (y < 1e-4) ? (1.0 - 0.5 * y) : (1.0 - exp(-y)) / y;
  50    return B0 * (1.0 + term) / (2.0 * (1.0 + y) * (1.0 + y));
  51  }
  52  
  53  // [E6] Hapke (1984) macroscopic roughness
  54  float hapkeE1(float x, float tt) { return exp(-2.0 / CT18B1_PI / tt / tan(x)); }
  55  float hapkeE2(float x, float tt) { float t = tan(x); return exp(-1.0 / CT18B1_PI / (tt * tt) / (t * t)); }
  56  
  57  void hapkeRoughness(float i, float e, float psi, float theta,
  58                      out float mu0e, out float mue, out float S)
  59  {
  60    float mu0 = cos(i);
  61    float mu = cos(e);
  62    if (theta <= 0.0) {
  63      mu0e = mu0; mue = mu; S = 1.0;
  64      return;
  65    }
  66    float tt = tan(theta);
  67    float chi = 1.0 / sqrt(1.0 + CT18B1_PI * tt * tt);
  68    float ic = max(i, 1e-5);
  69    float ec = max(e, 1e-5);
  70    float E1i = hapkeE1(ic, tt), E2i = hapkeE2(ic, tt);
  71    float E1e = hapkeE1(ec, tt), E2e = hapkeE2(ec, tt);
  72    float si = sin(ic), se = sin(ec);
  73    float etaI = chi * (mu0 + si * tt * E2i / (2.0 - E1i));
  74    float etaE = chi * (mu + se * tt * E2e / (2.0 - E1e));
  75    float f = exp(-2.0 * tan(min(psi, CT18B1_PI - 1e-4) * 0.5));
  76    float s2 = sin(0.5 * psi); s2 *= s2;
  77    float cpsi = cos(psi);
  78    float fpi = psi / CT18B1_PI;
  79    if (i <= e) {
  80      float d = 2.0 - E1e - fpi * E1i;
  81      mu0e = chi * (mu0 + si * tt * (cpsi * E2e + s2 * E2i) / d);
  82      mue  = chi * (mu  + se * tt * (E2e - s2 * E2i) / d);
  83      S = (mue / etaE) * (mu0 / etaI) * chi / (1.0 - f + f * chi * (mu0 / etaI));
  84    }
  85    else {
  86      float d = 2.0 - E1i - fpi * E1e;
  87      mu0e = chi * (mu0 + si * tt * (E2i - s2 * E2e) / d);
  88      mue  = chi * (mu  + se * tt * (cpsi * E2i + s2 * E2e) / d);
  89      S = (mue / etaE) * (mu0 / etaI) * chi / (1.0 - f + f * chi * (mu / etaE));
  90    }
  91  }
  92  
  93  // [E1] radiance factor RADF = I/F = pi * r
  94  float hapkeRadf(float i, float e, float g, float psi, HapkeParams p) {
  95    if (cos(i) <= 0.0 || cos(e) <= 0.0) return 0.0;
  96    float mu0e, mue, S;
  97    hapkeRoughness(i, e, psi, p.theta, mu0e, mue, S);
  98    if (mu0e <= 0.0 || mue <= 0.0) return 0.0;
  99    float ls = mu0e / (mu0e + mue);
 100    float single = hapkeP(g, p.b, p.c) * (1.0 + hapkeBS(g, p.bs0, p.hs));
 101    float multi = hapkeH(mu0e / p.K, p.w) * hapkeH(mue / p.K, p.w) - 1.0;
 102    float r = p.K * p.w / (4.0 * CT18B1_PI) * ls * (single + multi) *
 103              (1.0 + hapkeBC(g, p.bc0, p.hc)) * S;
 104    return CT18B1_PI * max(r, 0.0);
 105  }
 106  
 107  // Same, from unit vectors: n = surface normal, l = toward light, v = toward viewer.
 108  // Also returns the phase angle g (radians) for diagnostics.
 109  float hapkeRadfVec(vec3 n, vec3 l, vec3 v, HapkeParams p, out float g) {
 110    float mu0 = dot(n, l);
 111    float mu = dot(n, v);
 112    g = acos(clamp(dot(l, v), -1.0, 1.0));
 113    if (mu0 <= 0.0 || mu <= 0.0) return 0.0;
 114    float i = acos(min(mu0, 1.0));
 115    float e = acos(min(mu, 1.0));
 116    vec3 lp = l - n * mu0;
 117    vec3 vp = v - n * mu;
 118    float nl = length(lp), nv = length(vp);
 119    float cpsi = (nl > 1e-6 && nv > 1e-6) ? dot(lp, vp) / (nl * nv) : 1.0;
 120    float psi = acos(clamp(cpsi, -1.0, 1.0));
 121    return hapkeRadf(i, e, g, psi, p);
 122  }
 123  
 124  #endif // CT18B1_HAPKE_GLSL

`shaders/ct18b1_lighting.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - light sources, operating modes and exposure.
   3   * Shared by the terrain and globe fragment shaders.
   4   *
   5   * RenderableModel gives a custom shader only a fixed set of uniforms, so three of
   6   * them are re-purposed (documented in AUDIT.md section 5.2 and README.md):
   7   *   ambientIntensity   -> camera exposure:  gain = 10^(CT18B1_EXPOSURE_DECADES * a)
   8   *   diffuseIntensity   -> rover lamp dimmer (0 = off, 1 = rated CT18B1_LAMP_CD)
   9   *   specularIntensity  -> mode = round(10 * s):
  10   *        10 physical (default)       0 opposition effect OFF (B_S0 = B_C0 = 0)
  11   *         1 shadow hiding only        2 coherent backscatter only
  12   *         3 "computer-graphics playbook" Lambert, same albedo at normal geometry
  13   *         4 false colour: lamp phase angle   5 false colour: Sun phase angle
  14   *         6 Sun visibility (shadow/penumbra) 7 heightfield level used
  15   *         8 rock mask over physical image    9 macroscopic roughness OFF
  16   *   lightIntensities[0] -> Sun on/off scale,  lightIntensities[1] -> Earth on/off scale
  17   * Light directions come from SceneGraphLightSource nodes CT18B1_SunHelper and
  18   * CT18B1_EarthHelper (see lua/ct18b1_light_helper_*.lua for why helpers are needed).
  19   *
  20   * Output: frag.color = gain * sum_k (E_k / E_sun) * vis_k * RADF_k   (grey, V band)
  21   * i.e. the radiance factor relative to 1 AU sunlight, times the camera gain, handed
  22   * to OpenSpace's HDR pipeline with disableLDR2HDR = true.
  23   ****************************************************************************************/
  24  #ifndef CT18B1_LIGHTING_GLSL
  25  #define CT18B1_LIGHTING_GLSL
  26  
  27  uniform mat4 modelViewTransform;
  28  uniform float ambientIntensity;
  29  uniform float diffuseIntensity;
  30  uniform float specularIntensity;
  31  uniform int nLightSources;
  32  uniform vec3 lightDirectionsViewSpace[8];
  33  uniform float lightIntensities[8];
  34  uniform float opacity;
  35  
  36  struct Ct18b1Frame {
  37    mat4 viewToModel;
  38    vec3 camera;      // model frame
  39    vec3 lampPos;     // model frame
  40    vec3 lampAxis;    // model frame, unit
  41    vec3 sunDir;      // model frame, unit, toward Sun
  42    vec3 earthDir;    // model frame, unit, toward Earth
  43    float sunScale;
  44    float earthScale;
  45    float lampScale;
  46    float gain;
  47    int mode;
  48  };
  49  
  50  Ct18b1Frame ct18b1Frame() {
  51    Ct18b1Frame f;
  52    f.viewToModel = inverse(modelViewTransform);
  53    mat3 v2m = mat3(f.viewToModel);
  54    f.camera = (f.viewToModel * vec4(0.0, 0.0, 0.0, 1.0)).xyz;
  55    f.lampPos = (f.viewToModel * vec4(0.0, CT18B1_LAMP_OFFSET, 0.0, 1.0)).xyz;
  56    f.lampAxis = normalize(v2m * vec3(0.0, 0.0, -1.0));
  57    f.sunDir = (nLightSources > 0) ? normalize(v2m * lightDirectionsViewSpace[0]) : vec3(0.0, 0.0, 1.0);
  58    f.earthDir = (nLightSources > 1) ? normalize(v2m * lightDirectionsViewSpace[1]) : vec3(0.0, 0.0, 1.0);
  59    f.sunScale = (nLightSources > 0) ? lightIntensities[0] : 0.0;
  60    f.earthScale = (nLightSources > 1) ? lightIntensities[1] : 0.0;
  61    f.lampScale = clamp(diffuseIntensity, 0.0, 1.0);
  62    f.gain = pow(10.0, CT18B1_EXPOSURE_DECADES * clamp(ambientIntensity, 0.0, 1.0));
  63    f.mode = int(floor(specularIntensity * 10.0 + 0.5));
  64    return f;
  65  }
  66  
  67  HapkeParams ct18b1Params(int mode) {
  68    HapkeParams p;
  69    p.w = CT18B1_W; p.b = CT18B1_B; p.c = CT18B1_C;
  70    p.bs0 = CT18B1_BS0; p.hs = CT18B1_HS;
  71    p.bc0 = CT18B1_BC0; p.hc = CT18B1_HC;
  72    p.theta = CT18B1_THETA_BAR; p.K = CT18B1_K;
  73    if (mode == 0) { p.bs0 = 0.0; p.bc0 = 0.0; }
  74    if (mode == 1) { p.bc0 = 0.0; }
  75    if (mode == 2) { p.bs0 = 0.0; }
  76    if (mode == 9) { p.theta = 0.0; }
  77    return p;
  78  }
  79  
  80  // Irradiance of each source relative to the Sun at 1 AU (dimensionless)
  81  float ct18b1EarthIrradiance(vec3 sunDir, vec3 earthDir) {
  82    // Earth's phase angle seen from the Moon: angle Sun-Earth-Moon.
  83    float a = acos(clamp(-dot(sunDir, earthDir), -1.0, 1.0));
  84    float lambertPhase = (sin(a) + (CT18B1_PI - a) * cos(a)) / CT18B1_PI;
  85    return CT18B1_EARTH_ALBEDO * CT18B1_EARTH_SOLID * lambertPhase;
  86  }
  87  
  88  float ct18b1LampIrradiance(Ct18b1Frame f, vec3 P, out vec3 l, out float dist) {
  89    vec3 d = f.lampPos - P;
  90    dist = length(d);
  91    l = d / dist;
  92    float th = acos(clamp(dot(-l, f.lampAxis), -1.0, 1.0));
  93    float beam = exp(-4.0 * log(2.0) * th * th / (CT18B1_LAMP_FWHM * CT18B1_LAMP_FWHM));
  94    return f.lampScale * CT18B1_LAMP_CD * beam / (dist * dist * CT18B1_E_SUN_LUX);
  95  }
  96  
  97  // Reflected radiance factor from one source (Hapke, or Lambert in mode 3)
  98  float ct18b1Reflect(vec3 N, vec3 l, vec3 v, HapkeParams p, int mode, out float g) {
  99    if (mode == 3) {
 100      g = acos(clamp(dot(l, v), -1.0, 1.0));
 101      return CT18B1_NORMAL_ALBEDO_LAMBERT * max(dot(N, l), 0.0) * (dot(N, v) > 0.0 ? 1.0 : 0.0);
 102    }
 103    return hapkeRadfVec(N, l, v, p, g);
 104  }
 105  
 106  // Inverse of OpenSpace's tone map (hdrAndFiltering.frag) at default settings, so
 107  // diagnostic false colours display exactly as authored.
 108  vec3 ct18b1DisplayToLinear(vec3 c) {
 109    c = pow(clamp(c, 0.0, 0.999), vec3(CT18B1_OS_GAMMA));
 110    return -log2(1.0 - c) / CT18B1_OS_EXPOSURE;
 111  }
 112  
 113  vec3 ct18b1AngleColour(float gdeg) {
 114    // bands: <0.25, <0.5, <1, <2, <4, <8, >=8 degrees
 115    if (gdeg < 0.25) return vec3(1.0, 1.0, 1.0);
 116    if (gdeg < 0.5)  return vec3(1.0, 0.1, 0.1);
 117    if (gdeg < 1.0)  return vec3(1.0, 0.6, 0.0);
 118    if (gdeg < 2.0)  return vec3(1.0, 1.0, 0.0);
 119    if (gdeg < 4.0)  return vec3(0.0, 0.9, 0.0);
 120    if (gdeg < 8.0)  return vec3(0.0, 0.6, 1.0);
 121    return vec3(0.25, 0.0, 0.5);
 122  }
 123  
 124  #endif // CT18B1_LIGHTING_GLSL

`shaders/ct18b1_terrain_fs.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - terrain fragment shader (RenderableModel custom FragmentShader).
   3   *
   4   * Retroreflective lunar regolith:  Hapke (2012) with shadow-hiding (SHOE) and coherent
   5   * backscatter (CBOE) opposition effects, macroscopic roughness, resolved-terrain cast
   6   * shadows (finite solar and terrestrial disks), inverse-square camera-mounted lamp.
   7   * No ambient term: an unlit lunar surface receives only earthshine.
   8   ****************************************************************************************/
   9  #include "fragment.glsl"
  10  #include "ct18b1_params.glsl"
  11  #include "ct18b1_hapke.glsl"
  12  #include "ct18b1_terrain_common.glsl"
  13  #include "ct18b1_lighting.glsl"
  14  
  15  in vec3 vs_positionModel;
  16  in float vs_spacing;
  17  in vec4 vs_positionCameraSpace;
  18  in float vs_screenSpaceDepth;
  19  
  20  Fragment getFragment() {
  21    Fragment frag;
  22    frag.depth = vs_screenSpaceDepth;
  23    frag.gPosition = vs_positionCameraSpace;
  24    frag.disableLDR2HDR = true;
  25    frag.color.a = opacity;
  26  
  27    Ct18b1Frame F = ct18b1Frame();
  28    HapkeParams hp = ct18b1Params(F.mode);
  29  
  30    // pixel footprint on the ground (metres), used to pick the heightfield level
  31    vec2 xy = vs_positionModel.xy;
  32    float fp = max(length(dFdx(xy)), length(dFdy(xy)));
  33    fp = max(fp, 1e-4);
  34  
  35    vec3 P = vec3(xy, ct18b1TerrainZ(xy, fp));
  36    if (vs_spacing > 2.0 * fp) {
  37      // mesh coarser than the pixel: intersect the view ray with the real heightfield
  38      P = ct18b1RefineHit(F.camera, vs_positionModel, vs_spacing, fp);
  39      xy = P.xy;
  40      frag.gPosition = modelViewTransform * vec4(P, 1.0);
  41    }
  42    vec3 N = ct18b1TerrainNormal(xy, fp);
  43    vec3 V = normalize(F.camera - P);
  44    frag.gNormal = vec4(normalize(mat3(modelViewTransform) * N), 0.0);
  45  
  46    float total = 0.0;
  47    float gSun = 0.0, gLamp = 0.0, gEarth = 0.0;
  48    float visSun = 0.0;
  49  
  50    // ---- Sun (finite disk, penumbra)
  51    if (F.sunScale > 0.0 && F.sunDir.z > -0.3) {
  52      visSun = ct18b1SourceVisibility(P, F.sunDir, CT18B1_SUN_ANG_RADIUS, fp);
  53      if (visSun > 0.0) {
  54        total += F.sunScale * visSun * ct18b1Reflect(N, F.sunDir, V, hp, F.mode, gSun);
  55      }
  56      else {
  57        gSun = acos(clamp(dot(F.sunDir, V), -1.0, 1.0));
  58      }
  59    }
  60  
  61    // ---- Earth (earthshine; disk-integrated, Lambert-sphere phase law)
  62    if (F.earthScale > 0.0 && F.earthDir.z > -0.3) {
  63      float eIrr = ct18b1EarthIrradiance(F.sunDir, F.earthDir);
  64      // earthshine is < 2e-4 of sunlight; its shadows are only traced when the Sun is
  65      // down or when this pixel is in sun shadow (otherwise it cannot change the pixel)
  66      float visE = 1.0;
  67      if (visSun < 1.0) {
  68        visE = ct18b1SourceVisibility(P, F.earthDir, CT18B1_EARTH_ANG_RADIUS, fp);
  69      }
  70      total += F.earthScale * eIrr * visE * ct18b1Reflect(N, F.earthDir, V, hp, F.mode, gEarth);
  71    }
  72  
  73    // ---- rover lamp: point source LAMP_OFFSET above the eye, inverse-square
  74    if (F.lampScale > 0.0) {
  75      vec3 l;
  76      float dist;
  77      float eL = ct18b1LampIrradiance(F, P, l, dist);
  78      if (eL > 1e-12) {
  79        float visL = ct18b1PointVisibility(P, F.lampPos, fp);
  80        total += eL * visL * ct18b1Reflect(N, l, V, hp, F.mode, gLamp);
  81      }
  82      else {
  83        gLamp = acos(clamp(dot(l, V), -1.0, 1.0));
  84      }
  85    }
  86  
  87    vec3 colour = vec3(total * F.gain);
  88  
  89    // ---- diagnostic modes
  90    if (F.mode == 4) {
  91      vec3 l; float dist;
  92      ct18b1LampIrradiance(F, P, l, dist);
  93      colour = ct18b1DisplayToLinear(ct18b1AngleColour(degrees(acos(clamp(dot(l, V), -1.0, 1.0)))));
  94    }
  95    else if (F.mode == 5) {
  96      colour = ct18b1DisplayToLinear(ct18b1AngleColour(degrees(acos(clamp(dot(F.sunDir, V), -1.0, 1.0)))));
  97    }
  98    else if (F.mode == 6) {
  99      float vs = ct18b1SourceVisibility(P, F.sunDir, CT18B1_SUN_ANG_RADIUS, fp);
 100      colour = ct18b1DisplayToLinear(vec3(vs, vs, 0.25 + 0.75 * vs));
 101    }
 102    else if (F.mode == 7) {
 103      float w0 = ct18b1LevelWeight(0, xy, fp);
 104      float w1 = ct18b1LevelWeight(1, xy, fp);
 105      vec3 c = mix(mix(vec3(0.2, 0.2, 1.0), vec3(0.1, 0.9, 0.1), w1), vec3(1.0, 0.2, 0.2), w0);
 106      colour = ct18b1DisplayToLinear(c * (0.35 + 0.65 * max(dot(N, vec3(0.0, 0.0, 1.0)), 0.0)));
 107    }
 108    else if (F.mode == 8) {
 109      int L = ct18b1MarchLevel(xy, fp);
 110      float rock = ct18b1LevelRock(L, xy);
 111      colour = mix(colour, ct18b1DisplayToLinear(vec3(1.0, 0.0, 0.0)), 0.6 * rock);
 112    }
 113  
 114    frag.color.rgb = colour;
 115    return frag;
 116  }

`shaders/ct18b1_terrain_common.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - terrain heightfield access, shared by vertex and fragment shader.
   3   *
   4   * The heightfield is a 3-level pyramid packed into one RGB8 atlas (the model's diffuse
   5   * texture):  R = high byte, G = low byte of a 16-bit height, B = rock mask.
   6   * Level 0 near (2.5 cm), level 1 mid (12.5 cm), level 2 far (4 m).  Model frame is the
   7   * site frame: x = East, y = North, z = Up, origin on the sphere at the site.
   8   *
   9   * Python mirror: tools/generate_terrain.py (bilinear, composite_height, curvature_drop).
  10   * Ghoul flips images vertically on load, so atlas row r of the PNG is GL row H-1-r.
  11   ****************************************************************************************/
  12  #ifndef CT18B1_TERRAIN_COMMON_GLSL
  13  #define CT18B1_TERRAIN_COMMON_GLSL
  14  
  15  uniform sampler2D texture_diffuse;
  16  
  17  vec4 ct18b1LevelRect(int L) {
  18    return (L == 0) ? CT18B1_L0_RECT : ((L == 1) ? CT18B1_L1_RECT : CT18B1_L2_RECT);
  19  }
  20  float ct18b1LevelTexel(int L) {
  21    return (L == 0) ? CT18B1_L0_TEXEL : ((L == 1) ? CT18B1_L1_TEXEL : CT18B1_L2_TEXEL);
  22  }
  23  ivec2 ct18b1LevelDims(int L) {
  24    return (L == 0) ? CT18B1_L0_DIMS : ((L == 1) ? CT18B1_L1_DIMS : CT18B1_L2_DIMS);
  25  }
  26  ivec2 ct18b1LevelAtlas(int L) {
  27    return (L == 0) ? CT18B1_L0_ATLAS : ((L == 1) ? CT18B1_L1_ATLAS : CT18B1_L2_ATLAS);
  28  }
  29  float ct18b1LevelHOff(int L) {
  30    return (L == 0) ? CT18B1_L0_HOFF : ((L == 1) ? CT18B1_L1_HOFF : CT18B1_L2_HOFF);
  31  }
  32  float ct18b1LevelHScale(int L) {
  33    return (L == 0) ? CT18B1_L0_HSCALE : ((L == 1) ? CT18B1_L1_HSCALE : CT18B1_L2_HSCALE);
  34  }
  35  float ct18b1LevelBlend(int L) {
  36    return (L == 0) ? CT18B1_L0_BLEND : ((L == 1) ? CT18B1_L1_BLEND : CT18B1_L2_BLEND);
  37  }
  38  
  39  // raw texel (i along x, j along y) of level L
  40  vec3 ct18b1Texel(int L, ivec2 ij) {
  41    ivec2 d = ct18b1LevelDims(L);
  42    ij = clamp(ij, ivec2(0), d - 1);
  43    ivec2 a = ct18b1LevelAtlas(L) + ij;
  44    a.y = CT18B1_ATLAS_H - 1 - a.y;       // undo Ghoul's vertical flip
  45    return texelFetch(texture_diffuse, a, 0).rgb;
  46  }
  47  
  48  float ct18b1DecodeHeight(int L, vec3 t) {
  49    float q = floor(t.r * 255.0 + 0.5) * 256.0 + floor(t.g * 255.0 + 0.5);
  50    return ct18b1LevelHOff(L) + q * ct18b1LevelHScale(L);
  51  }
  52  
  53  // bilinear height of one level, texel-centre convention, clamp to edge
  54  float ct18b1LevelHeight(int L, vec2 xy) {
  55    vec4 r = ct18b1LevelRect(L);
  56    float tx = ct18b1LevelTexel(L);
  57    vec2 f = (xy - r.xy) / tx - 0.5;
  58    ivec2 i0 = ivec2(floor(f));
  59    vec2 t = f - vec2(i0);
  60    float h00 = ct18b1DecodeHeight(L, ct18b1Texel(L, i0));
  61    float h10 = ct18b1DecodeHeight(L, ct18b1Texel(L, i0 + ivec2(1, 0)));
  62    float h01 = ct18b1DecodeHeight(L, ct18b1Texel(L, i0 + ivec2(0, 1)));
  63    float h11 = ct18b1DecodeHeight(L, ct18b1Texel(L, i0 + ivec2(1, 1)));
  64    return mix(mix(h00, h10, t.x), mix(h01, h11, t.x), t.y);
  65  }
  66  
  67  // nearest-texel height (used by the shadow ray march)
  68  float ct18b1LevelHeightNearest(int L, vec2 xy) {
  69    vec4 r = ct18b1LevelRect(L);
  70    ivec2 ij = ivec2(floor((xy - r.xy) / ct18b1LevelTexel(L)));
  71    return ct18b1DecodeHeight(L, ct18b1Texel(L, ij));
  72  }
  73  
  74  float ct18b1LevelRock(int L, vec2 xy) {
  75    vec4 r = ct18b1LevelRect(L);
  76    ivec2 ij = ivec2(floor((xy - r.xy) / ct18b1LevelTexel(L)));
  77    return ct18b1Texel(L, ij).b;
  78  }
  79  
  80  float ct18b1BorderDist(vec4 r, vec2 xy) {
  81    return min(min(xy.x - r.x, r.z - xy.x), min(xy.y - r.y, r.w - xy.y));
  82  }
  83  
  84  // Weight of level L relative to the next coarser one at footprint fp (metres/pixel):
  85  // full weight while a pixel covers <= 1 texel, zero at >= 2 texels (prefiltering, no
  86  // aliasing), and a smooth fade over BLEND metres inside the level's border.
  87  float ct18b1LevelWeight(int L, vec2 xy, float fp) {
  88    float wb = smoothstep(0.0, ct18b1LevelBlend(L), ct18b1BorderDist(ct18b1LevelRect(L), xy));
  89    float wr = 1.0 - smoothstep(1.0, 2.0, fp / ct18b1LevelTexel(L));
  90    return wb * wr;
  91  }
  92  
  93  // Height above the tangent plane's sphere datum, blended over levels.
  94  // Mirrors generate_terrain.composite_height().
  95  float ct18b1TerrainHeight(vec2 xy, float fp) {
  96    float h = ct18b1LevelHeight(2, xy);
  97    float w1 = ct18b1LevelWeight(1, xy, fp);
  98    if (w1 > 0.0) {
  99      vec4 r = CT18B1_L1_RECT;
 100      h = mix(h, ct18b1LevelHeight(1, clamp(xy, r.xy, r.zw)), w1);
 101    }
 102    float w0 = ct18b1LevelWeight(0, xy, fp);
 103    if (w0 > 0.0) {
 104      vec4 r = CT18B1_L0_RECT;
 105      h = mix(h, ct18b1LevelHeight(0, clamp(xy, r.xy, r.zw)), w0);
 106    }
 107    return h;
 108  }
 109  
 110  // Sphere: surface lies below the tangent plane by r^2 / (R + sqrt(R^2 - r^2)).
 111  float ct18b1CurvatureDrop(vec2 xy) {
 112    float r2 = dot(xy, xy);
 113    return r2 / (CT18B1_R_MOON + sqrt(CT18B1_R_MOON * CT18B1_R_MOON - r2));
 114  }
 115  
 116  // z coordinate of the surface in the model (site) frame
 117  float ct18b1TerrainZ(vec2 xy, float fp) {
 118    return ct18b1TerrainHeight(xy, fp) - ct18b1CurvatureDrop(xy);
 119  }
 120  
 121  // Outward surface normal, model frame, from central differences of the same height
 122  // function (the sphere's curvature is included through ct18b1TerrainZ).
 123  vec3 ct18b1TerrainNormal(vec2 xy, float fp) {
 124    float eps = max(fp, CT18B1_L0_TEXEL);
 125    float hx = ct18b1TerrainZ(xy + vec2(eps, 0.0), fp) - ct18b1TerrainZ(xy - vec2(eps, 0.0), fp);
 126    float hy = ct18b1TerrainZ(xy + vec2(0.0, eps), fp) - ct18b1TerrainZ(xy - vec2(0.0, eps), fp);
 127    return normalize(vec3(-hx / (2.0 * eps), -hy / (2.0 * eps), 1.0));
 128  }
 129  
 130  // Finest level usable for footprint fp at xy (for the shadow march)
 131  int ct18b1MarchLevel(vec2 xy, float fp) {
 132    if (fp <= 2.0 * CT18B1_L0_TEXEL && ct18b1BorderDist(CT18B1_L0_RECT, xy) > 0.0) return 0;
 133    if (fp <= 2.0 * CT18B1_L1_TEXEL && ct18b1BorderDist(CT18B1_L1_RECT, xy) > 0.0) return 1;
 134    return 2;
 135  }
 136  
 137  // Bilinear (not nearest-texel) so the horizon angle, and therefore shadow edges and
 138  // penumbrae, vary continuously instead of in 2.5 cm steps.
 139  float ct18b1MarchZ(vec2 xy, float fp) {
 140    return ct18b1LevelHeight(ct18b1MarchLevel(xy, fp), xy) - ct18b1CurvatureDrop(xy);
 141  }
 142  
 143  // Fraction of a circular source disk (angular radius rho) that lies above the terrain
 144  // horizon seen from P in direction dir.  Exact for a heightfield in the vertical
 145  // plane of the source; the disk-above-chord area fraction is analytic.
 146  float ct18b1SourceVisibility(vec3 P, vec3 dir, float rho, float fp) {
 147    float lh = length(dir.xy);
 148    float elev = asin(clamp(dir.z, -1.0, 1.0));
 149    if (lh < 1e-6) return 1.0;
 150    vec2 dh = dir.xy / lh;
 151    float maxTan = -1e9;
 152    float t = max(2.0 * fp, 1.5 * CT18B1_L0_TEXEL);
 153    float tanLow = tan(elev - rho);
 154    for (int k = 0; k < CT18B1_SHADOW_STEPS; k++) {
 155      vec2 q = P.xy + dh * t;
 156      if (ct18b1BorderDist(CT18B1_L2_RECT, q) < 0.0) break;
 157      float z = ct18b1MarchZ(q, max(fp, t * CT18B1_MARCH_FP_PER_M));
 158      maxTan = max(maxTan, (z - P.z) / t);
 159      if ((CT18B1_ZMAX - P.z) / t < tanLow) break;   // nothing further can block
 160      t *= CT18B1_MARCH_GROWTH;
 161    }
 162    float hor = atan(maxTan);
 163    float x = clamp((elev - hor) / rho, -1.0, 1.0);
 164    return 0.5 + (asin(x) + x * sqrt(1.0 - x * x)) / CT18B1_PI;
 165  }
 166  
 167  // Visibility of a point light at Lp from surface point P (hard test with a 1 cm
 168  // soft band, straight segment against the heightfield).
 169  float ct18b1PointVisibility(vec3 P, vec3 Lp, float fp) {
 170    vec3 d = Lp - P;
 171    float D = length(d);
 172    float vis = 1.0;
 173    float t0 = max(2.0 * fp, 1.5 * CT18B1_L0_TEXEL);
 174    for (int k = 1; k <= CT18B1_LAMP_STEPS; k++) {
 175      float t = mix(t0, D - 0.05, float(k) / float(CT18B1_LAMP_STEPS));
 176      if (t <= t0) continue;
 177      vec3 q = P + d * (t / D);
 178      float z = ct18b1MarchZ(q.xy, max(fp, t * CT18B1_MARCH_FP_PER_M));
 179      vis = min(vis, smoothstep(-0.01, 0.01, q.z - z));
 180    }
 181    return vis;
 182  }
 183  
 184  // The triangle mesh is only a bounding proxy: its vertices are spaced up to a few
 185  // times coarser than the 2.5 cm heightfield, so a small rock can be a low mound in the
 186  // mesh.  This finds the true first intersection of the view ray with the heightfield
 187  // within +/- 3 vertex spacings of the rasterised point (fixed-step search followed by
 188  // bisection), so shading, normals and shadows are evaluated on the real surface.
 189  // Returns the refined point, or the rasterised point if no crossing is found.
 190  vec3 ct18b1RefineHit(vec3 C, vec3 Pm, float spacing, float fp) {
 191    vec3 d = Pm - C;
 192    float tm = length(d);
 193    vec3 rd = d / tm;
 194    float D = 3.0 * spacing + CT18B1_REFINE_EXTRA;
 195    float t0 = max(tm - D, 0.02);
 196    float t1 = tm + D;
 197    vec3 q = C + rd * t0;
 198    if (q.z < ct18b1TerrainZ(q.xy, fp)) return Pm;          // started below: give up
 199    float tPrev = t0;
 200    for (int k = 1; k <= CT18B1_REFINE_STEPS; k++) {
 201      float t = mix(t0, t1, float(k) / float(CT18B1_REFINE_STEPS));
 202      q = C + rd * t;
 203      if (q.z < ct18b1TerrainZ(q.xy, fp)) {
 204        float a = tPrev, b = t;
 205        for (int j = 0; j < 8; j++) {
 206          float m = 0.5 * (a + b);
 207          vec3 qm = C + rd * m;
 208          if (qm.z < ct18b1TerrainZ(qm.xy, fp)) b = m; else a = m;
 209        }
 210        vec3 hit = C + rd * b;
 211        return vec3(hit.xy, ct18b1TerrainZ(hit.xy, fp));
 212      }
 213      tPrev = t;
 214    }
 215    return Pm;
 216  }
 217  
 218  #endif // CT18B1_TERRAIN_COMMON_GLSL

`shaders/ct18b1_terrain_vs.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - terrain vertex shader (RenderableModel custom VertexShader).
   3   * The mesh is flat (z = 0); heights come from the same heightfield the fragment shader
   4   * uses, prefiltered to the local vertex spacing, so geometry and shading agree.
   5   *   in_st.x = local vertex spacing in metres (negative = skirt vertex)
   6   ****************************************************************************************/
   7  #version __CONTEXT__
   8  
   9  #include "PowerScaling/powerScaling_vs.hglsl"
  10  #include "ct18b1_params.glsl"
  11  #include "ct18b1_terrain_common.glsl"
  12  
  13  layout(location = 0) in vec4 in_position;
  14  layout(location = 1) in vec2 in_st;
  15  
  16  out vec3 vs_positionModel;
  17  out float vs_spacing;
  18  out vec4 vs_positionCameraSpace;
  19  out float vs_screenSpaceDepth;
  20  
  21  uniform mat4 modelViewTransform;
  22  uniform mat4 projectionTransform;
  23  uniform mat4 meshTransform;
  24  
  25  void main() {
  26    vec2 xy = in_position.xy;
  27    float spacing = abs(in_st.x);
  28    float z = ct18b1TerrainZ(xy, spacing);
  29    if (in_st.x < 0.0) {
  30      z -= CT18B1_SKIRT_DEPTH;
  31    }
  32    vec4 pm = vec4(xy, z, 1.0);
  33    vs_positionModel = pm.xyz;
  34    vs_spacing = spacing;
  35    vs_positionCameraSpace = modelViewTransform * (meshTransform * pm);
  36    vec4 positionClipSpace = projectionTransform * vs_positionCameraSpace;
  37    vec4 positionScreenSpace = z_normalization(positionClipSpace);
  38    gl_Position = positionScreenSpace;
  39    vs_screenSpaceDepth = positionScreenSpace.w;
  40  }

`shaders/ct18b1_globe_fs.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - globe fragment shader.  Same Hapke model, Sun + earthshine + lamp.
   3   * The globe's own relief is too coarse to cast resolved shadows; sub-resolution
   4   * shadowing is carried by Hapke's macroscopic roughness term.
   5   ****************************************************************************************/
   6  #include "fragment.glsl"
   7  #include "ct18b1_params.glsl"
   8  #include "ct18b1_hapke.glsl"
   9  #include "ct18b1_lighting.glsl"
  10  
  11  in vec3 vs_positionModel;
  12  in vec3 vs_normalModel;
  13  in vec4 vs_positionCameraSpace;
  14  in float vs_screenSpaceDepth;
  15  
  16  Fragment getFragment() {
  17    Fragment frag;
  18    frag.depth = vs_screenSpaceDepth;
  19    frag.gPosition = vs_positionCameraSpace;
  20    frag.disableLDR2HDR = true;
  21    frag.color.a = opacity;
  22  
  23    Ct18b1Frame F = ct18b1Frame();
  24    HapkeParams hp = ct18b1Params(F.mode);
  25    vec3 P = vs_positionModel;
  26    vec3 N = normalize(vs_normalModel);
  27    vec3 V = normalize(F.camera - P);
  28    frag.gNormal = vec4(normalize(mat3(modelViewTransform) * N), 0.0);
  29  
  30    float total = 0.0;
  31    float gS = 0.0, gE = 0.0, gL = 0.0;
  32    if (F.sunScale > 0.0) {
  33      total += F.sunScale * ct18b1Reflect(N, F.sunDir, V, hp, F.mode, gS);
  34    }
  35    if (F.earthScale > 0.0) {
  36      total += F.earthScale * ct18b1EarthIrradiance(F.sunDir, F.earthDir) *
  37               ct18b1Reflect(N, F.earthDir, V, hp, F.mode, gE);
  38    }
  39    if (F.lampScale > 0.0) {
  40      vec3 l;
  41      float dist;
  42      float eL = ct18b1LampIrradiance(F, P, l, dist);
  43      if (eL > 1e-12) {
  44        total += eL * ct18b1Reflect(N, l, V, hp, F.mode, gL);
  45      }
  46    }
  47    vec3 colour = vec3(total * F.gain);
  48    if (F.mode == 5) {
  49      colour = ct18b1DisplayToLinear(ct18b1AngleColour(degrees(acos(clamp(dot(F.sunDir, V), -1.0, 1.0)))));
  50    }
  51    frag.color.rgb = colour;
  52    return frag;
  53  }

`shaders/ct18b1_globe_vs.glsl`

   1  /*****************************************************************************************
   2   * CT18.B1.Openspace - globe vertex shader: the whole test body (lunar radius, synthetic
   3   * relief), positions and normals baked in the mesh, model frame = IAU_MOON body-fixed.
   4   ****************************************************************************************/
   5  #version __CONTEXT__
   6  
   7  #include "PowerScaling/powerScaling_vs.hglsl"
   8  
   9  layout(location = 0) in vec4 in_position;
  10  layout(location = 1) in vec2 in_st;
  11  layout(location = 2) in vec3 in_normal;
  12  
  13  out vec3 vs_positionModel;
  14  out vec3 vs_normalModel;
  15  out vec4 vs_positionCameraSpace;
  16  out float vs_screenSpaceDepth;
  17  
  18  uniform mat4 modelViewTransform;
  19  uniform mat4 projectionTransform;
  20  uniform mat4 meshTransform;
  21  uniform mat4 meshNormalTransform;
  22  
  23  void main() {
  24    vec4 pm = meshTransform * in_position;
  25    vs_positionModel = pm.xyz;
  26    vs_normalModel = normalize(mat3(meshNormalTransform) * in_normal);
  27    vs_positionCameraSpace = modelViewTransform * pm;
  28    vec4 positionClipSpace = projectionTransform * vs_positionCameraSpace;
  29    vec4 positionScreenSpace = z_normalization(positionClipSpace);
  30    gl_Position = positionScreenSpace;
  31    vs_screenSpaceDepth = positionScreenSpace.w;
  32  }

`shaders/ct18b1_params.glsl`

   1  // GENERATED by tools/generate_terrain.py from tools/ct18b1_config.py - do not edit.
   2  // Every constant below is traceable in AUDIT.md section 4 (parameter table).
   3  #ifndef CT18B1_PARAMS_GLSL
   4  #define CT18B1_PARAMS_GLSL
   5  const float CT18B1_PI = 3.14159265358979324;
   6  const float CT18B1_R_MOON = 1737400.0;
   7  const float CT18B1_W = 0.12;
   8  const float CT18B1_B = 0.41;
   9  const float CT18B1_C = 0.1;
  10  const float CT18B1_BS0 = 2.81456954;
  11  const float CT18B1_HS = 0.066877;
  12  const float CT18B1_BC0 = 0.08;
  13  const float CT18B1_HC = 0.048915;
  14  const float CT18B1_THETA_BAR = 0.141371669;
  15  const float CT18B1_K = 1.0;
  16  const float CT18B1_E_SUN_LUX = 133800.0;
  17  const float CT18B1_SUN_ANG_RADIUS = 0.00465241753;
  18  const float CT18B1_EARTH_ALBEDO = 0.226;
  19  const float CT18B1_EARTH_SOLID = 0.000274693544;
  20  const float CT18B1_EARTH_ANG_RADIUS = 0.0165738814;
  21  const float CT18B1_LAMP_CD = 20000.0;
  22  const float CT18B1_LAMP_FWHM = 0.698131701;
  23  const float CT18B1_LAMP_OFFSET = 0.1016;
  24  const float CT18B1_NORMAL_ALBEDO_LAMBERT = 0.0947780753;
  25  const float CT18B1_SKIRT_DEPTH = 40.0;
  26  const float CT18B1_ZMAX = 33.6398947;   // highest representable terrain height (m)
  27  const int   CT18B1_ATLAS_H = 3072;
  28  const int   CT18B1_SHADOW_STEPS = 110;
  29  const float CT18B1_MARCH_GROWTH = 1.12;
  30  const float CT18B1_MARCH_FP_PER_M = 0.02;
  31  const int   CT18B1_LAMP_STEPS = 24;
  32  const int   CT18B1_REFINE_STEPS = 32;
  33  const float CT18B1_REFINE_EXTRA = 0.3;
  34  const float CT18B1_EXPOSURE_DECADES = 6.0;
  35  const float CT18B1_OS_EXPOSURE = 3.7;
  36  const float CT18B1_OS_GAMMA = 0.95;
  37  // level 0: near  texel 0.025 m  quantisation 0.1111 mm
  38  const vec4  CT18B1_L0_RECT = vec4(-64.0, -12.8, 64.0, 12.8);
  39  const float CT18B1_L0_TEXEL = 0.025;
  40  const ivec2 CT18B1_L0_DIMS = ivec2(5120, 1024);
  41  const ivec2 CT18B1_L0_ATLAS = ivec2(0, 0);
  42  const float CT18B1_L0_HOFF = -13.9109172;
  43  const float CT18B1_L0_HSCALE = 0.000111071132;
  44  const float CT18B1_L0_BLEND = 0.5;
  45  // level 1: mid  texel 0.125 m  quantisation 0.1418 mm
  46  const vec4  CT18B1_L1_RECT = vec4(-128.0, -128.0, 128.0, 128.0);
  47  const float CT18B1_L1_TEXEL = 0.125;
  48  const ivec2 CT18B1_L1_DIMS = ivec2(2048, 2048);
  49  const ivec2 CT18B1_L1_ATLAS = ivec2(0, 1024);
  50  const float CT18B1_L1_HOFF = -14.5560591;
  51  const float CT18B1_L1_HSCALE = 0.000141783527;
  52  const float CT18B1_L1_BLEND = 8.0;
  53  // level 2: far  texel 4.0 m  quantisation 1.6804 mm
  54  const vec4  CT18B1_L2_RECT = vec4(-3200.0, -3200.0, 3200.0, 3200.0);
  55  const float CT18B1_L2_TEXEL = 4.0;
  56  const ivec2 CT18B1_L2_DIMS = ivec2(1600, 1600);
  57  const ivec2 CT18B1_L2_ATLAS = ivec2(2048, 1024);
  58  const float CT18B1_L2_HOFF = -76.487997;
  59  const float CT18B1_L2_HSCALE = 0.00168044391;
  60  const float CT18B1_L2_BLEND = 0.0;
  61  #endif

`lua/ct18b1_light_helper_sun.lua`

   1  -- GENERATED by tools/make_assets.py (CT18.B1.Openspace).  Light-direction helper: Sun.
   2  --
   3  -- Why this node exists (AUDIT.md 5.3): in OpenSpace 0.21.4,
   4  -- SceneGraphLightSource::directionViewSpace computes
   5  --     combinedViewMatrix * vec4(lightPosition - nodePosition, 1.0)
   6  -- The w = 1 applies the camera translation to a DIRECTION, so the result is
   7  -- proportional to R (L - N - C), C = camera world position.  For the Sun the error is
   8  -- ~0.2 deg; for the Earth the direction would point roughly at the Sun.
   9  -- Placing the light-source node at P = L + C cancels the error exactly:
  10  --     R (P - N - C) = R (L - N).
  11  -- C is read from the navigation state; if that fails, C = position of CT18B1_Site
  12  -- (the camera is within a few km of it in all rover tests: error < 0.001 deg).
  13  
  14  local TARGET = "Sun"
  15  local FALLBACK = "CT18B1_Site"
  16  
  17  local function v3(t)
  18    if t == nil then return 0.0, 0.0, 0.0 end
  19    if t.x ~= nil then return t.x, t.y, t.z end
  20    return t[1], t[2], t[3]
  21  end
  22  
  23  function cameraWorld()
  24    local ns = openspace.navigation.getNavigationState("Root")
  25    local ax, ay, az = v3(openspace.worldPosition(ns.Anchor))
  26    local px, py, pz = v3(ns.Position)
  27    return ax + px, ay + py, az + pz
  28  end
  29  
  30  function translation(simTime, prevSimTime, wallTime)
  31    local lx, ly, lz = v3(openspace.worldPosition(TARGET))
  32    local ok, cx, cy, cz = pcall(cameraWorld)
  33    if not ok then
  34      cx, cy, cz = v3(openspace.worldPosition(FALLBACK))
  35    end
  36    return { lx + cx, ly + cy, lz + cz }
  37  end