Task Specification: Real-time Solar System Model (Single HTML File)

This document is an implementation specification. Copy the data tables and formulas exactly; do not rely on memory or guess the values. Write the rendering code, shaders, and UI code yourself. Read all of Section 7, “Pitfalls You Must Avoid,” before you begin—every item describes a problem you will encounter in practice.

---

0. Goal

Build a real-time model of the Solar System as a single HTML file that runs when opened directly in a browser.

There are two core requirements:

1. Be as realistic as possible—celestial positions must be calculated in real time from genuine orbital elements, not from a preset animation or approximate circular motion. Size and distance scales may be compressed (otherwise they cannot be seen together), but angular positions, orbital phases, periods, eccentricities, and inclinations must be real.
2. Max out the visuals—this is not merely a diagram; it should be something people want to keep looking at.

---

1. Hard Constraints

About “no Three.js”: this means you must write the matrix library, sphere-mesh generation, camera controls, and framebuffer management yourself. Inlining a 600 KB library would make the file size balloon, and features such as the bloom pipeline and mutual shadowing between planets and rings require custom shaders anyway.

Texture sources: external images are forbidden, so every celestial surface must be generated on the fly with procedural noise in fragment shaders. This is one of the main challenges of the task.

---

2. Astronomical Data (copy exactly; do not invent values)

2.1 Planetary orbital elements (NASA/JPL, “Approximate Positions of the Major Planets”)

Six elements plus their rates of change per Julian century. Nominal validity interval: 1800–2050.

Units: a is in AU; all other values are in degrees.

Rates of change (per Julian century, or 36,525 days):

Strictly speaking, the Earth row describes the Earth–Moon barycenter. The visual difference is imperceptible, so use it directly as Earth.

2.2 Rotation periods and axial tilts

rotH = rotation period in hours; a negative value denotes retrograde rotation.

tilt = angle in degrees between the rotation axis and the north ecliptic pole. poleLon = ecliptic longitude in degrees toward which the rotation axis projects on the ecliptic plane.

Note: textbook “axial tilt” is measured relative to a body’s own orbital plane, whereas the values here are relative to the ecliptic plane; they differ by several degrees. Use the values below for rendering because the world coordinate system is ecliptic, while the information panel may show textbook values. These values have already been converted from IAU pole right ascension and declination; copy them directly.

Direction vector of the rotation axis in ecliptic Cartesian coordinates:

axis_ecl = ( sin(tilt)·cos(poleLon),  sin(tilt)·sin(poleLon),  cos(tilt) )

The rotation axis must remain fixed in the inertial frame rather than changing with orbital motion. This lets Saturn’s ring opening angle naturally widen and narrow across its 29.5-year orbit, and makes Uranus genuinely “roll on its side.”

2.3 Planetary physical parameters (for the information panel and display radii)

The information panel should also show orbital eccentricity, orbital inclination, axial tilt, atmospheric composition, and a two- or three-sentence educational description.

2.4 Moons

Uniform circular motion with the real period is sufficient; you do not need to solve Kepler’s equation. The periods themselves must be real, however, so the 1:2:4 orbital resonance of Io, Europa, and Ganymede emerges naturally.

• The Moon’s reference plane is the ecliptic plane—its 5.145° inclination is relative to the ecliptic. Every other moon uses its parent planet’s equatorial plane. This reflects reality; do not force them all onto one reference plane.
• Triton’s negative period denotes retrograde orbital motion, which must be represented.

2.5 Halley’s Comet, 1P/Halley

Compute the mean anomaly directly from the perihelion epoch: M = 2π·(JD − Tp) / P_days.

2.6 Asteroid belt / Kuiper belt (procedural generation)

Solve Kepler’s equation for every small body in the vertex shader. Each particle stores eight floating-point values: a, e, I, Ω, ω, M₀, n, and size. The CPU generates them only once during initialization. Suggested counts: 9,000 main-belt objects, 900 Hildas, 2,000 Trojans, and 6,000 Kuiper-belt objects.

Main belt: a ∈ [2.06, 3.28] AU. Use rejection sampling to carve out the Kirkwood gaps. These are genuine resonance gaps and implementing them adds a great deal of credibility:

Retention probability = Π(1 − strength·exp(−((a−position)/σ)²)), multiplied by an overall density envelope centered at 2.72 AU with σ=0.55. Eccentricity e follows a half-normal distribution with a mean of about 0.14 and a maximum of 0.34; inclination I follows a half-normal distribution with a mean of about 8°.

Hilda group: a ≈ 3.90 ± 0.09 AU, e ∈ [0.10, 0.30].

Jupiter Trojans: this is the standout detail. Use a ≈ 5.202, but the mean anomaly must remain locked 60° ahead of or behind Jupiter:

