Hal Yeager's Solar-Electroc Home
Hal Yeager's Solar-Electric Home

Code Listing (C-language)
Copy and compile with your own C-compiler.
Sorry, I-m an old C-hack who has been out of the business for 14 years now, and so I don't
know how to put a fancy user interface on this thing.

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

};

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

};

}