#define PI 3.14159265359
#define MAX_STEPS 1000
#define MAX_DIST 100.
#define SURFACE_DIST 0.001
#define EPSILON 0.000001
#define UNDERSTEP 0.7

mat3 rotateX(float theta) {
    float c = cos(theta);
    float s = sin(theta);
    return mat3(
        vec3(1, 0, 0),
        vec3(0, c, -s),
        vec3(0, s, c)
    );
}

mat3 rotateY(float theta) {
    float c = cos(theta);
    float s = sin(theta);
    return mat3(
        vec3(c, 0, s),
        vec3(0, 1, 0),
        vec3(-s, 0, c)
    );
}

mat3 rotateZ(float theta) {
    float c = cos(theta);
    float s = sin(theta);
    return mat3(
        vec3(c, -s, 0),
        vec3(s, c, 0),
        vec3(0, 0, 1)
    );
}

vec4 smoothUnionSDF(vec4 d1, vec4 d2, float k) {
    float h = clamp( 0.5 + 0.5*(d2.x-d1.x)/k, 0.0, 1.0 );
    return vec4( mix( d2.x, d1.x, h ) - k*h*(1.0-h), mix(d2.yzw, d1.yzw, h));
}

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

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

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

float sdEllipsoid(vec3 p, vec3 r) {
  float k0 = length(p/r);
  float k1 = length(p/(r*r));
  return k0*(k0-1.0)/k1;
}

vec4 cp(vec4 sd1, vec4 sd2) {
    return (sd1.x < sd2.x) ? sd1 : sd2;
}

mat4 lookAt(vec3 eye, vec3 target, vec3 up) {
	vec3 z = normalize(eye - target);
    vec3 x = normalize(cross(up, z));
    vec3 y = normalize(cross(z, x));
    mat4 translation = mat4(
        vec4(1.0, 0.0, 0.0, eye.x),
    	vec4(0.0, 1.0, 0.0, eye.y),
        vec4(0.0, 0.0, 1.0, eye.z),
        vec4(0.0, 0.0, 0.0, 1.0)
    );
    mat4 rotation = transpose(
        mat4(
            vec4(x, 0.0),
            vec4(y, 0.0),
            vec4(z, 0.0),
            vec4(0.0, 0.0, 0.0, 1.0)
        )
    );
    return translation * rotation;
}

vec3 getRayDir(mat4 camera, float fov, vec2 pixel) {
    vec2 p = 2.0 * ((pixel / iResolution.xy) - vec2(0.5, 0.5));
    vec2 h = vec2(
        tan(fov / 2.0) * (iResolution.x / iResolution.y), 
        tan(fov / 2.0)
    );
    vec3 pCam = vec3(p * h, -1.0);
	return normalize((inverse(camera) * vec4(pCam, 0.0)).xyz);
}

vec4 sdFlower(vec3 p, float phaseShift) {
    float phase = abs(sin((iTime + phaseShift)*.2)) *  3.;
    vec4 dStem = vec4(sdCappedCylinder(p - vec3(sin(p.y*10.) * .1, 0., cos(p.y*10.) * .1), clamp(.2 * phase / 3., .04, 1.), phase), .1, .7, .1);
    
    float leave1Phase = clamp(phase, 1., 2.) - 1.;
    float dLeave1 = min(
            sdEllipsoid(rotateX(.8) * (p - vec3(0, 1., .5 * leave1Phase)), vec3(sin(p.y) * .6, .6, sin(p.y) * .2) * leave1Phase),
            sdEllipsoid(rotateX(-.8) * (p - vec3(0, 1., -.5 * leave1Phase)), vec3(sin(p.y) * .6, .6, sin(p.y) * .2) * leave1Phase));
    
    float leave2Phase = clamp(phase, 2., 3.) - 2.;
    float dLeave2 = min(
            sdEllipsoid(rotateZ(.8) * (p - vec3(-.5 * leave2Phase, 2., 0)), vec3(sin(p.y-1.) * .2, .6, sin(p.y - 1.) * .6) * leave2Phase),
            sdEllipsoid(rotateZ(-.8) * (p - vec3(.5 * leave2Phase, 2., 0)), vec3(sin(p.y-1.) * .2, .6, sin(p.y - 1.) * .6) * leave2Phase));
    
    float flowerR = .9;
    float flowerPhase = clamp(phase, 3. - flowerR, 3.) -3. + flowerR;
    float dFlower = sdSphere(p - vec3(sin(p.y*10.)* .1, 3. - flowerR + flowerPhase * 1.2, cos(p.y*10.) * .1), abs(sin(p.y*20.))*.1*flowerPhase+ flowerR * flowerPhase);
    
    vec4 hit = dStem;
    hit = cp(smoothUnionSDF(dStem, vec4(dLeave1, .5, .7, .1), .05), hit);
    hit = cp(smoothUnionSDF(dStem, vec4(dLeave2, .5, .7, .1), .05), hit);
    hit = cp(smoothUnionSDF(dStem, vec4(dFlower, .7, .2, .3), .1), hit);
    return hit;
}

