// inspired by https://www.shadertoy.com/view/ldfXzn

#ifdef GL_ES
precision highp float;
#endif

#define DELTA               0.01
#define RAY_LENGTH_MAX      80.0
#define MAX_STEPS           120
#define MAX_RAYS            7 // start ray, [reflection, refraction]...
#define MAX_SUBREFRACTIONS  20
#define REFRACTION          1.762

#define reflectionI         (ray * 2 + 1)
#define backupI             (ray + 1)
#define refractionI         (ray * 2 + 2)
#define REPEAT_PATTERN      vec3(10.0, 10.0, 10.0) * 3.2
#define MAX_REPEAT          1.0
#define USE_ADDITIONAL_RAYS false

// gem definition
#define TOP_CONE_ANGLE      vec2(0.5038710255240861, 0.8637789008984333)
#define BOT_CONE_ANGLE      vec2(0.7808688094430303, -0.6246950475544243)
#define GEM_SCALE           21.0
#define GEM_TOP_OFFSET      0.2
#define GEM_BOT_OFFSET      0.2
#define TOP_CONE_HEIGHT     0.05
#define MAGICSTEPS          10.0

#define REFLECTION_POWER     10.0

#define M_PI              3.1415926535897932384626433832795

mat3 mRotate (in vec3 angle) {
    float c = cos (angle.x);
    float s = sin (angle.x);
    mat3 rx = mat3 (1.0, 0.0, 0.0, 0.0, c, s, 0.0, -s, c);

    c = cos (angle.y);
    s = sin (angle.y);
    mat3 ry = mat3 (c, 0.0, -s, 0.0, 1.0, 0.0, s, 0.0, c);

    c = cos (angle.z);
    s = sin (angle.z);
    mat3 rz = mat3 (c, s, 0.0, -s, c, 0.0, 0.0, 0.0, 1.0);

    return rz * ry * rx;
}


float cone(in vec3 p, in vec3 c) {
    vec2 q = vec2(length(p.xz), p.y);
    float d1 = -p.y-c.z;
    float d2 = max(dot(q, c.xy), p.y);
    return length(max(vec2(d1, d2),0.0)) + min(max(d1, d2), 0.0);
}

vec3 rayForRepeat(in vec3 ray, in vec3 repeatPattern) {
    vec3 q = ray + repeatPattern * 0.5;
    vec3 k = floor (q / repeatPattern);

    q -= repeatPattern * (k + 0.5);

    if(max(max(abs(k.x), abs(k.y)), abs(k.z)) > MAX_REPEAT) q = ray;

    //vec3 k =vec3(1.0,1.0,1.0);
    mat3 rotateMatrix = mRotate (k + iTime * 0.1);
    //revertBy = inverse(rotateMatrix);
    return rotateMatrix * q;
}

vec3 rayForRepeat(in vec3 ray, in vec3 repeatPattern, out mat3 revertBy) {
    //vec3 overflowedRay   = mod(ray,repeatPattern);
    //vec3 patternCenter   = -0.5 * repeatPattern;
    //return patternCenter + overflowedRay;

    vec3 q = ray + repeatPattern * 0.5;
    vec3 k = floor (q / repeatPattern);
    //if(max(max(k.x, k.y), k.z) < 2.0)
    q -= repeatPattern * (k + 0.5);

    if(max(max(abs(k.x), abs(k.y)), abs(k.z)) > MAX_REPEAT) q = ray;
    //vec3 k =vec3(1.0,1.0,1.0);
    mat3 rotateMatrix = mRotate (k + iTime * 0.1);
    revertBy = rotateMatrix;//inverse(rotateMatrix);
    return rotateMatrix * q;
}

float conebination( vec3 rayDot )
{
    float scale = GEM_SCALE;
    //vec2 top = normalize(vec2(0.35, 0.6));
    vec2 top = TOP_CONE_ANGLE;
    //vec2 bot = normalize(vec2(0.5, -0.4));
    vec2 bot = BOT_CONE_ANGLE;

    float horizontal = length(vec2(rayDot.x,rayDot.z)) / scale;
    float vertical   = rayDot.y / scale;

    float gemEdgeY = abs(top.y) * (GEM_TOP_OFFSET + GEM_BOT_OFFSET) / (abs(top.y) + abs(top.x*bot.y/bot.x)) -GEM_TOP_OFFSET;

    float dotTop = top.x * horizontal + top.y * (vertical - GEM_TOP_OFFSET);
    float dotBot = bot.x * horizontal + bot.y * (vertical + GEM_BOT_OFFSET);

    return max(max(dotTop, dotBot), vertical - gemEdgeY - TOP_CONE_HEIGHT) * scale;
}

