Go Down

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

Austin Duff

About the compile error, I think Arduino requires a setup() and loop() function to initialize. You can leave the setup() empty and use loop() as your main. For example your code may look like this:

#define ....
#define ....

setup()
{
}

loop
{

(Put your main code here)

}

mowcius

Quote
Would some kind of inclinometer be able to tell your UGV if it were driving up a slope?

Yes, but to be honest, the small amount of slope I may be stopping on will not be a lot more than the innaccuracies of removing some decimal places from the calculation (because it is running on an arduino it is not as accurate using floating point values instead of doubles).
I am not too worried about the slight slope. If I stick in an accelerometer and gyroscope which is likely but it might be long after I get it finished so I might then fiddle with that if I get bored!

Quote

incidently, How fast do you think a standard arduino would be able to carry out the entire calculation?

Umm, no idea but best guess of less than 100ms. Faster than required :P

Quote
About the compile error, I think Arduino requires a setup() and loop() function to initialize. You can leave the setup() empty and use loop() as your main. For example your code may look like this:
#define ....
#define ....
setup()
{
}
loop
{
(Put your main code here)
}

It does indeed and I have tried that  :P
First thing I did but thanks for the reply.
That is why I am stumped.

Mowcius

mowcius

#32
Jun 30, 2010, 02:42 pm Last Edit: Jun 30, 2010, 02:51 pm by mowcius Reason: 1
Right.

I had a bit more time so I took another look at that code.

I still wasn't sure what was the issue with it so I removed it from the library and stuck it all in one sketch.

It seems to compile but I have not looked into the finer details yet (such as whether it works) as I have not got an arduino to test it with until I get home.

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

Anyway here it is:
Code: [Select]
#define pi    3.14159265358979323846
#define twopi (2*pi)
#define rad   (pi/180)
#define dEarthMeanRadius     6371.01      // In km
#define dAstronomicalUnit    149597890      // In km

     int iYear = 2010; //year
     int iMonth = 6; //month
     int iDay = 30; //day
     float dHours = 12; //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;


void setup() {
 //do setup stuff
}

void sunPos(){
     // Main variables
     float dElapsedJulianDays;
     float dDecimalHours;
     float dEclipticLongitude;
     float dEclipticObliquity;

     // 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;
           float 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;
     
}

void loop() {
 sunPos();
 Serial.print("Right Ascension:  ");
 Serial.println(dRightAscension);
 Serial.print("Declination:  ");
 Serial.println(dDeclination);
 delay(500);
 
}

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 :)

Mowcius

breaksbassbleeps

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

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?  ;D

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.

Austin Duff

My apologies :) Good luck mow

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

breaksbassbleeps

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


I think you're allowed to be after all the work you've put into this code lol!  ;D

BreaksBassBleeps

mowcius

#41
Jul 02, 2010, 09:29 pm Last Edit: Jul 03, 2010, 09:03 pm by mowcius Reason: 1
Right. After all that, there is some code guaranteed to work now.

Code: [Select]
//Sun Position Calculation by Mowcius (mowcius.co.uk)
//Provides sun position (relative) from static variables

#include <math.h>
#define pi    3.14159265358979323846
#define twopi (2*pi)
#define rad   (pi/180)
#define EarthMeanRadius     6371.01      // In km
#define AstronomicalUnit    149597890      // In km

//Input Variables --------------------- TIME HAS TO BE IN UT (UNIVERSAL TIME)! NO TIME ZONES OR SUMMER TIMES --------
//My last modifications were probably at this time on this date!
     int Year = 2010; //year
     int Month = 7; //month
     int Day = 3; //day
     float Hours = 16; //hour
     float Minutes = 38; //minutes

     float Longitude = 1.2967; //enter longitude here
     float Latitude = 1.5465; //enter latitude here
//--------

//Program Variables
     float ZenithAngle;
     float Azimuth;
       float RightAscension;
     float Declination;
       float Parallax;
       float ElevationAngle;

     float ElapsedJulianDays;
     float DecimalHours;
     float EclipticLongitude;
     float EclipticObliquity;
