source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/SunColor.cpp @ 2957

#include "SunColor.h"
#include "Vector3.h"
#include "common.h"


using namespace CHCDemoEngine;
using namespace std;


inline float CBQ(float x)
{
        return x * x * x;
}


inline float SQR(float x)
{
        return x * x ;
}


void SunColor::Compute(const Vector3 &sunDir, Vector3 &ambient, Vector3 &diffuse) const
{
        const float turbidity = 3.0f;

        // sunDir is sun direction
        // ambient color: shadow color
        // diffuse color: sun color

        float sun_theta[2];
        
        sun_theta[0] = acos(sunDir.z);

        const float cos_theta = cos(sun_theta[0]);
        sun_theta[1] = cos_theta * cos_theta;

        Vector3 zenithColor;

        zenithColor.x = ( 0.00165f * CBQ(sun_theta[0]) - 0.00374f * SQR(sun_theta[0]) + 0.00208f * sun_theta[0] + 0.0f)     * SQR(turbidity) +  
                                    (-0.02902f * CBQ(sun_theta[0]) + 0.06377f * SQR(sun_theta[0]) - 0.03202f * sun_theta[0] + 0.00394f) * turbidity +
                                        ( 0.11693f * CBQ(sun_theta[0]) - 0.21196f * SQR(sun_theta[0]) + 0.06052f * sun_theta[0] + 0.25885f);

        zenithColor.y = (0.00275f  * CBQ(sun_theta[0]) - 0.00610f * SQR(sun_theta[0]) + 0.00316f * sun_theta[0] + 0.0f)     * SQR(turbidity) +
                        (-0.04214f * CBQ(sun_theta[0]) + 0.08970f * SQR(sun_theta[0]) - 0.04153f * sun_theta[0] + 0.00515f) * turbidity +
                            ( 0.15346f * CBQ(sun_theta[0]) - 0.26756f * SQR(sun_theta[0]) + 0.06669f * sun_theta[0] + 0.26688f);
        
        zenithColor.z  = (float)((4.0453f * turbidity - 4.9710f) *      tan((4.0f / 9.0f - turbidity / 120.0f) * 
                                                         (M_PI - 2.0f * sun_theta[0])) - 0.2155f * turbidity + 2.4192f);

        // convert kcd/m² to cd/m²
        zenithColor.z *= 1000.0f;

        ambient = zenithColor;

        Vector3 zenith_XYZ;                                                                                             

        zenith_XYZ.x = (ambient.x / ambient.y) * ambient.z;                                                                        
        zenith_XYZ.y = ambient.z;                                                                                          
        zenith_XYZ.z = ((1.0f - ambient.x - ambient.y) / ambient.y) * ambient.z;                                                                

        ambient.x =  3.240479f * zenith_XYZ.x - 1.537150f * zenith_XYZ.y - 0.498535f * zenith_XYZ.z;
        ambient.y = -0.969256f * zenith_XYZ.x + 1.875992f * zenith_XYZ.y + 0.041556f * zenith_XYZ.z;    
        ambient.z =  0.055648f * zenith_XYZ.x - 0.204043f * zenith_XYZ.y + 1.057311f * zenith_XYZ.z;

        // downscale ambient color
        ambient *= 5e-5f;

        // simulate the sun intensity by modulating the ambient term.
        ambient *= (1.0f - 0.9f * DotProd(sunDir, Vector3::UNIT_Z()));
        ambient += Vector3(0.2f);

        float ABCDE_x[5], ABCDE_y[5], ABCDE_Y[5];

        ABCDE_x[0] = -0.01925f * turbidity - 0.25922f;
        ABCDE_x[1] = -0.06651f * turbidity + 0.00081f;
        ABCDE_x[2] = -0.00041f * turbidity + 0.21247f;
        ABCDE_x[3] = -0.06409f * turbidity - 0.89887f;
        ABCDE_x[4] = -0.00325f * turbidity + 0.04517f;

        ABCDE_y[0] = -0.01669f * turbidity - 0.26078f;
        ABCDE_y[1] = -0.09495f * turbidity + 0.00921f;
        ABCDE_y[2] = -0.00792f * turbidity + 0.21023f;
        ABCDE_y[3] = -0.04405f * turbidity - 1.65369f;
        ABCDE_y[4] = -0.01092f * turbidity + 0.05291f;

        ABCDE_Y[0] =  0.17872f * turbidity - 1.46303f;
        ABCDE_Y[1] = -0.35540f * turbidity + 0.42749f;
        ABCDE_Y[2] = -0.02266f * turbidity + 5.32505f;
        ABCDE_Y[3] =  0.12064f * turbidity - 2.57705f;
        ABCDE_Y[4] = -0.06696f * turbidity + 0.37027f;

        Vector3 num;
          
        num.x = (1.0f + ABCDE_x[0] * exp( ABCDE_x[1] / sunDir.z )) * (1.0f + ABCDE_x[2])+ ABCDE_x[4];    
        num.y = (1.0f + ABCDE_y[0] * exp( ABCDE_y[1] / sunDir.z )) * (1.0f + ABCDE_y[2])+ ABCDE_y[4];   
        num.z = (1.0f + ABCDE_Y[0] * exp( ABCDE_Y[1] / sunDir.z )) * (1.0f + ABCDE_Y[2])+ ABCDE_Y[4];   


        Vector3 den;

        den.x = (1.0f + ABCDE_x[0] * exp( ABCDE_x[1] )) * (1.0f + ABCDE_x[2] * exp( ABCDE_x[3] * sun_theta[0] ) + ABCDE_x[4] * sun_theta[1]);  
        den.y = (1.0f + ABCDE_y[0] * exp( ABCDE_y[1] )) * (1.0f + ABCDE_y[2] * exp( ABCDE_y[3] * sun_theta[0]) + ABCDE_y[4] * sun_theta[1]);   
        den.z = (1.0f + ABCDE_Y[0] * exp( ABCDE_Y[1] )) * (1.0f + ABCDE_Y[2] * exp( ABCDE_Y[3] * sun_theta[0] ) + ABCDE_Y[4] * sun_theta[1]);   


        Vector3 xyY = (num / den) * zenithColor;    

        Vector3 XYZ;                                                                                             

        XYZ.x = (xyY.x / xyY.y) * xyY.z;                                                                        
        XYZ.y = xyY.z;                                                                                          
        XYZ.z = ((1.0f - xyY.x - xyY.y) / xyY.y) * xyY.z;                                                                


        Vector3 color;  
        color.x = 3.240479f*XYZ.x -1.537150f*XYZ.y-0.498535f*XYZ.z;
        color.y = -0.969256f*XYZ.x +1.875992f*XYZ.y +0.041556f*XYZ.z;   
        color.z = 0.055648f*XYZ.x -0.204043f*XYZ.y +1.057311f*XYZ.z;

        // Calculate final sun diffuse color.
        diffuse = color * 0.00017f;
}