////////////////////////////////////////////////////////////////////////
//
//  H32_gendaylit.CPP - Perez Sky Model Calculation
//
//  Version:    1.00A
//
//  History:    09/10/01 - Created.
//				11/10/08 - Modified for Unix compilation.
//				11/10/12 - Fixed conditional __max directive.
//
//  Compilers:  Microsoft Visual C/C++ Professional V10.0    
//
//  Author:     Ian Ashdown, P.Eng.
//              byHeart Consultants Limited
//              620 Ballantree Road
//              West Vancouver, B.C.
//              Canada V7S 1W3
//              e-mail: ian_ashdown@helios32.com
//
//  References:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289.
//
//				Perez, R., R. Seals, and J. Michalsky. 1993.
//				“All-Weather Model for Sky Luminance Distribution -
//				Preliminary Configuration and Validation,” Solar Energy
//				50(3):235-245.
//
//				Perez, R., R. Seals, and J. Michalsky. 1993. "ERRATUM to
//				All-Weather Model for Sky Luminance Distribution -
//				Preliminary Configuration and Validation,” Solar Energy
//				51(5):423.
//
//  NOTE:		This program is a completely rewritten version of
//				gendaylit.c written by Jean-Jacques Delaunay (1994).
//
//  Copyright 2009-2011 byHeart Consultants Limited. All rights
//  reserved.
//
//  Redistribution and use in source and binary forms, with or without
//  modification, are permitted for personal and commercial purposes
//  provided that redistribution of source code must retain the above
//  copyright notice, this list of conditions and the following
//  disclaimer:
//
//    THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESSED OR IMPLIED
//    WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
//    OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//    DISCLAIMED. IN NO EVENT SHALL byHeart Consultants Limited OR
//    ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
//    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
//    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
//    USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
//    AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
//    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
//    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//    POSSIBILITY OF SUCH DAMAGE.
//
////////////////////////////////////////////////////////////////////////

// Zenith is along the Z-axis
// X-axis points east
// Y-axis points north

// Include files
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#if defined(__cplusplus)
  #undef __max
  #define __max(a, b) (((a) > (b)) ? (a) : (b))
#else
  #define __max    max
#endif

struct SkyPatch		// Sky patch
{
	double altitude;	// Patch center altitude (radians)
	double azimuth;		// Patch center azimuth (radians)
	double luminance;	// Patch luminance
};

char *progname;								// Program name
char errmsg[128];							// Error message buffer
const double GN_Pi = 3.141592654;			// Pi
const double DC_SolarConstantE = 1367.0;	// Solar constant W/m^2
const double DC_SolarConstantL = 127.5;		// Solar constant klux

double altitude;			// Solar altitude (radians)
double azimuth;				// Solar azimuth (radians)
double apwc;				// Atmospheric precipitable water content
double dew_point = 11.0;	// Surface dew point temperature (deg. C)
double diff_illum;			// Diffuse illuminance
double diff_irrad;			// Diffuse irradiance
double dir_illum;			// Direct illuminance
double dir_irrad;			// Direct irradiance
double julian_date = 150.0;	// Julian date
double norm_diff_illum;		// Normalized diffuse illuimnance
double perez_param[5];		// Perez sky model parameters
double sky_brightness;		// Sky brightness
double sky_clearness;		// Sky clearness
double solar_rad;			// Solar radiance
double sun_zenith;			// Sun zenith angle (radians)
double zlumin;				// Zenith luminance
int	input = 0;				// Input type

double CalcAirMass();
double CalcDiffuseIllumRatio( int );
double CalcDiffuseIrradiance();
double CalcDirectIllumRatio( int );
double CalcDirectIrradiance();
double CalcEccentricity();
double CalcPrecipWater( double );
double CalcRelHorzIllum( SkyPatch * );
double CalcRelLuminance( double, double );
double CalcSkyBrightness();
double CalcSkyClearness();
double DegToRad( double );
double RadToDeg( double );
int CalcSkyParamFromIllum();
int GetCategoryIndex();
void CalcPerezParam( double, double, double, int );
void CalcSkyPatchLumin( SkyPatch * );
void ComputeSky();
void UsageError( const char * );

