Shadertoy

float sdPlane( vec3 p ) {
	return p.y;
}

float sdSphere( in vec3 p, in vec4 s ) {
    return length(p-s.xyz) - s.w;
}

float sdBox( vec3 p, vec3 b )
{
    vec3 d = abs(p) - b;
    return min(max(d.x,max(d.y,d.z)),0.0) + length(max(d,0.0));
}

float length2( vec2 p )
{
	return sqrt( p.x*p.x + p.y*p.y );
}

float length8( vec2 p )
{
	p = p*p; p = p*p; p = p*p;
	return pow( p.x + p.y, 1.0/8.0 );
}

float sdTorus( vec3 p, vec2 t )
{
    return length( vec2(length(p.xz)-t.x,p.y) )-t.y;
}

float sdTorus82( vec3 p, vec2 t )
{
    vec2 q = vec2(length2(p.xz)-t.x,p.y);
    return length8(q)-t.y;
}

vec3 background() {
    return vec3(0.7, 0.9, 1.0);
}

vec2 opU( vec2 d1, vec2 d2 ) {
	return (d1.x<d2.x) ? d1 : d2;
}

float opS( float d1, float d2 )
{
    return max(-d2,d1);
}

void pR(inout vec2 p, float a) {
	p = cos(a)*p + sin(a)*vec2(p.y, -p.x);
}

void pR45(inout vec2 p) {
	p = (p + vec2(p.y, -p.x))*sqrt(0.5);
}

float fOpUnionSoft(float a, float b, float r) {
	float e = max(r - abs(a - b), 0.0);
	return min(a, b) - e*e*0.25/r;
}

vec2 wheel(vec3 p) {
    vec3 origin = p-vec3(-2.0,0.5,2.0);
    pR(origin.yz, 1.57);
    float d = sdTorus(origin, vec2(0.25,0.1));
    d = fOpUnionSoft(d, sdTorus82(origin, vec2(0.25,0.1)),0.5);
    return vec2(d, 20);
}

vec2 mapOpaque( vec3 p )
{
    vec2 dm = vec2(sdSphere(p, vec4(0.0, 0.25, 2.0, 0.25)), 40);
    dm = opU(dm, vec2(sdPlane(p), 1));
    dm = opU(dm, vec2(sdBox(p-vec3( -1.0,0.25,2.0), vec3(0.25)), 30));
    dm = opU(dm, wheel(p));
    return dm;
}

vec2 intersectOpaque( in vec3 ro, in vec3 rd,  float mindist,  float maxdist, out vec4 matInfo )
{
    vec2 res = vec2(-1.0);
    
    /*
    // bounding volume
    float tp1 = (0.0-ro.y)/rd.y; if( tp1>0.0 ) maxdist = min( maxdist, tp1 );
    float tp2 = (1.6-ro.y)/rd.y; if( tp2>0.0 ) { if( ro.y>1.6 ) mindist = max( mindist, tp2 );
                                                 else           maxdist = min( maxdist, tp2 ); }
    */

    float t = mindist;
    for( int i=0; i<64; i++ )
    {
        vec3 p = ro + t*rd;
        vec2 h = mapOpaque( p );
        res = vec2(t,h.y);

        if( h.x<(0.001*t) ||  t>maxdist ) break;
        
        t += h.x*0.9;
    }
	return res;
}

vec3 calcNormalOpaque( in vec3 pos, in float eps )
{
    vec2 e = vec2(1.0,-1.0)*0.5773*eps;
    return normalize( e.xyy*mapOpaque( pos + e.xyy ).x + 
					  e.yyx*mapOpaque( pos + e.yyx ).x + 
					  e.yxy*mapOpaque( pos + e.yxy ).x + 
					  e.xxx*mapOpaque( pos + e.xxx ).x );
}

float calcAO( in vec3 pos, in vec3 nor )
{
	float occ = 0.0;
    float sca = 1.0;
    for( int i=0; i<5; i++ )
    {
        float hr = 0.01 + 0.12*float(i)/4.0;
        vec3 aopos =  nor * hr + pos;
        float dd = mapOpaque( aopos ).x;
        occ += -(dd-hr)*sca;
        sca *= 0.95;
    }
    return clamp( 1.0 - 3.0*occ, 0.0, 1.0 );    
}

