#define T iTime
#define pix 2.0/iResolution.y

//Terrain Noise
#define FREQUENCY 0.20
#define OCTAVES 5.0
#define AMPLITUDE 3.0
#define GAIN 0.6
#define LACUNARITY 2.0

//Raymarch
#define EPSILON 1e-4
#define MAX_STEPS 64
#define MAX_DIST 90.0

//FabriceNeyret2 + IQ hash
//https://www.shadertoy.com/view/fsKBzw
float hash( vec2 f ) {   
    uvec2 x = uvec2( floatBitsToUint(f.x), floatBitsToUint(f.y) ),
          q = 1103515245U * ( x>>1U ^ x.yx    );
    return float( 1103515245U * (q.x ^ q.y>>3U) ) / float(0xffffffffU);
}

//Hash without Sine by "Dave_Hoskins"
//https://www.shadertoy.com/view/4djSRW
vec2 hash21(float p)
{
	vec3 p3 = fract(vec3(p) * vec3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.xx+p3.yz)*p3.zy);
}

float vNoise(vec2 uv)
{
    vec2 index = floor(uv);
    uv = fract(uv);
    uv = smoothstep(0., 1., uv);
    
    float x1 = mix(hash(index), hash(index + vec2(1., 0.)), uv.x);
    float x2 = mix(hash(index + vec2(0., 1.)), hash(index + vec2(1)), uv.x);
    return mix(x1, x2, uv.y);
}


float fbm(vec2 uv) 
{
    float frequency = FREQUENCY; 
    float amplitude = AMPLITUDE;
    float value = 0.;
    for (float i = 0.; i < OCTAVES; i++)
    {
        value += vNoise(uv * frequency) * amplitude;
        amplitude *= GAIN;
        frequency *= LACUNARITY;
    }
    return value;
}

vec3 normal(vec3 p)
{
    float eps = 0.001;
    return normalize(vec3(fbm(vec2(p.x - eps, p.z)) - fbm(vec2(p.x + eps, p.z)),
                          2. * eps,
                          fbm(vec2(p.x, p.z - eps)) - fbm(vec2(p.x, p.z + eps))));
}


//Thanks alro! :)
//Binary search for 0 crossing given two points on either side of the surface
float bisection(vec3 start, vec3 rayDir, float near_, float far_){
    float midpoint = (far_ + near_) * 0.5;
    //Sample point
    vec3 p = vec3(0);
    float near = near_;
    float far = far_;
    float height = 0.0;
    //Difference between sample point and terrain heights
    float diff = 0.0;
    
    for(int i = 0; i < 8; i++){
        p = start + rayDir * midpoint;
        height = fbm(vec2(p.x, p.z));
        diff = p.y - height;
        
        if(abs(diff) < EPSILON){
        	break;
        }else{
            
            if(diff < EPSILON){
                //Point is below terrain
                //Search first half
                far = midpoint;
            }else{
                //Point is above terrain
                //Search second half
                near = midpoint;
            }
            midpoint = (far + near) * 0.5;
        }
    }
    return midpoint;
}

float getIntersection(vec3 start, vec3 rayDir, float maxDist){
	//Distance between sample points. Set according to previous sample
    float stepSize = 0.0;
    //Height of the terrain
    float height = 0.0;
    //Length of the ray
    float dist = 0.0;
    //Difference between sample point and terrain heights
    float diff = 0.0;
    
    for(int i = 0; i < MAX_STEPS; i++){
        //Sample point
        vec3 p = start + rayDir * dist;
        
        //The height of the terrain at the xz coordinates of the sample point
        height = fbm(vec2(p.x, p.z));
        diff = abs(p.y - height);
        //If sample point is close enough to the terrain, return distance
        if(diff < EPSILON){
            return dist;
        }
        //If height of sample point is less than the height of the terrain,
        //the ray has hit the terrain. Use bisection to find the 0 crossing
        if(p.y < height){
        	dist = bisection(start, rayDir, dist - stepSize, dist);
            return dist;
        }
        
        //Static step size misses features and leads to banding. 
        //Set the step size to a fraction of the distance above the terrain.
        //Could also have a small step size which increases with distance, giving 
        //detailed results close to the camera and reaching far. However, 
        //this approach is used in many shaders and seems to give best results
        stepSize = diff * 0.5;
        
        //Increment ray
        dist += stepSize;
        
        if(dist > MAX_DIST){
        	return MAX_DIST;
        }
    }
    return dist;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord - 0.5 * iResolution.xy)/iResolution.x;
    
    //camera tutorial by BigWIngs (ArtOfCode)
    //https://www.youtube.com/watch?v=PBxuVlp7nuM&t=130s
    float mX = iMouse.x/iResolution.x * 2.0 * 3.141592;
    float mY = iMouse.y/iResolution.y * 30.0;
    vec3 offset = vec3(T*2., 0., 0.);
    vec3 ro = vec3(sin(mX) * (30. -mY), mY + 10., cos(mX) * (30. - mY)) + offset;
    vec3 lp = vec3(0., 1., 0.) + offset;
    
    vec3 f = normalize(lp - ro);
    vec3 r = normalize(cross(vec3(0., 1., 0.), f));
    vec3 u = cross(f, r);
    vec3 rd = f + r * uv.x + u * uv.y;
    
    float dist = getIntersection(ro, rd, MAX_DIST);
    vec3 pp = ro + rd * dist;
    vec3 n = normal(pp);
    
    //scan distance relative to pulse origin
    float sD = distance(floor((lp.xz)/8.0)*8.0, vec2(pp.x, pp.z)); 
    
    vec3 col = vec3(0);
    //pulse
    float pR = mod(T*30., 120.); //pulseRadius
    float m1 = smoothstep(pR + 2., pR, sD); //mask
    col += smoothstep(pR - 10., pR, sD) * vec3(0.0, 0.7, 1.) * m1;
    
    //stripes
    col += mix(vec3(0.4, 0.0, 0.), vec3(0., 0.25, 0.55), smoothstep(0.4, 0.8, dot(n, vec3(0.,1.,0.)))) //color
         * smoothstep(0.85, 1., abs(fract(sD * 3.0) - 0.5) * 2.0) //stripes shape
         * m1 * smoothstep(pR - 80., pR - 40., sD); //masks
         
    //point
    //"The simplest 3D" by BigWIngs https://www.youtube.com/watch?v=dKA5ZVALOhs
    vec3 p = lp + vec3(0., 5., 0.);
    col += vec3(smoothstep(0.3, 0.3 - pix*20.0, length(cross(p - ro, rd))/length(rd)));
    
    //terrain
    col += vec3(fbm(vec2(pp.x, pp.z))) * GAIN/AMPLITUDE * smoothstep(8., 0., distance(p.xz, vec2(pp.x, pp.z))) * vec3(0.0, 0.1, 0.4);
    
    //fade scan with distance
    col = mix(col, vec3(0.0, -0.2, 0.), smoothstep(MAX_DIST - 50.0, MAX_DIST, dist));
    
    //background
    col +=vec3(0.0, 0.0, 0.15);
    fragColor = vec4(col, 1.);
}