vec4 getDist(vec3 p) {
    vec4 dPlane = vec4(sdPlane(p, 0.), .2, .2, .2);
    //vec4 dFlower = sdFlower(rotateY(mod(p.x, 6.) / mod(p.z, 6.)) * vec3(mod(p.x, 20.) - 3., p.y, mod(p.z, 20.) - 3.), sin(p.z/p.x)*5.);
    vec4 dFlower = sdFlower(p, 0.);
    
    vec4 hit = dPlane;
    hit = cp(dFlower, hit);
    return hit;
}

vec4 rayMarch(vec3 ro, vec3 rd) {
 float dO = 0.;
 vec4 dS = vec4(0);
 for(int i = 0; i < MAX_STEPS; i++) {
     vec3 p = ro + dO * rd;
     dS = getDist(p);
     dO += dS.x * UNDERSTEP;
     if(dS.x < SURFACE_DIST || dO > MAX_DIST) break;
 }
 return vec4(dO, dS.y, dS.z, dS.w);
}

vec3 getNormal(vec3 p) {
    float d = getDist(p).x;
    vec2 e = vec2(EPSILON, 0);
    vec3 n = d - vec3(
        getDist(p - e.xyy).x,
        getDist(p - e.yxy).x,
        getDist(p - e.yyx).x);
    return normalize(n);
}

vec3 getPointLight(vec3 p, vec3 lightPos, vec3 n, vec3 fragToLight, vec3 fragToCam, vec3 texColor, vec3 lightColor) {
    float valDiffuse = max(0.0, dot(n, fragToLight));
    
    float d = rayMarch(p + n * SURFACE_DIST, fragToLight).x;
    if(d < length(lightPos-p)) {
        return valDiffuse * texColor * lightColor * .1;
    }
    
    vec3 blinnH = normalize(fragToLight + fragToCam);
    float valSpecular = pow(max(0.0, dot(n, blinnH)), 50.);
    return valDiffuse*texColor*lightColor + valSpecular*lightColor;
}

vec3 getLight(vec3 p, vec3 cam, vec3 col) {
    vec3 lightPos = vec3(-5, 2.+abs(sin(iTime*.1))*2., 0);
    vec3 lightPos1 = vec3(5, 2.+abs(sin(iTime*.1))*2., 0);
    vec3 lightPos2 = vec3(0, 2.+abs(sin(iTime*.1))*2., -5);
    vec3 lightPos3 = vec3(0, 2.+abs(sin(iTime*.1))*2., 5);
    vec3 ln = normalize(lightPos - p);
    vec3 ln1 = normalize(lightPos1 - p);
    vec3 ln2 = normalize(lightPos2 - p);
    vec3 ln3 = normalize(lightPos3 - p);
    
    
    vec3 lightCol = vec3(.7, .7, .7);
    vec3 n = getNormal(p);
    vec3 cn = normalize(cam - p);
    
    return 
        getPointLight(p, lightPos, n, ln, cn, col, lightCol) + 
        getPointLight(p, lightPos1, n, ln1, cn, col, lightCol) +
        getPointLight(p, lightPos2, n, ln2, cn, col, lightCol) +
        getPointLight(p, lightPos3, n, ln3, cn, col, lightCol);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    //vec3 ro = vec3(0., 2.5, 10.);
    vec3 ro = vec3(sin(iTime*.5) * 10., 4.5, cos(iTime*.5) * 10.);
    vec3 target = vec3(0.0, 2.0, 0.0);
    vec3 up = vec3(0.0, 1.0, 0.0);
    float fov = radians(45.0 + sin(iTime*.04) * 10.);
    
    mat4 camera = lookAt(ro, target, up); 
    vec3 rd = getRayDir(camera, fov, fragCoord);
    vec4 d = rayMarch(ro, rd);
    if(d.x <= MAX_DIST) {//sqrt(pow(MAX_DIST, 2.) + pow(length(ro - target - rd), 2.))) {
        vec3 p = ro + rd * d.x;
        vec3 dif = getLight(p, ro, d.yzw);
        fragColor = vec4(dif, 1.0);
        //fragColor = vec4(vec3(d.x * .05), 1.);
    } else {
        vec2 uv = (fragCoord - 0.5 * iResolution.xy) / iResolution.y;
        vec3 col =  -vec3(.5, .5, 0.) * uv.y + 0.5;
        fragColor = vec4(col, 1.0);
    }
}