// Sky patch angles
static const double SkyPatchAngles[145][2] =
{
	{ 84.0, 0.0 },
	{ 84.0, 12.0 },
	{ 84.0, 24.0 },
	{ 84.0, 36.0 },
	{ 84.0, 48.0 },
	{ 84.0, 60.0 },
	{ 84.0, 72.0 },
	{ 84.0, 84.0 },
	{ 84.0, 96.0 },
	{ 84.0, 108.0 },
	{ 84.0, 120.0 },
	{ 84.0, 132.0 },
	{ 84.0, 144.0 },
	{ 84.0, 156.0 },
	{ 84.0, 168.0 },
	{ 84.0, 180.0 },
	{ 84.0, 192.0 },
	{ 84.0, 204.0 },
	{ 84.0, 216.0 },
	{ 84.0, 228.0 },
	{ 84.0, 240.0 },
	{ 84.0, 252.0 },
	{ 84.0, 264.0 },
	{ 84.0, 276.0 },
	{ 84.0, 288.0 },
	{ 84.0, 300.0 },
	{ 84.0, 312.0 },
	{ 84.0, 324.0 },
	{ 84.0, 336.0 },
	{ 84.0, 348.0 },
	{ 72.0, 0.0 },
	{ 72.0, 12.0 },
	{ 72.0, 24.0 },
	{ 72.0, 36.0 },
	{ 72.0, 48.0 },
	{ 72.0, 60.0 },
	{ 72.0, 72.0 },
	{ 72.0, 84.0 },
	{ 72.0, 96.0 },
	{ 72.0, 108.0 },
	{ 72.0, 120.0 },
	{ 72.0, 132.0 },
	{ 72.0, 144.0 },
	{ 72.0, 156.0 },
	{ 72.0, 168.0 },
	{ 72.0, 180.0 },
	{ 72.0, 192.0 },
	{ 72.0, 204.0 },
	{ 72.0, 216.0 },
	{ 72.0, 228.0 },
	{ 72.0, 240.0 },
	{ 72.0, 252.0 },
	{ 72.0, 264.0 },
	{ 72.0, 276.0 },
	{ 72.0, 288.0 },
	{ 72.0, 300.0 },
	{ 72.0, 312.0 },
	{ 72.0, 324.0 },
	{ 72.0, 336.0 },
	{ 72.0, 348.0 },
	{ 60.0, 0.0 },
	{ 60.0, 15.0 },
	{ 60.0, 30.0 },
	{ 60.0, 45.0 },
	{ 60.0, 60.0 },
	{ 60.0, 75.0 },
	{ 60.0, 90.0 },
	{ 60.0, 105.0 },
	{ 60.0, 120.0 },
	{ 60.0, 135.0 },
	{ 60.0, 150.0 },
	{ 60.0, 165.0 },
	{ 60.0, 180.0 },
	{ 60.0, 195.0 },
	{ 60.0, 210.0 },
	{ 60.0, 225.0 },
	{ 60.0, 240.0 },
	{ 60.0, 255.0 },
	{ 60.0, 270.0 },
	{ 60.0, 285.0 },
	{ 60.0, 300.0 },
	{ 60.0, 315.0 },
	{ 60.0, 330.0 },
	{ 60.0, 345.0 },
	{ 48.0, 0.0 },
	{ 48.0, 15.0 },
	{ 48.0, 30.0 },
	{ 48.0, 45.0 },
	{ 48.0, 60.0 },
	{ 48.0, 75.0 },
	{ 48.0, 90.0 },
	{ 48.0, 105.0 },
	{ 48.0, 120.0 },
	{ 48.0, 135.0 },
	{ 48.0, 150.0 },
	{ 48.0, 165.0 },
	{ 48.0, 180.0 },
	{ 48.0, 195.0 },
	{ 48.0, 210.0 },
	{ 48.0, 225.0 },
	{ 48.0, 240.0 },
	{ 48.0, 255.0 },
	{ 48.0, 270.0 },
	{ 48.0, 285.0 },
	{ 48.0, 300.0 },
	{ 48.0, 315.0 },
	{ 48.0, 330.0 },
	{ 48.0, 345.0 },
	{ 36.0, 0.0 },
	{ 36.0, 20.0 },
	{ 36.0, 40.0 },
	{ 36.0, 60.0 },
	{ 36.0, 80.0 },
	{ 36.0, 100.0 },
	{ 36.0, 120.0 },
	{ 36.0, 140.0 },
	{ 36.0, 160.0 },
	{ 36.0, 180.0 },
	{ 36.0, 200.0 },
	{ 36.0, 220.0 },
	{ 36.0, 240.0 },
	{ 36.0, 260.0 },
	{ 36.0, 280.0 },
	{ 36.0, 300.0 },
	{ 36.0, 320.0 },
	{ 36.0, 340.0 },
	{ 24.0, 0.0 },
	{ 24.0, 30.0 },
	{ 24.0, 60.0 },
	{ 24.0, 90.0 },
	{ 24.0, 120.0 },
	{ 24.0, 150.0 },
	{ 24.0, 180.0 },
	{ 24.0, 210.0 },
	{ 24.0, 240.0 },
	{ 24.0, 270.0 },
	{ 24.0, 300.0 },
	{ 24.0, 330.0 },
	{ 12.0, 0.0 },
	{ 12.0, 60.0 },
	{ 12.0, 120.0 },
	{ 12.0, 180.0 },
	{ 12.0, 240.0 },
	{ 12.0, 300.0 },
	{ 0.0, 0.0 }
};

