EQUATIO: Sidereal & Solar Clock

Meet EQUATIO: Sidereal & Solar Clock, my latest investigation of the measurement and perceptions of time. As so much of our life revolves around the 24-hour per day clock, it seems only obvious to question whether such precision and regularity is really natural? Or is our common measurement of time a man-made contrivance to insulate us from the chaos of the natural world?

Peering behind the curtain, we find strange things afoot!

So, I give you (roll of drums)… EQUATIO: Sidereal & Solar Clock… a bold attempt to re-appropriate “real time” and expose those “meanies” who control us through the imposition of order, symmetry and regularity (although at my age, the latter can be pretty important!).

EQUATIO Graphic

EQUATIO Graphic: A revision to EQUATIO to display Local Mean Solar Time, Local Sidereal Time and Apparent Solar Time on a 192×64 pixel graphic display module. The difference in time between the Local Mean Solar Time and the Apparent Solar Time is the Equation of Time for 29 June 2015.

Sk

EQUATIO: Sidereal & Solar Clock showing current Eastern Daylight Time and today’s date (based on Mean Solar Time), Local Sidereal Time (LST) (based on time relative to the stars) and Apparent Solar Time (SOL) and this day’s Equation of Time (EOT) offset based on the actual observation of the sun.

Equatio: Sidereal & Solar Clock

EQUATIO: Sidereal & Solar Clock showing that the Apparent Solar Time currently lags Mean Solar Time (Eastern Daylight Time)  by 0.6 minutes (or 36 seconds).

“Equatio” was a word widely used in early astronomy as the difference between observed and averaged values. It seems quite apropos for this project.

Read on for a steamy exposé .  I NAME NAMES!

Now, it is commonly held that our 24 hour clock is based on the actual time it takes for Earth to rotate once about its axis such that the sun appears in the same position in the sky; commonly  referred to as a solar day.  OK, so far so good. However, the length of a solar day varies through the year as a result of two principal affects – the elliptical nature of our orbit around the sun, and that the Earth’s axis is not perpendicular to the plane of its orbit around the sun (referred to as the Obliguity of the Ecliptic).  The accumulated affect is that the length of solar days (an apparent or true solar day) varies according to the season by up to 16 minutes ! compared to an “average” solar day. In other words, our regularly ticking clocks run faster or slower by up to a quarter of an hour compared to our observations of the sun. But I don’t see anyone adjust their watch each day to accommodate these seasonal shifts. So, what gives?

Apart from British Rail’s timetables and my partner Lise’s punctuality, many of us would be flummoxed by dealing with a seasonal variation of quarter of an hour. So bright minds invented the notion of an averaged time, called “mean solar time” such that a “mean solar day” consists of 24 neatly and identically divided hours.

Ahh. Order is returned to the universe and our passage through it!

However, back to the not-so-pretty reality. The combination of the two effects (elliptical orbit and Obliquity of the Ecliptic) on time variance has been graphed in what is termed the “Equation of Time“. This, according to Wikipedia, is defined as “the discrepancy between … apparent solar time, which directly tracks the motion of the sun, and mean solar time” and was used historically to set mechanical clocks based on the time seen on a sundial.

Equation of Time describing the seasonal variation of time between apparent and mean solar days

Equation of Time:  describing seasonal variations of the length of apparent solar days compared to a mean solar day. Above the line means that my sundial is “faster” than your wristwatch. Image shamelessly stolen from Wikipedia.

The elliptical nature of our orbit around the sun means that the Earth travels faster (traveling up to 21,726 km more each day) when it is nearest the Sun (perihelion) and slower when it is farthest from the Sun (aphelion) and its effect manifests itself in the Equation of Time as a sine wave with a period of a year.

As the Earth’s rotational axis with respect to the sun is tilted to about 23.5°, the angle of the sun’s rays (declination) varies through the year. This effect  – the Obliguity of the Ecliptic (OoE) – creates changes in the position of the sun such that it is at a maximum at the solstices (when the movement of the Sun is parallel to the equator) and at a minimum at the equinoxes. The OoE moves from +23.5° at the summer solstice to -23.5° at the winter solstice, and is 0° during the equinoxes.  In the Equation of Time this time variance is a sine wave with a period of a half of a year.