M₀ = L_Jupiter(J2000) ± 60° + Gaussian perturbation (σ≈26°)     // L_Jupiter(J2000) = 34.396°
n  = Jupiter’s mean motion (do not derive n from each object’s own a, or the clusters will disperse after a few decades)

This keeps the L4 and L5 clusters moving with Jupiter indefinitely.

Kuiper belt: 22% plutinos (a ≈ 39.45 ± 0.55, in a 3:2 resonance with Neptune), 64% classical-belt objects (a ∈ [42.5, 47.7], e < 0.11, mean I of 7°), and 14% scattered-disc objects (a ∈ [48, 80], e ∈ [0.15, 0.5], mean I of 22°).

Rendering: use gl.POINTS. Point size falls with distance and is clamped to [0.75, 4.5] pixels. Use additive blending and multiply brightness by a factor that falls with heliocentric distance so distant objects are dimmer.

2.7 Ring systems

Describe each ring system with a radial density function, ringDensity(r), where r is measured in parent-planet radii. This directly models optical depth rather than using a texture.

Saturn (1.235–2.34 parent-planet radii):

Finally, multiply the whole profile by high-frequency noise, 0.70 + 0.60·fbm(r·230), to create fine ring banding.

Uranus: a set of extremely narrow, dark rings built by adding narrow Gaussians. The ring geometry should span 1.60–2.00; the actual rings lie within it. The ε ring is at 1.951 (τ≈0.55, σ≈0.0055), with another six at 1.909 / 1.863 / 1.823 / 1.786 / 1.751 / 1.690. Their τ decreases from 0.20 to 0.07 and σ ranges from 0.0015 to 0.0025. Their color is nearly charcoal black because the real albedo is only 0.03. Uranus has 13 known rings in reality; implementing the seven principal rings is sufficient.

Jupiter: an extremely tenuous dust ring spanning 1.40–1.81, with τ ≈ 0.016. It must remain very faint; see Section 7.4 for the reason.

---

3. Core Algorithms

3.1 Time

JD = unix_ms / 86400000 + 2440587.5
T  = (JD − 2451545.0) / 36525          // Julian centuries since J2000

3.2 Kepler’s equation

Given mean anomaly M and eccentricity e, find eccentric anomaly E by solving E − e·sin(E) = M:

First normalize M to [−π, π]
Initial value: E = M + e·sin(M)·(1 + e·cos(M))
Perform 8–10 Newton iterations: E −= (E − e·sin(E) − M) / (1 − e·cos(E))
Convergence criterion: |ΔE| < 1e-13

For orbits with e < 0.97, including Halley’s Comet, this initial estimate converges in three to five iterations.

3.3 Orbital elements → ecliptic Cartesian coordinates

1. Update all six elements with T: a = a₀ + (da/dt)·T, and likewise for the others
2. Argument of perihelion  ω = ϖ − Ω
   Mean anomaly            M = L − ϖ         (normalize to [−π, π])
3. Solve Kepler’s equation for E
4. Coordinates in the orbital plane:
   x' = a·(cos E − e)
   y' = a·√(1 − e²)·sin E
5. Rotate into ecliptic coordinates:
   x = (cosω·cosΩ − sinω·sinΩ·cosI)·x' + (−sinω·cosΩ − cosω·sinΩ·cosI)·y'
   y = (cosω·sinΩ + sinω·cosΩ·cosI)·x' + (−sinω·sinΩ + cosω·cosΩ·cosI)·y'
   z = (sinω·sinI)·x'              + (cosω·sinI)·y'

3.4 Coordinate-system mapping

Ecliptic coordinates form a right-handed system with the z-axis pointing toward the north ecliptic pole, while WebGL conventionally uses an upward y-axis. Map them as follows:

world = ( x_ecl,  z_ecl,  −y_ecl )

This mapping preserves handedness.

Milky Way background orientation: the north galactic pole is approximately (−0.8678, 0.4972, 0.0003) in ecliptic coordinates; the vector shown is already in world coordinates. Use it to orient the Milky Way band so the relative orientation of the star field and ecliptic plane is correct.

3.5 Display-scale compression

Neptune is at 30 AU and Mercury at 0.39 AU, so the inner planets collapse to a point at a linear scale. Use radial power-law compression:

r_display = SCALE · r_AU ^ COMPRESS        // Suggested: SCALE = 34, COMPRESS = 0.62
Keep direction unchanged; compress radial distance only

Key property: the power law changes only radial distance, leaving angular position completely unchanged. Planetary longitudes, relative positions, and conjunction/opposition relationships therefore remain exact. Explain this in the UI.

Enlarge celestial display radii with another power law:

R_display = 0.72 · (radius_km / 6371) ^ 0.42        // Use 4.6 separately for the Sun