float atan2(in float y, in float x)
{
    float s = (abs(x) > abs(y)) ? 1.0 : 0.0;
    return mix(M_PI/2.0 - atan(x,y), atan(y,x), s);
}

vec3 conebinationNormal( vec3 rayDot )
{

    float scale = GEM_SCALE;
    vec2 top = TOP_CONE_ANGLE;
    vec2 bot = BOT_CONE_ANGLE;


    mat3 revertBy;
    vec3 modifiedRay = rayForRepeat(rayDot, REPEAT_PATTERN, revertBy) / scale;

    vec2 horizontalNormal  = normalize(vec2(modifiedRay.x,modifiedRay.z));

    // fake normals
    vec2 topAdd = vec2(0.0, 0.0);

    // fake normals end

    float horizontal = length(vec2(modifiedRay.x,modifiedRay.z));
    float vertical   = modifiedRay.y;

    float gemEdgeY = abs(top.y) * (GEM_TOP_OFFSET + GEM_BOT_OFFSET) / (abs(top.y) + abs(top.x*bot.y/bot.x)) -GEM_TOP_OFFSET;

    float dotTop = top.x * horizontal + top.y * (vertical - GEM_TOP_OFFSET);
    float dotBot = bot.x * horizontal + bot.y * (vertical + GEM_BOT_OFFSET);
    float dotCut = vertical - gemEdgeY - TOP_CONE_HEIGHT;

    float maxDot = max(dotTop, dotBot);
    maxDot = max(maxDot, dotCut);

    vec2 verticalNormal = (maxDot == dotCut) ?
        vec2(0.0, 1.0) :
        ((maxDot == dotTop) ? top + topAdd: bot)  ;

    float angle = atan2(horizontalNormal.x, horizontalNormal.y);

    float sideStepB = (2.0 * M_PI / 14.0);
    float stepAngleB = (floor(angle / sideStepB) + 0.5) * sideStepB;
    float subAngleB = angle - stepAngleB;
    vec2 horizontalNormalB = normalize(vec2(sin(stepAngleB), cos(stepAngleB)));

    if (maxDot == dotTop) {
        float sideStepT = (2.0 * M_PI / 7.0);
        float stepAngleT = (floor(angle / sideStepT) + 0.5) * sideStepT;
        float subAngleT = angle - stepAngleT;

        vec2 horizontalNormalT = normalize(vec2(sin(stepAngleT), cos(stepAngleT)));

        float fullConeHeight = (gemEdgeY - GEM_TOP_OFFSET);
        float hitConeHeight = (vertical);
        float maxConeHeight = (TOP_CONE_HEIGHT);



        float ratio = abs(gemEdgeY - vertical)/TOP_CONE_HEIGHT;
        bool outside = (abs(fullConeHeight * sin(subAngleT) /(vertical)) > mix(0.0, sin(sideStepT / 2.0), ratio));

        horizontalNormal = outside ? horizontalNormalB : horizontalNormalT;
    } else  {
        horizontalNormal = horizontalNormalB;
    }


    return normalize(vec3(
        horizontalNormal.x * verticalNormal.x,
        verticalNormal.y,
        horizontalNormal.y * verticalNormal.x
    )) * revertBy;
}




float cube(in vec3 ray)
{
    vec3 sides = vec3 (2.5, 2.5, 2.5) * 1.5;
    vec3 distances = abs(ray) - sides/2.0;
    float maxDistance = max(distances.x, max(distances.y, distances.z));


    return min(maxDistance, 0.0) + length(max(maxDistance, 0.0));
}

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

float gem(in vec3 ray)
{
    float heightFactor = 1.3;
    vec3 coneBase = vec3(0.35, 0.4, 1.64 * heightFactor);
    vec3 coneBase2 = vec3(0.4, 0.3, 2.5 * heightFactor);


    // your magical distance function
    float cone1 = cone(vec3(ray.x, ray.y - coneBase.z , ray.z), coneBase);
    float cone2 = cone(vec3(ray.x, -ray.y - coneBase2.z, ray.z), coneBase2);
    return max(ray.y - coneBase.z * 0.3, min(cone1, cone2));

}



float getDistance (in vec3 ray) {
    //return sdSphere(rayForRepeat(ray, REPEAT_PATTERN), 2.0);
    ray = rayForRepeat(ray, REPEAT_PATTERN);
    return conebination(ray);
}