In the graph of the Equation of Time (above) positive values indicate that the time on a sundial (apparent solar time) is faster than a regular clock (mean solar time). Similarly, negative values indicate that apparent solar time is slower than mean solar time.

EQUATIO shows both the Equation of Time variance (in minutes) for the day of the year, and the calculated Apparent Solar Time (SOL).

What about the “sidereal” element to EQUATIO?

Well, we are not only spinning on our own axis, but orbiting around the sun as well. So, in a single solar day, we not only revolve our own 360º but we have also moved nearly 1º degree (we complete 360 degrees in 365.2422 days, or 0.9856 degrees per day) in our orbit around the sun. In other words, we have to travel all the way around and then a little further (equivalent of just under 4 minutes) in order for us to observe the sun in the same position in the sky – an apparent solar day.

Now, to all intents and purposes, stars are so far away that the Earth’s orbit around the sun makes no difference to their apparent direction. So, if we wanted to calculate time relative to seeing the same stars in the same position of the sky, we need to remove this additional “orbital” time.

Enter “Sidereal time“, a timekeeper based on the Earth’s rotation measured relative to the fixed stars. such that from a given location, a star found at one location in the sky will be found at the same location on another night at the same sidereal time. Therefore, a “sidereal day” is the time taken for the Earth to make one full revolution (or the rotational period) and is 23 hours, 56 minutes, 4.0916 seconds, or 23.9344699 hours or 0.99726958 mean solar days.  Now, as the Earth orbits the sun once a year, the sidereal time at any given place and time will gain about four minutes every day relative to mean solar time, until, after a year has passed, there is one additional sidereal “day” compared to the number of elapsed solar days. Neat huh?

EQUATIO shows the Local Sidereal Time (LST) which is based on calculating the Greenwich Apparent Sidereal Time, which according to Wikipedia is “… is the hour angle of the vernal equinox at the prime meridian at Greenwich, England”. The GAST is then adjusted for the longitude of the user for the calculation of Local Apparent Sidereal Time.

EQUATIO is currently built on an Arduino Nano. There is also a real-time clock and a 4×20 character LCD, both of which are connected using an I2C interface. To model the daily time variance to calculate real solar time, I average the results from three different algorithms that approximate the Equation of Time. I used EXCEL to model the three algorithms as shown here. The average of the algorithms closely tracks more accurate EoT online calculators such as this and correctly calculates the difference between apparent solar time and mean solar time to be about -14 minutes in February and around +16 minutes in October.

Modeling of Equation of Time Algorithms

Modeling of Equation of Time Algorithms

The algorithm for calculating Sidereal Time was based on a paper entitled “ASTR 310 – Observational Astronomy: Formula for Greenwich Sidereal Time (GST)“. I modeled the algorithm and created a simple way to extrapolate the calculation of the principal variant (G) to account for any year beyond 1 January 2000. The modeling for the extrapolation is seen here:

Greenwich

Greenwich Sidereal Time algorithm: extrapolation of one of the variables.

The algorithm also includes factors for the day of the year and the current time (in UTC). Once GST is calculated, I then add a time adjustment for my current longitude equal to (my_longitude/360)*length_of_sidereal_day. Similar to EoT, this algorithm tracks very closely to calculations of sidereal time found in online calculators, such as this and this.

For those of you star gazers, here’s the Arduino code…

[codesyntax lang=”php” title=”EQUATIO: Sidereal & Solar Clock: v1.1 : Arduino Code” blockstate=”collapsed”]

//**************************************************************************************************
//                               EQUATIO: SIDEREAL & SOLAR CLOCK
//                                   Adrian Jones, June 2015
//
//**************************************************************************************************