// Perez sky model coefficients
//
// Reference:	Perez, R., R. Seals, and J. Michalsky, 1993. "All-
//				Weather Model for Sky Luminance Distribution -
//				Preliminary Configuration and Validation," Solar
//				Energy 50(3):235-245, Table 1.
//
static const double PerezCoeff[8][20] =
{
	// Sky clearness (epsilon): 1.000 to 1.065
	{   1.3525,  -0.2576,  -0.2690,  -1.4366,   -0.7670,
	    0.0007,   1.2734,  -0.1233,   2.8000,    0.6004,
	    1.2375,   1.0000,   1.8734,   0.6297,    0.9738,
	    0.2809,   0.0356,  -0.1246,  -0.5718,    0.9938 },
    // Sky clearness (epsilon): 1.065 to 1.230
	{  -1.2219,  -0.7730,   1.4148,   1.1016,   -0.2054,
	    0.0367,  -3.9128,   0.9156,   6.9750,    0.1774,
		6.4477,  -0.1239,  -1.5798,  -0.5081,   -1.7812,
		0.1080,   0.2624,   0.0672,  -0.2190,   -0.4285 },
    // Sky clearness (epsilon): 1.230 to 1.500
	{  -1.1000,  -0.2515,   0.8952,   0.0156,    0.2782,
	   -0.1812, - 4.5000,   1.1766,  24.7219,  -13.0812,
	  -37.7000,  34.8438,  -5.0000,   1.5218,    3.9229,
	   -2.6204,  -0.0156,   0.1597,   0.4199,   -0.5562 },
    // Sky clearness (epsilon): 1.500 to 1.950
	{  -0.5484,  -0.6654,  -0.2672,   0.7117,   0.7234,
	   -0.6219,  -5.6812,   2.6297,  33.3389, -18.3000,
	  -62.2500,  52.0781,  -3.5000,   0.0016,   1.1477,
	    0.1062,   0.4659,  -0.3296,  -0.0876,  -0.0329 },
    // Sky clearness (epsilon): 1.950 to 2.800
	{  -0.6000,  -0.3566,  -2.5000,   2.3250,   0.2937,
	    0.0496,  -5.6812,   1.8415,  21.0000,  -4.7656 ,
	  -21.5906,   7.2492,  -3.5000,  -0.1554,   1.4062,
	    0.3988,   0.0032,   0.0766,  -0.0656,  -0.1294 },
    // Sky clearness (epsilon): 2.800 to 4.500
	{  -1.0156,  -0.3670,   1.0078,   1.4051,   0.2875,
	   -0.5328,  -3.8500,   3.3750,  14.0000,  -0.9999,
	   -7.1406,   7.5469,  -3.4000,  -0.1078,  -1.0750,
	    1.5702,  -0.0672,   0.4016,   0.3017,  -0.4844 },
    // Sky clearness (epsilon): 4.500 to 6.200
	{  -1.0000,   0.0211,   0.5025,  -0.5119,  -0.3000,
	    0.1922,   0.7023,  -1.6317,  19.0000,  -5.0000,
		1.2438,  -1.9094,  -4.0000,   0.0250,   0.3844,
		0.2656,   1.0468,  -0.3788,  -2.4517,   1.4656 },
    // Sky clearness (epsilon): 6.200 to ...
	{  -1.0500,   0.0289,   0.4260,   0.3590,  -0.3250,
	    0.1156,   0.7781,   0.0025,  31.0625, -14.5000,
	  -46.1148,  55.3750,  -7.2312,   0.4050,  13.3500,
	    0.6234,   1.5000,  -0.6426,   1.8564,   0.5636 }
};

