next topic »

[Reply #30 on:]

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)

}

[Reply #31 on:]

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!

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

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  
First thing I did but thanks for the reply.
That is why I am stumped.

Mowcius

[Reply #32 on:]

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:

#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

[Reply #33 on:]

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

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

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.

[Reply #34 on:]

My apologies Good luck mow

[Reply #35 on:]

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

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

#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

[Reply #36 on:]

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   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..

[Reply #37 on:]

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

[Reply #38 on:]

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:

       float dElevationAngle;

&

Code:

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

[Reply #39 on:]

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  

Mowcius

[Reply #40 on:]

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

[Reply #41 on:]

Right. After all that, there is some code guaranteed to work now.

Code:

//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

[Reply #42 on:]

Mowcius,
Nice work! Thank you.
I'm sure this will have value for a lot of projects.

[Reply #43 on:]

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

[Reply #44 on:]

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:

//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){}
}
}