# Why is my sky color calculation in Mathematica incorrect?

I'm trying to implement an algorithm to calculate sky color based on this paper (Perez' model). Before I start programming the shader I wanted to test the concept in Mathematica. There are already some problems I cannot get rid off. Maybe someone already has implemented the algorithm.

I started with equations for the absolute zenital luminances `Yz`, `xz` and `yz` as proposed in the paper (page 22). The values for `Yz` seem to be reasonable. The following diagram shows `Yz` as a function of the zenital distance of the sun for a turbidity `T` of 5:

The function gamma(zenith,azimuth,solarzenith,solarazimuth) calculates the angle between a point with the given zenital distance and azimuth and the sun at the given position. This function seems to work, too. The following diagram shows this angle for `solarzenith=0.5` and `solarazimuth=0`. `zenith` grows from top down (0 to Pi/2), `azimuth` grows from left to right (-Pi to Pi). You can clearly see the position of the sun (the bright spot, angle becomes zero):

The Perez function (F) and coefficients have been implemented as given in the paper. Then the color values Yxy should be `absolute value * F(z, gamma) / F(0, solarzenith)`. I expect those values to be within the range [0,1]. However, this is not the case for the Y component (see update below for details). Here are some sample values:

``````{Y, x, y}
{19.1548, 0.25984, 0.270379}
{10.1932, 0.248629, 0.267739]
{20.0393, 0.268119, 0.280024}
``````

Here is the current result:

The Mathematica Notebook with all calculations can be found here and the PDF version here.

Does anyone have an idea what I have to change to get the same results as in the paper?

## C like code

``````// this function returns the zenital Y component for
// a given solar zenital distance z and turbidity T
float Yz(float z, float T)
{
return (4.0453 * T - 4.9710)*tan( (4.0f/9-T/120)*(Pi-2*z) ) - 0.2155 * T + 2.4192
}

// returns zenital x component
float xz(float z, float T)
{
return //matrix calculation, see paper
}

// returns zenital y component
float yz(float z, float T)
{
return //matrix calculation, see paper
}

// returns the rgb color of a Yxy color
Color RGB(float Y, float x, float y)
{
Matrix m; //this is a CIE XYZ -> RGB conversion matrix
Vector v;
v.x = x/y*Y;
v.y = Y;
v.z = (1-x-y)/y*Y;
v = M * v; //matrix-vector multiplication;
return Color ( v.x, v.y, v.z );
}

// returns the 5 coefficients (A-E) for the given turbidity T
float[5] CoeffY(float T)
{
float[5] result;
result[0] = 0.1787 * T - 1.4630;
result[1] = -0.3554 * T + 0.4275;
...
return result;
}

//same for Coeffx and Coeffy

// returns the angle between an observed point and the sun
float PerezGamma(float zenith, float azimuth, float solarzenith, float solarazimuth)
{
return acos(sin(solarzenith)*sin(zenith)*cos(azimuth-solarazimuth)+cos(solarzenith)*cos(zenith));
}

// evalutes Perez' function F
// the last parameter is a function
float Perez(float zenith, float gamma, float T, t->float[5] coeffs)
{
return (1+coeffs(T)[0] * exp(coeffs(T)[1]/cos(zenith)) *
(1+coeffs(T)[2] * exp(coeffs(T)[3]*gamma) +
coeffs(T)[4]*pow(cos(gamma),2))
}

// calculates the color for a given point
YxyColor calculateColor(float zenith, float azimuth, float solarzenith, float solarazimuth, float T)
{
YxyColor c;
float gamma = PerezGamma(zenith, azimuth, solarzenith, solarazimuth);
c.Y = Yz(solarzenith, T) * Perez(zenith, gamma, T, CoeffY) / Perez(0, solarzenith, T, CoeffY);
c.x = xz(solarzenith, T) * Perez(zenith, gamma, T, Coeffx) / Perez(0, solarzenith, T, Coeffx);
c.y = yz(solarzenith, T) * Perez(zenith, gamma, T, Coeffy) / Perez(0, solarzenith, T, Coeffy);
return c;
}

// draws an image of the sky
void DrawImage()
{
for(float z from 0 to Pi/2) //zenithal distance
{
for(float a from -Pi to Pi) //azimuth
{
YxyColor c = calculateColor(zenith, azimuth, 1, 0, 5);
Color rgb = RGB(c.Y, c.x, c.y);
setNextColor(rgb);
}
newline();
}
}
``````

## Solution

As promised I wrote a blog article about rendering the sky. You can find it here.

-
I suspect that more people here would be able to help you if you were to try to implement the algorithm in actual code (shader or otherwise) instead of in Mathematica. – Tetrad Mar 27 '13 at 22:16
There is a Mathematica SE, though you would have to check their FAQ to see if your question is on topic over there. – Byte56 Mar 27 '13 at 22:20
Well, the question is not really about Mathematica, but about the algorithm. I added the PDF version of the notebook, so everyone can read it. I'm sure that the syntax is comprehensible for a common programmer and probably more comprehensible than shader code. – Nico Schertler Mar 27 '13 at 22:58
@NicoSchertler: The problem is that I don't think many people in here understand Mathematica syntax. You will probably have more luck if you rewrite your code in a C-like or Python-like language, at least for the purposes of this question. – Panda Pajama Mar 29 '13 at 9:44
The question is really too localised and might get closed, but thanks for the paper link, it's interesting. – sam hocevar Mar 29 '13 at 10:01

There are two errors in the matrix used for `xz`: 1.00166 should be 0.00166, and 0.6052 should be 0.06052.