// Perez irradiance component model coefficients
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289.

struct CategoryBounds
{
	double lower;	// Lower bound
	double upper;	// Upper bound
};

// Perez sky clearness (epsilon) categories (Table 1)
static const CategoryBounds SkyClearCat[8] =
{
	{ 1.000, 1.065 },	// Overcast
	{ 1.065, 1.230 },
	{ 1.230, 1.500 },
	{ 1.500, 1.950 },
	{ 1.950, 2.800 },
	{ 2.800, 4.500 },
	{ 4.500, 6.200 },
	{ 6.200, 12.00 }	// Clear
};

// Luminous efficacy model coefficients
struct ModelCoeff
{
	double a;
	double b;
	double c;
	double d;
};

// Diffuse luminous efficacy model coefficients (Table 4, Eqn. 7)
static const ModelCoeff DiffuseLumEff[8] =
{
	{  97.24, -0.46,  12.00,  -8.91 },
	{ 107.22,  1.15,   0.59,  -3.95 },
	{ 104.97,  2.96,  -5.53,  -8.77 },
	{ 102.39,  5.59, -13.95, -13.90 },
	{ 100.71,  5.94, -22.75, -23.74 },
	{ 106.42,  3.83, -36.15, -28.83 },
	{ 141.88,  1.90, -53.24, -14.03 },
	{ 152.23,  0.35, -45.27,  -7.98 }
};

// Direct luminous efficacy model coefficients (Table 4, Eqn. 8)
static const ModelCoeff DirectLumEff[8] =
{
	{  57.20, -4.55, -2.98, 117.12 },
	{  98.99, -3.46, -1.21,  12.38 },
	{ 109.83, -4.90, -1.71,  -8.81 },
	{ 110.34, -5.84, -1.99,  -4.56 },
	{ 106.36, -3.97, -1.75,  -6.16 },
	{ 107.19, -1.25, -1.51, -26.73 },
	{ 105.75,  0.77, -1.26, -34.44 },
	{ 101.18,  1.58, -1.10,  -8.29 }
};

int main( int argc, char *argv[])
{
	int  i;	// Loop index

	progname = argv[0];

	if (argc < 6)
		UsageError("arg count");

	altitude = DegToRad(atof(argv[1]));
	azimuth = DegToRad(atof(argv[2]));

	for (i = 3; i < argc; i++)
	{
		if (argv[i][0] == '-' || argv[i][0] == '+')
		{
			switch (argv[i][1])
			{
				case 'P':
					input = 0;	// Perez parameters: epsilon, delta
					sky_clearness = atof(argv[++i]);
					sky_brightness = atof(argv[++i]);
					break;
				case 'W':
					// Direct normal irradiance [W/m^2]
					// Diffuse horizontal irradiance [W/m^2]
					input = 1;
					dir_irrad = atof(argv[++i]);
					diff_irrad = atof(argv[++i]);
					break;
				case 'L':
					// Direct normal illuminance [lux]
					// Diffuse horizontal illuminance [lux]
					input = 2;
					dir_illum = atof(argv[++i]);
					diff_illum = atof(argv[++i]);
					break;
				default:
					sprintf(errmsg, "Unknown option: %s", argv[i]);
					UsageError(errmsg);
					break;
				}
		}
		else
			UsageError("bad option");
	}

	ComputeSky();	// Compute sky parameters

	exit(0);
}

