// https://www.shadertoy.com/view/tlcGWX
//  ROOF OF THE WORLD
#pragma iChannel0 rgba_noise.png
#pragma iChannel1 BufferA

//Ray marching terrain based on several examples: 
//https://iquilezles.org/www/articles/terrainmarching/terrainmarching.htm
//http://iquilezles.org/www/articles/normalsSDF/normalsSDF.htm
//https://iquilezles.org/www/articles/outdoorslighting/outdoorslighting.htm
//https://iquilezles.org/www/articles/fog/fog.htm
//https://www.shadertoy.com/view/lslfRN
//https://www.shadertoy.com/view/MdX3Rr
//http://www.kevs3d.co.uk/dev/shaders/
//https://www.decarpentier.nl/scape-procedural-basics

//Terrain:	FBM of noise with some octaves inverted. Amplitude and frequency changes are from 
//			experimenting and examples. A rotation each iteration gives more variety
//			as seen in multiple examples. Normals are found by sampling around the 
//			intersection point. The Y component is set to 2.0*eps. 
//			The normal sampling uses more FBM iterations to get sharper looking results and 
//			the sampling offset depends on the ray distance
//Marching:	Steps are limited by number and maximum distance. The step size is a fraction
//			of the difference between the last sample point height and the terrain height
//			at those coordinates. The fraction can be changed to better capture ridges 
//			depending on the terrain parameters. Once a sample point is under the terrain,
//			a bisection method is used to move the sample backwards and forwards to find
//			the 0 crossing in a binary search. Return value within EPSILON of surface or
//			best approximation after a set number of iterations.
//Textures:	The rock texture is noise FBM with offset height-based stripes. Close range rock
//			also has fine noise texture. Snow is a plain white colour on slopes facing 
//			upward. The snow boundary is offset by noise to avoid straight lines.
//Lighting:	Three lights based on iq article. Strong white light from the sun. Soft shadows
//			using ray marching. Sky light is a blueish light directly from above. A third
//			white light is used to simulate indirect lighting in a very bright environment.
//			Includes tinted fog to give a sense of scale.



const int MAX_STEPS = 250;
const float MAX_DIST = 3000.0;
const float EPSILON = 1e-4;
const float PI = 3.1415;
const int SHADOW_STEPS = 64;
//e.g. 8.0: soft, 128: hard
const float SHADOW_SHARPNESS = 16.0;
const float HEIGHT = 200.0;
const float SCALE = 0.005;

const int terrainLimit = 5;
const int normalLimit = 9;
const int cameraLimit = 2;

const vec3 skyColour = 0.3 * vec3(0.09, 0.33, 0.81);
const vec3 sunLightColour = vec3(1);
const vec3 sunColour = vec3(1.0, 1.0, 0.8);
const vec3 specularColour = vec3(1);

const float speed = 10.0;

//In a circle of 2*PI
const float sunLocation = 1.633;
//0: horizon, 1: zenith
const float sunHeight = 0.2;

vec3 rayDirection(float fieldOfView, vec2 fragCoord) {
    vec2 xy = fragCoord - iResolution.xy / 2.0;
    float z = iResolution.y / tan(radians(fieldOfView) / 2.0);
    return normalize(vec3(xy, -z));
}

//https://www.geertarien.com/blog/2017/07/30/breakdown-of-the-lookAt-function-in-OpenGL/
mat3 lookAt(vec3 camera, vec3 targetDir, vec3 up){
  vec3 zaxis = normalize(targetDir);    
  vec3 xaxis = normalize(cross(zaxis, up));
  vec3 yaxis = cross(xaxis, zaxis);

  return mat3(xaxis, yaxis, -zaxis);
}

//Darken sky when looking up
vec3 getSkyColour(vec3 rayDir){
    return mix(skyColour, 0.2*skyColour, rayDir.y);
}

//By iq
float noised( in vec2 x ){
    vec2 f = fract(x);
    vec2 u = f*f*(3.0-2.0*f);
  
    vec2 p = floor(x);
	float a = textureLod( iChannel0, (p+vec2(0.5,0.5))/256.0, 0.0 ).x;
	float b = textureLod( iChannel0, (p+vec2(1.5,0.5))/256.0, 0.0 ).x;
	float c = textureLod( iChannel0, (p+vec2(0.5,1.5))/256.0, 0.0 ).x;
	float d = textureLod( iChannel0, (p+vec2(1.5,1.5))/256.0, 0.0 ).x;
    
	float res = (a+(b-a)*u.x+(c-a)*u.y+(a-b-c+d)*u.x*u.y);
    res = res - 0.5;
    return res;
}

