Lesson 6 of 10 · The Solar System in GLSL

Jupiter

Bands that run against each other. And a storm wider than Earth.

The code

precision highp float;
// The exact sphere: how far along the ray it is hit, or -1 if the ray misses
float sphere(vec3 ro, vec3 rd, float r){
  float b = dot(ro, rd), c = dot(ro, ro) - r*r, h = b*b - c;
  return h < 0.0 ? -1.0 : -b - sqrt(h);
}
// WoGL close: knee at 0.70, sRGB, chroma 1.3
vec3 wogl(vec3 c){
  float l = dot(c, vec3(0.2126, 0.7152, 0.0722)), k = max(l - 0.7, 0.0);
  c *= l > 0.0 ? (l - k + k/(1.0 + k))/l : 0.0;
  c = pow(clamp(c, 0.0, 1.0), vec3(0.4545));
  l = dot(c, vec3(0.2126, 0.7152, 0.0722));
  return clamp(mix(vec3(l), c, 1.3), 0.0, 1.0);
}
// Noise in 3D: the usual fb with one more dimension. On a globe it leaves no seam.
float h3(vec3 p){ return h2(p.xy + p.z*vec2(37.0, 17.0)); }
float n3(vec3 p){
  vec3 i = floor(p), f = fract(p);
  f = f*f*(3.0 - 2.0*f);
  return mix(mix(mix(h3(i), h3(i + vec3(1,0,0)), f.x), mix(h3(i + vec3(0,1,0)), h3(i + vec3(1,1,0)), f.x), f.y),
             mix(mix(h3(i + vec3(0,0,1)), h3(i + vec3(1,0,1)), f.x), mix(h3(i + vec3(0,1,1)), h3(i + vec3(1,1,1)), f.x), f.y), f.z);
}
float f3(vec3 p){
  float v = 0.0, a = 0.5;
  for (int i = 0; i < 5; i++){ v += a*n3(p); p = p*2.03 + vec3(1.7, 9.2, 4.1); a *= 0.5; }
  return v;
}
// A 2D turn (to spin the globe about its axis)
mat2 turn(float a){ float c = cos(a), s = sin(a); return mat2(c, -s, s, c); }
// The background: black, with a star now and then
vec3 sky(vec2 f){
  float s = h2(floor(f));
  return vec3(0.0002, 0.0003, 0.0006) + vec3(0.7, 0.7, 0.8)*step(0.9982, s)*h2(floor(f) + 3.1);
}
float bands = 7.0;     // pairs of zones and belts
float wind = 0.6;      // how hard each band runs against the next
float storm = 1.0;     // the Great Red Spot's spin
// The clouds at a point of the globe: bands by latitude, with turbulence
vec3 cloudsAt(vec3 q){
  float lat = asin(clamp(q.y, -1.0, 1.0));
  vec3 w = q; w.xz = turn(T*0.03*wind*sin(lat*bands*2.0))*w.xz;   // shear: each band turns at its own speed
  float z = sin(lat*bands*2.0 + f3(w*vec3(3.0, 9.0, 3.0) + f3(w*5.0)*1.5)*2.6);
  vec3 zone = vec3(0.92, 0.85, 0.72), belt = vec3(0.55, 0.36, 0.22);
  return mix(belt, zone, smoothstep(-0.3, 0.3, z))*(0.85 + 0.3*f3(w*vec3(6.0, 20.0, 6.0)))*0.55;
}
void main(){
  vec2 p = (gl_FragCoord.xy*2.0 - R)/min(R.x, R.y);
  vec3 ro = vec3(0.0, 0.0, 3.2), rd = normalize(vec3(p, -2.4));   // the camera looks at the center
  vec3 sun = normalize(vec3(M.x*4.0 - 3.8, M.y*3.0 - 1.0, 1.0));   // you move the sun: it is the pointer
  vec3 col = sky(gl_FragCoord.xy);
  // One turn every 9 h 56 min: the globe flattens 6.5%. Shrink space in y, trace a sphere.
  vec3 k = vec3(1.0, 1.0/0.935, 1.0);
  vec3 ros = ro*k, rds = normalize(rd*k);
  float t = sphere(ros, rds, 1.0);
  if (t > 0.0){
    vec3 ps = ros + rds*t, n = normalize(ps*k);
    vec3 q = ps; q.xz = turn(T*0.05)*q.xz;
    // The vortex: near the spot's center the point is turned; far away, not
    vec3 c0 = vec3(cos(-0.39)*sin(0.25), sin(-0.39), cos(-0.39)*cos(0.25));   // 22° south
    vec3 e1 = normalize(cross(vec3(0.0, 1.0, 0.0), c0)), e2 = cross(c0, e1);
    vec2 d = vec2(dot(q - c0, e1)/0.18, dot(q - c0, e2)/0.09);              // twice as wide as it is tall
    float r = length(d), inside = 1.0 - smoothstep(0.6, 1.0, r);
    d = turn(storm*(1.0 - smoothstep(0.0, 1.0, r))*(3.0 + T*0.5))*d;
    vec3 qs = mix(q, c0 + e1*d.x*0.18 + e2*d.y*0.09, step(r, 1.0));
    vec3 albedo = cloudsAt(qs);
    albedo = mix(albedo, vec3(0.62, 0.24, 0.12)*(0.8 + 0.4*f3(qs*12.0)), inside*0.85);
    col = albedo*max((dot(n, sun) + 0.1)/1.1, 0.0)*3.0;
  }
  gl_FragColor = vec4(wogl(col), 1.0);
}