ArduinoBeam - Suntracking skylight

Hey everyone, I've been working on this project for quite a while now but have finally got some time together and built the first prototype. As I get things together I'll be putting the designs as they stand on my blog http://arduinobeam.blogspot.com/. At this stage I just have a video on there which gives a bit of an idea of what I'm going for.

As for an overview of the hardware, there are two axes driven by maxon dc motors with optical encoders. I use a kind of bastardised PI system to set the speed and drive to the calculated positions. The sun position is calculated using time signals from a DS1307 real timeclock, and using the sun position code generously supplied by mowcius http://www.mowcius.co.uk/suntrackercode.html.

Let me know what you guys think!

Mowcius silently hates the person who is doing his website hosting for not sorting out the ftp yet - broken website - >:(

Need a host mowcius? Got a working server here ;)

Need a host mowcius? Got a working server here

Thanks for the offer but I think I should be ok.

I just need to nudge someone a bit harder :D

Let me know what you guys think!

Nice start to the project, I did a Tinker Tech blog entry on it at http://www.musheen.com/arduino-driven-arduinobeam-sun-tracker-1312.html

Let me know if you would like any updates or would like to share any further details.

Keep up the good work.

Boz www.musheen.com

I've run through a little accuracy test of the sun position code I'm using right now, as provided by Mowcius. To check the accuracy I compared results calculated by Mowcius' code and the NREL online solar position calculator solpos (uncertainty of +/- 0.01 degrees). You can see the results here:http://arduinobeam.blogspot.com/

Mowcius' code does seem to suffer from some pretty serious error especially in azimuth angle around midday (~10deg).

Mowcius' code does seem to suffer from some pretty serious error especially in azimuth angle around midday (~10deg).

Hmm, what's your earth position?

I checked mine against stellarium and it seemed pretty accurate at all times of the day.

Mowcius

I did the calculation every 10 minutes for the 11th of November using location data - Longitude = 172, Latitude = -39 (New Zealand).

The error does seem to follow a predictable pattern. Perhaps this is the result of additional effects not accounted for in the simplified calculation? Mowcius is your code based on the NREL code?

Hmm well my code was pretty much taken from someone else's code. No idea what that was based on.

I will take another look at what it gives me here at some point but like I said, it seemed to work fine when I used it before. You are completely over the other side of the world and rather lower than me...

If you get the time, do you want to try it with roughly my location and see what you get from the test algorithm you're using? 54.14, -1.5 ish ;)

Mowcius

A really great thing would be a sun position library built on the NREL code. It sounds like Emdee over at heliostat.sourceforge.net has one but is yet to release it. Perhaps some other programming guru will want to take up the challenge? ;) -- I think it may be beyond me at this point..

If anyone's curious here's the link to the NREL sources http://rredc.nrel.gov/solar/codesandalgorithms/solpos/ I think I might have a crack, C is C right?

I'll take a look. Got a bit of time this morning. No promises though. ;)

That code is very very complex :p

Right, a quote from the code I basically copied:

Blanco-Muriel from the Plataforma Solar de Almerýa (PSA) review the accuracy of all the algorithms. Further they develop a simplified algorithm that is accurate to within 0.5 minutes of arc for the year 1999-2015.

Well like I say, I have not had any major issues. Have you converted your lat and lon right to decimal degrees?

Mowcius

Hello, It would appear that Mowcius' implementation of the PSA 'sunpos' code or how it is being called that is causing the 'accuracy' problem.

The following code based on the same PSA algorithm generates data that is within 0.05 of a degree of the NREL algorithm.

Note this code uses (abuses?) the Arduino Time library for passing dates and times - this makes it easy to use a RTC and ask for the suns position for now()

Hope this helps - Andrew

#include <Time.h>                       // http://www.arduino.cc/playground/Code/Time

// DEFAULT TEST DATA VALUES
double Latitude = -37.7926;
double Longitude = 144.932;
int testDay = 18;
int testMonth = 11;
int testYear = 2010;

double ZenithAngle;
double Azimuth;

void setup() {
  Serial.begin(115200);
  Serial.print("Time,");
  Serial.print("Date,");
  Serial.print("Elevation,");
  Serial.print("Azimuth");
  Serial.println();
}

