code to test gendaylit against the python version

gendaylit_test.c

#define  _USE_MATH_DEFINES
#include  <stdio.h>
#include  <string.h>
#include  <math.h>
#include <stdlib.h>

float coeff_perez[] = {
	1.3525,-0.2576,-0.2690,-1.4366,-0.7670,0.0007,1.2734,-0.1233,2.8000,0.6004,1.2375,1.000,1.8734,0.6297,
	0.9738,0.2809,0.0356,-0.1246,-0.5718,0.9938,-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,-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,-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,-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,-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,-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,-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};


float defangle_theta[] = {
	84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84,
	84, 84, 84, 84, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72,
	72, 72, 72, 72, 72, 72, 72, 72, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60,
	60, 60, 60, 60, 60, 60, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48,
	48, 48, 48, 48, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 24, 24, 24, 24,
	24, 24, 24, 24, 24, 24, 24, 24, 12, 12, 12, 12, 12, 12, 0};

float defangle_phi[] = {
	0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264,
	276, 288, 300, 312, 324, 336, 348, 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180,
	192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 0, 15, 30, 45, 60, 75, 90, 105,
	120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 15, 30, 45, 60, 75,
	90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 0, 20, 40, 60,
	80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 320, 340, 0, 30, 60, 90, 120, 150, 180, 210,
240, 270, 300, 330, 0, 60, 120, 180, 240, 300, 0};


/* Perez sky parametrization: epsilon and delta calculations from the direct and diffuse irradiances */
double sky_brightness();
double sky_clearness();


double	sunzenith;
double  calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[]);
double	radians(double degres);
void	theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z);

double integ_lv(float *lv,float *theta);

static int get_numlin(float epsilon);

void computesky();
void check_parametrization();


const double	AU = 149597890E3;
const double 	solar_constant_e = 1367;    /* solar constant W/m^2 */
const double  	solar_constant_l = 127500;   /* solar constant lux */

const double	half_sun_angle = 0.2665;
const double	half_direct_angle = 2.85;

const double	skyclearinf = 1.0;	    /* limitations for the variation of the Perez parameters */
const double	skyclearsup = 12.01;
const double	skybriginf = 0.01;
const double	skybrigsup = 0.6;
const double	WHTEFFICACY = 179.;

double skyclearness = 0;
double skybrightness = 0;
double	sunzenith, daynumber, atm_preci_water=2;
double	diffuseilluminance, directilluminance, diffuseirradiance, directirradiance, globalirradiance;

double 	glob_h_diffuse_effi_PEREZ();
double direct_n_effi_PEREZ();

/* astronomy and geometry*/
double 	get_eccentricity();
double air_mass();

int main() {
   computesky();
}


void computesky(){

	float 	*lv_mod;  /* 145 luminance values */
	float 	*theta_o, *phi_o;
    double dzeta, gamma;
	
	if ( (lv_mod = malloc(145*sizeof(float))) == NULL)
	{
		fprintf(stderr,"Out of memory in function main");
		exit(1);
	}

	/* read the angles */
	theta_o = defangle_theta;
	phi_o = defangle_phi;
	int j;
    double	solarradiance;
    sunzenith = 88.2;
    double diffnormalization;

    daynumber = 1;
    directirradiance = 22;
    diffuseirradiance = 10;
    /* parameters for the perez model */
	
	skybrightness = sky_brightness();
	skyclearness =  sky_clearness();
	check_parametrization();
	
	diffuseilluminance = diffuseirradiance*glob_h_diffuse_effi_PEREZ();/*diffuse horizontal illuminance*/
	directilluminance = directirradiance*direct_n_effi_PEREZ();
    // check_illuminances();
		
    
	
    /*calculation of the modelled luminance */
	for (j=0;j<145;j++)
	{
		theta_phi_to_dzeta_gamma(radians(*(theta_o+j)),radians(*(phi_o+j)),&dzeta,&gamma,radians(sunzenith));
        
		*(lv_mod+j) = calc_rel_lum_perez(dzeta,gamma,radians(sunzenith),skyclearness,skybrightness,coeff_perez);
		
		// fprintf(stderr,"theta, phi, lv_mod %f\t %f\t %f\n", *(theta_o+j),*(phi_o+j),*(lv_mod+j));
    }
    
    diffnormalization = integ_lv(lv_mod, theta_o);
    
    diffnormalization = diffuseilluminance/diffnormalization/WHTEFFICACY;
    
    solarradiance = directilluminance/(2*M_PI*(1-cos(half_sun_angle*M_PI/180)))/WHTEFFICACY;
    fprintf(stderr,"solarradiance,: %f", solarradiance);
    
}