// Degrees into radians
inline double DegToRad( double degrees )
{ return degrees * GN_Pi / 180.0; }

// Radiuans into degrees
inline double RadToDeg( double radians )
{ return radians * 180.0 / GN_Pi; }

// Compute sky parameters
void ComputeSky()
{
	SkyPatch sky_patch[145];	// Sky patch array
	int index;					// Category index

	// Calculate atmospheric precipitable water content
	apwc = CalcPrecipWater(dew_point);

	// Limit minimum altitude
	if (altitude > DegToRad(87.0))
		altitude = DegToRad(87.0);

	// Calculate sun zenith angle
	sun_zenith = DegToRad(90.0) - altitude;

	// Compute the inputs for the calculation of the light distribution
	// over the sky
	//
	if (input == 0)
	{
		// Calculate irradiance
		diff_irrad = CalcDiffuseIrradiance();
		dir_irrad = CalcDirectIrradiance();
		
		// Calculate illuminance
		index = GetCategoryIndex();
		diff_illum = diff_irrad * CalcDiffuseIllumRatio(index);
		dir_illum = dir_irrad * CalcDirectIllumRatio(index);
	}
	else if (input == 1)
	{
		sky_brightness = CalcSkyBrightness();
		sky_clearness =  CalcSkyClearness();

		// Calculate illuminance
		index = GetCategoryIndex();
		diff_illum = diff_irrad * CalcDiffuseIllumRatio(index);
		dir_illum = dir_irrad * CalcDirectIllumRatio(index);
	}
	else if (input == 2)
	{
		// Calculate sky brightness and clearness from illuminance values
		index = CalcSkyParamFromIllum();
	}

	// Calculate Perez sky model parameters
	CalcPerezParam(sun_zenith, sky_clearness, sky_brightness, index);

	// Calculate sky patch luminance values
	printf("DEBUG: Patch luminance values\n");
	CalcSkyPatchLumin(sky_patch);

	// Calculate relative horizontal illuminance
	norm_diff_illum = CalcRelHorzIllum(sky_patch);

	// Normalization coefficient
	norm_diff_illum = diff_illum / norm_diff_illum;

	// Calculate relative zenith luminance
	zlumin = CalcRelLuminance(sun_zenith, 0.0);

	// Calculate absolute zenith illuminance
	zlumin *= norm_diff_illum;

	printf("gendaylit - zenith luminance: %.0lf\n", zlumin);
}

// Print usage error message and quit
void UsageError( const char *pmsg )
{
	if (pmsg != NULL)
		fprintf(stderr, "%s: Use error - %s\n", progname, pmsg);
	fprintf(stderr, "Usage: %s altitude azimuth [-P|-W|-L] "
			"direct_value diffuse_value [options]\n", progname);
	fprintf(stderr, "	-P epsilon delta  (these are the Perez "
			"parameters)\n");
	fprintf(stderr, "	-W direct-normal-irradiance diffuse-horizontal-"
			"irradiance (W/m^2)\n");
	fprintf(stderr, "	-L direct-normal-illuminance diffuse-"
			"horizontal-illuminance (lux)\n");

	exit(1);
}

// Determine category index
int GetCategoryIndex()
{
	int index;	// Loop index

	for (index = 0; index < 8; index++)
		if ((sky_clearness >= SkyClearCat[index].lower) &&
				(sky_clearness < SkyClearCat[index].upper))
			break;

	return index;
}

// Calculate diffuse illuminance to diffuse irradiance ratio
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 7.
//
double CalcDiffuseIllumRatio( int index )
{
	ModelCoeff const *pnle;	// Category coefficient pointer
	
	// Get category coefficient pointer
	pnle = &(DiffuseLumEff[index]);

	return pnle->a + pnle->b * apwc + pnle->c * cos(sun_zenith) +
			pnle->d * log(sky_brightness);
}

