// Persistence of Vision 3.7 Scene Include File
// File: SkySim.inc
// Desc: This macro creates a realistic sky_sphere based          
//       on the given sun position and sky conditions.  
// Date: 7th June 2013
// Auth: Scott Boham
//
// Further details:
// ----------------
// Algorithm based on that used in Stellarium (www.stellarium.org),
// see the source files Skylight.cpp/.hpp for details.
// For the algorithm explanation see "A Practical Analytic Model
// for Daylight" by A. J. Preetham, Peter Shirley and Brian Smits.
//
#ifndef(SKYSIM_TEMP)
#declare SKYSIM_TEMP = version;
#version 3.7;

#include "transforms.inc"

// The SkySim macro will create a sky_sphere object with accurately
// coloured pigment and correct brightness
// Parameters:
//  SunPos         Vector specifying the direction towards the sun,
//                 magnitude of the vector doesn't matter
//  UpVector       Vector specifying the vertical up direction
//  Turbidity      Value in the range 1.25 to 32 that specifies the
//                 "haziness" of the sky, values around 4-6 are normal
//  ExposureFactor A conversion factor between the luminance of the
//                 sky calculated (in cd/m2) to POV values. As the sky
//                 luminance can be very high (around 10000 cd/m2)
//                 then an ExposureFactor around 1e-5 will be needed to
//                 get a correct exposure when looking at the sky
#macro SkySim(SunPos , UpVector, Turbidity ,ExposureFactor)

// Check T is in the correct range
#if(Turbidity<1.25 | Turbidity>32)
 #error "The Turbidity value T should be in the range 1.25 to 32"
#end

// ensure SunPos is normalized
#local SunPos = vnormalize(SunPos);

// Rotate the SunPos so it's relative to z rather than UpVector
#local SunPos = vtransform( SunPos , Reorient_Trans(UpVector,z) );

#if(SunPos.z < 0 )
 #error "The sun should be above the horizon"
#end

// calculate Thetas, the angle between "up" and the sun
#local Thetas = pi/2 - asin(SunPos.z);

// Get the luminance and colour at the zenith
// note the return vector is in Yxy colour space, not rgb
#local zenithColor = computeZenithColor(Turbidity,Thetas); 

// Setup constant factors
#local AY = 0.2787*Turbidity- 1.0630;
#local BY =-0.3554*Turbidity+ 0.4275;
#local CY =-0.0227*Turbidity+ 6.3251;
#local DY = 0.1206*Turbidity- 2.5771;
#local EY =-0.0670*Turbidity+ 0.3703;
#local Ax =-0.0148*Turbidity- 0.1703;
#local Bx =-0.0664*Turbidity+ 0.0011;
#local Cx =-0.0005*Turbidity+ 0.2127;
#local Dx =-0.0641*Turbidity- 0.8992;
#local Ex =-0.0035*Turbidity+ 0.0453;
#local Ay =-0.0131*Turbidity- 0.2498;
#local By =-0.0951*Turbidity+ 0.0092;
#local Cy =-0.0082*Turbidity+ 0.2404;
#local Dy =-0.0438*Turbidity- 1.0539;
#local Ey =-0.0109*Turbidity+ 0.0531;
#local cos_thetas = SunPos.z;
#local term_x = zenithColor.x / ((1 + Ax * exp(Bx)) * (1 + Cx * exp(Dx*Thetas) + Ex * cos_thetas * cos_thetas));
#local term_y = zenithColor.y / ((1 + Ay * exp(By)) * (1 + Cy * exp(Dy*Thetas) + Ey * cos_thetas * cos_thetas));
#local term_Y = zenithColor.z / ((1 + AY * exp(BY)) * (1 + CY * exp(DY*Thetas) + EY * cos_thetas * cos_thetas));

// Helper functions needed
#local cosDistSun = function(x,y,z,sx,sy,sz ){x*sx+y*sy+z*sz}
#local distSun = function(x,y,z,sx,sy,sz){acos(x*sx+y*sy+z*sz)}
#local cosDistSun_q = function(x,y,z,sx,sy,sz){cosDistSun(x,y,z,sx,sy,sz)*cosDistSun(x,y,z,sx,sy,sz)}
#local oneOverCosZenithAngle = function(x,y,z){1/z}

// This is needed as you can't write SunPos.x in a function
#local sx = SunPos.x;
#local sy = SunPos.y;
#local sz = SunPos.z;