You must provide a “true distance scale” toggle that moves COMPRESS between 0.62 and 1.0 with an animated transition over about 0.4 seconds. At 1.0, the inner planets collapse into a cluster and the asteroid belt becomes a dense narrow ring. This highly persuasive comparison is one of the work’s highlights. Rebuild the orbit curves every frame during the transition.

3.6 Rotation

Model matrix = translation(world position) × orthonormal basis with the rotation axis as Y × rotation about Y(rotation phase) × scale(display radius).

To construct the orthonormal basis, use axis as the Y column, take the cross product with any nonparallel reference vector to obtain the X column, then take another cross product for the Z column.

The rotation phase must be rate-limited; see Section 7.11:

Per-frame Δphase = clamp(2π·Δt_days / rotation_period_days,  −0.22,  +0.22)

3.7 Moon positions

θ = 2π·(JD − J2000)/period + initial phase
Local coordinates = (cos θ · d, 0, sin θ · d), then tilt by inc about the X axis
Then:
  Reference plane = ecliptic  → add directly to the parent planet’s world coordinates
  Reference plane = equatorial → first transform by the parent planet’s rotation-axis basis, then add

Display distance—there is a pitfall here; see Section 7.6:

a₀ = actual orbital radius of the planet’s innermost implemented moon
base = max( 2.35,  outer ring edge + innermost moon display radius / parent display radius + 0.55 )   // Units: parent-planet radii
d = parent display radius × base × (a_moon / a₀) ^ 0.55

---

4. Rendering Architecture

4.1 Pipeline overview

[MSAA framebuffer RGBA16F + DEPTH24]
   1. Background nebula       depth test off, depth writes off
   2. Star points             depth test off, additive blending
   3. Sun sphere              depth test on, depth writes on, back-face culling
   4. Planet/moon spheres     same as above
   5. Orbit ribbons           depth test on, depth writes off, additive blending
   6. Asteroid belts          same as above
   7. Atmosphere shells       same as above
   8. Ring systems            depth test on, depth writes off, normal alpha blending, no culling
   9. Comet tails + coma      additive blending
  10. Solar corona            additive blending
        ↓ resolve MSAA with blitFramebuffer
[Scene texture RGBA16F]
        ↓ brightness-threshold extraction + 6 downsample levels (13-tap)
        ↓ 6 upsample levels (3×3 tent, accumulated with additive blending; see 7.3)
[Bloom texture]
        ↓
[Composite: chromatic aberration → add bloom → lens flare → exposure → ACES tone mapping → gamma → vignette → grain]
[Default framebuffer]

The EXT_color_buffer_float extension is required for HDR. If unavailable, fall back to RGBA8 and lower the bloom threshold.

Implement 4× MSAA with renderbufferStorageMultisample and blitFramebuffer. If checkFramebufferStatus reports an incomplete framebuffer, fall back gracefully to a non-MSAA path.

4.2 Key points for each pass

Background: render a full-screen triangle and reconstruct the view direction from the inverse VP matrix. Draw the Milky Way band with 3D fbm noise, using exp(−(view direction·north galactic pole)²·24) for its band profile, then add a layer of dark dust. The deep-space base color must be extremely dark, around 0.002 in linear space, or the entire scene will look gray; see Section 7.7.

Stars: approximately 11,000 points. Generate directions with rejection sampling biased toward the galactic plane. Use pow(random, 3.6) for brightness, yielding many dim stars and a few bright ones. Sample color from a seven-step color-temperature palette according to realistic spectral-type proportions, from orange-leaning M-type stars to blue-leaning B-type stars. Place them at camera position + direction×4000, with depth testing disabled, and add subtle twinkling.

Sun: see Section 5.1.

Solar corona: use a camera-facing billboard quad with additive blending and depth testing enabled, so planets correctly occlude it when they pass in front. Draw two layers: an inner layer at 3.1 solar radii and an outer layer at 9 radii. Use fbm over polar angle to create radial coronal streamers.

Orbit ribbons: do not use gl.LINES, because line width is fixed at 1 px in most implementations and cannot be antialiased. Expand a triangle strip in screen space: emit two vertices for each sample point (side = ±1), project the tangent into screen space in the vertex shader, take its normal, and offset by the pixel width. Use 1−|side| in the fragment shader for soft edges.

Give opacity a gradient along the orbit: brightest at the body’s current position, with a trailing decay opposite its direction of motion, using alpha = base + peak·pow(max(0, 1−d·k), 3.2), where d is the lag measured around the orbit parameter. This “comet-tail” orbit looks far better than a uniformly bright ellipse.