float softshadow( in vec3 ro, in vec3 rd, in float mint, in float tmax )
{
	float res = 1.0;
    float t = mint;
    for( int i=0; i<16; i++ )
    {
		float h = mapOpaque( ro + rd*t ).x;
        res = min( res, 8.0*h/t );
        t += clamp( h, 0.02, 0.10 );
        if( h<0.001 || t>tmax ) break;
    }
    return clamp( res, 0.0, 1.0 );
}

vec3 shadeOpaque( in vec3 ro, in vec3 rd, in float t, in float m, in vec4 matInfo )
{
    float eps = 0.0005;
    
    vec3 pos = ro + t*rd;
    vec3 nor = calcNormalOpaque( pos, eps );
    vec3 ref = reflect( rd, nor );
    
    vec3 col = 0.45 + 0.35*sin( vec3(0.05,0.08,0.10)*(m-1.0) );

    float occ = calcAO( pos, nor );
    vec3  lig = normalize( vec3(-0.4, 0.7, -0.6) );
	float amb = clamp( 0.5+0.5*nor.y, 0.0, 1.0 );
    float dif = clamp( dot( nor, lig ), 0.0, 1.0 );
    float bac = clamp( dot( nor, normalize(vec3(-lig.x,0.0,-lig.z))), 0.0, 1.0 )*clamp( 1.0-pos.y,0.0,1.0);
    float dom = smoothstep( -0.1, 0.1, ref.y );
    float fre = pow( clamp(1.0+dot(nor,rd),0.0,1.0), 2.0 );
    float spe = pow(clamp( dot( ref, lig ), 0.0, 1.0 ),16.0);
    
    dif *= softshadow( pos, lig, 0.02, 2.5 );
    dom *= softshadow( pos, ref, 0.02, 2.5 );
    
	vec3 lin = vec3(0.0);
    lin += 1.30*dif*vec3(1.00,0.80,0.55);
    lin += 2.00*spe*vec3(1.00,0.90,0.70)*dif;
    lin += 0.40*amb*vec3(0.40,0.60,1.00)*occ;
    lin += 0.50*dom*vec3(0.40,0.60,1.00)*occ;
    lin += 0.50*bac*vec3(0.25,0.25,0.25)*occ;
    lin += 0.25*fre*vec3(1.00,1.00,1.00)*occ;
    
    col = col*lin;
    
    // Some sort of color bleed?!
    col = mix( col, vec3(0.8,0.9,1.0), 1.0-exp( -0.0005*t*t*t ) );

    return col;
}

vec3 render(in vec3 ro, in vec3 rd) {
    vec3 col = background();

    float mindist = 1.0;
    float maxdist = 20.0;

    vec4 matInfo;
    vec2 tm = intersectOpaque( ro, rd, mindist, maxdist, matInfo );
    if( tm.y>-0.5 && tm.x < maxdist )
    {
        col = shadeOpaque( ro, rd, tm.x, tm.y, matInfo );
        maxdist = tm.x;
    }
    
    // Gamma correct
    col = pow(col, vec3(0.45));

    return clamp(col, 0.0, 1.0);
}

mat3 setCamera( in vec3 ro, in vec3 ta, float cr )
{
	vec3 cw = normalize(ta-ro);
	vec3 cp = vec3(sin(cr), cos(cr),0.0);
	vec3 cu = normalize( cross(cw,cp) );
	vec3 cv = normalize( cross(cu,cw) );
    return mat3( cu, cv, cw );
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
    vec2 p = (-iResolution.xy + 2.0*fragCoord)/iResolution.y;

    // camera
    float time = 15.0 + iGlobalTime;
    vec3 ro = vec3( -0.5+3.5*cos(0.1*time), 1.0, 0.5 + 4.0*sin(0.1*time) );
    vec3 ta = vec3( -0.5, -0.4, 0.5 );
    // camera-to-world transformation
    mat3 ca = setCamera( ro, ta, 0.0 );
    // ray direction
    vec3 rd = ca * normalize( vec3(p.xy,2.0) );

    vec3 col = render(ro, rd);
	fragColor = vec4(col,1.0);
}