float fbm(vec3 pos, int limit){
    float res = 0.0;
    float freq = 1.72;
    float amp = 0.5;
    for(int i = 0; i < 9; i++){ 
        if(i == limit){break;}
        if(i < 3){
        	res += (1.0-abs(noised(freq*0.5*pos.xz)))*amp*2.5;
        }else{ 
        	res += noised(freq*pos.xz)*amp;
        }
        freq *= 2.5;
        amp *= 0.25;
        pos.xz *= mat2(1, -.75, .75, 1);
    }
	return res;
}

//Get height of terrain at xz coordinates
float getHeight(vec3 pos, int limit){
    return HEIGHT*fbm(SCALE*pos, limit);
}

//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 = getHeight(p, terrainLimit);
        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 = getHeight(p, terrainLimit);
        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;
}

//http://iquilezles.org/www/articles/normalsSDF/normalsSDF.htm
//Making the normal sample distance depend on the ray length leads to less noise
vec3 getNormal(vec3 p, float t){
    float eps = 0.001 * t;

    return normalize(vec3( 
        getHeight(vec3(p.x-eps, p.y, p.z), normalLimit) 
        - getHeight(vec3(p.x+eps, p.y, p.z), normalLimit),
        
        2.0*eps,
        
        getHeight(vec3(p.x, p.y, p.z-eps), normalLimit) 
        - getHeight(vec3(p.x, p.y, p.z+eps), normalLimit) 
    ));
}

//https://www.iquilezles.org/www/articles/rmshadows/rmshadows.htm
float softShadow(vec3 pos, vec3 rayDir){
    float res = 1.0;
    float t = 1.0;
    //float ph = 1e10;
    //Start some small distance away from the surface to avoid artifacts
    pos += rayDir * 5.0 * t;
	for(int i = 0; i < SHADOW_STEPS; i++){
	    vec3 p = pos + t * rayDir;
        if(p.y > 2.0*HEIGHT){
        	break;
        }
        float h = p.y - getHeight(p, terrainLimit);
		res = min(res, SHADOW_SHARPNESS * h / t );
        /*
		//An improved shadow approach that didn't quite work
        float y = h*h/(2.0*ph);
        float d = sqrt(h*h-y*y);
        res = min( res, SHADOW_SHARPNESS*d/max(0.0,t-y) );
        ph = h;
		*/
		t += h;
        if(res < EPSILON){
            break;
        }
	}
	return clamp(res, 0.0, 1.0);
}

float rockFbm(vec3 p, vec3 normal, float dist){
    
    float res = 0.0;
    float freq = 0.005;
    float amp = 0.5;
    
    //When close to the rock, add fine detail. 
    //Showing it at all distances introduces noise and shimmering.
    if(dist < 80.0){
    	float detail = smoothstep(75.0, 0.0, dist);
    	res += detail*noised(40.0*p.xz);
    }
    
    //Offset position height in a large wavy pattern
    p.y += noised(0.02*p.xz) * 150.0;
    
    //Add several horizontal lines (offset by the above) for layers
    for(int i = 0; i < 4; i++){
    	res += texture(iChannel0, freq * p.yy).r * amp;
        freq *= 1.5;
        amp *= 0.5;
    }    
    
    return clamp(res, 0.0, 1.0);
}
//Return colour of surface fragment based on light information
vec3 shading(vec3 position, vec3 normal, vec3 rayDir, float dist, vec3 lightDirection){
    
    float ambientStrength = 0.15;
    
    float specularStrength = 0.0;
    float shininess = 32.0;
    
    vec3 ambientColour = vec3(0);
    vec3 diffuseColour = vec3(0);
    
    //Offset the height of the position to give snow line a more organic border.
    float height = position.y + noised( 0.2 * position.xz) * 14.0;
	bool snow = false;
    //Snow occurs on horizontal enough areas
    if(((normal.y > 0.65)) || ((height > HEIGHT*1.55) && (normal.y > 0.5))){
        snow = true;
    	ambientColour = vec3(1);
        diffuseColour = vec3(1);
        //Make snow subtly specular
        specularStrength = 0.2;
    }
    
    if(!snow){
        //Colour rock by layering noise. 
    	//Add dark brown and grey lines and finer detail in close proximity
    	vec3 rockColour = mix(0.25*vec3(0.3, 0.2, 0.1), vec3(0.8), 
                              0.5 * rockFbm(position, normal, dist));
    	rockColour = mix(0.5*rockColour, vec3(0.3), 0.5);
    
        //Add darker blotches 
        if(normal.y < 0.5){
            rockColour *= 0.8;
        } 

        //Ambient colour slightly darker than diffuse
        ambientColour = 0.75*rockColour;
        diffuseColour = rockColour;
    }
    
	vec3 halfwayDir = normalize(lightDirection - rayDir);  
	float spec = pow(max(dot(normal, halfwayDir), 0.0), shininess);

	//Colour of light sharply reflected into the camera
	vec3 specular = spec * specularColour * sunLightColour; 
	
  	//Path to light blocked     
    float shadow = softShadow(position, lightDirection);
    
    //No specular reflections in shadow
    if(shadow < 0.75){
    	specularStrength = 0.0;
    }
    
	vec3 result = vec3(0.0); 
    
	//How much a fragment faces the sun
	float sun = max(dot(normal, lightDirection), 0.0);
    //Main sun
    vec3 sunLight = shadow * sun * sunLightColour;
    
    //How much the fragment faces up
    float sky = max(dot(normal, vec3(0,1,0)), 0.0);
    //Sky light. A blue light from directly above. A lighter blue than the sky.
	vec3 skyLight = sky * vec3(0.12, 0.29, 0.55);
    
    //Indirect light. A white light opposite to the sun direction simulates global 
    //illumination. Snow is very diffusely reflective and so valleys are not as dark. 
    //Also gives a nice scattering look to large snow banks.
    float indirect = max(dot(normal, normalize(vec3(-1.0, 0.0, -1.0) * lightDirection)), 0.0);
   	vec3 indirectLight = indirect * vec3(1.0);
    
    //Combine light
    result += 1.2 * sunLight;
    result += 0.2 * skyLight;
	result += 0.3 * indirectLight;
    
    //Light and material interaction
    result *= diffuseColour;
    result += ambientStrength * ambientColour + specularStrength * specular;
    
    
    return  result;
}