Comet tails: render two. The ion tail points straight away from the Sun and is blue; the dust tail is wider, warmer yellow, and curves backward along the velocity direction. Scale length and brightness by (3/r_AU)², and do not draw them when r > 6.5 AU. Build them as billboard ribbons facing the camera around the tail axis, and rebuild them on the CPU each frame; with only dozens of vertices, the cost is negligible.

4.3 Post-processing

• Bloom: extract with a soft-knee threshold (threshold ≈ 1.25, knee width ≈ 0.6), downsample six levels with a 13-tap filter, upsample and accumulate with a 3×3 tent filter, then add the result back at a strength of 0.34.
• Tone mapping: use the Narkowicz approximation of ACES, followed by gamma 2.2.
• Lens flare: requires the Sun’s screen position and a visibility factor. Compute visibility on the CPU: determine whether the Sun is in front of the camera and whether any planet occludes it by projecting each planet, comparing screen-space distance with apparent radius, and comparing depth. The flare should contain a core glow, six starburst rays, a horizontal anamorphic streak, and three or four ghost discs along the line from the Sun to screen center. Keep the overall intensity restrained; see Section 7.8.
• Add subtle chromatic aberration by slightly offsetting RGB samples according to distance from center, along with vignette and animated grain.

4.4 Depth precision

The scene spans enormous scales: planet radii around 0.5, Neptune’s compressed orbit around 274, and its true-scale orbit around 1020. Use dynamic near and far planes:

near = max(0.012, distance from camera to target × 0.0035)
far  = max(distance from camera to target × 14, 4200)

---

5. Procedural Celestial Surfaces

Sample all surface noise in object space using normalized position p = normalize(object-space coordinates). This makes the texture rotate with the model matrix instead of “floating” over the surface.

Required noise tools: hash13 / hash33, 3D value noise, fbm with a rotation matrix between octaves to avoid axial alignment, Worley F1, and dual-distance Worley (F1, F2).

Bump mapping: do not evaluate the height function three times just to obtain a gradient. Screen-space derivatives do it in one evaluation:

vec3 dpx = dFdx(worldPos), dpy = dFdy(worldPos);
float dhx = dFdx(height), dhy = dFdy(height);
vec3 r1 = cross(dpy, N), r2 = cross(N, dpx);
float det = dot(dpx, r1);
vec3 grad = (r1*dhx + r2*dhy) / det;
N = normalize(N − strength * grad);

5.1 Sun

• Granulation: use Worley F2−F1 to draw the dark cellular boundaries; see Section 7.9. Combine two frequencies, approximately 76 and 31.
• Supergranulation: use low-frequency fbm, around 6.5, for large-scale brightness variation.
• Sunspots: threshold the high end of low-frequency fbm, with separate umbra and penumbra levels. Set the threshold high enough that spots cover only about 5% of the surface, or the whole hemisphere will look dark. In addition to lowering color temperature, multiply sunspots by a separate darkening factor so the umbra appears dark brown rather than orange.
• Limb darkening: limb = 0.24 + 0.76·pow(max(dot(N,V),0), 0.38). Apply it to both brightness and color temperature; the limb should become darker and redder, not merely darker.
• Map the final color through a five-step ramp from deep red → orange → yellow → yellow-white → white. Use an overall brightness around 2.0 in linear HDR values so bloom can lift it.
• Let the surface flow slowly by translating the noise coordinates over time.

5.2 Rocky bodies (Mercury, the Moon, Callisto, and Charon)

Impact-crater function—this is essential; a simple Worley threshold produces only a field of equally sized black dots:

In the Worley loop, track both the distance d1 to the nearest feature point and that cell’s random value id
Independent crater radius per cell: rad = 0.13 + 0.22·frac(id·7.31)
t = d1 / rad
Bowl depression: bowl = (t < 1) ? −(1 − t²) : 0
Raised rim: rim  = exp(−((t−1)/0.26)²)
Height = (bowl·0.90 + rim·0.50) · depth · existence flag
Albedo modulation = (bowl·0.55 + rim·0.32) · existence flag       // dark floor, bright rim
Existence flag = step(id, 0.74) · smoothstep(0.60, 0.40, d1)       // see 7.5 for the latter term

Combine four octaves—frequencies around 5 / 12 / 28 / 62 and depths around 1.0 / 0.52 / 0.26 / 0.13—to obtain a realistic crater-size distribution.

Add the following to the Moon: maria, using a low-frequency fbm threshold for dark basalt plains with fewer and shallower internal craters, plus bright ray systems.

5.3 Earth

