/////////////////////////////////////////////////
//
// Sky.sl v0.1
// 
// written by Mike Andrews <m.c.andrews@rdg.ac.uk>
//
// 26th May 2003
//
// for use with PovMan v1.0 by Vahur Krouverk <vkrouverk@starman.ee>
//
// This file implements a model for clear sky colour and luminance.
//
// It is based on Skylight.inc v0.6b by Philippe Debar <phdebar@yahoo.fr>
// which is based on:
//
// Sources : Preetham
//           "A practical analytical model for daylight"
//           http://www.cs.utah.edu/vissim/papers/sunsky/index.html

//           Charles Poynton
//           Color and Gamma FAQ
//           http://www.inforamp.net/~poynton/ColorFAQ.html
//           http://www.inforamp.net/~poynton/GammaFAQ.html 

//           Adrian Ford & Alan Roberts
//           Colour Space Conversions
//           Can be found at Poynton's website : http://www.inforamp.net/~poynton/PDFs/coloureq.pdf

//           Colourware
//           Color Physic Faq
//           http://www.colourware.co.uk/cpfaq.htm
//
//
// notes:
//
//    north is fixed to -z
//
// to do:
//
//    need to re-integrate with Jaimes Vives Piqueres' light macros
//    add a model for aerial perspective
//    add a model for cloudy sky

// Preetham's uses :
// Perez, Seals, J.M. and Ineichen model for sky luminance distribution

  // T: theta, angle in radians
  // G: gamma, angle in radians
  // C[]: distribution coefficients
  // C[0]: darkening or brightening of the horizon?
  // C[1]: luminance gradient near the horizon?
  // C[2]: realtive intensity of the circumsolar region?
  // C[3]: width of the circumsolar region?
  // C[4]: relative backscattered light?
float PSJI_F(float T, G, C[5])
{
  return (1 + C[0]*exp(C[1]/cos(T)) ) * (1+ C[2]*exp(C[3]*G) + C[4]*pow(cos(G),2) );
}

  // Yz: zenith luminance (or chromacity)
  // T:  theta, angle (in radians) from vertical to observed direction
  // G:  gamma, angle (in radians) between observed direction and sun
  // Ts: theta_s, angle (in radians) from vertical to sun direction
  // C[]: distribution coefficients
float PSJI_Y(float Yz, T, G, Ts, C[5])
{
  return Yz * PSJI_F(T, G, C) / PSJI_F(0, Ts, C);
}

// Preetham's distribution coefficients

// For Y luminance
void Preetham_Y(uniform float CT; output uniform float C[5])
{
  C[0] = 0.1787 * CT - 1.4630;
  C[1] = -0.3554 * CT + 0.4275;
  C[2] = -0.0227 * CT + 5.3251;
  C[3] = 0.1206 * CT - 2.5771;
  C[4] = -0.0670 * CT + 0.3703;
}

// For x chromacity
void Preetham_x(uniform float CT; output uniform float C[5])
{
  C[0] = -0.0193 * CT - 0.2592;
  C[1] = -0.0665 * CT + 0.0008;
  C[2] = -0.0004 * CT + 0.2125;
  C[3] = -0.0641 * CT - 0.8989;
  C[4] = -0.0033 * CT + 0.0452;
}

// For y chromacity
void Preetham_y(uniform float CT; output uniform float C[5])
{
  C[0] = -0.0167 * CT - 0.2608;
  C[1] = -0.0950 * CT + 0.0092;
  C[2] = -0.0079 * CT + 0.2102;
  C[3] = -0.0441 * CT - 1.6537;
  C[4] = -0.0109 * CT + 0.0529;
}

// Preetham's zenith values

// Zenith luminance
// T turbidity, Ts angle vertical-sun
float Preetham_Yz(uniform float T, Ts)
{ 
  return (4.0453*T-4.9710)*tan((4/9-T/120)*(PI-2*Ts)) - 0.2455*T + 2.4192;
}

// Zenith x chromacity
// T turbidity, Ts angle vertical-sun
float Preetham_xz(uniform float T, Ts)
{
  float xz, T2 = T*T;
  
  xz = Ts*Ts*Ts *( 0.00166*T2 -0.02903*T +0.11693 ) 
      +Ts*Ts    *(-0.00375*T2 +0.06377*T -0.21196 )
      +Ts       *( 0.00209*T2 -0.03202*T +0.06052 )
      +          (             0.00394*T +0.25886 );
      
  return xz;
}

// Zenith y chromacity
// T turbidity, Ts angle vertical-sun
float Preetham_yz(uniform float T, Ts)
{ 
  float yz, T2 = T*T;
  
  yz = Ts*Ts*Ts *( 0.00275*T2 -0.04214*T +0.15346 ) 
      +Ts*Ts    *(-0.00610*T2 +0.08970*T -0.26756 )
      +Ts       *( 0.00317*T2 -0.04153*T +0.06670 )
      +          (             0.00516*T +0.26688 );
      
  return yz;
}