void loop() {

  for(int hr=0;hr<24;hr++)
    for(int min=0;min<60;min=min+10)
    {
      setTime(hr,min,0,testDay,testMonth,testYear);  // this is kinda corrupting the proper use of set time and now(), but hey its a test...
      sunpos(now(), Latitude, Longitude, &ZenithAngle, &Azimuth);

      printDateTime(now());
      Serial.print(",");
      Serial.print((90.0-ZenithAngle));
      Serial.print(",");
      Serial.print(Azimuth);
      Serial.println();
    }
    
  while(1); // wait here forever
}
//////////////////////////////////////////////
// SUN POSITION CALCULATION FUNCTION

/*
This file is based on an algorithm and code at http://www.pveducation.org/pvcdrom/properties-of-sunlight/sun-position-high-accuracy

PSA algorithm for High Accuracy Tracking of the Sun

Blanco-Muriel et al. from the Plataforma Solar de Almerýa (PSA) review the accuracy of all the algorithms. 
Further they develop a simplified algorithm that is accurate to within 0.5 minutes of arc for the year 1999-2015. 
The PSA algorithm has been specially optimised in C++ code for microcontrollers and is available at http://www.psa.es/sdg/sunpos.htm. 

The PSA algorithm uses Universal Time (UT) to remove the uncertainty caused by local time zones.
The location is entered using longitude and latitude with the minutes and seconds converted to fractions of a degree. 
The azimuth angle is measured from true north not magnetic north and the zenith angle is measured from the vertical. 
The elevation angle is measured from the horizontal.
*/

// Declaration of some constants 
#define pi    3.14159265358979323846
#define twopi (2*pi)
#define rad   (pi/180)
#define dEarthMeanRadius     6371.01      // In km
#define dAstronomicalUnit    149597890      // In km

void sunpos(time_t t, double dLatitude, double dLongitude, double *dZenithAngle, double *dAzimuth)
{
  // Main variables
  double dElapsedJulianDays;
  double dDecimalHours;
  double dEclipticLongitude;
  double dEclipticObliquity;
  double dRightAscension;
  double dDeclination;
  // Auxiliary variables
  double dY;
  double dX;

  // Calculate difference in days between the current Julian Day 
  // and JD 2451545.0, which is noon 1 January 2000 Universal Time
  {
    double dJulianDate;
    long int liAux1;
    long int liAux2;
    // Calculate time of the day in UT decimal hours
    dDecimalHours = hour(t) + (minute(t) + second(t) / 60.0 ) / 60.0;
    // Calculate current Julian Day
    liAux1 =(month(t)-14)/12;
    liAux2=(1461*(year(t) + 4800 + liAux1))/4 + (367*(month(t) - 2-12*liAux1))/12- (3*((year(t) + 4900 + liAux1)/100))/4+day(t)-32075;
    dJulianDate=(double)(liAux2)-0.5+dDecimalHours/24.0;
    // Calculate difference between current Julian Day and JD 2451545.0 
    dElapsedJulianDays = dJulianDate-2451545.0;
  }

  // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the 
  // ecliptic in radians but without limiting the angle to be less than 2*Pi 
  // (i.e., the result may be greater than 2*Pi)
  {
    double dMeanLongitude;
    double dMeanAnomaly;
    double dOmega;
    dOmega=2.1429-0.0010394594*dElapsedJulianDays;
    dMeanLongitude = 4.8950630+ 0.017202791698*dElapsedJulianDays; // Radians
    dMeanAnomaly = 6.2400600+ 0.0172019699*dElapsedJulianDays;
    dEclipticLongitude = dMeanLongitude + 0.03341607*sin( dMeanAnomaly ) + 0.00034894*sin( 2*dMeanAnomaly )-0.0001134-0.0000203*sin(dOmega);
    dEclipticObliquity = 0.4090928 - 6.2140e-9*dElapsedJulianDays+0.0000396*cos(dOmega);
  }

  // Calculate celestial coordinates ( right ascension and declination ) in radians 
  // but without limiting the angle to be less than 2*Pi (i.e., the result may be 
  // greater than 2*Pi)
  {
    double dSin_EclipticLongitude;
    dSin_EclipticLongitude= sin( dEclipticLongitude );
    dY = cos( dEclipticObliquity ) * dSin_EclipticLongitude;
    dX = cos( dEclipticLongitude );
    dRightAscension = atan2( dY,dX );
    if( dRightAscension < 0.0 ) dRightAscension = dRightAscension + twopi;
    dDeclination = asin( sin( dEclipticObliquity )*dSin_EclipticLongitude );
  }

  // Calculate local coordinates ( azimuth and zenith angle ) in degrees
  {
    double dGreenwichMeanSiderealTime;
    double dLocalMeanSiderealTime;
    double dLatitudeInRadians;
    double dHourAngle;
    double dCos_Latitude;
    double dSin_Latitude;
    double dCos_HourAngle;
    double dParallax;
    dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283*dElapsedJulianDays + dDecimalHours;
    dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime*15 + dLongitude)*rad;
    dHourAngle = dLocalMeanSiderealTime - dRightAscension;
    dLatitudeInRadians = dLatitude*rad;
    dCos_Latitude = cos( dLatitudeInRadians );
    dSin_Latitude = sin( dLatitudeInRadians );
    dCos_HourAngle= cos( dHourAngle );
    *dZenithAngle = (acos( dCos_Latitude*dCos_HourAngle*cos(dDeclination) + sin( dDeclination )*dSin_Latitude));
    dY = -sin( dHourAngle );
    dX = tan( dDeclination )*dCos_Latitude - dSin_Latitude*dCos_HourAngle;
    *dAzimuth = atan2( dY, dX );
    if ( *dAzimuth < 0.0 ) 
      *dAzimuth = *dAzimuth + twopi;
    *dAzimuth = *dAzimuth/rad;
    // Parallax Correction
    dParallax=(dEarthMeanRadius/dAstronomicalUnit)*sin(*dZenithAngle);
    *dZenithAngle=(*dZenithAngle + dParallax)/rad;
  }
}