vec3 getBackground(in vec3 start, in vec3 direction, in vec3 origin, in vec2 fragCoord) {
    //return texture(iChannel0, origin).rgb;
    //return texture(iChannel0, direction).rgb;

    vec3 coord = normalize(direction);

    float elevation = abs(asin(coord.x));
    float modifier = (1.0 - 0.5 * length(coord.xz));
    //elevation *= modifier;
    float rotation = atan2(coord.y, coord.z) + iTime * 0.5;
    bool rail;
    bool row = mod(rotation, M_PI / 20.0) < 0.04;
    float beg = 0.045;
    float off = 0.008;
    bool col = elevation < M_PI * (beg + off) && elevation > M_PI * (beg - off);
    vec3 ret;
    if(elevation < M_PI * 0.07  && (row || col) ){
        ret = vec3(1.0,0.8,1.0);
    } else {
        ret = vec3(0.4,0.1,0.05) * (M_PI / 2.0 - abs(elevation - M_PI * beg));
    }


    vec3 grain = normalize(vec3(
        mod(coord.x * coord.y + (1.0 - coord.z), 0.03 * mod(iTime, 0.016)),
        mod(coord.x * coord.z + (1.0 - coord.y), 0.07 * mod(iTime, 0.007)),
        mod(coord.z * coord.y + (1.0 - coord.x), 0.05 * mod(iTime, 0.011))
    ));

    return  ret + mod(mod(vec3(10000000.0 * iTime), abs(grain)), 0.0013) * 100.0;
}