//--------

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 JulianDate;
           long int liAux1;
           long int liAux2;
           // Calculate time of the day in UT decimal hours
           DecimalHours = Hours + (Minutes / 60.0);
           // Calculate current Julian Day
           liAux1 =(Month-14)/12;
           liAux2=(1461*(Year + 4800 + liAux1))/4 + (367*(Month
                 - 2-12*liAux1))/12- (3*((Year + 4900
           + liAux1)/100))/4+Day-32075;
           JulianDate=(float)(liAux2)-0.5+DecimalHours/24.0;
           // Calculate difference between current Julian Day and JD 2451545.0
           ElapsedJulianDays = JulianDate-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 MeanLongitude;
           float MeanAnomaly;
           float Omega;
           Omega=2.1429-0.0010394594*ElapsedJulianDays;
           MeanLongitude = 4.8950630+ 0.017202791698*ElapsedJulianDays; // Radians
           MeanAnomaly = 6.2400600+ 0.0172019699*ElapsedJulianDays;
           EclipticLongitude = MeanLongitude + 0.03341607*sin( MeanAnomaly )
                 + 0.00034894*sin( 2*MeanAnomaly )-0.0001134
                 -0.0000203*sin(Omega);
           EclipticObliquity = 0.4090928 - 6.2140e-9*ElapsedJulianDays
                 +0.0000396*cos(Omega);
     
     // 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 Sin_EclipticLongitude;
           Sin_EclipticLongitude= sin( EclipticLongitude );
           dY = cos( EclipticObliquity ) * Sin_EclipticLongitude;
           dX = cos( EclipticLongitude );
           RightAscension = atan2( dY,dX );
           if( RightAscension < 0.0 ) RightAscension = RightAscension + twopi;
           Declination = asin( sin( EclipticObliquity )*Sin_EclipticLongitude );
     
     // Calculate local coordinates ( azimuth and zenith angle ) in degrees
     
           float GreenwichMeanSiderealTime;
           float LocalMeanSiderealTime;
           float LatitudeInRadians;
           float HourAngle;
           float Cos_Latitude;
           float Sin_Latitude;
           float Cos_HourAngle;
           GreenwichMeanSiderealTime = 6.6974243242 +
                 0.0657098283*ElapsedJulianDays
                 + DecimalHours;
           LocalMeanSiderealTime = (GreenwichMeanSiderealTime*15
                 + Longitude)*rad;
           HourAngle = LocalMeanSiderealTime - RightAscension;
           LatitudeInRadians = Latitude*rad;
           Cos_Latitude = cos( LatitudeInRadians );
           Sin_Latitude = sin( LatitudeInRadians );
           Cos_HourAngle= cos( HourAngle );
           ZenithAngle = (acos( Cos_Latitude*Cos_HourAngle
                 *cos(Declination) + sin( Declination )*Sin_Latitude));
           dY = -sin( HourAngle );
           dX = tan( Declination )*Cos_Latitude - Sin_Latitude*Cos_HourAngle;
           Azimuth = atan2( dY, dX );
           if ( Azimuth < 0.0 )
             Azimuth = Azimuth + twopi;
           Azimuth = Azimuth/rad;
           // Parallax Correction
           Parallax=(EarthMeanRadius/AstronomicalUnit)
                 *sin(ZenithAngle);
           ZenithAngle=(ZenithAngle //Zenith angle is from the top of the visible sky (thanks breaksbassbleeps)
                 + Parallax)/rad;
               ElevationAngle = (90-ZenithAngle); //Retrieve useful elevation angle from Zenith angle
}

void loop(){
 sunPos(); //Run sun position calculations
 Serial.print("Elevation Angle:  ");
 Serial.println(ElevationAngle, 0); //Print Elevation (Vertical) with no decimal places as accuracy is not really great enough
 Serial.print("Azimuth:  ");
 Serial.println(Azimuth, 0); //Print Azimuth (Horizontal) with no decimal places
 if(ElevationAngle < 0)
 Serial.println("The sun has set. Get some sleep!");
 delay(10000); //Delay 10 seconds - Values aren't going to have changed anyway as they are currently static variables!
}