• Land and sea: use seven-octave fbm as a height field, with a sea-level threshold of 0.552. The mean of fbm is 0.5, so a 0.5 threshold creates 50% land; the real oceans cover 71%.
• Ocean: interpolate between deep and shallow blue by elevation and add a turquoise coastal shelf. For specular highlights, use Blinn–Phong with a shininess around 420; lower values smear into a broad white patch. Apply it only to cloud-free ocean.
• Land: use a second noise layer to divide vegetation and desert, transition high elevations to rocky gray, and turn high latitudes ice-white.
• Clouds: use an independent fbm layer with latitude-dependent zonal drift over time—faster at low latitudes, slower at high ones. Modulate coverage by latitude: more cloud in the equatorial convergence zone and mid-latitudes, less in the subtropics.
• Night-side lights: over land, without clouds, and where dot(N,L) < 0, raise high-frequency noise to a high power to create point-like, warm-yellow emission.
• Apply bump mapping only to land; perturbing the ocean normal breaks up its specular highlight.

5.4 Venus

Only the thick clouds are visible; the surface is not. Use medium- to high-frequency fbm for swirling bands, warped by latitude, with colors ranging from creamy yellow to off-white and a high albedo of 0.75.

5.5 Mars

• Use a low-albedo red-brown base color; the real albedo is only 0.17. Separate the linear RGB channel ratios to roughly 1 : 0.37 : 0.16, or the result will look beige rather than rust red.
• Use the same craters as Section 5.2, but at half the depth.
• North–south dichotomy: the northern hemisphere is young plains, so reduce craters there by about 78%; the southern hemisphere is ancient highland with dense cratering. This is Mars’s most prominent geological feature.
• Add polar caps using latitude plus noise, and Valles Marineris as a long, dark Gaussian trench near the equator.

5.6 Gas giants (Jupiter and Saturn)

• Zonal bands: use latitude in degrees as the independent variable and combine four sine waves at different frequencies, corresponding to roughly 4 / 9 / 17 / 30 bands. Use two layers of fbm to warp latitude laterally and create turbulence.
• Make distortion strength adjustable per planet: strong for Jupiter’s violent turbulence and weak for Saturn, at about 0.42 times Jupiter’s strength, because Saturn’s bands are much softer.
• Fade bands toward the poles, beginning around 46° latitude and reaching a 94% reduction, or the poles will show concentric “target rings.”
• For oval storms, flatten Worley noise along latitude to create ellipses stretched in longitude.
• Great Red Spot: place an elliptical region at a specified longitude and latitude with a major-to-minor axis ratio around 2.15:1. Sample noise from polar angle plus radius to form a spiral, mixing radial and time terms into the angle, and add a bright rim.
• Saturn’s north-polar hexagon: where latitude > 0.80, use the hexagonal distance function hexEdge = 0.115 / cos(mod(longitude, 60°) − 30°) and draw a dark line along its boundary. It is a particularly rewarding detail.
• Add limb darkening.

5.7 Ice giants (Uranus and Neptune)

Use mostly featureless cyan/blue surfaces with only faint banding and a few bright methane-cloud streaks. Give Neptune an additional dark spot. Important: do not put the bright cloud streaks into the height field; see Section 7.10.

5.8 Other bodies

• Europa: a high-albedo ice surface plus a network of reddish-brown lineae. Draw fine lines in the low-value regions of abs(fbm − 0.5), combining three frequencies.
• Io: sulfur-yellow and orange, scattered with dark volcanic vents; add faint emission at the center of each vent.
• Ganymede: dark ancient regions plus bright grooved terrain. The grooves can likewise use the low-value regions of abs(fbm−0.5).
• Titan: dense orange haze that almost completely erases surface detail, with an especially soft terminator.
• Pluto: a dark equatorial band, Cthulhu Macula, plus a heart-shaped bright region, Sputnik Planitia, built from two Gaussian spots.

5.9 Atmospheres

Do not use “a slightly larger sphere plus Fresnel rim light”; that creates a glass shell with a hard boundary.

The correct approach integrates optical depth by view-ray impact parameter:

vec3 rd = normalize(fragment world position − camera);
vec3 oc = planet center − camera;
float tca = dot(oc, rd);
float b = sqrt(dot(oc,oc) − tca*tca) / planet radius;      // normalized impact parameter
float dens = exp(−max(b − 1.0, 0.0) / scale height);       // thickest at the limb, then decays exponentially outward
dens *= 1.0 − smoothstep(1.0, 0.960, b) * 0.86;            // sight lines over the disc traverse only half the atmosphere
// Determine day/night from the normal at the view ray’s closest approach to the center, not from the fragment normal
vec3 an = normalize(camera + rd*tca − planet center);

Then multiply by the illumination factor and a forward Mie-scattering term, pow(dot(rd, L)*0.5+0.5, 4). Set the shell geometry radius to planet radius × (1 + scale height×7.5), allowing the exponential term to decay to zero before the geometry ends. Use a scale height of 0.02–0.05 planet radii; the real value is around 0.0013 and is invisible to the eye, so it must be exaggerated.

