/*********************************************************************
 *
 *  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.
 *				1/11/13 - Tweaks and optimizations (G.Ward)
 *
 *  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 */
#define	_USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

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

char *progname;								/* Program name */
char errmsg[128];							/* Error message buffer */
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 */

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

/* Degrees into radians */
#define DegToRad(deg)	((deg)*(M_PI/180.))

/* Radiuans into degrees */
#define RadToDeg(rad)	((rad)*(180./M_PI))


/* 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. */

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

/* 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 */
typedef struct
{
	double a;
	double b;
	double c;
	double d;
} ModelCoeff;

/* 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);
}

/* Return maximum of two doubles */
double dmax( double a, double b )
{ return (a > b) ? a : b; }

/* 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 */
	
	return dmax((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 * M_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  = (julian_date - 1.0) * (2.0 * M_PI / 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 -
			RadToDeg(sun_zenith), -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 = M_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 = M_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 * M_PI / num_patch;
		}

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