// Calculate direct illuminance to direct irradiance ratio
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 8.
//
double CalcDirectIllumRatio( int index )
{
	ModelCoeff const *pnle;	// Category coefficient pointer

	// Get category coefficient pointer
	pnle = &(DirectLumEff[index]);

	// Calculate direct illuminance from direct irradiance
	//
	// NOTE: The C equivalent of the C++ function __max() is max()
	//
	return __max((pnle->a + pnle->b * apwc + pnle->c * exp(5.73 *
			sun_zenith - 5.0) + pnle->d * sky_brightness),
			0.0);
}

// Calculate sky brightness
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 2.
//
double CalcSkyBrightness()
{
	return diff_irrad * CalcAirMass() / (DC_SolarConstantE *
			CalcEccentricity());
}

// Calculate sky clearness
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 1.
//
double CalcSkyClearness()
{
	double sz_cubed;	// Sun zenith angle cubed

	// Calculate sun zenith angle cubed
	sz_cubed = pow(sun_zenith, 3.0);

	return ((diff_irrad + dir_irrad) / diff_irrad + 1.041 *
			sz_cubed) / (1.0 + 1.041 * sz_cubed);
}

// Calculate diffuse horizontal irradiance from Perez sky brightness
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 2
//				(inverse).
//
double CalcDiffuseIrradiance()
{
	return sky_brightness * DC_SolarConstantE * CalcEccentricity() /
			CalcAirMass();
}

// Calculate direct normal irradiance from Perez sky clearness
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 1
//				(inverse).
//
double CalcDirectIrradiance()
{
	return CalcDiffuseIrradiance() * ((sky_clearness - 1.0) * (1 + 1.041
			* pow(sun_zenith, 3.0)));
}

// Calculate sky brightness and clearness from illuminance values
int CalcSkyParamFromIllum()
{
	double test1 = 0.1;
	double test2 = 0.1;
	int	counter = 0;
	int index = 0;			// Category index

	// Convert illuminance to irradiance
	diff_irrad = diff_illum * DC_SolarConstantE /
			(DC_SolarConstantL * 1000.0);
	dir_irrad = dir_illum * DC_SolarConstantE /
			(DC_SolarConstantL * 1000.0);

	// Calculate sky brightness and clearness
	sky_brightness = CalcSkyBrightness();
	sky_clearness =  CalcSkyClearness(); 

	// Limit sky clearness
	if (sky_clearness > 12.0)
		sky_clearness = 12.0;

	// Limit sky brightness
	if (sky_brightness < 0.05)
			sky_brightness = 0.01; 

	while (((fabs(diff_irrad - test1) > 10.0) ||
			(fabs(dir_irrad - test2) > 10.0)) && !(counter == 5))
	{
		test1 = diff_irrad;
		test2 = dir_irrad;	
		counter++;
	
		// Convert illuminance to irradiance
		index = GetCategoryIndex();
		diff_irrad = diff_illum / CalcDiffuseIllumRatio(index);
		dir_irrad = dir_illum / CalcDirectIllumRatio(index);
	
		// Calculate sky brightness and clearness
		sky_brightness = CalcSkyBrightness();
		sky_clearness =  CalcSkyClearness();

		// Limit sky clearness
		if (sky_clearness > 12.0)
			sky_clearness = 12.0;
	
		// Limit sky brightness
		if (sky_brightness < 0.05)
			sky_brightness = 0.01; 
	}

	return GetCategoryIndex();
}		

// Calculate relative luminance
//
// Reference:	Perez, R., R. Seals, and J. Michalsky. 1993.
//				“All-Weather Model for Sky Luminance Distribution -
//				Preliminary Configuration and Validation,” Solar Energy
//				50(3):235-245, Eqn. 1.
//
double CalcRelLuminance( double gamma, double zeta )
{
	return (1.0 + perez_param[0] * exp(perez_param[1] / cos(zeta))) *
		    (1.0 + perez_param[2] * exp(perez_param[3] * gamma) +
			perez_param[4] * cos(gamma) * cos(gamma));
}

