Go Down

Topic: Problem with calculations (Sunrise/Sunset times) (Read 200 times) previous topic - next topic

corinreveck

Hey, I'm trying to make a sunrise/sunset calculator for my arduino project but I have a problem with calculations.
My c++ code works fine and is able to tell sunrise/sunset time with +-1 accuracy. But when trying to rewrite the code into arduino I have a problem. 
C++ code:
Code: [Select]
#include <iostream>


using namespace std;
int main(){
int Y=2021;
int M=1;
int D=22;

double J=367*Y-int(7*(Y+int((M+9)/12))/4)+int(275*M/9)+D-730531.5;
cout<<J;
return 0;
}

With that code the output is 7691.5 which is a correct and allows for further calculations to proceed.
Arduino code:
Code: [Select]
void setup(){
  Serial.begin(9600);
}
void loop(){
 
int Y=2021;
int M=1;
int D=22;

double J=367*Y-int(7*(Y+int((M+9)/12))/4)+int(275*M/9)+D-730531.5;
  Serial.println("J:");
  Serial.println(J);
}

With that code the output is -713204.50 which is incorrect. 
Do you have any ideas what could be happening here and how to fix that? 
Of course its only a part of the code but since the first line in the code is already wrong I only posted that.

TheMemberFormerlyKnownAsAWOL

You could break your calculation down into individual expressions and see exactly where you're going wrong.
Hint: 16 bit int.
Please don't PM technical questions - post them on the forum, then everyone benefits/suffers equally

aarg

Probably, your int expression overflowed because the "C++" version (but Arduino is C++!!) which I assume runs on a PC or server, uses a 32 bit int, and we know that the AVR implementation has only a 16 bit int.

For readability, you should break that expression down. Not just to find the problem.
  ... with a transistor and a large sum of money to spend ...
Please don't PM me with technical questions. Post them in the forum.

corinreveck

#3
Jan 22, 2021, 10:37 pm Last Edit: Jan 22, 2021, 10:38 pm by corinreveck
Code: [Select]

void setup(){
  Serial.begin(9600);
}
void loop(){
 
int Y=2021;
int M=1;
int D=22;

double J=367*Y-int(7*(Y+int((M+9)/12))/4)+int(275*M/9)+D-730531.5;
double a=367*Y;
double b=Y+int((M+9)/12);
double c=int(7*b/4);
double d=int(275*M/9);
double e=a-c+d+D-730531.5;
Serial.println("a:");
  Serial.println(a);
  Serial.println("b:");
  Serial.println(b);
  Serial.println("c:");
  Serial.println(c);
  Serial.println("d:");
  Serial.println(d);
  Serial.println("e:");
  Serial.println(e);
  Serial.println("J:");
  Serial.println(J);
}

It fails at calculating 'a' which is 367*2021(year). Since int/double is capped at 32.767 and the expected output is bigger. I tried using long(swapped double->long for 'a' and 'e') for it since it can hold up to 2,147,483,647. It didn't help me much and the output is still the same.

TheMemberFormerlyKnownAsAWOL

Please don't PM technical questions - post them on the forum, then everyone benefits/suffers equally

jremington

#5
Jan 22, 2021, 10:43 pm Last Edit: Jan 22, 2021, 10:44 pm by jremington
One AVR based Arduinos, double is the same as float (32 bits), which is not accurate enough for most sorts of astronomical calculations, unless special steps are taken.

Which Arduino do you have?

aarg

Quote
able to tell sunrise/sunset time with +-1 accuracy
+/- 1 what? Year? Microsecond? :)
  ... with a transistor and a large sum of money to spend ...
Please don't PM me with technical questions. Post them in the forum.

jremington

#7
Jan 22, 2021, 10:55 pm Last Edit: Jan 22, 2021, 10:58 pm by jremington
The posted code just calculates the Julian day, since 1/1/2000.

For accurate solar position calculations on any Arduino, for any time of day, consider using the Solar Position Library.

corinreveck

I'm on arduino Uno,
+/- 1 minute so that kind of accuracy is fine for me. Thanks for the help I didn't change Year value to double.
Here is my code:
 
Code: [Select]
  void loop(){
  double R=2021;
  double M=1;
  double D=22;
  double Req=-0.833;
  double J=367*R-int(7*(R+int((M+9)/12))/4)+int(275*M/9)+D-730531.5;
  double Cent = J/36525;
  double L=fmod((4.8949504201433+628.331969753199*Cent),6.28318530718);
  double G=fmod((6.2400408+628.3019501*Cent),6.28318530718);
  double O=0.409093-0.0002269*Cent;
  double F=0.033423*sin(G)+0.00034907*sin(2*G);
  double E=(0.0430398*sin(2*L+F))-0.00092502*sin(4*(L+F))-F;
  double A=asin(sin(O)*sin(L+F));
  double C=(sin(0.017453293*Req)-sin(0.017453293*Lat)*sin(A))/(cos(0.017453293*Lat)*cos(A));
  double Wsch=(3.14-(E+0.017453293*Long+1*acos(C)))*57.29577951/15;
  double Zach=(3.14-(E+0.017453293*Long-1*acos(C)))*57.29577951/15;
  double WschA=int(Wsch);
  double WschB=round((Wsch-WschA)*60);
  double ZachA=int(Zach);
  double ZachB=round((Zach-ZachA)*60);
}