Let me know if you use the code for anything. I think I will connect the horizontal (Azimuth) values to a stepper motor this weekend and figure out how to use the RTC on my rDuino LEDHead!

The Finish Line (the song by Snow Patrol) just started playing from random... Seems rather appropriate. I can see the finish from here!
Hopefully this code will help other people too.

Mowcius

bHogan

Mowcius,
Nice work! Thank you.
I'm sure this will have value for a lot of projects.
"Data is not information, information is not knowledge, knowledge is not understanding, understanding is not wisdom."
~ Clifford Stoll

mowcius

Well I am watching it now and it seems to be working perfectly.
I am comparing it to stellarium and the sun. The sun's just about to set and it looks about right :)
It doesn't take altitude into consideration but to be honest, for the kind of things people might use it for, that degree of accuracy is not required. That's also why I removed the seconds value. It made hardly any difference and you might as well have the extra RAM   ;)

Mowcius

mowcius

#44
Jul 03, 2010, 11:10 am Last Edit: Jul 03, 2010, 09:03 pm by mowcius Reason: 1
Had a bit of a play around and wrote some code to calculate the sunrise and sunset for the next week.

Put in the date and it'll throw out the sunrise and sunset.

It runs through every minute until the elevation gets to ~0 then records the value. Takes about 1 second to calculate each day.

Code: [Select]
//Sun Position Calculation by Mowcius (mowcius.co.uk)
//Provides sunrise and sunset times for a week after a given date

#include <math.h>
#define pi    3.14159265358979323846
#define twopi (2*pi)
#define rad   (pi/180)
#define EarthMeanRadius     6371.01      // In km
#define AstronomicalUnit    149597890      // In km

//Input Variables --------------------- TIME HAS TO BE IN UT (UNIVERSAL TIME)! NO TIME ZONES OR SUMMER TIMES --------
//My last modifications were probably at this time on this date!
     int Year = 2010; //year
     int Month = 7; //month
     int Day = 3; //day
     float Hours = 10; //hour
     float Minutes = 10; //minutes

     float Longitude = 1.5458; //enter longitude here
     float Latitude = 1.5464; //enter latitude here
//--------

   int SunRiseHour = 0;
   int SunRiseMinute = 0;
   int SunSetHour = 0;
   int SunSetMinute = 0;
   int d = 0;

//Program Variables
     float ZenithAngle;
     float Azimuth;
       float RightAscension;
     float Declination;
       float Parallax;
       float ElevationAngle;

     float ElapsedJulianDays;
     float DecimalHours;
     float EclipticLongitude;
     float EclipticObliquity;