//////////////////////////////////////////////
// TIME & DATE OUTPUT FUNCTIONS

void leadingZero(int digits){ 
  if(digits < 10) Serial.print('0');
  Serial.print(digits);
}  

void printDateTime(time_t t){
  printTime(t);
  Serial.print(",");
  printDate(t);
}

void printDate(time_t t){
  leadingZero(day(t));
  Serial.print("/");
  leadingZero(month(t));
  Serial.print("/");
  leadingZero(year(t));

}
void printTime(time_t t){
  leadingZero(hour(t));
  Serial.print(":");
  leadingZero(minute(t));
  Serial.print(":");
  leadingZero(second(t));
}

Here is the link to the NREL data corresponding to the default position and time in the code above.
http://www.nrel.gov/midc/apps/solpos.pl?syear=2010&smonth=11&sday=18&eyear=2010&emonth=11&eday=18&step=10&stepunit=1&latitude=-37.7926&longitude=144.932&timezone=0&press=1013.0&temp=15&aspect=180&tilt=0&solcon=1367&sbwid=7.6&sbrad=31.7&sbsky=0.04&interval=0&field=2&field=10&zip=0

Hmm well it will be less accurate due to the use of floats instead of doubles but I can't see much of a difference apart from that.

Any ideas as to what it is which is throwing it all out?

Mowcius

I can't put my finger on why there is the difference. I came to your code via looking at the PSA site and my friend Google. I saw you had removed / changed a few things like dropping off calculating per second and changing to floats. I decided I wanted to start again with the raw PSA code and this is what I got. I'm tempted to wrap this up in a library with some sunrise /sunset and time zone calculations but this will have to wait for another day. Its well after pumpkin time here in Oz -Andrew

I wonder how important it is to adjust for altitude? I was testing that code a few months ago and found that I had to add an offset to get it to match up with the charts. However, the difference was pretty linear, and once I did that, I got pretty good results. (I'm at ~1800M). It doesn't sound like what Chojin and agmatthews are seeing, but it's just a thought.

I wonder how important it is to adjust for altitude?

Not very I wouldn't have thought.

Yea I read that, but I couldn’t help wondering why my azimuth was always high by ~ 11° compared to the NREL data.