// Calculate Perez sky model parameters
//
// Reference:	Perez, R., R. Seals, and J. Michalsky. 1993.
//				“All-Weather Model for Sky Luminance Distribution -
//				Preliminary Configuration and Validation,” Solar Energy
//				50(3):235-245, Eqns. 6 - 8.
//
void CalcPerezParam( double sz, double epsilon, double delta,
		int index )
{
	double x[5][4];		// Coefficents a, b, c, d, e
	int i, j;			// Loop indices

	// Limit sky brightness
	if (epsilon > 1.065 && epsilon < 2.8)
	{
		if (delta < 0.2)
			delta = 0.2;
	}

	// Get Perez coefficients
	for (i = 0; i < 5; i++)
		for (j = 0; j < 4; j++)
			x[i][j] = PerezCoeff[index][4 * i + j];

	if (index != 0)
	{
		// Calculate parameter a, b, c, d and e (Eqn. 6)
		for (i = 0; i < 5; i++)
			perez_param[i] = x[i][0] + x[i][1] * sz + delta * (x[i][2] +
					x[i][3] * sz);
	}
	else
	{
		// Parameters a, b and e (Eqn. 6)
		perez_param[0] = x[0][0] + x[0][1] * sz + delta * (x[0][2] +
				x[0][3] * sz);
		perez_param[1] = x[1][0] + x[1][1] * sz + delta * (x[1][2] +
				x[1][3] * sz);
		perez_param[4] = x[4][0] + x[4][1] * sz + delta * (x[4][2] +
				x[4][3] * sz);

		// Parameter c (Eqn. 7)
		perez_param[2] = exp(pow(delta * (x[2][0] + x[2][1] * sz),
				x[2][2])) - x[2][3];

		// Parameter d (Eqn. 8)
		perez_param[3] = -exp(delta * (x[3][0] + x[3][1] * sz)) + 
				x[3][2] + delta * x[3][3];
	}

	printf("Perez parameters:\n\ta = %f\n\tb = %f\n\tc = %f\n\td = %f\n"
			"\te = %f\n", perez_param[0], perez_param[1],
			perez_param[2], perez_param[3], perez_param[4]);
}

// Calculate relative horizontal illuminance
//
// Reference:	Perez, R., R. Seals, and J. Michalsky. 1993.
//				“All-Weather Model for Sky Luminance Distribution -
//				Preliminary Configuration and Validation,” Solar Energy
//				50(3):235-245, Eqn. 3.
//
double CalcRelHorzIllum( SkyPatch *sp )
{
	int i;
	double rh_illum = 0.0;	// Relative horizontal illuminance

	for (i = 0; i < 145; i++)
		rh_illum += sp[i].luminance * sin(sp[i].altitude);

	return rh_illum * 2.0 * GN_Pi / 145.0;
}

// Calculate earth orbit eccentricity correction factor
//
// Reference:	Sen, Z. 2008. Solar Energy Fundamental and Modeling 
//				Techniques. Springer, p. 72.
//
double CalcEccentricity()
{
	double day_angle;	// Day angle (radians)
	double E0;			// Eccentricity

	// Calculate day angle
	day_angle  = 2.0 * GN_Pi * (julian_date - 1.0) / 365.0;

	// Calculate eccentricity
	E0 = 1.00011 + 0.034221 * cos(day_angle) + 0.00128 * sin(day_angle)
			+ 0.000719 * cos(2.0 * day_angle) + 0.000077 * sin(2.0 *
			day_angle);

	return E0;
}