Wsch=Sunrise, A=hour,B=minute
Zach=Sunset, A=hour,B=minute

TheMemberFormerlyKnownAsAWOL

Quote
I'm on arduino Uno,
You could've saved some keystrokes - "float" is shorter than "double"
Please don't PM technical questions - post them on the forum, then everyone benefits/suffers equally

aarg

See reply #1. You didn't break it down at all.
  ... with a transistor and a large sum of money to spend ...
Please don't PM me with technical questions. Post them in the forum.

jremington

#11
Jan 23, 2021, 04:50 am Last Edit: Jan 23, 2021, 05:45 pm by jremington
The calculation for J below is correct only for integer multiplication and division, although it will sometimes give the correct answer when floats are used for dates.

Code: [Select]
  double R=2021;
  double M=1;
  double D=22;
  double J=367*R-int(7*(R+int((M+9)/12))/4)+int(275*M/9)+D-730531.5;

aarg

So, on the use of the library that was mentioned for solar position, if your existing method compares rise and set times, it would have to be changed because the library only yields elevation data and doesn't have reverse functions to calculate rise and set times. However, in the internal logic of your program, as I have done with my own use of it, you could directly compare the current sun position to the current time, instead of a calculated time of rise/set.

I wrote the library, paying close attention to the calculations, but as not completely immersed in astronomical calculations, I didn't understand them well. In the re-write (the math was written by a far more able person than I) I paid the most attention to calculation and algorithmic integrity, i.e. making sure the authors intent was accurately carried out. There was really, just a lot of re-formatting and re-naming. Then it was wrapped as a class.

I have heard wind of some other code that yields rise/set times. At least, consider, the Solar Position library can be used in a polling mode to determine sunrise/sunset. I actually used it for a clock, the display and backlight were controlled by the time, instead of by some sensor. So for example, I would turn the lights on strongly blue when the sun is 10 degrees below the horizon because that is twilight. Then some other colour until 20 degrees above, and so on. Perhaps your application can be changed so that it polls the sun position rather than comparing the time with some calculated rise/set times.
  ... with a transistor and a large sum of money to spend ...
Please don't PM me with technical questions. Post them in the forum.

jremington

#13
Jan 23, 2021, 05:49 pm Last Edit: Jan 23, 2021, 09:17 pm by jremington
The code posted in reply #8 will not work properly.

Interestingly, essentially the same code was posted a number of years ago on a Polish forum, because it gave the wrong results (no fix was proposed). It is amazing how long such useless bits are passed around the web. They have a life of their own!

Arduino forum member jurs posted sunrise/sunset code at https://forum.arduino.cc/index.php?topic=313942.0 (.zip file attached to reply #11 and modified version here, which works).
Code: [Select]
// Arduino library version of sunrise/sunset calculations by 'jurs' based on:
// A simple C++ program calculating the sunrise and sunset for
// the current date and a set location (latitude,longitude)
// e-mail: jjlammi@yahoo.com
// C++ program calculating the sunrise and sunset for
// the current date and a constant location(latitude,longitude)
// Jarmo Lammi 1999 - 2004
// Last update
// Jan 3rd, 2004 - change to compile with gcc-3.3.2 by Stephan Wynhoff
// http://wynhoff.home.cern.ch/wynhoff/weather/gex.html
// It compiles well also with old gcc version 2.95.4 20011002 (Debian prerelease)
// Try compiling in Linux with command g++ -o rscalc rscalc.cc
// Also the time formats are tidy now.

#define SUNDIAMETER 0.53  // Sun diameter in degrees
#define AIRREFRACTION (34.0/60.0) // athmospheric refraction in degrees

double FNday (int y, int m, int d, float h)
{
  //   Get the days to J2000
  //   h is UT in decimal hours
  //   FNday only works between 1901 to 2099 - see Meeus chapter 7
  long int count = - 7 * (y + (m + 9)/12)/4 + 275*m/9 + d;
  // type casting necessary on PC DOS and TClite to avoid overflow
  count+= (long int)y*367;
  return (double)count - 730531.5 + h/24.0;
}

double FNrange (double x)
{
  //  the function returns an angle in the range 0 to 2*PI
  double b = x / TWO_PI;
  double a = TWO_PI * (b - (long)(b));
  if (a < 0) a = TWO_PI + a;
  return a;
}

