So I can see that many people are becoming interested in the calculation of sunrise/sunset and moon phase for automating lighting in our systems. In particular I think with respect to PWM control over LED dimming. I have adopted a very nice piece of code and loaded it onto my PWM expansion controller and its working great, but the existing code (Deckoz2302 Progressive Sunphase thread) is locked into a particular lat/lon coordinate. I have no earthly idea why I care, but I was really curious about enabling the simple input of Lat/Lon coords and getting back correct sunrise/sunset/moon phase for that coordinate set. In the end you probably want to select your own longitude and then the latitude of the area of the world you want to emulate (i.e. great barrier reef etc etc).
So, I did a bunch of reading and working in Excel and figuring out the math (not overly complex but not trivial) when I found that Michael Rice of www.swfltek.com had basically allready done the entire thing, optimized the calculations to run on Arduino processors (floating point can be an issue) and posted the entire thing as a third generation of working code. So I have been corresponding with him and its available and he asks that this be posted regarding the source code... so here goes his adviso... we can do anything we want with it no charge... but if were going to use it we should acknowledge his work as its significant to get something like this working well and on arduino.
The code provided on this web site is provided in a free and open manner, and may be used in open or proprietary applications
under the terms of the license below.
© 2010 Michael Rice, SWFLTEK
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the
distribution.
* Neither the name of the copyright holders nor the names of
contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Here are the .h and .cpp files I generated from his code.
you can call it whatever you want but just remember to reference the correct file name and move the files into the library location in your set up.
.cpp file
Code: Select all
//==================================================================================================
#include <stdlib.h>
#include <math.h>
// arc-seconds per radian
#define _sec_rad 206264.806247096370813
// axial tilt of earth at epoch, in radians
#define _tilt 0.409092804222329
// tropical year in seconds... rounding error accumulates to 26 seconds by the year 2136
#define _tropical_year 31556925
// 'zenith' of rising (setting) sun in radians (360 - 2 * 90.833 degrees)
#define _zenith 3.11250383272322
//#include "SWFLTEK_Ephemera.h"
/*------------------------------------------------------------------------------------------------
convert degrees to seconds of arc
*/
// decimal degrees
long ddToSeconds(float dd){
return dd * 3600.0;
}
//Degrees, minutes, seconds
long dmsToSeconds(int d, unsigned char m, unsigned char s){
long ret;
ret = labs((long)d);
ret = ret * 3600L + 60L * m + s;
ret = (d<0L) ? -ret : ret;
return ret;
}
/* ------------------------------------------------------------------------------------------------
'Equation of Time'
We use the 'short for' equation, which has a theoretical accuracy of about 40 seconds.
The returned value is in seconds.
*/
int equation_of_time(unsigned long dt){
double t;
dt -= 192540UL; // refer to Jan 3 2000 05:29 (first periapsis)
dt %= _tropical_year;
t = dt;
t /= _tropical_year;
t *= 6.283185307179586;
t = -459.27672 * sin(t) + 575.333472 * sin(2.0 * t + 3.588414);
return t;
}
/*
'Solar Noon' adjusts the passed time stamp to the time (GMT) of local solar noon.
The accuracy is about 40 seconds (set by the equation of time).
*/
void SolarNoon(unsigned long * dt){
long r;
// Set stamp to noon GMT
*dt /= 86400UL;
*dt *= 86400UL;
*dt += 43200UL;
// adjust for equation of time, at noon GMT
*dt -= equation_of_time(*dt);
// rotate to our longitude
r = longitude / 15L;
*dt -= r;
}
/* -----------------------------------------------------------------------------------------------
'Solar Declination'
Returns declination in radians
Accurate to within 50 arc-seconds
*/
double SolarDeclination(unsigned long dt){
double y;
dt %= _tropical_year;
y = dt;
y /= _tropical_year; // fractional year
y *= 6.283185307179586;
y=0.006918-0.399912*cos(y)+0.070257*sin(y)-0.006758*cos(y*2)+0.000907*sin(y*2)-0.002697*cos(y*3)+0.00148*sin(y*3);
return y;
}
/* ------------------------------------------------------------------------------------------------
Return the period between sunrise and sunset, in seconds.
At high latitudes around the time of the solstices, this could be zero, or all day.
*/
unsigned long daylightseconds(unsigned long dt){
float l, d, e;
long n;
d = -SolarDeclination(dt); // will be positive in Northern winter
l = latitude / _sec_rad; // latitude in radians
e += 60.0 * l * tan(l + d); // latitudinal error
d = tan(l) * tan(d); //
if(d>1.0) return 86400UL;
if(d < -1.0) return 0UL;
d = acos(d);
d /= _zenith;
n = 86400UL * d;
n += e;
return n;
}
/* ------------------------------------------------------------------------------------------------
Modify the passed time stamp to the time of sunrise (or sunset if 'set' is non-zero).
Returns 0 to signal 'normal' completion. If the position is in a polar circle, 1 will be
returned if the sun is above the horizon all day, and -1 if the sun is below the horizon
all day.
*/
char SunRiseSet(unsigned long * dt, char set){
unsigned long daylen;
daylen = daylightseconds(*dt);
if(daylen == 86400UL) return 1; // there is no 'night' today (midnight sun)
if(daylen == 0UL) return -1; // there is no 'day' today
*dt /= 86400UL;
*dt *= 86400UL;
*dt += 43200UL; // set the time stamp to 12:00:00 GMT
*dt -= daylen / 2; // sunrise at the prime meridian
if(set) *dt += daylen; // sunset at the prime meridian
*dt -= equation_of_time(*dt);
*dt -= longitude / 15.0; // rotate to our own meridian
return 0;
}
// 'short' forms of SunRiseSet
char SunRise(unsigned long* when){
return SunRiseSet(when, 0);
}
char SunSet(unsigned long* when){
return SunRiseSet(when, 1);
}
/* ------------------------------------------------------------------------------------------------
Mean Sidereal Time.
Accurate to within 1 sidereal second.
*/
unsigned long GMSidereal(unsigned long gmt){
unsigned long long sidereal;
sidereal = gmt * 10027379094LL; // multiply by 10 Billion times the ratio
sidereal /= 10000000000LL; // and divide by 10 billion
sidereal += 23992LL; // add sidereal time at the epoch
return sidereal;
}
unsigned long LMSidereal(unsigned long gmt){
return GMSidereal(gmt ) + longitude / 15L;
}
/* ------------------------------------------------------------------------------------------------
An approximation of the moons phase.
Magnitude of the result approximates the percentage of the moons illuminated surface.
Sign of the result indicates a Waxing (+) or Waning (-) moon.
It uses the mean lunar cycle, which may differ from the actual by many hours.
As such, it may vary from the correct value by 20%.
*/
char MoonPhase(unsigned long d){
long r;
// the first full moon of the epoch was Jan 21 at 04:40
// but the 'mean' full moon was 4 hours later
d -= 1759365UL; // subtract time of first full moon
d %= 2551443L; // mod by the mean lunar cycle
r = d - 1275721L;
r /= 12757L;
return r;
}
// Season of the year
// 0 = winter, 1 = spring, 2 = summer, 3 = fall
unsigned char season(unsigned long dt){
dt += 838800UL;// refer to prior winter solstice
dt %= _tropical_year;
dt /= 7889400UL; // 91.3125 days
if(latitude<0) dt += 2UL;
return dt % 4;
}
Code: Select all
/*
SWFLTEK_Ephemera copyright 2010 by Michael Rice
Swfltek Ephemera is a set of astronomical functions often of interest.
Though written for avr-gcc, it should convert easily to other compilers.
Most functions require an unsigned long time stamp as an argument. This 32 bit value represents
the number of seconds elapsed since midnight Jan 1 2000 GMT.
*/
#ifndef SWFLTEK_EPHEMERA_h
#define SWFLTEK_EPHEMERA_h
// geographic coordinates, in seconds of arc North and East
long latitude, longitude;
// conveniences to convert two typical representations into seconds of arc
long ddToSeconds(float);
long dmsToSeconds(int, unsigned char, unsigned char);
// return equation of time in seconds
int equation_of_time(unsigned long);
// adjust stamp to Solar noon
void SolarNoon(unsigned long * );
// return solar declination in radians
double SolarDeclination(unsigned long dt);
// return seconds between sunrise and sunset
unsigned long daylightseconds(unsigned long);
// compute time of sun rise or sunset
char SunRiseSet(unsigned long*, char);
// shorthand form
char SunRise(unsigned long*);
char SunSet(unsigned long*);
// convert from GMT to Mean Sidereal Time
unsigned long GMSidereal(unsigned long);
unsigned long LMSidereal(unsigned long);
// approximate phase of the moon
char MoonPhase(unsigned long);
// season
unsigned char season(unsigned long);
#endif // SWFLTEK_EPHEMERA_h
The 'set' parameter is required for the 'worker' SunRiseSet() function. The 'helper' functions SunRise() and SunSet will pass the correct value to the worker function.
Since your code is using the helper functions, you need not worry about it.
I have noted several things I must document more clearly however.
The Ephemera Library is straight 'C', not a C++ class. All the functions are global, so there is no need to instantiate any SWFLTEK_EPHEMERA object. Just #include the library and call the functions.
The timestamp is an unsigned 32 bit value, using an 'int' will not work.
Coordinates are in seconds, requiring a 'long' to represent them. For your application you will want to set the latitude to the same latitude as the reef you are simulating. You can also set the longitude same as the reef, but you'll probably want to set it to your local longitude to minimize disturbances from natural light and human activities.
The Ephemera functions use Greenwich Mean Time. Passing the local time will return a value which may have twice the 'daily' error, which probably isn't significant in your application but you might want to consider it.
The sun functions arguments are passed by reference. The returned value indicates if the call was successful. Since the reef in question is (probably ?) not inside either polar circle, you can ignore the returned value in your application.
Here is a framework for what you want to do...
============================================================================
#include SWFLTEK_EPHEMERA_h
unsigned long sunrise, sunset; // 32 bit unsigned integer representing seconds since year 2000
char moonphase; // 8 bit signed value, -100 ... 0 ... 100 %, sign indicates waning or waxing
void setup(){
// we are mimicking a coral reef at the Bahama Islands but 'rotated' to our longitude
latitude = dmsToSeconds(22,0,0);
longitude = dmsToSeconds(20,0,0);
...
...
}
void CalSun{
unsigned long SecInput;
// compute the current timestamp as SecInput
...
...
// copy to sunrise and sunset
sunrise = sunset = SecInput;
// compute actual times
// we know we are not in the Arctic or Antarctic, so we can ignore the returned value
SunRise( & sunrise);
SunSet( & sunset);
// compute moon phase
moonphase = MoonPhase(SecInput);
}
============================================================================
Hope this helps.
if (statement) seems a common error! I read it and read it and never noticed that I started using it correctly then fell back in to VB! 
) that the .h file I download from Michaels site which is this