/*
    [References]:
        Kuzminykh, A. (2021). Physically Based Real-Time Rendering of the Moon. https://elib.dlr.de/203152/1/Bachelorarbeit_Alexander_Kuzminykh_20210818.pdf
*/

float Hapke_2012(vec3 normal, vec3 viewDirection, vec3 lightDirection, float w) {
    const float b       = 0.23955;
    const float c       = 0.30452;
    const float h_S     = 0.07145;
    const float h_C     = 1.0;
    const float B_S0    = 1.80238;
    const float B_C0    = 0.0;
    const float theta_p = 23.4;
    const float phi     = 0.3;

    float cos_g = dot(viewDirection, lightDirection);
    float cos_i = saturate(dot(normal, lightDirection));
    float cos_e = saturate(dot(normal, viewDirection ));

    float g = acos(cos_g), i = acos(cos_i), e = acos(cos_e);

    float sin_i = sin(i), sin_e = sin(e);

    // Calculate Psi by projecting L and V on a plane with N as the normal and getting
    // the cosine of the angle between the projections.
    //float Psi = acos((cos_g - cos_e * cos_i) / (sin_e * sin_i));
    float cos_Psi = dot(normalize(lightDirection - cos_i * normal), normalize(viewDirection - cos_e * normal));
          cos_Psi = clamp(cos_Psi, -1.0, 1.0);
    
    float Psi     = acos(cos_Psi);
    float PsiHalf = Psi * 0.5;

    float f_Psi        = exp(-2.0 * tan(PsiHalf));
    float sin_PsiHalf  = sin(PsiHalf);
    float sin2_PsiHalf = sin_PsiHalf * sin_PsiHalf;
    float PsiPerPI     = Psi * RCP_PI;

    float tan_theta_p = tan(theta_p * deg_to_rad);

    // The tan_theta_p is zero when theta_p is zero.
    // A zero for theta_p is somewhat unrealistic.
    float cot_theta_p  = 1.0 / tan_theta_p;
    float cot2_theta_p = cot_theta_p * cot_theta_p;

    float tan_i = tan(i);
    float tan_e = tan(e);

    // Here, tan_i and tan_e are zero each when i and e are zero each.
    // Because 1 / 0.0 is positive infinity we should check the usages.
    // The results are only used as factors for exponential functions,
    // where the whole expressions for the arguments are negated and
    // therefore result in negative infinity. Thus overall, these
    // exponential functions simply result in zero. No fix needed.
    float cot_i = 1.0 / tan_i;
    float cot_e = 1.0 / tan_e;

    float cot2_i = cot_i * cot_i;
    float cot2_e = cot_e * cot_e;

    float E_1_i = exp(-2.0 / PI * cot_theta_p  * cot_i );
    float E_1_e = exp(-2.0 / PI * cot_theta_p  * cot_e );
    float E_2_i = exp(-1.0 / PI * cot2_theta_p * cot2_i);
    float E_2_e = exp(-1.0 / PI * cot2_theta_p * cot2_e);

    float chi_theta_p = 1.0 / sqrt(1.0 + PI * tan_theta_p * tan_theta_p);
    float eta_i = chi_theta_p * (cos_i + sin_i * tan_theta_p * E_2_i / (2.0 - E_1_i));
    float eta_e = chi_theta_p * (cos_e + sin_e * tan_theta_p * E_2_e / (2.0 - E_1_e));

    float my_0e = chi_theta_p;
    float my_e  = chi_theta_p;

    float S;
    if (i <= e) {
        my_0e *= cos_i + sin_i * tan_theta_p * (cos_Psi * E_2_e + sin2_PsiHalf * E_2_i) / (2.0 - E_1_e - PsiPerPI * E_1_i);
        my_e  *= cos_e + sin_e * tan_theta_p * (E_2_e - sin2_PsiHalf * E_2_i) / (2.0 - E_1_e - PsiPerPI * E_1_i);
        S      = my_e / eta_e * cos_i / eta_i * chi_theta_p / (1.0 - f_Psi + f_Psi * chi_theta_p * (cos_i / eta_i));
    } else {
        my_0e *= cos_i + sin_i * tan_theta_p * (E_2_i - sin2_PsiHalf * E_2_e) / (2.0 - E_1_i - PsiPerPI * E_1_e);
        my_e  *= cos_e + sin_e * tan_theta_p * (cos_Psi * E_2_i + sin2_PsiHalf * E_2_e) / (2.0 - E_1_i - PsiPerPI * E_1_e);
        S      = my_e / eta_e * cos_i / eta_i * chi_theta_p / (1.0 - f_Psi + f_Psi * chi_theta_p * (cos_e / eta_e));
    }

    float KphiTerm = 1.209 * pow(phi, 2.0 / 3.0);
    // This goes into complex numbers already within [0, 1].
    float K = -log(1.0 - KphiTerm) / KphiTerm;

    float gHalf           = g * 0.5;
    float tan_gHalf       = tan(gHalf);
    float tan_gHalfPerh_C = tan_gHalf / h_C;
    float B_C             = (1.0 + (1.0 - exp(-tan_gHalfPerh_C)) / tan_gHalfPerh_C) / (2.0 * pow(1.0 + tan_gHalfPerh_C, 2.0));

    // When g = 0, the division in the upper part causes a NaN.
    if (isnan(B_C)) B_C = 1.0;

    float r_0Term = sqrt(1.0 - w);
    float r_0     = (1.0 - r_0Term) / (1.0 + r_0Term);

    float LS = my_0e / (my_0e + my_e);
    float b2 = b * b;

    // An approximation from the publication.
    //c = 3.29 * exp(-17.4 * b2) - 0.908;
    float oneMinusb2 = 1.0 - b2;
    float twobcos_g  = 2.0 * b * cos_g;
    float p_g = (1.0 + c) / 2.0 * oneMinusb2 / pow(1.0 - twobcos_g + b2, 1.5) + (1.0 - c) / 2.0 * oneMinusb2 / pow(1.0 + twobcos_g + b2, 1.5);

    float B_S = 1.0 / (1.0 + tan_gHalf / h_S);

    float x_i = my_0e / K;
    float x_e = my_e  / K;
    float H_i = 1.0 / (1.0 - w * x_i * (r_0 + (1.0 - 2.0 * r_0 * x_i) / 2.0 * log((1.0 + x_i) / x_i)));
    float H_e = 1.0 / (1.0 - w * x_e * (r_0 + (1.0 - 2.0 * r_0 * x_e) / 2.0 * log((1.0 + x_e) / x_e)));

    float M = H_i * H_e - 1.0;

    return LS * K * w / (2.0 * TAU) * (p_g * (1.0 + B_S0 * B_S) + M) * (1.0 + B_C0 * B_C) * S / cos_i;
}