Keep intensity low. Haze over the disc should be only 5%–12% of the surface brightness; only the limb should form a thin, bright arc. Too much intensity washes the entire day side white.

5.10 Ring systems (a major focus and visual highlight)

Use the density profiles in Section 2.7. During rendering:

μ_view = max(|dot(ring normal, view direction)|, 0.10)          // the 0.10 floor is essential; see 7.4
μ_sun  = max(|dot(ring normal, light direction)|, 0.10)
τ_view = d · 1.45 / μ_view
τ_sun  = d · 1.45 / μ_sun
alpha  = 1 − exp(−τ_view)

Front-lit (observer and Sun on the same side of the ring plane):
    brightness = (1 − exp(−τ_sun)) · 0.85
Back-lit (opposite sides):
    brightness = exp(−τ_sun·0.42) · 1.15 · (1 − exp(−τ_sun·0.75))
    // Result: the thin C Ring and Cassini Division glow through forward scattering, while the thick B Ring blocks light and darkens—this is a real phenomenon

Two-way shadowing, which looks excellent when implemented:

1. Planet shadow on rings: cast a ray from the ring point toward the Sun and perform an analytic intersection with the parent planet’s sphere. Let oc = ring point − planet center and bb = dot(oc, L). If bb < 0, compute the closest distance dmin = sqrt(dot(oc,oc) − bb²), then use smoothstep(R·0.97, R·1.09, dmin) for a shadow coefficient with a soft penumbra.
2. Ring shadow on planet: in the planet fragment shader, transform the Sun direction into object space. From the current fragment, find the intersection along that direction with the ring plane, object-space y = 0. If the intersection radius lies within the ring system, attenuate illumination by exp(−ringDensity(r)·1.45 / |sunDir.y|).

---

6. Interaction and Interface

Camera: use an orbital camera with yaw, pitch, and distance. Drag with the mouse to rotate, use the wheel to zoom, and support two-finger pinch zoom on touch devices. Clicking a body focuses it. Use screen-space projection for picking: project every body to the screen and select the nearest body within a threshold; GPU picking is unnecessary. When focus changes, smoothly transition both the camera target and distance over about one second rather than snapping.

Interface elements:

Compute real-time orbital speed with the vis-viva equation: v = 29.7847 · sqrt(2/r − 1/a) km/s, where r and a are in AU.

Time speed: use a logarithmic mapping from real time, one second per second, to about ten years per second, and support negative values for reverse playback. Suggested formula:

speed_days_per_second = sign(v) · 1.157407e-5 · 10^(|v|/100 · 8.55)      // v ∈ [−100, 100]

Automatically change the displayed unit by magnitude: seconds, minutes, hours, days, or years per second.

Honesty requirement: when the date falls outside 1800–2050, display “Outside the valid interval of the orbital elements (extrapolated).” Also explain somewhere that distance scale is compressed while angular positions remain exact.

Keyboard shortcuts: Space to pause, Left/Right Arrow to change speed, and H to hide the interface.

Visual style: dark glassmorphism panels using backdrop-filter: blur, hairline borders, small labels with generous letter spacing, and monospaced digits with font-variant-numeric: tabular-nums, preventing the interface from jittering as numbers change.

---

7. Pitfalls You Must Avoid

This is the most important section of the document. Every item below was encountered while building the project. Several can make the result look completely broken while producing no errors at all.

7.1 gl.clear(DEPTH_BUFFER_BIT) is constrained by depthMask ★★★

If depthMask is false when clear is called, the depth buffer is not cleared at all.

Typical symptoms: the first frame looks correct, but the depth buffer then remains frozen at its first-frame values. From a distant view, this may be invisible because the depth values happen to remain self-consistent. As soon as the camera moves closer or the scene scale changes, stale depth values cull all geometry—the Sun becomes a black disc and the planets disappear, while the star field and background remain because their depth tests are disabled. Meanwhile, gl.getError() returns 0 and the console remains clean.

Most pipelines disable depthMask during some pass, such as transparent geometry or post-processing, then forget to restore it. The next frame’s opening clear consequently does nothing.

Remedy: explicitly call gl.depthMask(true) before clear every frame.

7.2 Units for screen-space line width ★★

When expanding a line into a triangle strip, project the tangent into screen space, obtain its normal, then divide that normal by half the viewport size. The resulting vector’s length in NDC is exactly one pixel. The offset should therefore be multiplied directly by the half-width in pixels.

If you mistakenly divide by the viewport dimensions a second time, line width becomes roughly 0.001 pixels—the orbit lines become completely invisible, with no error reported.

7.3 Texture feedback loops during bloom upsampling ★★

