Pages: [1]   Go Down
 Author Topic: (SOLVED)Problems in "Mowcius Sun Position" Arduino Sketch  (Read 1380 times) 0 Members and 1 Guest are viewing this topic.
Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « on: September 12, 2011, 08:16:41 am » Bigger Smaller Reset

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

Code:
//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:
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 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 HourAngle;
float Cos_Latitude;
float Sin_Latitude;
float Cos_HourAngle;
GreenwichMeanSiderealTime = 6.6974243242 +
0.0657098283*ElapsedJulianDays
+ DecimalHours;
LocalMeanSiderealTime = (GreenwichMeanSiderealTime*15
HourAngle = LocalMeanSiderealTime - RightAscension;
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;
// Parallax Correction
*sin(ZenithAngle);
ZenithAngle=(ZenithAngle //Zenith angle is from the top of the visible sky (thanks breaksbassbleeps)
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

Whereas the sunposition calculators on the website gave me following results:
http://www.sunposition.info/sunposition/spc/locations.php#1

Date: 12th Sept, 2011 time 13:15   Azimuth:210°   Altitude:66°

Why their is so much difference?and How can i solve the problem?
 Solar Position.xls (277.5 KB - downloaded 4 times.) « Last Edit: September 13, 2011, 03:34:33 am by Khalid » Logged

Massachusetts, USA
Offline
Tesla Member
Karma: 168
Posts: 7711
 « Reply #1 on: September 12, 2011, 11:26:55 am » Bigger Smaller Reset

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

Send Bitcoin tips to: 1L3CTDoTgrXNA5WyF77uWqt4gUdye9mezN
Send Litecoin tips to : LVtpaq6JgJAZwvnVq3ftVeHafWkcpmuR1e

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #2 on: September 12, 2011, 08:51:21 pm » Bigger Smaller Reset

Hi John,
Reference following thread on which Page No.03 bottom Mowcius provided fully working code and many persons their used that sketch without any trouble:
http://www.arduino.cc/cgi-bin/yabb2/YaBB.pl?num=1274361815/30

However, I have no good experience with programming.. Can someone elaborate how the changes in the sketch can be made to work successfully? Any help?
 Logged

Massachusetts, USA
Offline
Tesla Member
Karma: 168
Posts: 7711
 « Reply #3 on: September 12, 2011, 09:58:11 pm » Bigger Smaller Reset

Karachi, Pakistan  N24.86,E67.21

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

http://www.sunposition.info/sunposition/spc/locations.php#1
Date: 12th Sept, 2011 time 13:15   Azimuth:210°   Altitude:66°

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
 Logged

Send Bitcoin tips to: 1L3CTDoTgrXNA5WyF77uWqt4gUdye9mezN
Send Litecoin tips to : LVtpaq6JgJAZwvnVq3ftVeHafWkcpmuR1e

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #4 on: September 12, 2011, 10:54:31 pm » Bigger Smaller Reset

Hi John ,
Intelligent you are:)...My Salute to you for tracking the problem..I want to confirm Should i add LocalTime+5 in the Sketch?
Code:
// Calculate time of the day in UT decimal hours
DecimalHours = Hours + 5+(Minutes / 60.0); //Where HOURS is from RTC

What is the better way?..Please need some more help from you in this:)
 « Last Edit: September 12, 2011, 11:05:39 pm by Khalid » Logged

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #5 on: September 13, 2011, 02:39:47 am » Bigger Smaller Reset

Result does match by subtracting 5 in following code.
Code:
// Calculate time of the day in UT decimal hours
DecimalHours = Hours - 5+(Minutes / 60.0); //Where HOURS is from RTC
Thank You John for all of your help:)
 Logged

North Yorkshire, UK
Offline
Karma: 104
Posts: 5531
 « Reply #6 on: September 16, 2011, 01:41:55 am » Bigger Smaller Reset

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

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #7 on: September 16, 2011, 01:55:10 pm » Bigger Smaller Reset

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

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #8 on: September 16, 2011, 02:34:46 pm » Bigger Smaller Reset

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

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #9 on: September 27, 2011, 01:05:00 pm » Bigger Smaller Reset

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...
Code:
/******************************
// //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++;
}
}

Code:
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...
 Logged

Global Moderator
Dallas
Offline
Shannon Member
Karma: 159
Posts: 11401
 « Reply #10 on: September 28, 2011, 02:18:46 pm » Bigger Smaller Reset

Thank you for the follow-up and posting your code.

Have you contacted mowcius about the changes?
 « Last Edit: September 28, 2011, 02:21:38 pm by Coding Badly » Logged

Pakistan
Offline
Sr. Member
Karma: 6
Posts: 325
Arduino rocks
 « Reply #11 on: June 08, 2012, 12:12:35 pm » Bigger Smaller Reset

Hi Coding Badly, I think Mowcius is much busy in other projects and shall amazed to see how badly i coded his sketch...
 Logged

 Pages: [1]   Go Up