Go Down

Topic: New sun tracker idea? (Read 6 times) previous topic - next topic

mowcius

Quote
An hour and a half until I can get home and try this out on my mega!! Oh the pain. lol

Yeah well I just got home, tried it and it's wrong. Correct code at the bottom :P

Quote
Quote:
You will need to then work out what the right ascension and declination mean for positioning but that shouldn't be too hard.


I thought it was the Azimuth & Zenith angle which gave the position of the sun in degrees.. Is Wikipedia lying to me again?  Grin

Yes indeed it is. I was a bit tired. I realised that when I tried it...

Quote
Quote:
It needs some work and cleaning up a bit but if you set it to take the time from an RTC then it should then track the sun


I was thinking of getting the datalogging shield to sit on my motor shield. Just so happend that it had a RTC on it and I can log various things such as light and power output to an SD card.

I have tidied it up a bit now.
That sounds good with the RTC. I might try it with one at the weekend.

New code:
Code: [Select]
#include <math.h>
#define pi    3.14159265358979323846
#define twopi (2*pi)
#define rad   (pi/180)
#define dEarthMeanRadius     6371.01      // In km
#define dAstronomicalUnit    149597890      // In km

int i = 0;

     int iYear = 2010; //year
     int iMonth = 6; //month
     int iDay = 1; //day
     float dHours = 1; //hour
     float dMinutes = 6; //minutes
     float dSeconds = 30; //seconds



     float dLongitude = 1.4356; //enter longitude here
     float dLatitude = 1.496; //enter latitude here


     float dZenithAngle;
     float dAzimuth;
       float dRightAscension;
     float dDeclination;
       float dParallax;

// Main variables
     float dElapsedJulianDays;
     float dDecimalHours;
     float dEclipticLongitude;
     float dEclipticObliquity;

void setup() {
Serial.begin(9600);
}

void sunPos(){
     

     // Auxiliary variables
     float dY;
     float dX;

     // Calculate difference in days between the current Julian Day
     // and JD 2451545.0, which is noon 1 January 2000 Universal Time
     
           float dJulianDate;
           long int liAux1;
           long int liAux2;
           // Calculate time of the day in UT decimal hours
           dDecimalHours = dHours + (dMinutes
                 + dSeconds / 60.0 ) / 60.0;
           // Calculate current Julian Day
           liAux1 =(iMonth-14)/12;
           liAux2=(1461*(iYear + 4800 + liAux1))/4 + (367*(iMonth
                 - 2-12*liAux1))/12- (3*((iYear + 4900
           + liAux1)/100))/4+iDay-32075;
           dJulianDate=(float)(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)
     
           float dMeanLongitude;
           float dMeanAnomaly;
           float 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)
     
           float 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
     
           float dGreenwichMeanSiderealTime;
           float dLocalMeanSiderealTime;
           float dLatitudeInRadians;
           float dHourAngle;
           float dCos_Latitude;
           float dSin_Latitude;
           float dCos_HourAngle;
           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;
     
}

void loop(){
 sunPos(); //Run sun position calculations
 Serial.print("Zenith Angle:  ");
 Serial.println(dZenithAngle);
 Serial.print("Azimuth:  ");
 Serial.println(dAzimuth);
 Serial.print("Julian Days:  ");
 Serial.println(dElapsedJulianDays);
 dHours = dHours++;
 delay(1000);
 
}


Mowcius

breaksbassbleeps

Like Richard Hammond when he said he was a "driving god" I believe you Mowcius can also say the same about yourself about Arduino code  :D hehe

As the time and date are hardcoded as there's no RTC, it should generate only one set of results. Yet after the first results are generated it starts creating incorrect results on a loop. I was just wondering what's causing this.

Right, time to get my RTC bought and wire up some steppers :)

Thanks again for all your help Mowcius..

mowcius

Quote
As the time and date are hardcoded as there's no RTC, it should generate only one set of results. Yet after the first results are generated it starts creating incorrect results on a loop. I was just wondering what's causing this.


 dHours = dHours++;

Can you see it yet?  :P
I just told it to keep increasing the hour by 1 each loop so I could see the value changing. Feel free to delete that so it works properly.

I also included the Julain date in the serial print just for fun.

The code seems to work for the Azimuth angle (I checked it using stellarium) but the Zenith angle seems to be rather off for my position.
Not sure why this is yet. I am looking into it. If you don't want to go too complicated and don't care about the up/down movement then it'll work fine currently but I will try to find out what the issue is. Maybe my position was entered in the wrong format or something. Hmm :-?

A quick comment though. The time has to be in UT (Universal time - no time zones or summer times etc) for it to work but that's not difficult to do...

Mowcius

breaksbassbleeps

I think I can see the problem.
Elevation is measured from the horizon up  to the sun, whereas the Zenith is measured from straight up, down to the sun.
I've added a tiny bit to your code and the results match what's generated on the NOAA site.

http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html

Code: [Select]
       float dElevationAngle;
&
Code: [Select]
           dZenithAngle=(dZenithAngle
                 + dParallax)/rad;
               dElevationAngle=(90-dZenithAngle);

}

void loop(){
 sunPos(); //Run sun position calculations
 Serial.print("Elevation Angle:  ");
 Serial.println(dElevationAngle);
 Serial.print("Azimuth:  ");
 Serial.println(dAzimuth);
 Serial.print("Julian Days:  ");
 Serial.println(dElapsedJulianDays);
 dMinutes = dMinutes++;
 delay(10000);

mowcius

Quote
I think I can see the problem.
Elevation is measured from the horizon up  to the sun, whereas the Zenith is measured from straight up, down to the sun.
I've added a tiny bit to your code and the results match what's generated on the NOAA site.

That seems to make sense with the values I was getting out.  :)

Thanks for the insight. I thought that the Zenith was bottom up but we are all wrong sometimes  :P

Mowcius

Go Up