// Functions to return the colour (xy) and luminance (L) depending on
// position in the sky (variable) and sun position (constant)
#local Col_x = function(x,y,z){ term_x * (1 + Ax * exp(Bx*oneOverCosZenithAngle(x,y,z)))
		* (1 + Cx * exp(Dx*distSun(x,y,z,sx,sy,sz)) + Ex * cosDistSun_q(x,y,z,sx,sy,sz) )}
#local Col_y = function(x,y,z){ term_y * (1 + Ay * exp(By*oneOverCosZenithAngle(x,y,z)))
		* (1 + Cy * exp(Dy*distSun(x,y,z,sx,sy,sz)) + Ey * cosDistSun_q(x,y,z,sx,sy,sz) )}
#local Col_Y = function(x,y,z){ term_Y * (1 + AY * exp(BY*oneOverCosZenithAngle(x,y,z)))
		* (1 + CY * exp(DY*distSun(x,y,z,sx,sy,sz)) + EY * cosDistSun_q(x,y,z,sx,sy,sz) )}

// Estimate the maximum luminance in the sky (to allow the color_map scale to work below)
// This is done by calculating the luminance at a point near the sun and then multiplying by 100
// If you see strange colour banding in the sky then try using a number greater than 100
#local ClosePos = vnormalize(<SunPos.x+0.05,SunPos.y+0.05,SunPos.z>);
#local SkyMaxLuminance = Col_Y(ClosePos.x,ClosePos.y,ClosePos.z)*100;

// Standard functions to convert Yxy values to linear RGB
#local xyYtoR = function (x,y,Y){ (x/y* 3.2406  -1.5372 + (1-x-y)/y * -0.4986)*Y/SkyMaxLuminance }
#local xyYtoG = function (x,y,Y){ (x/y*-0.9689  +1.8758 + (1-x-y)/y *  0.0415)*Y/SkyMaxLuminance }
#local xyYtoB = function (x,y,Y){ (x/y* 0.0557  -0.2040 + (1-x-y)/y *  1.0570)*Y/SkyMaxLuminance }
#local Pos = function(x){ select(x,0,x) }

// Define three pigments for the RGB contributions
#local Pr =  pigment{function{Pos(xyYtoR(Col_x(x,y,abs(z)),Col_y(x,y,abs(z)),Col_Y(x,y,abs(z))))}
                     color_map{[0 rgb <0,0,0>][1 rgb <3,0,0>]}};
#local Pg =  pigment{function{Pos(xyYtoG(Col_x(x,y,abs(z)),Col_y(x,y,abs(z)),Col_Y(x,y,abs(z))))}
                     color_map{[0 rgb <0,0,0>][1 rgb <0,3,0>]}};
#local Pb =  pigment{function{Pos(xyYtoB(Col_x(x,y,abs(z)),Col_y(x,y,abs(z)),Col_Y(x,y,abs(z))))}
                     color_map{[0 rgb <0,0,0>][1 rgb <0,0,3>]}};

// Combine into one pigment in the sky_sphere
sky_sphere
{ 
 pigment{
   average pigment_map{[1 Pr][1 Pg][1 Pb]}    // Average RGB pigments 
   Reorient_Trans(z,UpVector)                 // Reorient back to UpVector pointing up
  }
 emission rgb ExposureFactor*SkyMaxLuminance 
}
#end // end macro

//-------------------------------------------------------------------------------------
#macro computeZenithColor(T,thetas)
 // Calculate luminance
 #local zenithLuminance = 1000 * ((4.0453*T - 4.9710) * tan( (0.4444 - T/120) * (pi-2*Thetas) ) - 0.2155*T + 2.4192);
 #if(zenithLuminance<=0) 
  #local zenithLuminance=1e-11;
 #end
 // Calculate colour
 #local thetas2 = thetas*thetas;
 #local thetas3 = thetas*thetas2;
 #local T2 = T*T;
 #local zx= ( 0.00216*thetas3 - 0.00375*thetas2 + 0.00209*thetas) * T2 +
	               (-0.02903*thetas3 + 0.06377*thetas2 - 0.03202*thetas + 0.00394) * T +
	               ( 0.10169*thetas3 - 0.21196*thetas2 + 0.06052*thetas + 0.25886);
 #local zy= ( 0.00275*thetas3 - 0.00610*thetas2 + 0.00317*thetas) * T2 +
	               (-0.04214*thetas3 + 0.08970*thetas2 - 0.04153*thetas + 0.00516) * T +
	               ( 0.14535*thetas3 - 0.26756*thetas2 + 0.06670*thetas + 0.26688);
 <zx,zy,zenithLuminance>
#end
//-------------------------------------------------------------------------------------


#version SKYSIM_TEMP;
#end