800 600 1 52
// Algorithm
#define MAX_MARCHING_STEPS 200
#define MIN_DIST 0.5
#define MAX_DIST 8.0
#define EPSILON 0.01
#define MARCHING_STEP 0.25

// Times (part 1)
/*#define T0 0.0
#define T1 10.0
#define T2 20.0
#define T3 30.0
#define T4 35.0
#define T5 45.0
#define T6 50.0*/
#define T0 0.0
#define T1 1.0
#define T2 2.0
#define T3 3.0
#define T4 3.5
#define T5 4.0
#define T6 4.5
// part 2 starts here
#define T7 5.0

// Animation constants (part 2)
#define HEART_BEAT_T0 T6
#define HEART_BEAT_T1 T6 + 15.0
#define HEART_ROTATE_T0 T6 + 10.0
#define HEART_ROTATE_T1 T6 + 15.0
#define SPIKE_T0 HEART_ROTATE_T1 + 1.0
#define SPIKE_T1 SPIKE_T0 + 0.5
#define SPIKE_HEIGHT 2.0
#define SPIKE_YFINAL -2.4
#define SPIKE_UNION_Y 0.75
#define BLUE_COLOR_AHEAD 0.55
#define BLUE_COLOR_SPAN 10.0
#define BLUE_RADIUS 7.0
#define CAMERA_ANGLE0 25.
#define CAMERA_Y1 1.3
#define CAMERA_OMEGA 0.2
#define VIEW_ANGLE0 75.0
#define VIEW_ANGLE_DELTA 28.0
#define SPIKE_Y0 SPIKE_YFINAL - SPIKE_HEIGHT

// Operations
#define OP_U(X, Y) min(X, Y)
#define OP_S(X, Y) max( -Y, X )
#define OP_I(X, Y) max( X, Y )
#define SMOOTHTR(X0, X1) smoothstep(X0, X1, iGlobalTime)
mat3 Rot4Z(float a) { // angle in radians
    float c = cos( a );
    float s = sin( a );
    return mat3(
        c,-s, 0,
        s, c, 0,
        0, 0, 1
    );
}

float timeTransition(float t0, float t1) {
    return clamp((iGlobalTime - t0)/(t1-t0), 0., 1.);
}

// SDF primitives
float cubeSDF(vec3 p, vec3 h) {
    vec3 d = abs(p) - h;
    float insideDistance = min(max(d.x, max(d.y, d.z)), 0.0);
    float outsideDistance = length(max(d, 0.0));
    return insideDistance + outsideDistance;
}

float cylinderSDF( vec3 p, vec2 c, float r ) {
    return length(p.xy - c.xy) - r;
}
float vcylinderSDF( vec3 p, vec2 c, float r ) {
    return length(p.xz - c.xy) - r;
}

float sphereSDF(vec3 p, float r) {
    return length(p) - r;
}

float planeSDF( vec3 p, vec4 n ) {
    return dot( p, n.xyz ) + n.w;
}

float coneSDF( vec3 p, vec2 c ) {
    float q = length(p.xz);
    return dot(c,vec2(q,p.y));
}

