#version 400 core

#define S_PI 3.1415926535

out vec4 FragColor;

in vec3 rdi;
in vec3 rd;

uniform int samples;
uniform int samples_light;
uniform float light_strength;

uniform float planet_radius_L;
uniform float atmo_height_L;
uniform float height_offset_L;
uniform float hr_L;
uniform float hm_L;
uniform float mie_g_L;
uniform vec3 beta_R_L;
uniform vec3 beta_M_L;
uniform float atmo_height_1;
uniform float hr_1;
uniform float hm_1;
uniform float mie_g_1;
uniform vec3 beta_R_1;
uniform vec3 beta_M_1;

uniform float radius_1;
uniform vec3 pos_1;
uniform float rot_1;
uniform sampler2D tex_1;
uniform sampler2D norm_1;

uniform float radius_2;
uniform vec3 pos_2;
uniform float rot_2;
uniform sampler2D tex_2;
uniform sampler2D norm_2;

uniform float radius_3;
uniform vec3 pos_3;
uniform float rot_3;
uniform sampler2D tex_3;
uniform sampler2D norm_3;

uniform float radius_4;
uniform vec3 pos_4;
uniform float rot_4;
uniform sampler2D tex_4;
uniform sampler2D norm_4;

uniform float radius_5;
uniform vec3 pos_5;
uniform float rot_5;
uniform sampler2D tex_5;
uniform sampler2D norm_5;

uniform vec3 light_direction_A;
uniform vec3 light_direction_B;

vec3 renderPlanet(vec3 res, vec3 rd, vec3 pos, sampler2D tex, sampler2D norm, float rot, vec3 lightDir);
vec3 raySphereIntersection(vec3 pos, vec3 dir, float radius);
vec3 computeIncidentLight(vec3 a, vec3 b, vec3 ro, vec3 rd, vec3 lightDir, float lightStrength, vec3 planetPos, float planetRadius, float atmoRadius, float hr, float hm, float mieG, vec3 betaR, vec3 betaM);
vec2 xyzToLatlong(vec3 pos) {
	pos = normalize(pos);
	return vec2(
		acos(pos.y),
		atan(pos.x, pos.z)
	);
}

void main() {
    vec3 col = vec3(0, 0, 0);
    
    vec3 rda = normalize(rd);
    vec3 res1 = raySphereIntersection(pos_1, rda, radius_1);
    vec3 res2 = raySphereIntersection(pos_2, rda, radius_2);
    vec3 res3 = raySphereIntersection(pos_3, rda, radius_3);
    vec3 res4 = raySphereIntersection(pos_4, rda, radius_4);
    vec3 res5 = raySphereIntersection(pos_5, rda, radius_5);
    
    float closest = 10000000000000.0;
    int closestIdx = 0;
    if(res1.z > 0) {
		closest = res1.x;
		closestIdx = 1;
	}
	if(res2.z > 0 && res2.x < closest) {
		closest = res2.x;
		closestIdx = 2;
	}
	if(res3.z > 0 && res3.x < closest) {
		closest = res3.x;
		closestIdx = 3;
	}
	if(res4.z > 0 && res4.x < closest) {
		closest = res4.x;
		closestIdx = 4;
	}
	if(res5.z > 0 && res5.x < closest) {
		closest = res5.x;
		closestIdx = 5;
	}
	
	switch(closestIdx) {
		default:
			col = vec3(0,0,0);
			break;
		case 1:
			col = renderPlanet(res1, rda, pos_1, tex_1, norm_1, rot_1, -light_direction_B);
			break;
		case 2:
			col = renderPlanet(res2, rda, pos_2, tex_2, norm_2, rot_2, -light_direction_B);
			break;
		case 3:
			col = renderPlanet(res3, rda, pos_3, tex_3, norm_3, rot_3, -light_direction_B);
			break;
		case 4:
			col = renderPlanet(res4, rda, pos_4, tex_4, norm_4, rot_4, -light_direction_B);
			break;
		case 5:
			col = renderPlanet(res5, rda, pos_5, tex_5, norm_5, rot_5, -light_direction_B);
			break;
	}
	vec3 res = raySphereIntersection(pos_1, rda, radius_1 + atmo_height_1);
	if(res.z > 0 && res.x < closest) {
		vec3 atmoColor = computeIncidentLight(
			res1,
			res,
			vec3(0,0,0),
			rda,
			-light_direction_B,
			light_strength,
			pos_1,
			radius_1,
			radius_1 + atmo_height_1,
			hr_1,
			hm_1,
			mie_g_1,
			beta_R_1,
			beta_M_1
		);
		if(atmoColor.x < 10000) {
			col = (vec3(1, 1, 1) - col) * atmoColor + col;
		}
	}
    
    vec3 localAtmoColor = computeIncidentLight(
        vec3(0,0,0),
        vec3(0,0,0),
        vec3(0, 0, -planet_radius_L - height_offset_L * atmo_height_L),
        vec3(rdi.x, rdi.y, -abs(rdi.z)),
        -light_direction_A,
        light_strength * 2,
        vec3(0,0,0),
        planet_radius_L,
        planet_radius_L + atmo_height_L, //Atmo radius
        hr_L,
        hm_L,
        mie_g_L,
        beta_R_L,
        beta_M_L
    );
    if(localAtmoColor.x < 10000) {
        col = (vec3(1, 1, 1) - col) * localAtmoColor + col;
    }
    
    vec3 sRGBLo = col * 12.92;
    vec3 sRGBHi = (pow(max(abs(col), 1.192092896e-07), vec3(1.0 / 2.4, 1.0 / 2.4, 1.0 / 2.4)) * 1.055) - 0.055;
    col = vec3(col.r <= 0.0031308 ? sRGBLo.r : sRGBHi.r, col.g <= 0.0031308 ? sRGBLo.g : sRGBHi.g, col.b <= 0.0031308 ? sRGBLo.b : sRGBHi.b);
    FragColor = vec4(col, 1);
}