/*check the range of epsilon and delta indexes of the perez parametrization*/
void check_parametrization()
{

if (skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup)
		{

/*  limit sky clearness or sky brightness, 2009 11 13 by J. Wienold */
		
		if (skyclearness<skyclearinf){
			/* if (suppress_warnings==0)
			    fprintf(stderr,"Range warning: sky clearness too low (%lf)\n", skyclearness); */
			skyclearness=skyclearinf;
		}
		if (skyclearness>skyclearsup){
			/* if (suppress_warnings==0)
			    fprintf(stderr,"Range warning: sky clearness too high (%lf)\n", skyclearness); */
			skyclearness=skyclearsup-0.001;
		}
		if (skybrightness<skybriginf){
			/* if (suppress_warnings==0)
			    fprintf(stderr,"Range warning: sky brightness too low (%lf)\n", skybrightness); */
			skybrightness=skybriginf;
		}
		if (skybrightness>skybrigsup){
			/* if (suppress_warnings==0)
			    fprintf(stderr,"Range warning: sky brightness too high (%lf)\n", skybrightness); */
			skybrightness=skybrigsup;
		}

	return;	}
	else return;
}

/* sky luminance perez model */
double calc_rel_lum_perez(double dzeta,double gamma,double Z,double epsilon,double Delta,float coeff_perez[])
{
			
	float x[5][4];
	int i,j,num_lin;
	double c_perez[5];

	if ( (epsilon <  skyclearinf) || (epsilon >= skyclearsup) )
	{
		fprintf(stderr,"Error: epsilon out of range in function calc_rel_lum_perez!\n");
		exit(1);
	}

	/* correction de modele de Perez solar energy ...*/
	if ( (epsilon > 1.065) && (epsilon < 2.8) )
	{
		if ( Delta < 0.2 ) Delta = 0.2;
	}
	
	
	num_lin = get_numlin(epsilon);
	
	for (i=0;i<5;i++)
		for (j=0;j<4;j++)
		{
			x[i][j] = *(coeff_perez + 20*num_lin + 4*i +j);
			/* fprintf(stderr,"x %d %d vaut %f\n",i,j,x[i][j]); */
		}


	if (num_lin)
	{
		for (i=0;i<5;i++)
			c_perez[i] = x[i][0] + x[i][1]*Z + Delta * (x[i][2] + x[i][3]*Z);
	}
	else
	{
		c_perez[0] = x[0][0] + x[0][1]*Z + Delta * (x[0][2] + x[0][3]*Z);
		c_perez[1] = x[1][0] + x[1][1]*Z + Delta * (x[1][2] + x[1][3]*Z);
		c_perez[4] = x[4][0] + x[4][1]*Z + Delta * (x[4][2] + x[4][3]*Z);
		c_perez[2] = exp( pow(Delta*(x[2][0]+x[2][1]*Z),x[2][2])) - x[2][3];
		c_perez[3] = -exp( Delta*(x[3][0]+x[3][1]*Z) )+x[3][2]+Delta*x[3][3];
	}


	return (1 + c_perez[0]*exp(c_perez[1]/cos(dzeta)) ) *
	    (1 + c_perez[2]*exp(c_perez[3]*gamma) +
	    c_perez[4]*cos(gamma)*cos(gamma) );
}

/* degrees into radians */
double radians(double degres)
{
	return degres*M_PI/180.0;
}

/* calculation of the angles dzeta and gamma */
void theta_phi_to_dzeta_gamma(double theta,double phi,double *dzeta,double *gamma, double Z)
{
	*dzeta = theta; /* dzeta = phi */
	if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1 && (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi) < 1.1 ) )
		*gamma = 0;
	else if ( (cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi)) > 1.1 )
	{
		printf("error in calculation of gamma (angle between point and sun");
		exit(1);
	}
	else
		*gamma = acos(cos(Z)*cos(theta)+sin(Z)*sin(theta)*cos(phi));
}


double integ_lv(float *lv,float *theta)
{
	int i;
	double buffer=0.0;
	
	for (i=0;i<145;i++)
	{
		buffer += (*(lv+i))*cos(radians(*(theta+i)));
	}
			
	return buffer*2*M_PI/144;
}


static int get_numlin(float epsilon)
{
	if (epsilon < 1.065)
		return 0;
	else if (epsilon < 1.230) 
		return 1;
	else if (epsilon < 1.500)
		return 2;
	else if (epsilon < 1.950)
		return 3;
	else if (epsilon < 2.800)
		return 4;
	else if (epsilon < 4.500)
		return 5;
	else if (epsilon < 6.200)
		return 6;
	return 7;
}

/* Perez sky's brightness */
double sky_brightness()
{
double value;

value = diffuseirradiance * air_mass() / ( solar_constant_e*get_eccentricity());

return(value);
}

/* Perez sky's clearness */
double sky_clearness()
{
	double value;

	value = ( (diffuseirradiance + directirradiance)/(diffuseirradiance) + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180 ) / (1 + 1.041*sunzenith*M_PI/180*sunzenith*M_PI/180*sunzenith*M_PI/180) ;

	return(value);
}


