This code was corrected on June 26, 2001, to reflect a more accuract TOU multiplier.
#include <math.h>
#include <stdio.h>
#include <alloc.h>
#include <process.h>
#include <stdlib.h>
/******************************************************************************/
/**********
PROGRAM SET-UP AND ORGANIZATION
**********/
/******************************************************************************/
/******************************************************************************/
/**********
Global Data Structures
**********/
/******************************************************************************/
struct XYZvector {double XX, YY, ZZ; };
/******************************************************************************/
/************
MATH SUPPORT ROUTINES
***********/
/******************************************************************************/
double sindeg(degrees)
double degrees;
{ return(sin(2*3.14159*(degrees/360)));
};
double arcsindeg(sin_value)
double sin_value;
{ return( (360/(2*3.14159))*asin(sin_value) );
};
double cosdeg(degrees)
double degrees;
{ return(cos(2*3.14159*(degrees/360)));
};
double tandeg(degrees)
double degrees;
{ return(tan(2*3.14159*(degrees/360)));
};
/******************************************************************************/
/************
END
***********/
/******************************************************************************/
/******************************************************************************/
/*** The following routine, Phi_from_day(), computes the
angle of the ***/
/*** Earth's axis given the day of the year. A value
of zero for Phi ***/
/*** means that the Earth's axis is perpendicular to the
orbital plane of ***/
/*** the earth about the sun. This condition corresponds
to two days a ***/
/*** year, roughly March 21 and September 21 (the equinoxes).
A positive ***/
/*** value of Phi corresponds to the North pole being tilted
toward the ***/
/*** sun (summar for the northern Hemisphere), and a negative
value of Phi ***/
/*** corresponds to the North pole being tilted away from
the sun (winter ***/
/*** for the northern Hemisphere). Phi ranges between
-23.5 and +23.5 degs.**/
/******************************************************************************/
double Phi_from_day(input_day)
/*-------------------------------------------------------------*/
int input_day; /** Day ranges between 0 and 364, with Day=0 corresponding
*/
/** to June 21st
*/
/*-------------------------------------------------------------*/
{ double day_in_radians;
day_in_radians = 2.0*3.14159*((double)input_day)/((double)365.0);
return( 23.5*cos(day_in_radians) );
};
/******************************************************************************/
/*** The following routines, Sunrise() and Sunset(), compute
the time ***/
/*** times of sumrise and sunset for a given day of the year
and the ***/
/*** latitude of the site. The time determine the start
and stop times ***/
/*** for the integration loops which integrate the amount
of sun exposure. ***/
/*** Time=0 is defined as "midnight", and corresponds to
the time when the ***/
/*** site is furtherest from the sun. Time=12 corresponds
to noon, where ***/
/*** the site is generally closest to the sun. Time
ranges between 0 and ***/
/*** 24, and generally corresponds to our time during the
winter. During ***/
/*** the summar when day-light savings time is in effect,
the value of ***/
/*** time is one hour less than our time becuase we have
"sprung ahead" ***/
/*** hour.
***/
/******************************************************************************/
double Sunrise(latitude, input_day)
double latitude; /* Northern latitude in degrees */
int input_day; /*range of day is 0-364, with June 21st
being day 0 */
{ double singamma, gamma, sin_Phi;
sin_Phi=sindeg(Phi_from_day(input_day));
/** First determine if the site is so far north that it is
always in sun **/
/** during summar (in which case sunrise is set to 0, corresponding
to **/
/** midnight), or in darkness during winter (in which case
sunrise is set **/
/** set to 12, corresponding to noon.
**/
if (sindeg(latitude)*sin_Phi > cosdeg(latitude)) return
(0.0);
if (sindeg(latitude)*sin_Phi < -cosdeg(latitude)) return
(12.0);
/** If neither of the above conditions exist, then compute
sunrise. This **/
/** done differentially with respect to the equinoxes (March
21 and **/
/** September 21) where sunrise=6 (6am) and some complicated
geometry. **/
singamma = (tandeg(latitude))*sin_Phi;
gamma = arcsindeg(singamma);
return(6.0 - gamma*6.0/90.0);
};
double Sunset(latitude, input_day)
double latitude; /* Northern latitude in degrees */
int input_day; /*range of day is 0-364, with June 21st being day
0 */
{ double singamma, gamma, sin_Phi;
sin_Phi=sindeg(Phi_from_day(input_day));
/** First determine if the site is so far north that it is
always in sun **/
/** during summar (in which case sunset is set to 24, corresponding
to **/
/** midnight), or in darkness during winter (in which case
sunset is set **/
/** set to 12, corresponding to noon.
**/
if (sindeg(latitude)*sin_Phi > cosdeg(latitude)) return
(24.0);
if (sindeg(latitude)*sin_Phi < -cosdeg(latitude)) return
(12.0);
/** If neither of the above conditions exist, then compute
sunset. This **/
/** done differentially with respect to the equinoxes (March
21 and **/
/** September 21) where sunset=18 (6pm) and some complicated
geometry. **/
singamma = (tandeg(latitude))*sin_Phi;
gamma = arcsindeg(singamma);
return(18.0 + gamma*6.0/90.0);
};
/******************************************************************************/
/****** ROUTINE EarthNormal()-- computes the normal vector of the
Earth's *****/
/******
Surface at the given latitude, time of *****/
/******
day, and tilt of the Earth's Axis
*****/
/****** The routine recieves a pointer to an XYZ vector,
which the rou- *****/
/****** tine modifies to reflect the computed normal vector.
For this *****/
/****** computation, the Z-axis is defined as being perpendicular
to the *****/
/****** orbital plane of the Earth, the Y-axis always points
to the sun, *****/
/****** and the X-axis is perpendicular to the Y- and Z-axes.
*****/
/******
*****/
/****** THIS ROUTINE IS NOT USED IN THE COMPUTATION PROCESS,
BUT WAS *****/
/****** DEVLOPED AS A TESTING ROUTINE TO AID ME IN THE DEVELOPMENT
OF *****/
/****** THE NEXT ROUTINE PanelNormal(). EarthNormal()
applies to the *****/
/****** SPECIAL CASE WHERE THE PANEL HAS SOUTH TILT=0 AND
WEST TILT=0. *****/
/******************************************************************************/
void EarthNormal(latitude, time, Phi, E_N)
/*-----------------------------------------------------------*/
double latitude, /* in degrees, 0 being the equator, 90 being north
pole */
time,
/* in hours, 0 being midnight, 12 being noon, 18 being 6 pm */
Phi;
/* in degrees, 23.5 degs is summar, -23.5 degs is winter
*/
XYZvector *E_N; /* pointer to an XYZ vector which updated
by this routine */
/*-----------------------------------------------------------*/
{ double time_angle, xtemp, ytemp, ztemp;
time_angle = -(360/24.0)*time;/*negative sign for counter-clockwise rotation*/
/* Compute Normal before rotation due to Earth's Axis (Phi)
*/
xtemp = cosdeg(latitude)*sindeg(time_angle);
ytemp = cosdeg(latitude)*cosdeg(time_angle);
ztemp = sindeg(latitude);
/* Now we account for the rotation due to Earth's Axis(Phi)
*/
/* This affects the values of the Y and Z coordinates only.
*/
E_N->XX = xtemp;
E_N->YY = ytemp*cosdeg(-Phi) + ztemp*sindeg(-Phi);
E_N->ZZ = ztemp*cosdeg(-Phi) - ytemp*sindeg(-Phi);
};
/******************************************************************************/
/****** ROUTINE PanelNormal()-- computes the integrated light exposure
*****/
/******
of a given day.
*****/
/******************************************************************************/
/******************************************************************************/
/****** ROUTINE PanelNormal()-- computes the normal vector of the
solar *****/
/******
Panel for the given latitude, time of *****/
/******
day, and tilt of the Earth's Axis
*****/
/****** The routine recieves a pointer to an XYZ vector,
which the rou- *****/
/****** tine modifies to reflect the computed normal vector.
For this *****/
/****** computation, the Z-axis is defined as being perpendicular
to the *****/
/****** orbital plane of the Earth, the Y-axis always points
to the sun, *****/
/****** and the X-axis is perpendicular to the Y- and Z-axes.
*****/
/******
*****/
/****** THIS ROUTINE WAS DEVELOPED FROM ROUTINE EarthNormal(),
see above.*****/
/******************************************************************************/
void PanelNormal(latitude, time, Phi, South_tilt, West_tilt, P_N)
double latitude, /* in degrees, 0 being the equator, 90 being
north pole */
time,
/* in hours, 0 being midnight, 12 being noon, 18 being 6 pm */
Phi,
/* in degs, 23.5 degrees is summar, -23.5 degrees is winter */
South_tilt, /* in degs, 0 being
flat, 90 being south, -90 being north */
West_tilt; /* in degs, 0 being
flat, 90 being west, -90 being east */
XYZvector *P_N; /* pointer to an XYZ vector which updated
by this routine */
{ double time_angle, xtemp, ytemp, ztemp;
time_angle = -(360/24.0)*time;/*negative sign for counter-clockwise rotation*/
/* Compute Normal without accouting for rotation due to Earth
Axis (Phi)*/
xtemp = cosdeg(latitude-South_tilt)*sindeg(time_angle+West_tilt);
ytemp = cosdeg(latitude-South_tilt)*cosdeg(time_angle+West_tilt);
ztemp = sindeg(latitude-South_tilt);
/*accouting for rotation due to Earth Axis (Phi)*/
P_N->XX = xtemp;
P_N->YY = ytemp*cosdeg(-Phi) + ztemp*sindeg(-Phi);
P_N->ZZ = ztemp*cosdeg(-Phi) - ytemp*sindeg(-Phi);
};
/******************************************************************************/
/****** ROUTINE Day_Exposure()-- compute the integrated light exposure
*****/
/******
of a given day.
*****/
/******************************************************************************/
double Day_exposure(day, latitude, South_tilt, West_tilt, TOU)
/**-------------------------------------------------------**/
int day;
/** range of day is 0-364, with June 21st being day 0
**/
double latitude, /** in degrees, 0 being the equator,
90 being north pole **/
South_tilt, /** in degs, 0
being flat, 90 being south, -90 being north**/
West_tilt, /** in degs,
0 being flat, 90 being west, -90 being east **/
*TOU;
/** pointer to address to fill with exposure with TOU metering*/
/**-------------------------------------------------------**/
{ double sunrise_time, sunset_time, time, exposure,
Phi,
TOU_exposure, TOU_multiplier;
XYZvector P_Normal;
exposure=0.0;
TOU_exposure=0.0;
TOU_multiplier=0.0;
sunrise_time=Sunrise(latitude, day);
sunset_time =Sunset(latitude, day);
Phi = Phi_from_day(day);
/* printf("sunrise=%4.2f sunset=%4.2f Phi=%4.2f\n", sunrise_time,
sunset_time, Phi); */
/* Compute time of use (TOU) Multiplier factor for TOU EXPOSURE
COMPUTATION */
/* Summer is defined as May 1st through October 31, which
corresponds
to days 313-364 and 0-128, Winter corresponds to days
129-313 */
if ((day>130)&&(day<312)) TOU_multiplier=1.22;
/* It's a winter day*/
else
TOU_multiplier=2.93; /* It's a summar day*/
/** TOU_Multiplier corrected from 3.7 to 2.93 on June 26. 2001 **/
/*******************/
/* Main Loop */
/*******************/
for (time=sunrise_time; time<sunset_time; time=time+0.1)
{ PanelNormal(latitude, time, Phi, South_tilt, West_tilt,
&P_Normal);
/* The next line computs the dot product of the sun's
ray vector with */
/* the normal vector of the panel. Sun's ray
vector is XYZ=(0, -1, 0) */
/* for the coordinate system defined in Panel_Normal().
*/
P_Normal.YY = -P_Normal.YY;
/* Assume that the system will not generate power unless
the angle */
/* between the sun's rays and the panel is between
-78 deg and 78 deg. */
/* Cos(+-78 deg) = 0.2. At +-78 deg, the array
puts out 20% of its */
/* maximum current, so this does not seem to be an
unreasonable asmpt. */
if (P_Normal.YY > 0.2)
{ exposure = exposure + 0.1*P_Normal.YY; /* Integrate
by 0.1 hour increments*/
/* TOU EXPOSURE COMPUTATION */
if ((day>130)&&(day<312)) /*
Winter Time - No day light savings time*/
{ if ((time>12.0)&&(time<18.0))
TOU_exposure = TOU_exposure + 0.1*P_Normal.YY*TOU_multiplier;
else
TOU_exposure = TOU_exposure + 0.1*P_Normal.YY;
}
else /* Summer Time - compensate
one hour for day light savings time*/
{ if ((time>11.0)&&(time<17.0))
TOU_exposure = TOU_exposure + 0.1*P_Normal.YY*TOU_multiplier;
else
TOU_exposure = TOU_exposure + 0.1*P_Normal.YY;
};
};
};
/* printf("\n Exposure=%4.2f TOU_Exposure=%4.2f\n", exposure, TOU_exposure);*/
*TOU=TOU_exposure;
return(exposure);
};
double Winter_Cloud_factor(day)
int day;
{
if (day<121) return(1.0); /* Summar days are assume
to be 0% cloudy. */
if (day<151) return(0.75); /* Fall days are assume to
be 25% cloudy. */
if (day<242) return(0.5); /* Winter days are assume
to be 50% cloudy. */
if (day<273) return(0.75); /* Spring days are assume
to be 25% cloudy. */
if (day<365) return(1.0); /* Summar days are assume
to be 0% cloudy. */
return(1.0);
};
/******************************************************************************/
/**** Year_Exposure compute the average light exposure over a year
****/
/******************************************************************************/
double Year_exposure(latitude, South_tilt, West_tilt, TOU)
double latitude, /* in degrees, 0 being the equator, 90 being north
pole */
South_tilt, /* in degs, 0 being flat,
90 being south, -90 being north */
West_tilt, /* in degs, 0 being flat,
90 being west, -90 being east */
*TOU;
/* t */
{
int day; /*range of day is 0-364,
with June 21st being day 0 */
double exposure, temp1, TOU_exposure, temp2;
exposure = 0.0;
TOU_exposure = 0.0;
for (day=0; day<365; day++)
{ temp1 = Day_exposure(day, latitude, South_tilt, West_tilt,&temp2);
exposure = exposure
+ temp1*Winter_Cloud_factor(day);
TOU_exposure = TOU_exposure + temp2*Winter_Cloud_factor(day);
/* printf("day=%d exposure=%4.2f\n", day, temp);
*/
};
*TOU = TOU_exposure/365.0;
return(exposure/365.0);
};
/******************************************************************************/
/**** Main program, collect input, call Year_Exposure, and output
results. ****/
/******************************************************************************/
void main ()
{ double mysin, latitude, Phi, time, length, Xdot, Ydot,
Zdot,
South_tilt, West_tilt, dotprod, exposure, TOU;
int ii, day;
XYZvector E_Normal, P_Normal;
time = 0.0;
latitude = 38.0;
Phi = 0.0;
E_Normal.XX = 0.0;
E_Normal.YY = 0.0;
E_Normal.ZZ = 0.0;
P_Normal.XX = 0.0;
P_Normal.YY = 0.0;
P_Normal.ZZ = 0.0;
South_tilt=0.0;
West_tilt=0.0;
printf("\n");
printf("South |
West Tilt\n");
printf("Tilt | ");
for (West_tilt= -20.0; West_tilt < 21.0; West_tilt=West_tilt+5.0)
{ printf(" %4.2f",West_tilt);
};
printf("\n");
printf("------------------------------------------------------------------------------");
for (South_tilt=-20.0; South_tilt < 71; South_tilt=South_tilt+5.0)
{
printf("%4.2f | ",South_tilt);
for (West_tilt= 25.0; West_tilt < 26.0; West_tilt=West_tilt+5.0)
{
exposure = Year_exposure(latitude, South_tilt,
West_tilt, &TOU);
printf(" %4.2f(%4.2f)",exposure, TOU);
};
printf("\n");
};
for (time=0.0; time < 0.0; time=time+1.0)
{ EarthNormal(latitude, time, Phi, &E_Normal);
PanelNormal(latitude, time, Phi, 20.0, 0.0, &P_Normal);
length = sqrt(P_Normal.XX*P_Normal.XX +
P_Normal.YY*P_Normal.YY +
P_Normal.ZZ*P_Normal.ZZ);
dotprod = P_Normal.XX*E_Normal.XX
+
P_Normal.YY*E_Normal.YY +
P_Normal.ZZ*E_Normal.ZZ;
printf("latitude=%4.2f time=%4.2f Phi=%4.2f",latitude,
time,Phi);
printf(" X=%4.2f, Y=%4.2f, Z=%4.2f l=%4.2f
dp=%4.2f\n",P_Normal.XX,
P_Normal.YY,
P_Normal.ZZ,
length,
dotprod);
};
}