//https://iquilezles.org/www/articles/fog/fog.htm
vec3 applyFog(vec3  rgb, float dist, vec3 rayOri, vec3 rayDir, vec3 sunDir){
    //Make horizon more hazy
    if(dist == MAX_DIST){dist = 4000.0;}
    if(abs(rayDir.y) < 0.0001){rayDir.y = 0.0001;}
    //Rate of fade
    float b = 0.014;
    float fogAmount = 1.0 * exp(-rayOri.y*b) * (1.0-exp(-dist*rayDir.y*b))/rayDir.y;
    float sunAmount = max( dot( rayDir, sunDir ), 0.0 );
    vec3  fogColor  = mix( vec3(0.5,0.6,0.7), vec3(1.0), pow(sunAmount, 8.0) );
    return mix(rgb, fogColor, fogAmount);
}

float getGlow(float dist, float radius, float intensity){
    dist = max(dist, 1e-6);
	return pow(radius/dist, intensity);	
}

//https://knarkowicz.wordpress.com/2016/01/06/aces-filmic-tone-mapping-curve/
//Thanks for the suggestion loicvdb
vec3 ACESFilm(vec3 x){
    float a = 2.51;
    float b = 0.03;
    float c = 2.43;
    float d = 0.59;
    float e = 0.14;
    return clamp((x*(a*x+b))/(x*(c*x+d)+e), 0.0, 1.0);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ){
    //Get the default direction of the ray (along the negative Z direction)
    vec3 rayDir = rayDirection(45.0, fragCoord);
    vec4  mouse  = (iMouse / iResolution.xyxy) * 2.0 - 1.0;
    if(iMouse.xy == vec2(0.0)){
		mouse.x += 1.0;
        mouse.y += 0.9;
    }
    mouse.x *= PI;
    //----------------- Define a camera -----------------
    
    vec3 cameraPos = vec3(120.0, 40.0, -400.0-iTime * speed);
    cameraPos.y = getHeight(cameraPos, cameraLimit) + 30.0;
    
    vec3 lightDirection = normalize(vec3(sin(sunLocation), sunHeight, cos(sunLocation)));
	
    vec3 targetDir = vec3(sin(mouse.x), mouse.y, -cos(mouse.x));
    vec3 up = vec3(0.0, 1.0, 0.0);
    
    //---------------------------------------------------
    
    //Get the view matrix from the camera orientation
    mat3 viewMatrix = lookAt(cameraPos, targetDir, up);
    //Transform the ray to point in the correct direction
    rayDir = normalize(viewMatrix * rayDir);
    
    //Find the distance to where the ray stops
    float dist = getIntersection(cameraPos, rayDir, MAX_DIST);
    vec3 col;
    
    bool skyVisible = false;
    if(dist == MAX_DIST){
        skyVisible = true;
    	col = getSkyColour(rayDir);
    }else{
    	vec3 position = cameraPos + rayDir * dist;
    	vec3 normal = getNormal(position, dist);
		col = shading(position, normal, rayDir, dist, lightDirection);
    }
    
    col = applyFog(col, dist, cameraPos, rayDir, lightDirection);
    
    //Display the sun as a glow in the light direction
    if(skyVisible){
    	float mu = dot(rayDir, lightDirection);
    	col += sunColour*getGlow(1.0-mu, 0.00015, 0.9);
    }
   
    //Tonemapping
    col = ACESFilm(col);

    //Gamma correction 1.0/2.2 = 0.4545...
    col = pow(col, vec3(0.4545));
    
    //Output to screen
    fragColor = vec4(col, 1.0);
}