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.
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).
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 ;)
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..
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?
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()
#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));
}
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.
