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

  };

}