// Build 1
//   r1 150601  initial build
//********************************************************************************************
#define build 1
#define revision 1
//********************************************************************************************

#include <avr/pgmspace.h>             // for memory storage in program space
#include <math.h>                     // math for solar & moon position calculations 

// sidereal calculation constants
#define dc 0.0657098244
#define tc 1.00273791
#define gc 6.648605
#define g2000 6.5988098
#define lc 0.0497958000000001
#define nc -0.0159140999999998
#define fudge -0.013922               // fudge factor (unnecessary?
#define LONGITUDE -75.64316            
#define LATITUDE  45.66167
#define siderealday 23.9344699        // length of sidereal day (23:56:04)
double GST,LST,utc;                   // greenwich & local apparent sidereal time
int dh,dm,ds;                         // local sidereal time
int sh,sm,ss;                         // solar time

// equation of time
double tv,tv1,tv2,tv3;                // total time variation (mins)
int tvm,tvs,nd;                       // EoT time mins and seconds, days from Jan 1.

// RTC
#include <Wire.h>
#include "RTClib.h"              
RTC_DS1307 RTC;                       // create instance of RTC
int hr,mn,se,osec=-1,dy,mo,yr,dw;     // time variables
int TZ;                               // time zone offset (hours)

// LCD
#include <LCD.h>                      // Standard lcd library
#include <LiquidCrystal_I2C.h>        // library for I2C interface

#define I2C_ADDR  0x27                // address found from I2C scanner
#define RS_pin    0                   // pin configuration for LCM1602 interface module
#define RW_pin    1
#define EN_pin    2
#define BACK_pin  3
#define D4_pin    4
#define D5_pin    5
#define D6_pin    6
#define D7_pin    7

LiquidCrystal_I2C lcd(I2C_ADDR, EN_pin, RW_pin, RS_pin, D4_pin, D5_pin, D6_pin, D7_pin, BACK_pin, POSITIVE);
//Pins for the LCD are SCL A5, SDA A4

// eeprom storage
#include <EEPROM.h>               // for alarm and mode data storage
int dst;                          // dst setting                      
int *set[] = { &dst };            // values to set/recover from EEPROM


//***************************
#define SERIAL_BAUDRATE 57600
#define serEn  true              // serial output
#define sp Serial.print   
#define spln Serial.println

#define lcdEn  true              // lcd output
#define lp lcd.print
#define ls lcd.setCursor
//***************************


//*****************************************************************************************//
//                                      Initial Setup
//*****************************************************************************************//
void setup() {
  
#if(serEn) 
  Serial.begin(SERIAL_BAUDRATE);
  spln(F("EQUATIO: SIDEREAL & SOLAR CLOCK by Adrian Jones"));
  sp(F("Build ")); sp(build); sp(F(".")); spln(revision);
  sp(F("Free RAM: "));  sp(freeRam()); spln(F("B"));
  spln(F("**************************************"));
#endif 
  
  Wire.begin();
  RTC.begin();                                         // start RTC
  getTime();                                           // get current time

  if (! RTC.isrunning() || (String(__DATE__).lastIndexOf(String(yr)) == -1)) {
    RTC.adjust(DateTime(__DATE__, __TIME__));          // adjust if year is off or RTC not running
    #if(serEn) 
      sp(F("RTC reset to ")); sp(__DATE__); sp(F(" ")); spln(__TIME__);
    #endif
  } 
  // RTC.adjust(DateTime(__DATE__, __TIME__));         // reset to system time
  // RTC.adjust(DateTime(2015, 1, 18, 1, 59, 50 ));    // test time 

#if(lcdEn)
  lcd.begin(20, 4);
  lcd.clear();
  ls(0, 0); lp(F("SIDEREAL SOLAR CLOCK")); 
  ls(5, 1); lp(F("b")); lp(build); lp(F(".")); lp(revision); lp(F(" "));  lp(freeRam()); lp(F("B"));
  delay(1000);
#endif

  restoreSettings();                                 // restore EEPROM settings

  getTime();
  TZ = (IsDST)?-4:-5;                                // determine time offset (hours) from GMT

}
 