/* global horizontal diffuse efficacy model, according to PEREZ */
double glob_h_diffuse_effi_PEREZ()
{
	double 	value;
	double 	category_bounds[10], a[10], b[10], c[10], d[10];
	int   	category_total_number, category_number, i;

	check_parametrization();
	
	
/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
 fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function glob_h_diffuse_PEREZ \n"); */
     
/* initialize category bounds (clearness index bounds) */

	category_total_number = 8;

//XXX:  category_bounds > 0.1
	category_bounds[1] = 1;
	category_bounds[2] = 1.065;
	category_bounds[3] = 1.230;
	category_bounds[4] = 1.500;
	category_bounds[5] = 1.950;
	category_bounds[6] = 2.800;
	category_bounds[7] = 4.500;
	category_bounds[8] = 6.200;
	category_bounds[9] = 12.01;


	/* initialize model coefficients */
	a[1] = 97.24;
	a[2] = 107.22;
	a[3] = 104.97;
	a[4] = 102.39;
	a[5] = 100.71;
	a[6] = 106.42;
	a[7] = 141.88;
	a[8] = 152.23;

	b[1] = -0.46;
	b[2] = 1.15;
	b[3] = 2.96;
	b[4] = 5.59;
	b[5] = 5.94;
	b[6] = 3.83;
	b[7] = 1.90;
	b[8] = 0.35;

	c[1] = 12.00;
	c[2] = 0.59;
	c[3] = -5.53;
	c[4] = -13.95;
	c[5] = -22.75;
	c[6] = -36.15;
	c[7] = -53.24;
	c[8] = -45.27;

	d[1] = -8.91;
	d[2] = -3.95;
	d[3] = -8.77;
	d[4] = -13.90;
	d[5] = -23.74;
	d[6] = -28.83;
	d[7] = -14.03;
	d[8] = -7.98;



	category_number = -1;
	for (i=1; i<=category_total_number; i++)
	{
		if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
			category_number = i;
	}

	if (category_number == -1) {
		fprintf(stderr, "Warning: sky clearness (= %.3f) too high, printing error sky\n", skyclearness);
		exit(0);
	}
		

	value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*cos(sunzenith*M_PI/180) + 
	    d[category_number]*log(skybrightness);

	return(value);

}


/* direct normal efficacy model, according to PEREZ */

double direct_n_effi_PEREZ()

{
double 	value;
double 	category_bounds[10], a[10], b[10], c[10], d[10];
int   	category_total_number, category_number, i;


/*if ((skyclearness<skyclearinf || skyclearness>skyclearsup || skybrightness<skybriginf || skybrightness>skybrigsup) && suppress_warnings==0)
   fprintf(stderr, "Warning: skyclearness or skybrightness out of range in function direct_n_effi_PEREZ \n");*/


/* initialize category bounds (clearness index bounds) */

category_total_number = 8;

category_bounds[1] = 1;
category_bounds[2] = 1.065;
category_bounds[3] = 1.230;
category_bounds[4] = 1.500;
category_bounds[5] = 1.950;
category_bounds[6] = 2.800;
category_bounds[7] = 4.500;
category_bounds[8] = 6.200;
category_bounds[9] = 12.1;


/* initialize model coefficients */
a[1] = 57.20;
a[2] = 98.99;
a[3] = 109.83;
a[4] = 110.34;
a[5] = 106.36;
a[6] = 107.19;
a[7] = 105.75;
a[8] = 101.18;

b[1] = -4.55;
b[2] = -3.46;
b[3] = -4.90;
b[4] = -5.84;
b[5] = -3.97;
b[6] = -1.25;
b[7] = 0.77;
b[8] = 1.58;

c[1] = -2.98;
c[2] = -1.21;
c[3] = -1.71;
c[4] = -1.99;
c[5] = -1.75;
c[6] = -1.51;
c[7] = -1.26;
c[8] = -1.10;

d[1] = 117.12;
d[2] = 12.38;
d[3] = -8.81;
d[4] = -4.56;
d[5] = -6.16;
d[6] = -26.73;
d[7] = -34.44;
d[8] = -8.29;



for (i=1; i<=category_total_number; i++)
{
if ( (skyclearness >= category_bounds[i]) && (skyclearness < category_bounds[i+1]) )
category_number = i;
}

value = a[category_number] + b[category_number]*atm_preci_water  + c[category_number]*exp(5.73*sunzenith*M_PI/180 - 5) +  d[category_number]*skybrightness;

if (value < 0) value = 0;

return(value);
}


/* enter sunzenith angle (degrees) return relative air mass (double) */
double 	air_mass()
{
double	m;
if (sunzenith>90)
	{
	fprintf(stderr, "Warning: air mass has reached the maximal value\n");
	sunzenith=90;
	}
m = 1/( cos(sunzenith*M_PI/180)+0.15*exp( log(93.885-sunzenith)*(-1.253) ) );
return(m);
}


double get_eccentricity()
{
	double day_angle;
	double E0;

	day_angle  = 2*M_PI*(daynumber -1)/365;
	E0         = 1.00011+0.034221*cos(day_angle)+0.00128*sin(day_angle)+
	    0.000719*cos(2*day_angle)+0.000077*sin(2*day_angle);

	return (E0);
}