Dear All,
I was experimenting with the great work posted by Mowcius on his Sun Position mathematical formulas implemented in Arduino Sketch. I have problem wrong Azimuth elevation calculations when i put my Location Latitude and Longitude in the code.
Following are the information related to my Locations:
Country: Pakistan
City : Karachi
Locating in : Northern Hemisphere
Time: +5 //Not used in the Mathematical Calculations
Longitude = 67.21
Latitude = 24.86
//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 = 2011; //year
int Month = 9; //month
int Day = 12; //day
float Hours = 13; //hour
float Minutes = 10; //minutes
float Longitude = 67.21; //enter longitude here "KARCHI PAKISTAN"
float Latitude = 24.86; //enter latitude here
Now The following is complete Sketch, I just made above changes only:
//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 = 2011; //year
int Month = 9; //month
int Day = 12; //day
float Hours = 13; //hour
float Minutes = 10; //minutes
float Longitude = 67.21; //enter longitude here "KARCHI PAKISTAN"
float Latitude = 24.86; //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!
}
PROBLEM:
After Uploading the sketch into Arduino and opening the Arduino Serial Monitor I am getting:
Elevation Angle: 6
Azimuth: 272
I think the problem is float precision. I read in the Arduino reference that 'float' and 'double' are the same and support '6 or 7' digits of precision. The julian date calculation uses the constant 2451545.0 which need 8 digits of precision.
You might be able to work around the problem by using 'long int' values and keeping date in days and time in minutes separately.
After Uploading the sketch into Arduino and opening the Arduino Serial Monitor I am getting:
12-September-2011 13:10 UCT
Elevation Angle: 6 Azimuth: 272
Why their is so much difference?and How can i solve the problem?
It looks like the website is using LOCAL TIME and the sketch is using COORDINATED UNIVERSAL TIME. If you look at the values the website produces for 18:00 and 18:15 (local time at 13:00 and 13:15 UTC) it looks like a match for the values you are getting:
18:00 Azimuth: 270 Altitude: 8
18:15 Azimuth: 272 Altitude: 5
I haven't looked at this code for a long time but you really want to have the input time as UT rather than subtracting 5 because it will go from -5:00 to 19:00 which will no doubt not play too happy in the code.
Thank You Mowcius...
I want your code to add LCD and Keypad...The program to be menu driven.. For example, user can input Longitude, Latitude ,UT, date and time..Also the user able to put Gearing ratio, backlash values etc..These values can be stored in EEPROM (It is difficult to place signed integer, and float values but can be sorted out).. The LCD in Idle position shows Current Azimuth , Altitude and the Sunset time etc..
I want to make a Universal open source Suntracker Arduino based..Currently searching the forum for LCD Menu sketch...so far not found a easy sketch...
mowcius:
I haven't looked at this code for a long time but you really want to have the input time as UT rather than subtracting 5 because it will go from -5:00 to 19:00 which will no doubt not play too happy in the code.
Hi Mowcius,
Thank you for your help and the nice well-commented sketch.. My Question:
Using UT 13:15Hrs I am getting:
Elevation Angle: 6
Azimuth: 272
But in my country at 13:15Hrs the elevation and Azimuth found to be : Azimuth:210° Altitude:66°
which is true for sun position..How can i get the actual Azimuth and Elevation of the sun without subtracting 5?The right values of Azimuth and Elevation is required to move both steppers to sun position..
just for information and record.. the Sunset and sunrise timing calculated by the Mowcius sketch is inaccurate to several minutes.. I have to add/subtract few constant to get the values accurate as per NOAA SOLAR CALCULATION excel sheet..
To work easily i converted the long sketch of Mowcius into two function for SUNSET and SUNRISE...
/******************************
// //SUNRISE HOUR CALCULATED HERE
******************************/
void SunRise(){ 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 for me to get UT to British Summer Time
SunRiseMinute = Minutes;
i = 1;
}
Minutes = 0;
Hours++;
}
}
/******************************
// //SUNSET HOUR CALCULATED HERE
******************************/
void SunSet(){
int 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 for me to get UT to British Summer Time
SunSetMinute = Minutes;
i = 1;
}
Minutes = 0;
Hours++;
}
}
Then i added
void loop(){
sunPos();
Serial.print("AZ: ");
Serial.print(Azimuth);
Serial.print(" EL: ");
Serial.print(ElevationAngle);
Serial.print(" Zenith: ");
Serial.print(ZenithAngle);
SunRise();
SunSet();
Serial.println("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-2); //The "-2" here correct the sunrise minutes as per NOAA
Serial.print("Sunset: ");
Serial.print(SunSetHour);
Serial.print(":");
if(SunSetMinute < 10)
Serial.print("0");
Serial.println(SunSetMinute+5); //The "+5" here correct the sunset minutes as per NOAA
Serial.println(SunSetMinute+5);
}
}
Above sketch now gives realistic SUNSET and SUNRISE time, i have checked all the months and upto several years and with different time zone i am getting accurate results...
]
]