mat3 rotXY(float rX, float rY){
 float cx = cos(rX),sx = sin(rX),cy = cos(rY),sy=sin(rY);
 return mat3(-sy, 0, -cy, -sx * cy, cx, sx * sy,cx * cy, sx, -cx * sy);
}
vec3 getdir(vec2 p,mat3 m){
  return normalize(vec3(p.x * m[0].x + p.y * m[1].x - m[2].x,
       p.x * m[0].y + p.y * m[1].y - m[2].y,
       p.x * m[0].z + p.y * m[1].z - m[2].z));
}
#define DELTA 0.025

float sam(vec2 uv){ return texture(iChannel0,uv).x; }
float fbm(vec2 p){
    return 0.5*sam(p)+ 0.35*sam(p*0.125)+ 0.125*sam(p*0.35);}

float hash(float n) { return fract(sin(n)*43758.5453123); }
float noise2(vec2 x){ vec2 p = floor(x),f = fract(x);
    f = f*f*(3.0-2.0*f);
    float n = p.x + p.y*17.0;
    return mix(mix( hash(n), hash(n+1.0),f.x),
               mix( hash(n+ 17.0), hash(n+ 18.0),f.x),f.y);
}
//inversed noise and texture together for heightmap
float map(vec3 p){
   float ns=fbm(p.xz*0.25)*0.15;
   float n=noise2(p.xz*0.3)*2.0-1.0;
   ns+=smoothstep(-0.65,0.6,n-0.28*ns);
   p.y+=ns*0.3;
   return p.y;
}
//large delta to avoid pixel flicking
vec3 normal(vec3 p){
    vec2 e=vec2(DELTA,-DELTA);
	return normalize(vec3(e.xyy*map(p+e.xyy)
         +e.yyx*map(p+e.yyx) +e.yxy*map(p+e.yxy) +e.x*map(p+e.x)));
}
#define NUM_STEPS 45
#define FAR 50.0
float trace(vec3 ori,vec3 dir, inout vec3 p) {
    float t = 0.0;
    float tot= 0.0;
    for(int i = 0; i < NUM_STEPS; i++) {
    	t = map(p);
		if(t < 0.02||t>FAR) break;
        tot+=t;
        p = ori + dir * tot;
    }
    return tot;
}

float n11(float p) { return fract(sin(p*154.101)*313.019);}
float n21(vec2 p) { return sin(dot(p, vec2(7., 157.)));}

float star(vec3 p) {
	vec2 gv = fract(p.xy) - 0.5;
	vec2 id = floor(p.xy);
	gv.x += sin(n21(id)*354.23) * 0.3;
	gv.y += sin(n11(n21(id))*914.19) * 0.3;
	float r = n11(n21(id));
	return 0.1*n11(r)*abs(sin(p.z+r*133.12))*0.2/length(gv)*0.1;
}

float stars(vec2 p) {
	float z = 1.0, m = 0.0,t=iTime*0.25;
	for(int i=0; i<5 ;i++){ z *= 2.0;m += star(vec3(p*z,1.0+t*2.0));}
	return m; }

float saturate(float x){ return clamp(x,0.0,1.0);}
//simplifyed sphere raytest
float raysphere(vec3 dir){
     float b = -10.0*dir.z;
     float discr = b * b - 99.0;
     if (discr < 0.0) return -1.0;
     return -b - sqrt(discr);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = fragCoord/iResolution.xy-0.5;
    uv.x*=iResolution.x/iResolution.y;
    float time= iTime*0.2;

    mat3 camrot=rotXY(-0.15,0.12);
    vec3 ro=camrot[2]*5.0;

    ro.y+=1.25;
    vec3 rd=getdir(uv,camrot);
    ro.x-=iTime*0.25;

    vec3 light=normalize(vec3(0.15,0.125,0.15));
    vec3 sky=vec3(0.05,0.09,0.1)*
          clamp(-(rd.y+0.85)/(ro.y-3.0),0.25,1.0);

    sky+=stars(uv);

    float d = length(uv-vec2(0.0,0.125));

    float land=smoothstep(0.25,0.75,fbm(2.2*uv+0.125));
    sky += exp(-5.0*d*d)*vec3(0.25,0.28,0.26);
    vec3 moon =vec3(0.2,0.4,0.6)+ vec3(0.9,0.6,0.3)*land;
    sky = mix(moon,sky, smoothstep(0.245,0.25,d) );

    vec3 raydir=normalize(vec3(uv.x,uv.y-0.125,2.465));
    float dist=raysphere(raydir);
   //earth
    if(dist>0.0){
      vec3 hit=(raydir*dist);
      hit.z-=10.0;

      float d2n=dot(hit,light)+0.65;
      float edge= 1.0-dot(hit,reflect(light,hit));
      edge=pow(edge,5.5);
      sky*=d2n;
      sky= mix(sky+land*(0.5-d2n)*vec3(1.5,1.2,0.65),vec3(0.72,0.825,0.9),edge*0.1);
    }

    vec3 col =vec3(0.0);
    vec3 pos=ro;
    //raystep matching is the slowest part...
    float t= trace(ro,rd,pos);

    if(t<FAR){
       vec3 n =normal(pos);
       light.xz=-light.xz;
       float d2l=max(dot(n,light),0.0);
       float fre = clamp(1.0+dot(light,n),0.0,1.0); // Fresnel reflection.
       float Schlick = pow(1.0- max(dot(rd, normalize(rd + light)),0.0),5.0);
       fre *= mix(0.2, 1.0, Schlick);//Hard clay.
       col+=d2l*vec3(1.25,1.22,1.2)+fre*vec3(0.65,0.725,0.8);
     }
     else col=sky;
  //blur the landscape edge to sky
   col=mix(col,sky,smoothstep(0.8,1.0,t/FAR));
   float uy=fragCoord.y/iResolution.y-0.5;
   float vig = (1.1-0.8*uy*uy)*(1.0-0.6*uv.x*uv.x);
   col*=col;
   fragColor = vec4(col*vig,1.0);
}