vec3 renderPlanet(vec3 res, vec3 rd, vec3 pos, sampler2D tex, sampler2D norm, float rot, vec3 lightDir) {
	vec3 isect = normalize((res.x * rd) - pos);
	
	vec2 uv = xyzToLatlong(isect).yx / vec2(S_PI * 2, S_PI);
	uv.x += rot;
	
	vec3 col = texture(tex, uv).rgb;
	vec3 normal = texture(norm, uv).xyz * 2 - 1;
		
	float alpha = rot * S_PI * 2;
	float sina = sin(alpha), cosa = cos(alpha);
	float wxn = normal.x * cosa + normal.z * sina;
	normal.z = -normal.x * sina + normal.z * cosa;
	normal.x = wxn;
	normal = normalize(normal);
	
	float lightLevel = dot(normal, lightDir);
	lightLevel = max(0.01, lightLevel);
	col *= lightLevel;
	return col;
}

vec3 raySphereIntersection(vec3 pos, vec3 dir, float radius) {
    float radiusSq = radius * radius;
    float tca = dot(dir, pos);
    
    float L_2 = dot(pos, pos);
    if(tca <= 0 && L_2 > radiusSq) return vec3(0, 0, -2);
    
    float d = L_2 - tca * tca;
    if(d >= radiusSq) return vec3(0, 0, -1);
    
    float thc = radiusSq - d;
    thc = sqrt(thc);
    
    return vec3(tca - thc, tca + thc, 1);
}

vec3 computeIncidentLight(vec3 a, vec3 b, vec3 ro, vec3 rd, vec3 lightDir, float lightStrength, vec3 planetPos, float planetRadius, float atmoRadius, float hr, float hm, float mieG, vec3 betaR, vec3 betaM) {
    float tmin = 0;
    float tmax = 10000000.0;
    
    rd = normalize(rd);
    vec3 temp3 = planetPos - ro;
    a = raySphereIntersection(temp3, rd, planetRadius);
    if(a.z > 0 && a.y > 0) tmax = max(0, a.x);
    
    if(b.z == 0) b = raySphereIntersection(temp3, rd, atmoRadius);
    if(b.z < 0 || b.y < 0) return vec3(10100, 10100, 10100);
    if(b.x > tmin && b.x > 0) tmin = b.x;
    if(b.y < tmax) tmax = b.y;
    
    float segmentLength = (tmax - tmin) / samples;
    float segmentLengthLight;
    float tCurrent = tmin;
    float tCurrentLight;
    vec3 sumR = vec3(0,0,0);
    vec3 sumM = vec3(0,0,0);
    float opticalDepthR = 0;
    float opticalDepthM = 0;
    
    float temp2 = mieG * mieG;
    float mu = dot(lightDir, rd);
    float temp = 1 + mu * mu;
    float phaseR = 3.0 / (16.0 * S_PI) * temp;
    float phaseM = 3.0 / (8.0 * S_PI) * ((1.0 - temp2) * temp) / ((2.0 + temp2) * pow(1.0 + temp2 - 2.0 * mieG * mu, 1.5));
    vec3 pos, posLight;
    float Hr_a = hr;
    float Hm_a = hm;
    betaR = betaR * 0.01;
    betaM = betaM * 0.01;
    
    int j;
    for(int i = 0; i < samples; i++) {
        pos = ro + (tCurrent + segmentLength * 0.5) * rd;
        
        a = raySphereIntersection(planetPos - pos, lightDir, atmoRadius);
        if(a.z < 0) {
            tCurrent += segmentLength;
            continue;
        }
        segmentLengthLight = a.y / samples_light, tCurrentLight = 0;
        float opticalDepthLightR = 0;
        float opticalDepthLightM = 0;
        vec3 startPosLight = pos;
        for(j = 0; j < samples_light; j++) {
            posLight = startPosLight + (tCurrentLight + segmentLengthLight * 0.5) * lightDir;
            float height = length(planetPos - posLight) - planetRadius;
            if(height < 0){ j = -5; break; }
            
            float hrLight = exp(-height / Hr_a) * segmentLengthLight;
            float hmLight = exp(-height / Hm_a) * segmentLengthLight;
            opticalDepthLightR += hrLight;
            opticalDepthLightM += hmLight;
            
            tCurrentLight += segmentLengthLight;
        }
        if(j != -5) {
            float height = length(planetPos - pos) - planetRadius;
            float hrr = exp(-height / Hr_a) * segmentLength;
            float hmm = exp(-height / Hm_a) * segmentLength;
            opticalDepthR += hrr;
            opticalDepthM += hmm;
            vec3 tau = betaR * (opticalDepthLightR * 700 + opticalDepthR) + betaM * (opticalDepthLightM * 700 + opticalDepthM);
            tau = exp(-tau);
            sumR += tau * hrr;
            sumM += tau * hmm;
        }
        tCurrent += segmentLength;
    }
    
    return vec3(sumR * betaR * phaseR + sumM * betaM * phaseM) * lightStrength * 10.0;
}
//what