float cappedCylinderSDF( vec3 p, vec2 h )
{
  vec2 d = abs(vec2(length(p.xz),p.y)) - h;
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

// Figures
float simpleEmblemSDF(vec3 p) { 
    // Heart
    vec3 q = p;
    q.z = q.z * (3. - q.y/15.);
    q.y = - 1.2*q.y - abs(q.x)*sqrt((20.-abs(q.x))/15.);
    float sphere = sphereSDF(q, 0.3);
    q = p;
    q.x = abs(q.x);
    float cylinder = cylinderSDF(q, vec2(0.25, 0.1), 0.18);
    float heart = OP_S(sphere, cylinder);
    
    // Spike
    vec3 spikeQ = p;
    spikeQ.y -= 0.5;
    q = abs(spikeQ);
    q.y = q.y*(1. + 2.*q.x);
    q.z = q.z*(3.5 + 0.2*q.x);
    sphere = cubeSDF(q, vec3(4., 4., 4.));
    float plane = planeSDF(q, normalize(
        vec4(-1., -1., -1., 0.4)));
    float diamond = OP_S(sphere, plane);
    q = abs(spikeQ);
    cylinder = cylinderSDF(q, vec2(0.18, 0.18), 0.13);
    float spike = OP_S(diamond, cylinder);
    
    return OP_U(heart, spike);
}

float altarSDF(vec3 p) {
    // Altar
    vec3 q = p;
    q.y -= 4.5;
    float cone = coneSDF(q, normalize(vec2(1., 0.5)));
    q = p;
    float p1 = planeSDF(q, vec4(0., -1., 0., 0.));
    float res = OP_S(cone, p1);
    
    // Stairs
    q.z -= 2.5;
    q.y += 0.8;
    float stair = cubeSDF(q,
       vec3(0.85, 0.8, 0.8));
    res = OP_U(res, stair);
    q.y = p.y - 0.6;
    for (int i = 0; i != 5; ++i) {
        stair = cubeSDF(q, vec3(0.95, 0.8, 0.2));
        res = OP_S(res, stair);
        q.y += 0.2;
        q.z -= 0.2;
    };
        
    // Stair sides
    q = p;
    q.x = abs(p.x);
    q.y += 0.7;
    q.z -= 2.8;
    q.x -= 0.8;
 	float stairSides = cubeSDF(q, vec3(0.1, 0.85, 1.));
    float stairPlane = planeSDF(q, vec4(0., -1., -1., 0.35));
    stairSides = OP_S(stairSides, stairPlane);
    res = OP_U(res, stairSides);
    
    // Stair blocks
    q = p;
    q.x = abs(q.x);
    q.z -= 2.1;
    q.x -= 0.8;
    q.y -= 0.2;
    float stairBlock = cubeSDF(q, vec3(0.2, 0.08, 0.25));
    res = OP_U(res, stairBlock);
    
    // Emblems
    q = p;
	q.x = abs(q.x);
    q.x -= 0.8;
    q.y -= 0.8;
    q.z -= 2.2;
	res = OP_U(res, simpleEmblemSDF(q));
    
    // Edge
    q = p;
    float cyl1 = cappedCylinderSDF(q, vec2(2.35, 0.15));
    float cyl2 = vcylinderSDF(q, vec2(0.0, 0.0), 2.0);
    q.z -= 1.4;
    float entry = cubeSDF(q, vec3(0.85, 1., 1.));
    float edge = OP_S(cyl1, cyl2);
    edge = OP_S(edge, entry);
    res = OP_U(res, edge);
   
    return res;
}

// Part 2
vec3 moveSpike( vec3 p ) {
    return vec3(
        p.x,
        p.y + SPIKE_Y0 + SPIKE_HEIGHT*
            smoothstep(SPIKE_T0, SPIKE_T1, iGlobalTime),
        p.z
    );
}

vec3 moveHeart (vec3 p) {
    //p.y -= 4.;
    float transition = 1.-timeTransition(
        HEART_BEAT_T0, HEART_BEAT_T1);
    float amplitude = mix(0., 0.14, transition);
    float omega = mix(2., 4., transition);
    float heartbeat = 1. +
        amplitude*pow(abs(sin(omega*iGlobalTime)), 4.);
    vec3 res = vec3(
        p.x,
        p.y,
        p.z
    ) * heartbeat;
    float angle = 180.0*(1.-smoothstep(
        HEART_ROTATE_T0, HEART_ROTATE_T1, iGlobalTime));
    mat3 rot = Rot4Z(radians(angle));
    return rot*res;
}

// Figures
float heartSDF( vec3 p ) {
    p = moveHeart(p);
    vec3 q = p;
    q.z = q.z * (3. - q.y/15.);
    q.y = 0.2 - 1.2*q.y - abs(q.x)*sqrt((20.-abs(q.x))/15.);
    float sphere = sphereSDF(q, 1.5);
    q = p;
    q.x = abs(q.x);
    float cylinder = cylinderSDF(q, vec2(0.6, 0.4), 0.42);
    return OP_S(sphere, cylinder);
}

float spikeSDF( vec3 samplePoint ) {
    samplePoint = moveSpike(samplePoint);
    vec3 q = abs(samplePoint);
    q.y = q.y*(1. + 2.*q.x);
    q.z = q.z*(3.5 + 0.2*q.x);
    float sph = cubeSDF(q, vec3(4., 4., 4.));
    vec3 norm = normalize(vec3(-1., -1., -1.));
    float plane = planeSDF(q, vec4(norm, 1.));
    float diamond = OP_S(sph, plane);
    q = samplePoint;
    q = abs(q);
    float cylinder = cylinderSDF(q, vec2(0.6, 0.6), 0.42);
    return OP_S(diamond, cylinder);
}

float sceneSDF(vec3 p) {
    if (iGlobalTime < T6)
        return altarSDF(p);
    else
        return OP_U(heartSDF(p), spikeSDF(p));
}

// Marching algorithm
float shortestDistanceToSurface(vec3 eye, vec3 marchingDirection, float start, float end) {
    float depth = start;
    vec3 p;
    float dist = 2.*EPSILON;
    int i = 0;
    for (; 
        i < MAX_MARCHING_STEPS &&
        abs(dist) > EPSILON &&
        depth < end; ++i) {
        depth += MARCHING_STEP*dist;
        dist = sceneSDF(eye + depth * marchingDirection);
    }
    if (depth >= end || i >= MAX_MARCHING_STEPS) {
		return end;
    }
    return depth;
}
            
// Get ray direction given frag coord
vec3 rayDirection(float fieldOfView, vec2 size, vec2 fragCoord) {
    vec2 xy = fragCoord - size / 2.0;
    float z = size.y / tan(radians(fieldOfView) / 2.0);
    return normalize(vec3(xy, -z));
}

// Transformation matrix from view to world
mat4 viewMatrix(vec3 eye, vec3 center, vec3 up) {
    vec3 f = normalize(center - eye);
    vec3 s = normalize(cross(f, up));
    vec3 u = cross(s, f);
    return mat4(
        vec4(s, 0.0),
        vec4(u, 0.0),
        vec4(-f, 0.0),
        vec4(0.0, 0.0, 0.0, 1)
    );
}

// Get normal
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        sceneSDF(vec3(p.x + EPSILON, p.y, p.z)) - sceneSDF(vec3(p.x - EPSILON, p.y, p.z)),
        sceneSDF(vec3(p.x, p.y + EPSILON, p.z)) - sceneSDF(vec3(p.x, p.y - EPSILON, p.z)),
        sceneSDF(vec3(p.x, p.y, p.z  + EPSILON)) - sceneSDF(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

/**
 * Lighting contribution of a single point light source via Phong illumination.
 *
 * The vec3 returned is the RGB color of the light's contribution.
 *
 * k_a: Ambient color
 * k_d: Diffuse color
 * k_s: Specular color
 * alpha: Shininess coefficient
 * p: position of point being lit
 * eye: the position of the camera
 * lightPos: the position of the light
 * lightIntensity: color/intensity of the light
 *
 * See https://en.wikipedia.org/wiki/Phong_reflection_model#Description
 */
vec3 phongContribForLight(vec3 k_d, vec3 k_s, float alpha, vec3 p, vec3 eye,
                          vec3 lightPos, vec3 lightIntensity) {

    vec3 N = estimateNormal(p);
    vec3 L = normalize(lightPos - p);
    vec3 V = normalize(eye - p);
    vec3 R = normalize(reflect(-L, N));
    
    float dotLN = dot(L, N);
    float dotRV = dot(R, V);
    
    if (dotLN < 0.0) {
        // Light not visible from this point on the surface
        return vec3(0.0, 0.0, 0.0);
    }
    
    if (dotRV < 0.0) {
        // Light reflection in opposite direction as viewer, apply only diffuse
        // component
        return lightIntensity * (k_d * dotLN);
    }
    return lightIntensity * (k_d * dotLN + k_s * pow(dotRV, alpha));
}

/**
 * Lighting via Phong illumination.
 *
 * The vec3 returned is the RGB color of that point after lighting is applied.
 * k_a: Ambient color
 * k_d: Diffuse color
 * k_s: Specular color
 * alpha: Shininess coefficient
 * p: position of point being lit
 * eye: the position of the camera
 *
 * See https://en.wikipedia.org/wiki/Phong_reflection_model#Description
 */
vec3 phongIllumination(vec3 k_a, vec3 k_d, vec3 k_s, float alpha, vec3 p, vec3 eye) {
    const vec3 ambientLight =  vec3(0.2, 0.2, 0.2);
    vec3 color = ambientLight * k_a;
    
    vec3 light1Pos = vec3(-4.0,
                          0.0,
                          4.0);
    vec3 light1Intensity = vec3(0.3, 0.3, 0.3);
    
    color += phongContribForLight(k_d, k_s, alpha, p, eye,
                                  light1Pos,
                                  light1Intensity);
    
    vec3 light2Pos = vec3(4.,
                          0.,
                          -4.);
    vec3 light2Intensity = vec3(0.3, 0.3, 0.3);
    
    color += phongContribForLight(k_d, k_s, alpha, p, eye,
                                  light2Pos,
                                  light2Intensity);
    vec3 light3Pos = vec3(0.0, 8.0, 0.0);
    vec3 light3Intensity = vec3(1., 1.0, 1.0);
    color += phongContribForLight(k_d, k_s, alpha, p, eye,
                                  light3Pos,
                                  light3Intensity);
    return color;
}

#define iterations 15
#define formuparam 0.53

#define stepsize 0.2

#define zoom   2.00
#define tile   0.850

#define brightness 0.0010
#define darkmatter 0.500
#define distfading 0.30
#define saturation 0.850

vec4 getBackground(vec2 fragCoord, in vec3 dir) {
	vec3 from=vec3(0.,0.0,0.0);
	
	//volumetric rendering
	float s=0.2,fade=1.2;
	vec3 v=vec3(0.);
	for (int r = 0; r < 18; r++) {
		vec3 p=from+s*dir;
		p = abs(vec3(tile)-mod(p,vec3(tile*2.))); // tiling fold
		float pa,a=pa=0.;
		for (int i=0; i<iterations; i++) { 
			p=abs(p)/dot(p,p)-formuparam; // the magic formula
			a+=abs(length(p)-pa); // absolute sum of average change
			pa=length(p);
		}
		float dm=max(0.,darkmatter-a*a); //dark matter
		a*=a*a; // add contrast
		v+=fade;
		v+=vec3(s,s*s,s*s)*a*brightness*fade; // coloring based on distance
		fade*=distfading; // distance fading
		s+=stepsize;
	}
	v=mix(vec3(length(v)),v,saturation); //color adjust
	return vec4(v*vec3(.01, 0.005, 0.01),1.);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    float angle = radians(
          25.*SMOOTHTR(T2, T3)
        - 45.*SMOOTHTR(T4, T5)
        + 20.*SMOOTHTR(T5, T6)
    );
    vec3 worldOrigin;
    worldOrigin.x =
          0.8*SMOOTHTR(T1, T2)
        - 1.6*SMOOTHTR(T3, T4)
        + 0.8*SMOOTHTR(T5, T6);
    worldOrigin.y =
          1.*SMOOTHTR(T0, T1)
        + 1.*SMOOTHTR(T1, T2)
        + 5.*SMOOTHTR(T5, T6)
         -1.;//-4.+1.;
    worldOrigin.z =
          2.2*SMOOTHTR(T1, T2);
    float eyeDist =
        - 1.4*SMOOTHTR(T1, T2)
        + 0.4*SMOOTHTR(T2, T3)
        + 5.; //+ 3.;
    float eyeY =
        - 0.5*SMOOTHTR(T1, T2)
        + 2.0*SMOOTHTR(T5, T6)
        + 1.5; //- 2.;

    vec3 eye = vec3(
        eyeDist*sin(angle),
        eyeY,
        eyeDist*cos(angle)
    );
    float viewAngle = 75.;

    vec3 viewDir = rayDirection(viewAngle, 
            iResolution.xy, fragCoord);
    
    mat4 viewToWorld = viewMatrix(
        eye,
        worldOrigin,
        vec3(0.0, 1.0, 0.0) // y axis in view ref system
    );
    vec3 worldDir = (viewToWorld * vec4(viewDir, 0.0)).xyz;
    float dist = shortestDistanceToSurface(eye, worldDir, MIN_DIST, MAX_DIST);
    
    if (dist > MAX_DIST - EPSILON) {
        // Didn't hit anything
        fragColor = getBackground(
            fragCoord,
            worldDir
        );
        return;
    }
    
    // The closest point on the surface to the eyepoint along the view ray
    vec3 p = eye + dist * worldDir;
    vec3 K_a = vec3(0.4, 0.4, 0.4);
    vec3 K_d = vec3(0.5, 0.5, 0.7);
    vec3 K_s = vec3(0.8, 0.8, 1.0);
    float shininess = 2.0;
    vec3 color = phongIllumination(K_a, K_d, K_s, shininess, p, eye);
    
    fragColor = vec4(color, 1.0);
}
uniform float iGlobalTime;
uniform vec2 iResolution;


// Algorithm
#define MIN_DIST 0.0
#define MAX_DIST 15.0
#define EPSILON 0.005

#define MAX_MARCHING_STEPS_1 120
#define MARCHING_STEP_1 0.4

#define MAX_MARCHING_STEPS_2 200
#define MARCHING_STEP_2 0.2

// Times (part 1)
#define T0 0.9
#define D1 4.0
#define D2 4.0
#define D3 4.0
#define D4 4.0
#define D5 4.0
#define D6 4.0

// part 2 starts here
#define D7 2.0 // sight on heart
#define D8 5. // heart starts spinning
#define D9 15.0 // heart ends spinning && zoom out
#define D10 0.5 // spike going down

// Animation constants (part 2)
#define HEART_BEAT_T0 T6
#define HEART_BEAT_T1 T9
#define HEART_BEAT_OMEGA 30.0
#define HEART_ROTATE_T0 T8
#define HEART_ROTATE_T1 T9
#define HEART_RADIUS 2.3 // for color
#define SPIKE_HEIGHT 2.0
#define SPIKE_YFINAL -2.4


#define CAMERA_ANGLE0 25.
#define CAMERA_Y1 1.3
#define CAMERA_OMEGA 0.2

// Computed
#define SPIKE_Y0 SPIKE_YFINAL - SPIKE_HEIGHT
#define T1 T0 + D1
#define T2 T1 + D2
#define T3 T2 + D3
#define T4 T3 + D4
#define T5 T4 + D5
#define T6 T5 + D6
#define T7 T6 + D7
#define T8 T7 + D8
#define T9 T8 + D9
#define T10 T9 + D10
#define TTRANS T6
#define SPIKE_T0 T9
#define SPIKE_T1 T10

// Operations
#define OP_U(X, Y) min(X, Y)
#define OP_S(X, Y) max( -Y, X )
#define OP_I(X, Y) max( X, Y )
#define SMOOTHTR(X0, X1) smoothstep(X0, X1, iGlobalTime)
mat3 Rot4Z(float a) { // angle in radians
    float c = cos( a );
    float s = sin( a );
    return mat3(
        c,-s, 0,
        s, c, 0,
        0, 0, 1
    );
}

float timeTransition(float t0, float t1) {
    return clamp((iGlobalTime - t0)/(t1-t0), 0., 1.);
}

// SDF primitives
float cubeSDF(vec3 p, vec3 h) {
    vec3 d = abs(p) - h;
    float insideDistance = min(max(d.x, max(d.y, d.z)), 0.0);
    float outsideDistance = length(max(d, 0.0));
    return insideDistance + outsideDistance;
}

float cylinderSDF( vec3 p, vec2 c, float r ) {
    return length(p.xy - c.xy) - r;
}
float vcylinderSDF( vec3 p, vec2 c, float r ) {
    return length(p.xz - c.xy) - r;
}

float sphereSDF(vec3 p, float r) {
    return length(p) - r;
}

float planeSDF( vec3 p, vec4 n ) {
    return dot( p, n.xyz ) + n.w;
}

float coneSDF( vec3 p, vec2 c ) {
    float q = length(p.xz);
    return dot(c,vec2(q,p.y));
}

float cappedCylinderSDF( vec3 p, vec2 h )
{
  vec2 d = abs(vec2(length(p.xz),p.y)) - h;
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

// Figures (Part 1)
float simpleEmblemSDF(vec3 p) { 
    // Heart
    vec3 q = p;
    q.z = q.z * (- q.y/15. + 3.);
    q.y = - 1.2*q.y - abs(q.x)*sqrt((20.-abs(q.x))/15.);
    float sphere = sphereSDF(q, 0.3);
    q = p;
    q.x = abs(q.x);
    float cylinder = cylinderSDF(q, vec2(0.25, 0.1), 0.18);
    float heart = OP_S(sphere, cylinder);
    
    // Spike
    vec3 spikeQ = p;
    spikeQ.y -= 0.5;
    q = abs(spikeQ);
    q.y = q.y*(2.*q.x + 1.);
    q.z = q.z*(0.2*q.x + 3.5);
    sphere = cubeSDF(q, vec3(4., 4., 4.));
    float plane = planeSDF(q, normalize(
        vec4(-1., -1., -1., 0.4)));
    float diamond = OP_S(sphere, plane);
    q = abs(spikeQ);
    cylinder = cylinderSDF(q, vec2(0.18, 0.18), 0.13);
    float spike = OP_S(diamond, cylinder);
    
    return OP_U(heart, spike);
}

float altarSDF(vec3 p) {
    // Altar
    vec3 q = p;
    q.y -= 4.5;
    float cone = coneSDF(q, normalize(vec2(1., 0.5)));
    q = p;
    float p1 = planeSDF(q, vec4(0., -1., 0., 0.));
    float res = OP_S(cone, p1);
    
    // Stairs
    q.z -= 2.5;
    q.y += 0.8;
    float stair = cubeSDF(q,
       vec3(0.85, 0.8, 0.8));
    res = OP_U(res, stair);
    q.y = p.y - 0.6;
    for (int i = 0; i != 5; ++i) {
        stair = cubeSDF(q, vec3(0.95, 0.8, 0.2));
        res = OP_S(res, stair);
        q.y += 0.2;
        q.z -= 0.2;
    };
        
    // Stair sides
    q = p;
    q.x = abs(p.x);
    q.y += 0.7;
    q.z -= 2.8;
    q.x -= 0.8;
 	float stairSides = cubeSDF(q, vec3(0.1, 0.85, 1.));
    float stairPlane = planeSDF(q, vec4(0., -1., -1., 0.35));
    stairSides = OP_S(stairSides, stairPlane);
    res = OP_U(res, stairSides);
    
    // Stair blocks
    q = p;
    q.x = abs(q.x);
    q.z -= 2.1;
    q.x -= 0.8;
    q.y -= 0.2;
    float stairBlock = cubeSDF(q, vec3(0.2, 0.08, 0.25));
    res = OP_U(res, stairBlock);
    
    // Emblems
    q = p;
	q.x = abs(q.x);
    q.x -= 0.8;
    q.y -= 0.8;
    q.z -= 2.2;
	res = OP_U(res, simpleEmblemSDF(q));
    
    // Edge
    q = p;
    float cyl1 = cappedCylinderSDF(q, vec2(2.35, 0.15));
    float cyl2 = vcylinderSDF(q, vec2(0.0, 0.0), 2.0);
    q.z -= 1.4;
    float entry = cubeSDF(q, vec3(0.85, 1., 1.));
    float edge = OP_S(cyl1, cyl2);
    edge = OP_S(edge, entry);
    res = OP_U(res, edge);
   
    return res;
}

// Part 2
vec3 moveSpike( vec3 p ) {
    return vec3(
        p.x,
        SPIKE_HEIGHT*smoothstep(SPIKE_T0, SPIKE_T1, iGlobalTime) + p.y + SPIKE_Y0,
        p.z
    );
}

vec3 moveHeart (vec3 p) {
    //p.y -= 4.;
    float transition = 1.-timeTransition(
        HEART_BEAT_T0, HEART_BEAT_T1);
    float amplitude = mix(0., 0.14, transition);
    float omega = mix(HEART_BEAT_OMEGA/2., HEART_BEAT_OMEGA, transition);
    float heartbeat = amplitude*pow(abs(sin(omega*iGlobalTime)), 4.) + 1.;
    vec3 res = vec3(
        p.x,
        p.y,
        p.z
    ) * heartbeat;
    float angle = 180.0*(1.-smoothstep(
        HEART_ROTATE_T0, HEART_ROTATE_T1, iGlobalTime));
    mat3 rot = Rot4Z(radians(angle));
    return rot*res;
}

// Figures
float heartSDF( vec3 p ) {
    p = moveHeart(p);
    vec3 q = p;
    q.z = q.z * (3. - q.y/15.);
    q.y = - 1.2*q.y - abs(q.x)*sqrt((20.-abs(q.x))/15.) + 0.2;
    float sphere = sphereSDF(q, 1.5);
    q = p;
    q.x = abs(q.x);
    float cylinder = cylinderSDF(q, vec2(0.6, 0.4), 0.42);
    return OP_S(sphere, cylinder);
}
const vec3 gSpikeNorm = normalize(vec3(-1., -1., -1.));

float spikeSDF( vec3 samplePoint ) {
    samplePoint = moveSpike(samplePoint);
    vec3 q = abs(samplePoint);
    q.y = q.y*(2.*q.x + 1.);
    q.z = q.z*(0.2*q.x + 3.5);
    float sph = cubeSDF(q, vec3(4., 4., 4.));
    float plane = planeSDF(q, vec4(gSpikeNorm, 1.));
    float diamond = OP_S(sph, plane);
    q = samplePoint;
    q = abs(q);
    float cylinder = cylinderSDF(q, vec2(0.6, 0.6), 0.42);
    return OP_S(diamond, cylinder);
}

float sceneSDF(vec3 p) {
    if (iGlobalTime < TTRANS)
        return altarSDF(p);
    else
        return OP_U(heartSDF(p), spikeSDF(p));
}

// Marching algorithm
float shortestDistanceToSurface(vec3 eye, vec3 marchingDirection, int max_steps, float march_step) {
    float depth = MIN_DIST;
    vec3 p;
    for (int i = 0; i < max_steps; i++) {
        p = depth * marchingDirection + eye;
        float dist = sceneSDF(p);
        if (abs(dist) < EPSILON) {
            return depth;
        }
        depth += march_step*dist;
        if (depth >= MAX_DIST) {
            return MAX_DIST;
        }
    }
    return MAX_DIST;
}
            
// Get ray direction given frag coord
vec3 rayDirection(float fieldOfView, vec2 size, vec2 fragCoord) {
    vec2 xy = - size / 2.0 + fragCoord;
    float z = size.y / tan(radians(fieldOfView) / 2.0);
    return normalize(vec3(xy, -z));
}

// Transformation matrix from view to world
mat4 viewMatrix(vec3 eye, vec3 center, vec3 up) {
    vec3 f = normalize(center - eye);
    vec3 s = normalize(cross(f, up));
    vec3 u = cross(s, f);
    return mat4(
        vec4(s, 0.0),
        vec4(u, 0.0),
        vec4(-f, 0.0),
        vec4(0.0, 0.0, 0.0, 1)
    );
}

// Get normal
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        sceneSDF(vec3(p.x + EPSILON, p.y, p.z)) - sceneSDF(vec3(p.x - EPSILON, p.y, p.z)),
        sceneSDF(vec3(p.x, p.y + EPSILON, p.z)) - sceneSDF(vec3(p.x, p.y - EPSILON, p.z)),
        sceneSDF(vec3(p.x, p.y, p.z  + EPSILON)) - sceneSDF(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

/**
 * Lighting contribution of a single point light source via Phong illumination.
 *
 * The vec3 returned is the RGB color of the light's contribution.
 *
 * k_a: Ambient color
 * k_d: Diffuse color
 * k_s: Specular color
 * alpha: Shininess coefficient
 * p: position of point being lit
 * eye: the position of the camera
 * lightPos: the position of the light
 * lightIntensity: color/intensity of the light
 *
 * See https://en.wikipedia.org/wiki/Phong_reflection_model#Description
 */
vec3 phongContribForLight(vec3 k_d, vec3 k_s, float alpha, vec3 p, vec3 eye,
                          vec3 lightPos, vec3 lightIntensity) {

    vec3 N = estimateNormal(p);
    vec3 L = normalize(lightPos - p);
    vec3 V = normalize(eye - p);
    vec3 R = normalize(reflect(-L, N));
    
    float dotLN = dot(L, N);
    float dotRV = dot(R, V);
    
    if (dotLN < 0.0) {
        // Light not visible from this point on the surface
        return vec3(0.0, 0.0, 0.0);
    }
    
    if (dotRV < 0.0) {
        // Light reflection in opposite direction as viewer, apply only diffuse
        // component
        return lightIntensity * (k_d * dotLN);
    }
    return lightIntensity * (k_d * dotLN + k_s * pow(dotRV, alpha));
}

/**
 * Lighting via Phong illumination.
 *
 * The vec3 returned is the RGB color of that point after lighting is applied.
 * k_a: Ambient color
 * k_d: Diffuse color
 * k_s: Specular color
 * alpha: Shininess coefficient
 * p: position of point being lit
 * eye: the position of the camera
 *
 * See https://en.wikipedia.org/wiki/Phong_reflection_model#Description
 */
vec3 phongIllumination(vec3 k_a, vec3 k_d, vec3 k_s, float alpha, vec3 p, vec3 eye) {
    const vec3 ambientLight =  vec3(0.2, 0.2, 0.2);
    vec3 color = ambientLight * k_a;
    
    vec3 light1Pos = vec3(-4.0,
                          0.0,
                          4.0);
    vec3 light1Intensity = vec3(0.3, 0.3, 0.3);
    
    color += phongContribForLight(k_d, k_s, alpha, p, eye,
                                  light1Pos,
                                  light1Intensity);
    
    vec3 light2Pos = vec3(4.,
                          0.,
                          -4.);
    vec3 light2Intensity = vec3(0.3, 0.3, 0.3);
    
    color += phongContribForLight(k_d, k_s, alpha, p, eye,
                                  light2Pos,
                                  light2Intensity);
    vec3 light3Pos = vec3(0.0, 8.0, 0.0);
    vec3 light3Intensity = vec3(1., 1.0, 1.0);
    color += phongContribForLight(k_d, k_s, alpha, p, eye,
                                  light3Pos,
                                  light3Intensity);
    return color;
}

#define iterations 15
#define formuparam 0.53

#define stepsize 0.2

#define zoom   2.00
#define tile   0.850

#define brightness 0.0010
#define darkmatter 0.500
#define distfading 0.30
#define saturation 0.850

vec4 getBackground(vec2 fragCoord, in vec3 dir) {
	vec3 from=vec3(0.,0.0,0.0);
	
	//volumetric rendering
	float s=0.2,fade=1.2;
	vec3 v=vec3(0.);
	for (int r = 0; r < 18; r++) {
		vec3 p=from+s*dir;
		p = abs(vec3(tile)-mod(p,vec3(tile*2.))); // tiling fold
		float pa,a=pa=0.;
		for (int i=0; i<iterations; i++) { 
			p=abs(p)/dot(p,p)-formuparam; // the magic formula
			a+=abs(length(p)-pa); // absolute sum of average change
			pa=length(p);
		}
		float dm=max(0.,darkmatter-a*a); //dark matter
		a*=a*a; // add contrast
		v+=fade;
		v+=vec3(s,s*s,s*s)*a*brightness*fade; // coloring based on distance
		fade*=distfading; // distance fading
		s+=stepsize;
	}
	v=mix(vec3(length(v)),v,saturation); //color adjust
	return vec4(v*vec3(.01, 0.005, 0.01),1.);
}

void main()
{
    vec2 fragCoord=gl_FragCoord.xy;
    vec4 fragColor=vec4(0.0,0.0,0.0,1.0);

    float angle = radians(
          25.*SMOOTHTR(T2, T3)
        - 45.*SMOOTHTR(T4, T5)
        + 20.*SMOOTHTR(T5, T6)
		+ 360.*SMOOTHTR(T8, T9)
    );
    vec3 worldOrigin;
    worldOrigin.x =
          0.8*SMOOTHTR(T1, T2)
        - 1.6*SMOOTHTR(T3, T4)
        + 0.8*SMOOTHTR(T5, T6);
    worldOrigin.y =
          1.*SMOOTHTR(T0, T1)
        + 1.*SMOOTHTR(T1, T2)
        + 5.*SMOOTHTR(T5, T6)
        - 4.5*SMOOTHTR(T6, T7)
		+ 0.5*SMOOTHTR(T8, T9)
         -1.;
    worldOrigin.z =
          2.2*SMOOTHTR(T1, T2);
    float eyeDist =
        - 1.4*SMOOTHTR(T1, T2)
        + 0.4*SMOOTHTR(T2, T3)
        + 1.0*SMOOTHTR(T6, T7)
		+ 4.0*SMOOTHTR(T8, T9)
        + 5.;
    float eyeY =
        - 0.5*SMOOTHTR(T1, T2)
        + 2.0*SMOOTHTR(T5, T6)
        + 1.5;

    vec3 eye = vec3(
        eyeDist*sin(angle),
        eyeY,
        eyeDist*cos(angle)
    );
    float viewAngle = 75.;

    vec3 viewDir = rayDirection(viewAngle, 
            iResolution.xy, fragCoord);
    
    mat4 viewToWorld = viewMatrix(
        eye,
        worldOrigin,
        vec3(0.0, 1.0, 0.0) // y axis in view ref system
    );
    vec3 worldDir = (viewToWorld * vec4(viewDir, 0.0)).xyz;
	int max_steps = (iGlobalTime < TTRANS) ? MAX_MARCHING_STEPS_1 : MAX_MARCHING_STEPS_2;
	float step_size = (iGlobalTime < TTRANS) ? MARCHING_STEP_1 : MARCHING_STEP_2;
    float dist = shortestDistanceToSurface(eye, worldDir, max_steps, step_size);
    
    if (dist > MAX_DIST - EPSILON) {
        // Didn't hit anything
        gl_FragColor= getBackground(
            fragCoord,
            worldDir
        );

        return;
    }
    
    // The closest point on the surface to the eyepoint along the view ray
    vec3 p = eye + dist * worldDir;
    vec3 K_a = vec3(0.4, 0.4, 0.4);
    vec3 K_d = vec3(0.5, 0.5, 0.7);
    vec3 K_s = vec3(0.8, 0.8, 1.0);
	if (iGlobalTime > TTRANS && length(p) < HEART_RADIUS) {
		float tr = timeTransition(T7, HEART_ROTATE_T0);
		K_d = mix(vec3(1., 0., 0.), K_d, tr);
		K_s = mix(vec3(1., 0., 0.), K_s, tr);
	}
    float shininess = 2.0;
    vec3 color = phongIllumination(K_a, K_d, K_s, shininess, p, eye);
    
    gl_FragColor= vec4(color, 1.0);
}