/*****************************************************************************************
 *                                                                                       *
 * OpenSpace                                                                             *
 *                                                                                       *
 * Copyright (c) 2014-2020                                                               *
 *                                                                                       *
 * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
 * software and associated documentation files (the "Software"), to deal in the Software *
 * without restriction, including without limitation the rights to use, copy, modify,    *
 * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
 * permit persons to whom the Software is furnished to do so, subject to the following   *
 * conditions:                                                                           *
 *                                                                                       *
 * The above copyright notice and this permission notice shall be included in all copies *
 * or substantial portions of the Software.                                              *
 *                                                                                       *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
 * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
 ****************************************************************************************/
#version __CONTEXT__

#include "atmosphere_common.glsl"

layout(location = 0) out vec4 renderTarget1;
layout(location = 1) out vec4 renderTarget2;

uniform float r;
uniform vec4 dhdH;

//uniform sampler2D transmittanceTexture;

void integrand(const float r, const float mu, const float muSun, const float nu, 
                const float y, out vec3 S_R, out vec3 S_M) {
  // The integral's integrand is the single inscattering radiance:
  // S[L0] = P_M*S_M[L0] + P_R*S_R[L0]
  // where S_M[L0] = T*(betaMScattering * exp(-h/H_M))*L0 and
  // S_R[L0] = T*(betaRScattering * exp(-h/H_R))*L0.
  // T = transmittance.
  // One must remember that because the occlusion on L0, the integrand
  // here will be equal to 0 in that cases.
  // Also it is important to remember that the phase function for the
  // Rayleigh and Mie scattering are added during the rendering time
  // to increase the angular precision
  S_R = vec3(0.0);
  S_M = vec3(0.0);
  
  // cosine law
  float ri = max(sqrt(r * r + y * y + 2.0 * r * mu * y), Rg);
  
  // Considering the Sun as a parallel light source,
  // thew vector s_i = s.
  // So muSun_i = (vec(y_i) dot vec(s))/r_i = ((vec(x) + vec(yi-x)) dot vec(s))/r_i
  // muSun_i = (vec(x) dot vec(s) + vec(yi-x) dot vec(s))/r_i = (r*muSun + yi*nu)/r_i
  float muSun_i = (nu * y + muSun * r) / ri;

  // If the muSun_i is smaller than the angle to horizon (no sun radiance
  // hitting the point y), we return S_R = S_M = 0.0f.
  if (muSun_i >= -sqrt(1.0 - Rg * Rg / (ri * ri))) {
    // It's the transmittance from the point y (ri) to the top of atmosphere
    // in direction of the sun (muSun_i) and the transmittance from the observer
    // at x (r) to y (ri).
    vec3 transmittanceY = transmittance(r, mu, y) * transmittanceLUT(ri, muSun_i);
    // exp(-h/H)*T(x,v)
    if (ozoneLayerEnabled) {
      S_R = (exp(-(ri - Rg) / HO) + exp( -(ri - Rg) / HR )) * transmittanceY;
      S_M = exp( -(ri - Rg) / HM ) * transmittanceY;
    }
    else {
      S_R = exp( -(ri - Rg) / HR ) * transmittanceY;
      S_M = exp( -(ri - Rg) / HM ) * transmittanceY;
    }
    // The L0 (sun radiance) is added in real-time.
  }
}

void inscatter(const float r, const float mu, const float muSun, const float nu,
               out vec3 S_R, out vec3 S_M) {
  // Let's calculate S_M and S_R by integration along the eye ray path inside
  // the atmosphere, given a position r, a view angle (cosine) mu, a sun
  // position angle (cosine) muSun, and the angle (cosine) between the sun position
  // and the view direction, nu.
  // Integrating using the Trapezoidal rule:
  // Integral(f(y)dy)(from a to b) = (b-a)/2n_steps*(Sum(f(y_i+1)+f(y_i)))
  S_R = vec3(0.0f);
  S_M = vec3(0.0f);
  float rayDist = rayDistance(r, mu);
  float dy =  rayDist / float(INSCATTER_INTEGRAL_SAMPLES);
  float yi = 0.0f;
  vec3 S_Ri;
  vec3 S_Mi;
  integrand(r, mu, muSun, nu, 0.0, S_Ri, S_Mi);
  for (int i = 1; i <= INSCATTER_INTEGRAL_SAMPLES; ++i) {
    float yj = float(i) * dy;
    vec3 S_Rj;
    vec3 S_Mj;
    integrand(r, mu, muSun, nu, yj, S_Rj, S_Mj);
    S_R += (S_Ri + S_Rj);
    S_M += (S_Mi + S_Mj);
    yi = yj;
    S_Ri = S_Rj;
    S_Mi = S_Mj;
  }
  S_R *= betaRayleigh * (rayDist / (2.0f * float(INSCATTER_INTEGRAL_SAMPLES)));
  S_M *= betaMieScattering * (rayDist / (2.0f * float(INSCATTER_INTEGRAL_SAMPLES)));
}

void main(void) {
  vec3 S_R; // First Order Rayleigh InScattering 
  vec3 S_M; // First Order Mie InScattering
  float mu, muSun, nu; // parametrization angles

  // From the layer interpolation (see C++ code for layer to r)
  // and the textures parameters (uv), we unmapping mu, muSun and nu.
  unmappingMuMuSunNu(r, dhdH, mu, muSun, nu);
  
  // Here we calculate the single inScattered light.
  // Because this is a single inscattering, the light
  // that arrives at a point y in the path from the
  // eye to the infinity (top of atmosphere or planet's
  // ground), comes only from the light source, i.e., the
  // sun. So, the there is no need to integrate over the
  // whole solid angle (4pi), we need only to consider
  // the Sun position (cosine of sun pos = muSun).
  // Then, following the paper notation:
  // S[L] = P_R*S_R[L0] + P_M*S_M[L0] + S[L*]
  // For single inscattering only:
  // S[L0] = P_R*S_R[L0] + P_M*S_M[L0]
  // In order to save memory, we just store the red component
  // of S_M[L0], and later we use the proportionality rule
  // to calcule the other components.
  inscatter(r, mu, muSun, nu, S_R, S_M);
  renderTarget1 = vec4(S_R, 1.0);
  renderTarget2 = vec4(S_M, 1.0);
}