//--------

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 JulianDate;
           long int liAux1;
           long int liAux2;
           // Calculate time of the day in UT decimal hours
           DecimalHours = Hours + (Minutes / 60.0);
           // Calculate current Julian Day
           liAux1 =(Month-14)/12;
           liAux2=(1461*(Year + 4800 + liAux1))/4 + (367*(Month
                 - 2-12*liAux1))/12- (3*((Year + 4900
           + liAux1)/100))/4+Day-32075;
           JulianDate=(float)(liAux2)-0.5+DecimalHours/24.0;
           // Calculate difference between current Julian Day and JD 2451545.0
           ElapsedJulianDays = JulianDate-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 MeanLongitude;
           float MeanAnomaly;
           float Omega;
           Omega=2.1429-0.0010394594*ElapsedJulianDays;
           MeanLongitude = 4.8950630+ 0.017202791698*ElapsedJulianDays; // Radians
           MeanAnomaly = 6.2400600+ 0.0172019699*ElapsedJulianDays;
           EclipticLongitude = MeanLongitude + 0.03341607*sin( MeanAnomaly )
                 + 0.00034894*sin( 2*MeanAnomaly )-0.0001134
                 -0.0000203*sin(Omega);
           EclipticObliquity = 0.4090928 - 6.2140e-9*ElapsedJulianDays
                 +0.0000396*cos(Omega);
     
     // 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 Sin_EclipticLongitude;
           Sin_EclipticLongitude= sin( EclipticLongitude );
           dY = cos( EclipticObliquity ) * Sin_EclipticLongitude;
           dX = cos( EclipticLongitude );
           RightAscension = atan2( dY,dX );
           if( RightAscension < 0.0 ) RightAscension = RightAscension + twopi;
           Declination = asin( sin( EclipticObliquity )*Sin_EclipticLongitude );
     
     // Calculate local coordinates ( azimuth and zenith angle ) in degrees
     
           float GreenwichMeanSiderealTime;
           float LocalMeanSiderealTime;
           float LatitudeInRadians;
           float HourAngle;
           float Cos_Latitude;
           float Sin_Latitude;
           float Cos_HourAngle;
           GreenwichMeanSiderealTime = 6.6974243242 +
                 0.0657098283*ElapsedJulianDays
                 + DecimalHours;
           LocalMeanSiderealTime = (GreenwichMeanSiderealTime*15
                 + Longitude)*rad;
           HourAngle = LocalMeanSiderealTime - RightAscension;
           LatitudeInRadians = Latitude*rad;
           Cos_Latitude = cos( LatitudeInRadians );
           Sin_Latitude = sin( LatitudeInRadians );
           Cos_HourAngle= cos( HourAngle );
           ZenithAngle = (acos( Cos_Latitude*Cos_HourAngle
                 *cos(Declination) + sin( Declination )*Sin_Latitude));
           dY = -sin( HourAngle );
           dX = tan( Declination )*Cos_Latitude - Sin_Latitude*Cos_HourAngle;
           Azimuth = atan2( dY, dX );
           if ( Azimuth < 0.0 )
             Azimuth = Azimuth + twopi;
           Azimuth = Azimuth/rad;
           // Parallax Correction
           Parallax=(EarthMeanRadius/AstronomicalUnit)
                 *sin(ZenithAngle);
           ZenithAngle=(ZenithAngle //Zenith angle is from the top of the visible sky (thanks breaksbassbleeps)
                 + Parallax)/rad;
               ElevationAngle = (90-ZenithAngle); //Retrieve useful elevation angle from Zenith angle
}

void loop(){
 Hours = 3;
 Minutes = 0;
 sunPos(); //Run sun position calculations
 int i = 0;
 while(i < 1){
 while(Minutes < 60 && ElevationAngle < -0.5){
   Minutes++;
   sunPos();
 }
 if(ElevationAngle < 0.5 && ElevationAngle > -0.5){
   SunRiseHour = Hours+1;                               //+1 for me to get UT to British Summer Time
   SunRiseMinute = Minutes;
   i = 1;
 }
 Minutes = 0;
 Hours++;
}
i = 0;
Hours = 16;
Minutes = 0;
sunPos();
while(i < 1){
while(Minutes < 60 && ElevationAngle > 0.5){
   Minutes++;
   sunPos();
 }
 if(ElevationAngle < 0.5 && ElevationAngle > -0.5){
   SunSetHour = Hours+1;                                //+1 for me to get UT to British Summer Time
   SunSetMinute = Minutes;
   i = 1;
 }  
 Minutes = 0;
 Hours++;
}
Serial.print("Date: ");
Serial.print(Day);
Serial.print("/");
Serial.print(Month);
Serial.print("/");
Serial.println(Year);
Serial.print("Sunrise: ");
Serial.print(SunRiseHour);
Serial.print(":");
if(SunRiseMinute < 10)
Serial.print("0");
Serial.println(SunRiseMinute);
Serial.print("Sunset: ");
Serial.print(SunSetHour);
Serial.print(":");
if(SunSetMinute < 10)
Serial.print("0");
Serial.println(SunSetMinute);
Serial.print("Sketch Time Taken: ");
Serial.println(millis());
Day++;
d++;
if(d == 7){
while(i == 1){}
}
}

Go Up