color Yxy_to_RGB(float Y, xx, yy)
{
  return color "XYZ" (max(Y*xx/yy,0), max(Y,0), max(Y*(1-xx-yy)/yy,0));
}

color SkyCol(
  point PP;
  uniform point SolarPosition;
  uniform float Current_Ts, UnderHorizonCorrection, Horizon_Epsilon, flIntensity_Final;
  uniform float P_Y, P_x, P_y, Current_Yz, Current_xz, Current_yz;
)
{
  uniform point y = point(0,1,0);
  float Theta, Y_Mult, Gamma;
  
  Theta = acos(min(1, y.normalize(PP)));
  
  if ( (Theta>(PI/2-Horizon_Epsilon)) && (UnderHorizonCorrection != 0) )
  {
    Y_Mult = pow(sin(Theta),10);  // set a light attenuation
    Theta = PI/2-Horizon_Epsilon; // keep horizon color
  }
  else
  {
    Y_Mult = 1; // else do nothing
  }
  Gamma = acos(min(1,(normalize(SolarPosition).normalize(PP))));
  
  return Yxy_to_RGB(
          // Y luminance
          Y_Mult* // Correction for what happens under the horizon
          flIntensity_Final* // Intensity correction to get vuable results
          PSJI_Y(
            Current_Yz,
            Theta,
            Gamma,
            Current_Ts,
            P_Y
          ),
          
          // x chrominance
          PSJI_Y(
            Current_xz,
            Theta,
            Gamma,
            Current_Ts,
            P_x
          ),
          
          // y chrominance
          PSJI_Y(
            Current_yz,
            Theta,
            Gamma,
            Current_Ts,
            P_y
          )
      );
}

surface Sky(
  output float Al = 35, Az = 40;
  float Intensity_Mult = .75, Int_Sun_weight = 2, Int_Zenith_weight = 1, Int_Horizon_weight = 1;
  float Current_Turbidity = 3.;
  float UnderHorizonCorrection = 1, Horizon_Epsilon = 1e-3;
  float UseSolarPosition = 0;
  output vector SolarPosition = vector(0,3/5,4/5);
  output color SunColor = color(1,1,1), ZenithColor = color(0,0,0), HorizonColor = color(0,0,0);
)
{
  uniform float P_Y[5], P_x[5], P_y[5];
  
  uniform float Current_Ts, Current_Azimuth, flIntensity_Final;
  uniform float Current_Yz, Current_xz, Current_yz;
  uniform point origin = point(0,0,0), x = point(1,0,0), y = point(0,1,0), z = point(0,0,1);
  uniform vector xz = vector(1,0,1);
  
  float Theta, Y_Mult, Gamma;
  
  Preetham_Y(Current_Turbidity, P_Y);
  Preetham_x(Current_Turbidity, P_x);
  Preetham_y(Current_Turbidity, P_y);
  
  if (UseSolarPosition != 0) 
  {
    Az = degrees(atan(xcomp(SolarPosition), zcomp(SolarPosition)));
    Al = degrees(atan(ycomp(SolarPosition), length(SolarPosition*xz)));
    Current_Ts = radians(Al);
    Current_Azimuth = radians(Az);
  }
  else
  {
    Current_Ts = radians(Al);
    Current_Azimuth = radians(Az);
    SolarPosition = rotate(z,-Current_Ts, origin, x);
    SolarPosition = rotate(SolarPosition, Current_Azimuth, origin, y);
  }
  
  Current_Yz=Preetham_Yz(Current_Turbidity,PI/2-Current_Ts);
  Current_xz=Preetham_xz(Current_Turbidity,PI/2-Current_Ts);
  Current_yz=Preetham_yz(Current_Turbidity,PI/2-Current_Ts);
  
  SunColor = SkyCol(SolarPosition, SolarPosition, Current_Ts, UnderHorizonCorrection, Horizon_Epsilon, 1,
                    P_Y, P_x, P_y, Current_Yz, Current_xz, Current_yz);
  ZenithColor = SkyCol(y, SolarPosition, Current_Ts, UnderHorizonCorrection, Horizon_Epsilon, 1,
                    P_Y, P_x, P_y, Current_Yz, Current_xz, Current_yz);
  HorizonColor = SkyCol(rotate(-z,-Current_Ts, origin, x), SolarPosition, Current_Ts, UnderHorizonCorrection, Horizon_Epsilon, 1,
                    P_Y, P_x, P_y, Current_Yz, Current_xz, Current_yz);

  flIntensity_Final= Intensity_Mult * ( Int_Sun_weight + Int_Zenith_weight + Int_Horizon_weight )/
    (Int_Sun_weight*comp(SunColor,1) + Int_Zenith_weight*comp(ZenithColor,1) + Int_Horizon_weight*comp(HorizonColor,1));

  Ci = SkyCol(normalize(P), SolarPosition, Current_Ts, UnderHorizonCorrection, Horizon_Epsilon, flIntensity_Final,
                    P_Y, P_x, P_y, Current_Yz, Current_xz, Current_yz);
}