double f0(double lat, double declin)
{
  // Calculating the hourangle
  double fo,dfo;
  // Correction: different sign at S HS
  dfo = DEG_TO_RAD*(0.5*SUNDIAMETER + AIRREFRACTION); if (lat < 0.0) dfo = -dfo;
  fo = tan(declin + dfo) * tan(lat*DEG_TO_RAD);
  if (fo>0.99999) fo=1.0; // to avoid overflow //
  fo = asin(fo) + PI/2.0;
  return fo;
}

double f1(double lat, double declin)
{
  // Calculating the hourangle for twilight times
  double fi,df1;
  // Correction: different sign at S HS
  df1 = DEG_TO_RAD * 6.0; if (lat < 0.0) df1 = -df1;
  fi = tan(declin + df1) * tan(lat*DEG_TO_RAD);
  if (fi>0.99999) fi=1.0; // to avoid overflow //
  fi = asin(fi) + PI/2.0;
  return fi;
}

double FNsun (double d, double &L)
{
  //   Find the ecliptic longitude of the Sun
  double g;
  //   mean longitude of the Sun
  L = FNrange(280.461 * DEG_TO_RAD + .9856474 * DEG_TO_RAD * d);

  //   mean anomaly of the Sun
  g = FNrange(357.528 * DEG_TO_RAD + .9856003 * DEG_TO_RAD * d);

  //   Ecliptic longitude of the Sun
  return FNrange(L + 1.915 * DEG_TO_RAD * sin(g) + .02 * DEG_TO_RAD * sin(2 * g));
};


void sunCalc(int year, int month, int day, int timezone, double latitude, double longitude, int &sunrise, int &sunset, int &dawn, int &dusk)
// calculates times of morning dawn, sunrise, sunset, evening dusk
{
  double L, daylen;
  double h = 12; // assume sun position at high noon
  double d = FNday(year, month, day, h);

  //   Use FNsun to find the ecliptic longitude of the Sun
  double lambda = FNsun(d, L);

  //   Obliquity of the ecliptic
  double obliq = 23.439 * DEG_TO_RAD - .0000004 * DEG_TO_RAD * d;

  //   Find the RA and DEC of the Sun
  double alpha = atan2(cos(obliq) * sin(lambda), cos(lambda));
  double delta = asin(sin(obliq) * sin(lambda));

  // Find the Equation of Time
  // in minutes
  // Correction suggested by David Smith
  double LL = L - alpha;
  if (L < PI) LL += TWO_PI;
  double equation = 1440.0 * (1.0 - LL / TWO_PI);
  double ha = f0(latitude,delta);
  double hb = f1(latitude,delta);
  double twx = hb - ha; // length of twilight in radians
  twx = 12.0*twx/PI;    // length of twilight in hours
  // Conversion of angle to hours and minutes //
  daylen = RAD_TO_DEG*ha/7.5;
  if (daylen<0.0001) {daylen = 0.0;}
  // arctic winter //

  double riset = 12.0 - 12.0 * ha/PI + timezone - longitude/15.0 + equation/60.0;
  double settm = 12.0 + 12.0 * ha/PI + timezone - longitude/15.0 + equation/60.0;
 // double noont = riset + 12.0 * ha/PI;
  double altmax = 90.0 + delta * RAD_TO_DEG - latitude;
  // Correction for S HS suggested by David Smith
  // to express altitude as degrees from the N horizon
  if (latitude < delta * RAD_TO_DEG) altmax = 180.0 - altmax;

  double twam = riset - twx;      // morning twilight begin
  double twpm = settm + twx;      // evening twilight end

  if (riset > 24.0) riset-= 24.0;
  if (settm > 24.0) settm-= 24.0;

  sunrise=round(riset*60); //minutes past midnight
  sunset=round(settm*60);
  dawn=round(twam*60);
  dusk=round(twpm*60);
}

void setup() {
Serial.begin(9600);
while(!Serial);
double myLat = ?.?;  //south latitude negative
double myLon = ?.?; //west longitude negative
int TZ = ?? ; //west negative
int sunrise,sunset,dawn,dusk;
sunCalc(2021, 1, 23, TZ, myLat, myLon, sunrise, sunset, dawn, dusk);
char buf[30];
snprintf(buf,sizeof(buf),"Sunrise %02d:%02d, Sunset %02d:%02d",sunrise/60,sunrise%60,sunset/60,sunset%60);
Serial.println(buf);
}

void loop() {}

6v6gt

I use the TimeLord library for obtaining sunrise/sunset at my location.
https://github.com/probonopd/TimeLord

Actually I use it for a time switch in an ESP8266 development. It claims an accuracy of +/- 5 minutes and I verified this by simulation, running through a couple of years and comparing the results from various on line sources. I use the library in UTC and use the native ESP8266 time functions for local time conversion.

In an earlier development, I simply put 12 sunrise and sunset time pairs in a table (22nd of each month) and did a linear interpolation for the intermediate results. These times don't change substantially from year to year. I seem to remember that the accuracy was around +/- 12 minutes per day which should be good enough for any non critical application.

 Incidentally, If you live near North/South pole you may not even see a sunrise or sunset every day.

Go Up