vec3 getFragmentColor (in vec3 origin, in vec3 direction, in vec2 fragCoord) {
    vec3 originalDirection = direction;
    vec3 startPoints[MAX_RAYS];
    vec3 directions[MAX_RAYS];
    vec3 luminocities[MAX_RAYS];
    vec2 data[MAX_RAYS];
    vec3 rgb = vec3(0.0, 0.0, 0.0);

    vec3 reflectFactor   = vec3(1.0, 0.9, 0.9) * 0.5;
    vec3 refractFactor   = vec3(1.0, 1.0, 1.0) * 1.0;
    vec3 airAbsorbFactor = vec3(0.003, 0.003, 0.001) * 0.05; // random values for now;
    vec3 gemAbsorbFactor = vec3(0.1, 0.6, 0.6) * 0.5; //random values for now;

    vec3 normal;
    vec3 reflection;
    vec3 refraction;
    vec3 luminocityA;
    vec3 luminocityB;

    float side = (getDistance(origin) <= 0.0) ? -1.0 : 1.0;
    float refractionIndex = 1.0 / REFRACTION;

    directions[0]   = direction;
    startPoints[0]  = origin;
    luminocities[0] = vec3(1.0, 1.0, 1.0) * 1.2;
    data[0] = vec2(1.0, refractionIndex);

    bool outOfSteps = false;
    bool outOfTime = false;

    for(int ray = 0; ray < MAX_RAYS; ray++){
        vec3 luminocity = luminocities[ray];

        vec3 startPoint = startPoints[ray];
        vec3 direction  = normalize(directions[ray]);


        vec2 delta      = vec2 (DELTA, 0.0);
        side            = data[ray].x;
        refractionIndex = data[ray].y;

        outOfSteps = false;
        outOfTime = false;
        bool modelOuted = (side == 1.0);

        if(max(luminocity.z, max(luminocity.x, luminocity.y)) >= 2.0 * DELTA){
            outOfSteps = true;

            int totalSteps = 0;
            float rayLength = 0.0;


            for(int step=0; step < MAX_SUBREFRACTIONS; step++){
                float totalMinimum = RAY_LENGTH_MAX;
                float prevRayLength = rayLength;
                // have not reached max refractions
                outOfSteps = false;

                // launch ray


                float distance = RAY_LENGTH_MAX;
                for (int rayStep = 0; rayStep < MAX_STEPS; ++rayStep) {
                    outOfSteps = false;
                    distance = side * getDistance (startPoint) * 0.8;
                    totalMinimum = min(distance, totalMinimum);
                    float distMin = max (distance, DELTA);
                    rayLength += distMin;
                    outOfTime = false;

                    if (distance <= 0.0 || ((rayLength > RAY_LENGTH_MAX || totalSteps >= MAX_STEPS))) {

                        outOfSteps = totalSteps > MAX_STEPS;
                        outOfTime = (rayLength > RAY_LENGTH_MAX || outOfSteps);
                        break;
                    }
                    startPoint += direction * distMin;
                    totalSteps += 1;
                    outOfSteps = true;
                }

                float subRayLength = rayLength - prevRayLength;


                // decrease luminocity removing absorb_value * traveled_distance
                //luminocity -= -rayLength * absorber;

                vec3 absorber = (modelOuted ? airAbsorbFactor : gemAbsorbFactor);
                //float boost = (1.5 - min(totalMinimum, 1.5)) * 0.04;
                //boost += (1.0 - min(totalMinimum, 1.0)) * 0.07;
                //boost += (0.5 - min(totalMinimum, 0.5)) * 0.1;
                //boost = boost * boost;
                luminocity = max(luminocity - absorber * subRayLength, 0.0);//  + ((side == 1.0 && distance > 0.0) ? boost : 0.0 );
                //luminocity = luminocity * (-rayLength * absorber + 1.0);

                if(outOfTime) break;

                normal = conebinationNormal(startPoint) * side;

                reflection = reflect (direction, normal);
                refraction = refract (direction, normal, refractionIndex);

                // next intersection is new surface
                if (modelOuted) break;

                // detect total internal reflection
                if (dot(refraction, refraction) <= 2.0 * DELTA) {
                    // ray is reflected
                    direction = reflection;
                    startPoint += DELTA * MAGICSTEPS * direction;
                } else {
                    // ray is refracted
                    direction = refraction;

                    // next distance is negative
                    // and normal is reverted
                    side = -side;

                    // next refraction index is reveted
                    refractionIndex = 1.0/refractionIndex;

                    // exiting model internal reflection
                    modelOuted = (side == 1.0);
                }

                // have reached max refractions
                outOfSteps = true;
            }
        }

        bool rayEnd =  outOfTime || max(luminocity.z, max(luminocity.x, luminocity.y)) < 2.0 * DELTA ;
        float reflectMultiplier = pow(max (0.0, REFLECTION_POWER * dot(reflection, direction)), 0.3);
        vec3 lumen = min((reflectFactor * reflectMultiplier), 1.0);
        luminocityA = rayEnd ? luminocity * (USE_ADDITIONAL_RAYS ? luminocity : luminocity * 0.0) : lumen * luminocity;
        luminocityB = rayEnd ? vec3(0.0, 0.0, 0.0) : luminocity * (refractFactor);

        if(rayEnd){
            // still traversing, didn't hit nothing
            reflection  = direction;
            refraction  = direction;
        }

        if(USE_ADDITIONAL_RAYS && backupI < MAX_RAYS && luminocities[backupI] == vec3(0.0, 0.0, 0.0)) {
            directions[backupI]   = reflection;
            startPoints[backupI]  = startPoint + reflection * DELTA;
            luminocities[backupI] = luminocity;
            data[backupI] = vec2(1.0, 1.0/REFRACTION);
        }

        if(refractionI < MAX_RAYS) {
            directions[reflectionI]   = reflection;
            startPoints[reflectionI]  = startPoint + reflection * DELTA * MAGICSTEPS;
            luminocities[reflectionI] = luminocityA;
            data[reflectionI] = vec2(1.0, 1.0/REFRACTION);
        }
        if(refractionI < MAX_RAYS) {
            directions[refractionI]   = refraction;
            startPoints[refractionI]  = startPoint; //+ refraction * DELTA;
            luminocities[refractionI] = luminocityB;
            data[refractionI] = vec2(-1.0, REFRACTION);
        }

        if(refractionI >= MAX_RAYS || (!USE_ADDITIONAL_RAYS && rayEnd)){
            vec3 finisher = (side == 1.0) ? airAbsorbFactor : gemAbsorbFactor;
            //vec3 col = (direction * 0.3 + vec3(1.0,1.0,1.0)*0.7) * (luminocity  - 500.0 * finisher);
            //vec3 col = (direction * 0.3 + vec3(1.0,1.0,1.0)*0.7);
            //rgb = vec3(max(rgb.x, col.x), max(rgb.y, col.y), max(rgb.z, col.z));

            //vec3 col = texture(iChannel0, direction).rgb;
            vec3 col;


            col = getBackground(startPoint, direction, originalDirection, fragCoord);
                //col = abs(direction - 0.5);
                //col = (direction * 0.3 + vec3(1.0,1.0,1.0)*0.7);
                //   col = (direction * 0.2 + vec3(1.0,1.0,1.0)*0.7);

            col *= (luminocity  - luminocity * 1000.0 * finisher);
            rgb = max(rgb, col);
        }
    }

    rgb = clamp(rgb, 0.18, 1.0);
    return  rgb;
}

void mainImage (out vec4 fragColor, in vec2 fragCoord) {

    // Define the ray corresponding to this fragment
    vec2 frag = (2.0 * fragCoord.xy - iResolution.xy) / iResolution.y;
    vec3 direction = normalize (vec3 (frag, 2.0));

    // Set the camera
    vec3 origin = vec3 ((15.0 * cos (iTime * 0.1)), 10.0 * sin (iTime * 0.2), 15.0 * sin (iTime * 0.1));
    vec3 forward = -origin;
    vec3 up = vec3 (sin (iTime * 0.5), 2.0, 0.0);
    mat3 rotation;
    rotation [2] = normalize (forward);
    rotation [0] = normalize (cross (up, forward));
    rotation [1] = cross (rotation [2], rotation [0]);
    direction = rotation * direction;

    // Set the fragment color
    fragColor = vec4 (getFragmentColor (origin, direction, frag), 1.0);
}