If the target mip is sampled as a texture while it is also being written during upsampling, perhaps to “read the previous value and add to it,” the same texture is both read and written in a feedback loop. Behavior is undefined; some drivers silently produce garbage, and no GL error is guaranteed.

Remedy: use additive blending, blendFunc(ONE, ONE), to accumulate the blurred result into the target. The shader should sample only the source mip.

7.4 Thin rings turn into black lines at grazing angles ★★

The optical-depth formula τ = d / μ diverges as μ→0. When the view direction approaches the ring plane, even Jupiter’s extremely tenuous dust ring with d ≈ 0.016 produces alpha ≈ 1, drawing a solid black horizontal line across the planet’s face that looks like a rendering bug.

Remedy: (1) clamp μ to a minimum of 0.10 rather than 0.045; (2) keep Jupiter’s ring density sufficiently low, around τ ≈ 0.016. Saturn’s ring really should become a solid line at a grazing angle; that behavior is correct.

7.5 Voronoi features are abruptly clipped at cell boundaries ★

If a crater’s influence—especially its ejecta blanket—extends beyond its Voronoi cell boundary, the “winning cell” changes abruptly at that boundary. Crater parameters such as radius and existence jump, producing segments of straight seams across the surface that look like scratches.

Remedy: make every crater contribution fade to zero before approaching the cell boundary. For example, multiply the whole result by smoothstep(0.60, 0.40, d1), and also limit the maximum crater radius.

7.6 Moon orbits collide with ring systems ★

If moon display distances use a fixed baseline, such as 2.35 planet radii, they nearly coincide with Saturn’s outer ring edge at 2.34. Titan is the only implemented Saturnian moon and therefore receives a scale factor of 1, leaving half of its sphere buried in the rings.

Remedy: the baseline distance must be outer ring edge + innermost moon’s own display radius + margin; see Section 3.7.

7.7 Deep-space base color ★★

It is very easy to set the background ambient floor too high. A linear value of 0.01 becomes roughly 13% gray after ACES and gamma, covering the entire image with gray haze and destroying contrast in both stars and planets.

Remedy: use a deep-space base color around 0.002 in linear space, and reduce the Milky Way’s brightness accordingly. The completed image should be nearly pure black so the stars can genuinely appear bright.

7.8 HDR budget and overexposure ★★★

This is among the most time-consuming categories of problem to debug. ACES tone mapping rapidly pushes linear values above 1.5 toward white while also stripping saturation. Therefore:

• Use realistic albedo values: Mars 0.17, the Moon 0.12, Mercury 0.14. These are all dark. An “attractive” value of 0.7 blows them out to off-white, turning red Mars cream-colored.
• Do not make sunlight excessively intense. A suggested formula is intensity = clamp(1.5 · (1/r_AU)^0.40, 0.26, 2.1). Note that this attenuation is intentionally compressed; a physically exact 1/r² law makes Neptune completely invisible.
• Every time you add an emissive contribution—atmosphere, rim scattering, specular reflection, or self-emission—recheck for overexposure.
• Debugging method: use gl.readPixels to read pixel values directly from the HDR scene buffer and inspect the linear values instead of judging screenshots by eye. Target typical planet-disc values of 0.3–0.7; only the Sun, cloud tops, ice caps, and rings should exceed 1.0.

7.9 Worley F1 does not produce cellular structure ★★

F1, the distance to the nearest feature point, is smallest at the feature points themselves. Thresholding F1 therefore creates a collection of dots, not cells. Solar granulation and any “cell-wall network” require F2 − F1, the second-nearest distance minus the nearest distance, which approaches zero at the cell boundaries.

Track the two smallest distances in the same loop; the performance cost is virtually unchanged.

7.10 Step-like features must not enter the height field ★

If a sharp feature produced by smoothstep, such as a cloud streak or line marking, enters the height field at substantial weight, screen-space derivative bump mapping magnifies its edge into a deep, dark groove that looks like a scratch.

Remedy: give such features a small height-field weight, around 0.1, or let them affect albedo only and not height.

7.11 Rotation flicker at high time speeds ★

At a time speed of one year per second, Earth rotates about six times per frame, turning into randomly jittering visual noise.

Remedy: clamp the rotation phase’s per-frame increment to roughly ±0.22 radians. Continue displaying the real rotation period in the information panel. Clamp the increment, not the speed, so low-speed motion remains exact.

7.12 Triangle winding for UV spheres ★

When generating a sphere mesh yourself, reversed winding makes back-face culling remove the front faces. The difficulty is that front and back look almost identical on a self-emissive sphere such as the Sun, so the mistake may not be immediately obvious; it becomes visible only when lit planets are added.