// Calculate atmospheric precipitable water content
//
// Reference:	Perez, R., P. Ineichen, R. Seals, J. Michalsky, and R.
//				Stewart. 1990. “Modeling Daylight Availability and
//				Irradiance Components from Direct and Global
//				Irradiance,” Solar Energy 44(5):271-289, Eqn. 3.
//
// Note:	The default surface dew point temperature is 11 deg. C
//			(52 deg. F). Typical values are:
//
//			Celsius 	Fahrenheit	 	Human Perception
//			> 24 		> 75 			Extremely uncomfortable
//			21 - 24 	70 - 74 		Very humid
//			18 - 21		65 - 69		 	Somewhat uncomfortable
//			16 - 18 	60 - 64 		OK for most people
//			13 - 16 	55 - 59		 	Comfortable
//			10 - 12 	50 - 54		 	Very comfortable
//			< 10 		< 49		 	A bit dry for some
//
double CalcPrecipWater( double dpt )
{ return exp(0.07 * dpt - 0.075); }

// Calculate relative air mass
//
// Reference:	Kasten, F. 1966. "A New Table and Approximation Formula
//				for the Relative Optical Air Mass," Arch. Meteorol.
//				Geophys. Bioklimataol. Ser. B14, pp. 206-233.
//
// Note:		More sophisticated relative air mass models are
//				available, but they differ significantly only for
//				sun zenith angles greater than 80 degrees.
//
double CalcAirMass()
{
	return (1.0 / (cos(sun_zenith) + 0.15 * pow((93.885 -
			sun_zenith * 180.0 / GN_Pi), -1.253)));
}

// Calculate Perez All-Weather sky patch luminances
//
// NOTE: The sky patches centers are determined in accordance with the
//       BRE-IDMP sky luminance measurement procedures. (See for example
//       Mardaljevic, J. 2001. "The BRE-IDMP Dataset: A New Benchmark
//       for the Validation of Illuminance Prediction Techniques,"
//       Lighting Research & Technology 33(2):117-136.)
//
void CalcSkyPatchLumin( SkyPatch *patch_array )
{
	int i, j;				// Loop indices
	int k = 0;				// Patch index
	int num_patch;			// Number of patches in horizontal band
	double aas;				// Sun-sky point azimuthal angle
	double pc_altitude;		// Patch center altitude (radians)
	double pc_azimuth;		// Patch center azimuth (radians)
	double sspa;			// Sun-sky point angle
	double zsa;				// Zenithal sun angle
	SkyPatch *ppatch;		// Sky patch pointer

	// Initialize patch array element pointer
	ppatch = patch_array;

	pc_altitude = GN_Pi * 0.5;	// Initialize altitude to zenith

	for (i = 0; i < 8; i++)
	{
		// Determine number of patches in horizontal band
		switch (i) 
		{
			case 0:
				num_patch = 1;
				break;
			case 1:
				num_patch = 6;
				break;
			case 2:
				num_patch = 12;
				break;
			case 3:
				num_patch = 18;
				break;
			case 4:
				num_patch = 24;
				break;
			case 5:
				num_patch = 24;
				break;
			case 6:
				num_patch = 30;
				break;
			case 7:
				num_patch = 30;
				break;
			default:
				break;
		}

		pc_azimuth = 0.0;	// Initialize azimuth

		for (j = 0; j < num_patch; j++)
		{
			// Save patch center coordinates
			ppatch->altitude = pc_altitude;
			ppatch->azimuth = pc_azimuth;

			// Calculate sun-sky point azimuthal angle
			aas = fabs(pc_azimuth - azimuth);

			// Calculate zenithal sun angle
			zsa = GN_Pi * 0.5 - pc_altitude;

			// Calculate sun-sky point angle (Equation 8-20)
			sspa = acos(cos(sun_zenith) * cos(zsa) + sin(sun_zenith) *
					sin(zsa) * cos(aas));

			// Calculate patch luminance
			ppatch->luminance = CalcRelLuminance(sspa, zsa);

			// Output patch luminance
			printf("  Patch %d theta = %f phi = %f luminance = %f\n", k,
					RadToDeg(zsa), RadToDeg(pc_azimuth),
					ppatch->luminance);

			ppatch++;	// Increment patch array element pointer
			k++;

			// Update azimuth
			pc_azimuth += 2.0 * GN_Pi / num_patch;
		}

		// Update altitude (12-degree increment)
		pc_altitude -= GN_Pi / 15;	
	}
}