Rendering convincing planets in real-time requires balancing physical accuracy with performance. This implementation combines Cook-Torrance PBR for realistic surface lighting, FBM noise for procedural terrain, Rayleigh-inspired atmospheric scattering, and subsurface scattering for translucent materials like water and lava. The result is planets that respond naturally to multiple dynamic light sources while maintaining high frame rates.
Sphere Reconstruction from Screen Space
Rather than rendering a 3D sphere mesh, we draw billboards (camera-facing quads) and reconstruct the sphere mathematically in the fragment shader. This approach avoids tessellation artifacts and provides perfect silhouettes at any zoom level. Given a point on the quad, we calculate whether it lies within the sphere's projected radius, then derive the surface normal from the implicit sphere equation.
// vUV is the quad coordinate relative to center [-1, 1] float d = length(vUV); float planetRadius = 0.6; // Normalized sphere radius on quad // Calculate Z coordinate on sphere surface float rotZSq = planetRadius * planetRadius - d * d; float rotZ = rotZSq > 0.0 ? sqrt(rotZSq) : 0.0; // 3D position on sphere vec3 spherePos = vec3(vUV, rotZ); // Derive surface normal (points outward from sphere center) vec3 baseN = d < planetRadius ? normalize(vec3(vUV, sqrt(planetRadius * planetRadius - d * d))) : vec3(0.0, 0.0, 1.0); // Mask for planet surface vs atmosphere float planetMask = 1.0 - smoothstep(planetRadius - 0.003, planetRadius + 0.001, d);
Planet Rotation with Rodrigues Formula
Each planet rotates around a unique axis derived from its index. We use Rodrigues' rotation formula to rotate the sphere position, giving each planet a distinct spin direction and speed. This rotation affects where terrain features appear and how they move over time.
// Generate unique rotation axis per planet from index float axisSeed = vIndex * 1.618033988749; // Golden ratio for distribution vec3 rotAxis = normalize(vec3( sin(axisSeed * 2.3 + 0.5), cos(axisSeed * 3.7 + 1.2), sin(axisSeed * 1.9 + 2.8) )); // Rotation speed varies per planet float rotSpeed = 0.03 + 0.05 * fract(axisSeed * 7.3); float rotAngle = uTime * rotSpeed; // Rodrigues' rotation formula: efficient axis-angle rotation // v_rot = v*cos(a) + (k x v)*sin(a) + k*(k.v)*(1-cos(a)) float cosA = cos(rotAngle); float sinA = sin(rotAngle); vec3 rotatedPos = spherePos * cosA + cross(rotAxis, spherePos) * sinA + rotAxis * dot(rotAxis, spherePos) * (1.0 - cosA);
Procedural Terrain with FBM Noise
Terrain is generated using Fractional Brownian Motion (FBM) - layered simplex noise with decreasing amplitude and increasing frequency per octave. We sample 5 octaves, with each octave adding finer detail. The noise is evaluated in 3D space using the rotated sphere position, creating seamless terrain that wraps around the planet.
float terrainHeight(vec3 pos, float scale) { float height = 0.0; float amplitude = 0.5; float frequency = scale; // uNoiseScale uniform // 5 octaves of simplex noise height += snoise3D(pos * frequency) * amplitude; frequency *= 2.0; amplitude *= 0.5; height += snoise3D(pos * frequency) * amplitude; frequency *= 2.0; amplitude *= 0.5; height += snoise3D(pos * frequency) * amplitude; frequency *= 2.0; amplitude *= 0.5; height += snoise3D(pos * frequency) * amplitude; frequency *= 2.0; amplitude *= 0.5; height += snoise3D(pos * frequency) * amplitude; return height; // Range approximately [-1, 1] } // Sample terrain at rotated position // Offset by planet index for unique terrain per planet vec3 planetOffset = vec3(vIndex * 10.0, vIndex * 7.3, vIndex * 5.1); vec3 terrainCoord = normalize(rotatedPos) + planetOffset; float heightC = terrainHeight(terrainCoord, uNoiseScale);
Normal Mapping from Height Field
To create the illusion of terrain geometry without actually displacing vertices, we compute a perturbed normal from the height field. By sampling the terrain at small offsets along the surface tangent and bitangent, we calculate height gradients that translate to normal perturbation.
float eps = 0.02; // Sample offset distance // Build tangent frame on sphere surface vec3 sphereUp = abs(rotatedPosNorm.y) < 0.999 ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0); vec3 sphereTangent = normalize(cross(sphereUp, rotatedPosNorm)); vec3 sphereBitangent = cross(rotatedPosNorm, sphereTangent); // Sample heights at tangent offsets float heightC = terrainHeight(terrainCoord, uNoiseScale); float heightX = terrainHeight(rotatedPosNorm + sphereTangent * eps + planetOffset, uNoiseScale); float heightY = terrainHeight(rotatedPosNorm + sphereBitangent * eps + planetOffset, uNoiseScale); // Convert height differences to normal perturbation float terrainGradX = -(heightX - heightC) / eps * uNormalStrength; float terrainGradY = -(heightY - heightC) / eps * uNormalStrength; // Perturb base normal using tangent frame vec3 N = normalize(baseN + tangent * terrainGradX + bitangent * terrainGradY);
Dual-Material System
Each planet supports two materials blended by terrain height. Material A (below sea level) handles water, lava, or low-lying terrain. Material B (above sea level) handles land, rock, or elevated terrain. Each material has independent base color, roughness, and subsurface scattering properties.
// Material parameters (uniforms) // Material A: below sea level (water, lava, etc.) uniform vec3 uMatABaseColor; // e.g., deep blue for ocean uniform float uMatARoughness; // e.g., 0.1 for smooth water uniform vec3 uMatASSSColor; // subsurface scatter color uniform float uMatASSSDistance; // scatter depth // Material B: above sea level (land, rock, etc.) uniform vec3 uMatBBaseColor; // e.g., green/brown for land uniform float uMatBRoughness; // e.g., 0.7 for rough terrain uniform vec3 uMatBSSSColor; uniform float uMatBSSSDistance; // Smooth blend based on terrain height vs sea level float materialMask = smoothstep(uSeaLevel - 0.02, uSeaLevel + 0.02, heightC); // Interpolate all material properties vec3 baseColor = mix(uMatABaseColor, uMatBBaseColor, materialMask); float roughness = mix(uMatARoughness, uMatBRoughness, materialMask); vec3 sssColor = mix(uMatASSSColor, uMatBSSSColor, materialMask); float sssDistance = mix(uMatASSSDistance, uMatBSSSDistance, materialMask);
Cook-Torrance PBR BRDF
Surface lighting uses the Cook-Torrance microfacet BRDF, the industry standard for physically-based rendering. This model combines three components: the GGX normal distribution function (D), the Smith geometry function (G), and the Fresnel-Schlick approximation (F). Together they create realistic specular highlights that respond correctly to roughness and viewing angle.
// GGX/Trowbridge-Reitz distribution // Models the statistical distribution of microfacet orientations float distributionGGX(float NdH, float roughness) { float a = roughness * roughness; float a2 = a * a; // Disney's roughness remapping float NdH2 = NdH * NdH; float denom = NdH2 * (a2 - 1.0) + 1.0; return a2 / (PI * denom * denom); } // The GGX distribution has a long tail, creating // realistic highlight falloff vs Blinn-Phong
// Schlick-GGX geometry term for a single direction // Accounts for microfacet self-shadowing float geometrySchlickGGX(float NdV, float roughness) { float r = roughness + 1.0; float k = (r * r) / 8.0; // Remapping for direct lighting return NdV / (NdV * (1.0 - k) + k); } // Smith's method: combine view and light shadowing float geometrySmith(float NdV, float NdL, float roughness) { float ggx1 = geometrySchlickGGX(NdV, roughness); float ggx2 = geometrySchlickGGX(NdL, roughness); return ggx1 * ggx2; // Joint masking-shadowing }
// Fresnel effect: reflectance increases at grazing angles vec3 fresnelSchlick(float cosTheta, vec3 F0) { return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0); } // Roughness-aware Fresnel for IBL vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness) { return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0); } // F0 for dielectrics (non-metals): typically 0.04 vec3 F0 = vec3(0.04);
Combining the BRDF
The complete Cook-Torrance BRDF multiplies D, G, and F terms together, divided by a normalization factor. The result gives us a specular contribution that we add to a diffuse term weighted by how much light wasn't reflected (energy conservation).
vec4 cookTorranceBRDF(vec3 N, vec3 V, vec3 L, float roughness, vec3 F0) { vec3 H = normalize(V + L); // Half-vector // Dot products (clamped to avoid artifacts) float NdV = max(dot(N, V), 0.001); float NdL = max(dot(N, L), 0.0); float NdH = max(dot(N, H), 0.0); float HdV = max(dot(H, V), 0.0); // BRDF terms float D = distributionGGX(NdH, roughness); float G = geometrySmith(NdV, NdL, roughness); vec3 F = fresnelSchlick(HdV, F0); // Cook-Torrance specular vec3 numerator = D * G * F; float denominator = 4.0 * NdV * NdL + 0.0001; // Avoid div by zero vec3 specular = numerator / denominator; // Return specular term and average Fresnel for diffuse weighting float avgF = (F.r + F.g + F.b) / 3.0; return vec4(specular * NdL, avgF); }
Multi-Light Accumulation
The shader supports 3 dynamic point lights plus an interactive mouse light. Each light contributes independently with its own color, intensity, and attenuation. We transform normals to world space to ensure correct lighting regardless of camera orientation.
// Transform normal to world space using camera basis vec3 cameraForward = vec3(sinRotY * cosRotX, -sinRotX, cosRotY * cosRotX); vec3 cameraRight = vec3(cosRotY, 0.0, -sinRotY); vec3 cameraUp = cross(cameraForward, cameraRight); vec3 N_world = normalize(cameraRight * N.x - cameraUp * N.y - cameraForward * N.z); vec3 V_world = -cameraForward; // View direction // World position of this surface point float radiusWorld = vOriginalRadius * worldScale; vec3 surfaceWorldPos = planetWorldPos + N_world * radiusWorld;
// Accumulate diffuse and specular from each light vec3 totalDiffuse = vec3(0.0); vec3 totalSpecular = vec3(0.0); // Light 0 contribution vec3 L0 = normalize(uLight0WorldPos - surfaceWorldPos); float dist0 = length(uLight0WorldPos - surfaceWorldPos); float atten0 = 1.0 / (1.0 + uLight0Atten * 200.0 * dist0 * dist0); // Quadratic float NdL0 = max(dot(N_world, L0), 0.0); vec4 pbr0 = cookTorranceBRDF(N_world, V_world, L0, roughness, F0); totalSpecular += uLightColor0 * pbr0.xyz * atten0 * uLight0Intensity; // Energy-conserving diffuse: weight by (1 - Fresnel) float diffuseWeight0 = (1.0 - pbr0.w) / PI; totalDiffuse += uLightColor0 * NdL0 * atten0 * diffuseWeight0 * uLight0Intensity; // Repeat for lights 1, 2, and mouse light...
Subsurface Scattering Approximation
For translucent materials like water or ice, light penetrates the surface and scatters internally before exiting. We approximate this with a wrap lighting term and backlight contribution, creating soft illumination that bleeds around the terminator edge.
// Wrap lighting: extends lit area past terminator // sssWrap of 0.5 means 50% wrap around the terminator float sssWrap = 0.5; float sssNdL0 = max(0.0, (dot(N_world, L0) + sssWrap) / (1.0 + sssWrap)); float sssNdL1 = max(0.0, (dot(N_world, L1) + sssWrap) / (1.0 + sssWrap)); float sssNdL2 = max(0.0, (dot(N_world, L2) + sssWrap) / (1.0 + sssWrap)); // Backlight: light passing through thin areas float backLight0 = pow(max(0.0, -dot(N_world, L0)), 2.0); float backLight1 = pow(max(0.0, -dot(N_world, L1)), 2.0); float backLight2 = pow(max(0.0, -dot(N_world, L2)), 2.0); // Combine SSS contribution vec3 sss = sssColor * sssDistance * ( uLightColor0 * (sssNdL0 + backLight0 * 0.5) * atten0 * uLight0Intensity + uLightColor1 * (sssNdL1 + backLight1 * 0.5) * atten1 * uLight1Intensity + uLightColor2 * (sssNdL2 + backLight2 * 0.5) * atten2 * uLight2Intensity );
Physically-Based Atmospheric Scattering
The atmosphere extends beyond the planet's solid surface and requires special handling. We ray-march through the atmospheric shell, computing optical depth based on density that falls off exponentially with altitude. This drives wavelength-dependent Rayleigh scattering for realistic sky coloring.
// Atmosphere extends beyond planet radius float atmosRadius = planetRadius + 0.4 * uAtmosThickness; float atmosThickness = atmosRadius - planetRadius; // Rayleigh scattering coefficients (wavelength-dependent) vec3 beta = vec3(uScatterR, uScatterG, uScatterB); float densityFalloff = uAtmosPower; // Exponential decay rate
float opticalDepth = 0.0; if (d < atmosRadius) { float atmosZ = sqrt(max(0.0, atmosRadius * atmosRadius - d * d)); if (d < planetRadius) { // Ray enters atmosphere, hits planet surface float planetZ = sqrt(max(0.0, planetRadius * planetRadius - d * d)); float pathLength = atmosZ - planetZ; // Integrate density along path (4 samples) for (float i = 0.0; i < 4.0; i++) { float t = (i + 0.5) / 4.0; float sampleZ = planetZ + t * (atmosZ - planetZ); float sampleR = sqrt(d * d + sampleZ * sampleZ); float sampleHeight01 = (sampleR - planetRadius) / atmosThickness; float sampleDensity = exp(-sampleHeight01 * densityFalloff) * (1.0 - sampleHeight01); opticalDepth += sampleDensity * (pathLength / 4.0); } } else { // Ray passes through atmosphere only (limb view) float pathLength = 2.0 * atmosZ; // 6 samples for longer path for (float i = 0.0; i < 6.0; i++) { float t = (i + 0.5) / 6.0; float sampleZ = atmosZ * (1.0 - 2.0 * t); float sampleR = sqrt(d * d + sampleZ * sampleZ); float sampleHeight01 = clamp((sampleR - planetRadius) / atmosThickness, 0.0, 1.0); float sampleDensity = exp(-sampleHeight01 * densityFalloff) * (1.0 - sampleHeight01); opticalDepth += sampleDensity * (pathLength / 6.0); } } } opticalDepth *= uScatterScale;
Rayleigh Color Shift
Rayleigh scattering depends on wavelength^-4, meaning blue light scatters much more than red. We compute transmittance through the atmosphere using Beer-Lambert law, then derive the scattered color. This creates blue skies where atmosphere is thin and warm sunset colors where the optical path is long.
// Beer-Lambert transmittance through atmosphere vec3 transmittance = exp(-beta * opticalDepth * uSunsetStrength); // Normalize transmittance to get color shift vec3 transmitColor = transmittance; float maxT = max(transmittance.r, max(transmittance.g, transmittance.b)); if (maxT > 0.001) { transmitColor = transmittance / maxT; } else { // Deep atmosphere: warm sunset colors transmitColor = vec3(1.0, 0.3, 0.1); } // Blue scatter for thin atmosphere vec3 blueScatter = vec3(1.0) - exp(-beta * 8.0); // Blend between blue scatter and transmit color based on depth float depthFactor = 1.0 - exp(-opticalDepth * 2.0 * uSunsetStrength); vec3 scatterColor = mix(blueScatter, transmitColor, depthFactor); // Final atmosphere contribution vec3 atmosColor = incomingLight * scatterColor * atmosDensity * uAtmosIntensity;
Atmosphere Self-Shadowing
The planet casts a shadow on its own atmosphere. We compute where each light's rays would intersect the planet, creating shadow cones that darken the atmosphere on the planet's far side from each light source.
// For pixels in the atmosphere (outside planet surface) float atmosShadow0 = 1.0; // Start fully lit if (d >= planetRadius && d < atmosRadius) { float r = planetRadius; float penumbraWidth = r * 0.2; // Soft shadow edge // 2D light direction from planet center vec2 D0 = normalize(uLight0ScreenPos - vCenter); // Distance from light ray to current point float perpDist0 = abs(vUV.x * D0.y - vUV.y * D0.x); float alongRay0 = dot(vUV, D0); // In shadow if: behind planet center AND close to light ray if (alongRay0 < 0.0 && perpDist0 < r + penumbraWidth) { float behindPlanet = -alongRay0 / r; float shadowStart = smoothstep(0.0, 1.0, behindPlanet); float shadowCore = smoothstep(r + penumbraWidth, r - penumbraWidth, perpDist0); atmosShadow0 = 1.0 - shadowCore * shadowStart; } atmosShadow0 = max(atmosShadow0, 0.05); // Never fully black }
Final Composition
The final color combines diffuse lighting, specular highlights, subsurface scattering, ambient light, Fresnel rim, and atmospheric contribution. We also apply depth fog based on camera distance for integration with the background.
vec3 col = vec3(0.0); // Surface lighting (planet only) col += baseColor * totalDiffuse * 2.5 * planetMask; // Diffuse col += totalSpecular * 3.0 * planetMask; // Specular col += sss * planetMask; // Subsurface // Ambient from average light color vec3 lightEnvColor = (uLightColor0 * atten0 * uLight0Intensity + uLightColor1 * atten1 * uLight1Intensity + uLightColor2 * atten2 * uLight2Intensity) / max(totalAttenuation, 0.001); col += baseColor * lightEnvColor * uAmbientIntensity * 0.1 * planetMask; // Fresnel rim highlighting vec3 fresnelRim = fresnelSchlickRoughness(NdV, F0, roughness) * planetMask; float rimIntensity = mix(0.04, 0.15, 1.0 - roughness); col += lightEnvColor * fresnelRim * totalAttenuation * rimIntensity; // Atmosphere blend over terrain at limb float limbFactor = pow(1.0 - NdV, 2.0); // Stronger at edges float terrainAtmosBlend = atmosAlpha * (0.3 + limbFactor * 0.7) * planetMask; col = mix(col, atmosColor, terrainAtmosBlend); // Depth fog for distance fade float distToCamera = length(surfaceWorldPos - cameraPos); float fogAmount = 1.0 - exp(-distToCamera * uFogIntensity); col = mix(col, fogColor, fogAmount * planetMask); // Atmosphere glow outside planet surface col += foggedAtmosColor * (1.0 - planetMask);
Performance Considerations
- Billboard rendering: No mesh tessellation needed, perfect silhouettes at any distance
- Noise evaluation: 5 FBM octaves is expensive; consider reducing for mobile or distant planets
- Atmosphere ray marching: Fixed sample counts (4-6) balance quality with performance
- Multiple lights: Cost scales linearly with light count; 3 lights + mouse is reasonable
- Early discard: Pixels outside the atmosphere shell are discarded early
- Depth fade: Distant planets can reduce quality without visible difference
Parameter Reference
// Oceanic planet (Type A) - water world planetParamsA = { noiseScale: 2.5, // Terrain frequency seaLevel: 0.0, // Water/land boundary normalStrength: 0.4, // Terrain bumpiness matABaseColor: [0.1, 0.3, 0.5], // Deep ocean matARoughness: 0.15, // Smooth water matASSSColor: [0.2, 0.5, 0.6], // Subsurface tint matASSSDistance: 0.3, matBBaseColor: [0.3, 0.45, 0.2], // Green land matBRoughness: 0.8, // Rough terrain atmosIntensity: 1.2, atmosThickness: 0.5, scatterR: 0.4, scatterG: 0.7, scatterB: 1.0 // Blue sky }; // Volcanic planet (Type B) - lava world planetParamsB = { noiseScale: 3.0, seaLevel: -0.2, matABaseColor: [0.8, 0.2, 0.05], // Molten lava matARoughness: 0.3, matBBaseColor: [0.15, 0.1, 0.08], // Dark rock atmosIntensity: 0.6, scatterR: 1.0, scatterG: 0.4, scatterB: 0.2 // Orange haze };