//*****************************************************************************************//
//                                      MAIN LOOP
//*****************************************************************************************//
void loop() {
  getTime();                                        // get current time
  doDSTAdjust();                                    // adjust for DST if required
  
  if(se != osec) {                                  // do every second
    doSiderealCalc();                               // do sidereal calculations
    doEoTCalc();                                    // do EoT calculations
    doSerialTime();                                 // write out time to serial port
    doLCDTime();                                    // write out time to LCD
    osec = se;                                      // old sec = curent
  }
  delay(10);
}



// ********************************************************************************** //
//                                 OUTPUT ROUTINES
// ********************************************************************************** //

void doLCDTime() {
#if(lcdEn) 
    ls(0, 1); lp(IsDST()?"EDT: ":"EST: ");  
    if(hr<10) lp(F("0")); lp(hr); lp(F(":"));
    if(mn<10) lp(F("0")); lp(mn); lp(F(":"));
    if(se<10) lp(F("0")); lp(se);

    lp(F(" ")); lp(yr-2000); lp(F(""));
    if(mo<10) lp(F("0")); lp(mo); lp(F(""));
    if(dy<10) lp(F("0")); lp(dy); 

    ls(0,2);  lp(F("LST: "));     
    if(dh<10) lp(F("0")); lp(dh); lp(F(":"));
    if(dm<10) lp(F("0")); lp(dm); lp(F(":"));
    if(ds<10) lp(F("0")); lp(ds); 

    ls(0,3);  lp(F("SOL: "));     
    if(sh<10) lp(F("0")); lp(sh); lp(F(":"));
    if(sm<10) lp(F("0")); lp(sm); lp(F(":"));
    if(ss<10) lp(F("0")); lp(ss); 
    
    lp(F(" "));
    lp(tv,1);  lp(F("m "));
#endif
}

void doSerialTime() {
#if(serEn) 
    sp(F("MEAN SOLAR (CLOCK) TIME: ")); sp(yr-2000); sp(F("/"));
    if(mo<10) sp(F("0")); sp(mo); sp(F("/"));
    if(dy<10) sp(F("0")); sp(dy); sp(F(" "));
    if(hr<10) sp(F("0")); sp(hr); sp(F(":"));
    if(mn<10) sp(F("0")); sp(mn); sp(F(":"));
    if(se<10) sp(F("0")); sp(se);
    
    sp(F("\tLOCAL SIDEREAL TIME: "));     
    if(dh<10) sp(F("0")); sp(dh); sp(F(":"));
    if(dm<10) sp(F("0")); sp(dm); sp(F(":"));
    if(ds<10) sp(F("0")); sp(ds);
    
    sp(F("\tEQUATION OF TIME: "));     
    sp(tvm);  sp(F("m ")); sp(tvs); sp(F("s ")); 
    sp((tv>0.0)?"s<c":"s<c");
    
    sp(F("\tAPPARENT SOLAR (SUNDIAL) TIME: "));     
    if(sh<10) sp(F("0")); sp(sh); sp(F(":"));
    if(sm<10) sp(F("0")); sp(sm); sp(F(":"));
    if(ss<10) sp(F("0")); sp(ss);

    spln(); 
#endif
}

// ********************************************************************************** //
//                             EQUATION OF TIME CALCULATIONS
// ********************************************************************************** //