Verification method: choose a point on the equator, calculate the cross product by hand, and confirm that the geometric normal points outward.

7.13 GLSL smoothstep(a, b, x) is undefined when a > b ★

When creating a “reversed gradient,” it is natural to write smoothstep(0.7, 0.3, x). The specification says the result is undefined when edge0 >= edge1. Although it works on most implementations, cross-platform consistency is not guaranteed.

Remedy: write your own sstep(a,b,x) { float t = clamp((x-a)/(b-a), 0.0, 1.0); return t*t*(3.0-2.0*t); } and use it throughout.

7.14 The mean of fbm is 0.5, not 0 ★

Multi-octave fbm with amplitudes 0.5, 0.25, 0.125, and so on has a range of roughly [0,1] and a mean of 0.5. If several fbm values are added and then fed into a color ramp, the mean drifts to 0.6–0.8 and saturates the upper end of the ramp, turning the entire body into a single color. The Sun is especially vulnerable and may become a pure white sphere.

Remedy: explicitly recenter and rescale values before sending them into a color ramp, and estimate on paper where typical values fall along that ramp.

7.15 Reading the default framebuffer with readPixels is unreliable ★

With preserveDrawingBuffer: false, reading the default framebuffer outside the rAF callback returns all black. For debugging, read from your own FBOs, which persist, or read immediately after drawing in the same frame.

---

8. Acceptance Criteria

8.1 Functionality

• [ ] Runs offline, produces no console errors, and keeps gl.getError() at 0 throughout
• [ ] Runs at ≥ 60 fps at 1080p with devicePixelRatio 1.8
• [ ] Clicking any celestial body focuses it with a smooth camera transition
• [ ] Every display toggle works, and “true distance scale” transitions with animation
• [ ] Time can run forward, run backward, pause, and jump back to the current moment
• [ ] Window resizing does not break the layout, and framebuffers are rebuilt correctly

8.2 Astronomical correctness (use these values for direct self-checks)

Set the date to 2026-08-07 and inspect Earth:

Additional checks:

If Earth’s longitude is off by 180°, you probably wrote M = L + ϖ instead of M = L − ϖ, or reversed the handedness of the coordinate mapping.

8.3 Visuals

• [ ] Deep space is nearly pure black, with a clearly visible Milky Way band in the star field
• [ ] The Sun’s cellular granulation network and sunspots are visible; it is not a plain white or yellow sphere
• [ ] Earth has blue oceans, white clouds, polar ice, night-side city lights, and a blue atmospheric limb
• [ ] Jupiter has multiple bands and a swirling Great Red Spot
• [ ] Saturn’s rings show the Cassini Division, the planet’s shadow on the rings, and the rings’ shadow on the planet
• [ ] Craters on rocky bodies vary in size and have raised rims rather than appearing as equally sized black dots
• [ ] Mars is rust red, not beige
• [ ] Bloom and lens flare are noticeable without obscuring the image

---

9. Recommended Implementation Order

Do not try to write everything at once. Build and run the project incrementally in this order, with a visible result at every step:

1. Foundation: WebGL2 context + matrix library + UV sphere + orbital camera + one solid-color sphere. Resolve the winding issue in Section 7.12 first.
2. Orbital mechanics: Kepler solver + positions of the nine planets + time controls. Validate with unlit, solid-color spheres. Match the astronomy against the Section 8.2 values before working on visuals.
3. Depth and framebuffers: HDR FBO + MSAA + dynamic near/far planes. The Section 7.1 trap appears here.
4. Basic lighting: point light + soft terminator. It should already begin to look planetary.
5. Post-processing: bloom + ACES + lens flare. Sections 7.3, 7.7, and 7.8 are concentrated here.
6. Procedural surfaces: write the noise library first, then implement one body type at a time. Take a close-up screenshot after each one.
7. Orbit lines + asteroid belts + star field. Section 7.2 applies here.
8. Ring systems and atmospheres, including both directions of shadowing.
9. Moons and comet. Section 7.6 applies here.
10. Interface and information panel.
11. Tuning: balance exposure, color saturation, and every intensity coefficient. This will take longer than expected; quantify with readPixels instead of relying on eyesight.

---

Appendix: An Easily Overlooked Overall Recommendation

The difficulty of this project lies not in any single technical element, but in the overall balance of exposure and color. The Sun, planet day sides, planet night sides, rings, star field, atmospheres, and bloom span several orders of magnitude in brightness. Overdriving any one of them ruins all the others.

A recommended approach is to render everything in neutral gray first and verify that all geometry and motion are correct. Then add color and tune brightness one element at a time, using readPixels after each addition to confirm that linear values remain in the expected range. Do not tune coefficients by feel while simultaneously writing shaders; that leads to endless back-and-forth adjustments.