void doEoTCalc() {
  nd = doNumDays(yr, mo, dy);                          // days from Jan 1 (inc. leap year)
  
// 1. based on "The Equation of Time" by Murray Bourne, 26 Aug 201
//    see http://www.intmath.com/blog/mathematics/the-equation-of-time-5039
  float dd = (2*PI*(float)nd/365.0);
  tv1 = -7.655 * sin(dd) + 9.873 * sin(2*dd + 3.588); // Effect of Orbit Eccentricity & Effect of Obliquit

// 2. based on http://www.susdesign.com/popups/sunangle/eot.php
  float B = 2*PI*((float)nd - 79.0)/365.0;               // B in radians (note should be -81)
  tv2 = 9.87*sin(2*B) - 7.53*cos(B) - 1.5*sin(B);

// 3. based on http://naturalfrequency.com/Tregenza_Sharples/Daylight_Algorithms/algorithm_1_12.htm
  tv3 = 60.0 * (0.170*sin(4*PI*(nd-80)/373.0) - 0.129*sin(2*PI*(nd-8)/355.0) );
  
  tv = (tv1+tv2+tv3)/3.0;
  tvm = int( abs(tv));                                  // whole mins
  tvs = int((abs(tv) - float(tvm))*60.0 );              // whole secs

  DateTime now = RTC.now();
  DateTime solTime (now.unixtime() + int(tv*60));       // solat time (if EoT is -ve, clock is ahead of solar time)
  sh = solTime.hour();
  sm = solTime.minute();
  ss = solTime.second();
}


// ********************************************************************************** //
//                       GREENWICH & LOCAL SIDEREAL TIME CALCULATIONS
// ********************************************************************************** //
// based on "ASTR 310 - Observational Astronomy: Formula for Greenwich Sidereal Time (GST)" 
// see http://www.astro.umd.edu/~jph/GST_eqn.pdf formulas

void doSiderealCalc() {
// UTC calculation
   DateTime now = RTC.now();                                  // current time
   DateTime utcTime (now.unixtime() - TZ*3600);               // adjust to GMT
   utc = (float)utcTime.hour() + (float)utcTime.minute()/60.0 + (float)utcTime.second()/3600.0; // decimal form
   
// calculate G (based on extrapolation)
   int g = (utcTime.year() - 2000);
   int leap = int((g+1.0)/4.0);                              // number of leap years since 2000
   int nleap = g-leap;                                       // number of non-leap years since 2000   
   double G = g2000 + (float)leap*lc + (float)nleap*nc;      // number of days
   
// calculate nd
  int nd = doNumDays(utcTime.year(), utcTime.month(), utcTime.day());
  
// calculate GST and Local Sidereal Time (LST)
  GST = G + (dc*nd) + (tc*utc) + fudge;                      // Grenwich Sidereal Tim
  LST = GST + 24.0 + (float)(LONGITUDE/360*siderealday);     // adjust for longitude (longitude portion of siderail day
  while(LST>24.0) {  LST -= 24.0; }                          // adjust to bring into 0-24 hours
  dh = int( LST );                                           // translate into hours, ...
  dm = int( (LST - (float)dh)*60.0 );                        //... mins and ...
  ds = int( (LST - (float)dh - (float)dm/60.0)*3600.0 );     //... seconds
 }

// number of days of month (m) and date (d) since beginning of year (y)
int doNumDays(int y, int m, int d) {
  int v=0;
  byte leapyear = ((y-2000)%4 == 0)? 1 : 0;
  switch(m) {
      case 12:  v += 30;        // Dec
      case 11:  v += 31;        // Nov
      case 10:  v += 30;        // Oct
      case 9:   v += 31;        // Sep
      case 8:   v += 31;        // Aug
      case 7:   v += 30;        // Jul
      case 6:   v += 31;        // Jun
      case 5:   v += 30;        // May
      case 4:   v += 31;        // Apr
      case 3:   v += 28 + leapyear;   // May (if year is leapyear, add extra day after February)
      case 2:   v += 31; break; // Feb
  }
  return v+d;                   // days from Jan 1 of given year
}
 

/* ********************************************************************************** //
//                                    TIME ROUTINES
// ********************************************************************************** //
   In most of Canada Daylight Saving Time begins at 2:00 a.m. local time on the second Sunday in March. 
   On the first Sunday in November areas on DST return to Standard Time at 2:00 a.m. local time. 
*/

// IsDST returns true if DST, false otherwise
boolean IsDST() {
  if (mo < 3 || mo > 11) { return false; }                // January, February, and December are out.
  if (mo > 3 && mo < 11) { return true;  }                // April to October are in
  int previousSunday = dy - dw;                               
  if (mo == 3) { return previousSunday >= 8; }            // In March, we are DST if our previous Sunday was on or after the 8th.
  return previousSunday <= 0;                             // In November we must be before the first Sunday to be DST. That means the previous Sunday must be before the 1st.
}

// getTime obtains current time from RTC
void getTime() {
  DateTime now = RTC.now();
  hr = now.hour(); if(hr==0) hr=24;                      // adjust to 1-24
  mn = now.minute();
  se = now.second();
  yr = now.year();
  mo = now.month();
  dy = now.day();
  dw = now.dayOfWeek();
}

// doDSTAdjust increments (or decrements) by one hour when entering (or leaving) DST
void doDSTAdjust() {
  if(dst == IsDST()) return;                    // if prior setting is same as DST setting, do nothing
  if(hr != 2) return;                           // do nothing until 2pm

  DateTime now = RTC.now();                     // get time
  if(IsDST() && !dst) {
    DateTime newTime (now.unixtime() + 3600);   // add one hour
    RTC.adjust(newTime);                        
  }
  if(!IsDST() && dst) {
    DateTime newTime (now.unixtime() - 3600);   // subtract one hour
    RTC.adjust(newTime);
  }
  dst = IsDST(); 
  saveSettings();                // save change to DST
}


// ********************************************************************************** //
//                                      OPERATION ROUTINES
// ********************************************************************************** //
// FREERAM: Returns the number of bytes currently free in RAM  
int freeRam(void) {
  extern int  __bss_end, *__brkval; 
  int free_memory; 
  if((int)__brkval == 0) {
    free_memory = ((int)&free_memory) - ((int)&__bss_end); 
  } 
  else {
    free_memory = ((int)&free_memory) - ((int)__brkval); 
  }
  return free_memory; 
}


// ********************************************************************************** //
//                                       EEPROM ROUTINES
// ********************************************************************************** //
// EEPROM save, restore and set/save defaults
void saveSettings()    { for(int addr = 0; addr < sizeof(set)/2; addr++) { EEPROM.write(addr, *set[addr]); } }    // save
void restoreSettings() { for(int addr = 0; addr < sizeof(set)/2; addr++) { *set[addr] = EEPROM.read(addr); } }   // recover EEPROM value

[/codesyntax]

9 thoughts on “EQUATIO: Sidereal & Solar Clock

  1. Alexander Lang

    I think your implementation is excellent – I’m trying to roll my own version of this. Well done. I’m trying to use your code to calculate my local sidereal time for the UK. As far as I can tell I need to update the latitude and longitude and everything should work! – Thanks

    Reply
      1. Adrian Post author

        One thing to note is that the software is on the hairy edge of both program and SRAM memory limits of the Nano (you would have lots more on a Mega). I needed a simple and small library so I used what i was familiar with.
        What do you think you could save using the time.h library?
        In addition, change the line in the doDSTAdjust() routine from:
        “if(hr != 2) return;”
        to
        “if(hr <= 2) return;" This will appropriately handle DST adjustment. Just email if you have any questions or issues.

        Reply
        1. Alexander Lang

          I don’t know if there is an improvement on using time.h over RTClib.h I was just curious as to why you were using RTClib.h as it’s not a standard library.

          Reply
    1. Adrian Post author

      Thanks for your encouraging words, Alexander.
      Yes… Just change these two variables and the timezone. Note that my timezone is EST so this has to be factored into the software (for you, set TZ=0).

      Reply
      1. Alexander Lang

        It all works and I’m using an arduino Uno. I really wanted to understand the code a little better though. How did you figure out how to calculate GMST and LMST from the formulae provided. Wikipedia has been less than helpful as I cannot seem to get their formula to concur with the online calculators – I must be doing something wrong and I expect it’s to do with the floor functions (I am terrible at Mathematics)….

        Reply

Come on... leave a comment...