NOAA Astronomical Calculations

Squirrel Code
1

Hello, I am attempting to perform the sunrise and sunset calculations found here (http://www.esrl.noaa.gov/gmd/grad/solcalc/solareqns.PDF) [Thanks Peter!] using just the agent. I have most of it going, however I am not sure if I am understanding the final arccos step correctly. I am getting NaN for the acos step…

`local longitude = -83;
local latitude = 30
local timezone = -5;
local datedata = date();
local dayofyear = datedata.yday;
local hour = datedata.hour;
local min = datedata.min;
local sec = datedata.sec;
server.log("Day of year = " + dayofyear);
server.log("Hour of day = " + hour);
server.log("Minute of day = " + min);
server.log("Second of day = " + sec);

local fryear = ((2PI)/365)(dayofyear-1+((hour-12)/24));
server.log("Fractional year = " + fryear);

local eqtime = 229.18*(0.000075+(0.001868math.cos(fryear))-(0.032077math.sin(fryear))-(0.014615math.cos(2fryear))-(0.040849math.sin(2fryear)));
server.log("Equation time = " + eqtime);

local declination = 0.006918-(0.399912math.cos(fryear))+(0.070257math.sin(fryear))-(0.006758math.cos(2fryear))+(0.000907math.sin(2fryear))-(0.002697math.cos(3fryear))+(0.00148math.sin(3fryear));
server.log("Declination = " + declination);

local timeoffset = eqtime-(4longitude)+(60timezone);
server.log("Time ofset = " + timeoffset);

local tst = hour*60+min+(sec/60)+timeoffset;
server.log("tst = " + tst);

local ha = (tst/4)-(180);
server.log("ha = " + ha);

local zenithang = math.sin(latitude)*math.sin(declination)+math.cos(latitude)*math.cos(declination)*math.cos(ha)
server.log("Zenith angle = " + zenithang);

local sunrise = 720+(4*(longitude-ha)-eqtime);
server.log("Sunrise = " + sunrise);

local ha3 = math.cos((90.8333)/(math.cos(latitude)*math.cos(declination))-math.tan(latitude)*math.tan(declination));

local newha = math.acos(math.cos(90.8333)/(math.cos(latitude)*math.cos(declination))-math.tan(latitude)*math.tan(declination));
server.log("Newha = " + newha);

local sunrise2 = 720+(4*(longitude-ha3)-eqtime);
server.log(sunrise2/60);

local sunset2 = 720-(4*(longitude-ha3)-eqtime);
server.log(sunset2/60);

`

I was wondering if someone could lend me a hand to understand how to correctly get the sunrise and sunset times from the solar equations provided by NOAA. I think I might be “skipping” steps, or doing something out of order (even though I am following the equations from top to bottom. I feel like i’m missing something… Any help is appreciated. Thanks!

2

Hello -

If you are doing this at the Agent level, I’m curious why not query a web API to do give the sunrise/sunset time?

When you run this is sunrise working right? Just trying to confirm that you know the problem is in the last portion of your code.

3

Yes, I am using wunderground already. But need this for worst case, like no Wi-Fi for say 2 months or something. Or if there is no memory chips to store tables, I can hardcode coordinates in and not have to worry.

I’m not too sure what to make of the results really… I don’t know if I’m interpreting the equations correctly…

4

To clarify, I’m only running the code in agent for testing purposes. This will be on the device in the end…

5

FWIW, I’d run it on the agent and send the result to the device. Agent has more processing power than the device

6

I agree, however have to be prepared for worst case for customers… The plan is to have the device be the “master” and agent is a “copy” of whatever runtime variables are created for web services… Worst case is cold boot with no Wi-Fi or memory relying solely on hardcoded install specifics…

7

No problem running it on the device (it’s not that much maths!).

I tried to get the NOAA one working but failed, so I found another algorithm and got that working in the end too. The numbers seem to match within a minute or two for me (at least in London and Palo Alto, where I checked them).

// Sunrise/sunset calculation in squirrel, hugo@electricimp.com
// Based on http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

// dayofyear is 1 based (Jan 1st = 1)
// longitude is positive for east, negative for west
// latitude is positive for north, negative for south
// localoffset is vs UTC (eg BST=1, PDT=-7)

function riseorset(rise, dayofyear, longitude, latitude, localoffset) {
    // "Official" zenith is 90.833 degrees
    local zenith = 90.833;
    
    // convert the longitude to hour value and calculate an approximate time
    local lngHour = longitude / 15.0;
    
    local t;
    if (rise) {
        // rising time
        t = dayofyear + ((6 - lngHour) / 24.0);
    } else {
        // setting time
        t = dayofyear + ((18 - lngHour) / 24.0);
    }
    
    // calculate the Sun's mean anomaly
    local M = ((0.9856 * t) - 3.289);
    
    // calculate the Sun's true longitude, and adjust to 0-360
    local L = M + (1.916 * math.sin((PI/180)*M)) + (0.020 * math.sin(2 * (PI/180)*M)) + 282.634;
    if (L<0) L+=360; else if (L>360) L-=360;
    
    // calculate the Sun's right ascension
    local RA = (180/PI)*math.atan(0.91764 * math.tan((PI/180)*L));
    if (RA<0) RA+=360; else if (RA>360) RA-=360;
    
    // right ascension value needs to be in the same quadrant as L
    local Lquadrant  = (math.floor( L/90)) * 90;
    local RAquadrant = (math.floor(RA/90)) * 90;
    RA = RA + (Lquadrant - RAquadrant);
    
    // right ascension value needs to be converted into hours
    RA = RA / 15;
    
    // calculate the Sun's declination
    local sinDec = 0.39782 * math.sin((PI/180)*L);
    local cosDec = math.cos(math.asin(sinDec));
    
    // calculate the Sun's local hour angle
    local cosH = (math.cos((PI/180)*zenith) - (sinDec * math.sin((PI/180)*latitude))) / (cosDec * math.cos((PI/180)*latitude))
    
    if (cosH > 1) {
    	server.log("The sun never rises on this location (on the specified date)");
    	return -1;
    } else if (cosH < -1) {
        server.log("The sun never sets on this location (on the specified date)");
        return -1;
    }
    
    // finish calculating H and convert into hours
    local H = (180/PI)*math.acos(cosH);
    if (rise) H = 360-H;
    H = H / 15.0;
    
    // calculate local mean time of rising/setting
    local T = H + RA - (0.06571 * t) - 6.622;
    local UT = T - lngHour;
    
    // Add local offset and normalize
    UT += localoffset;
    if (UT < 0) UT += 24;
    if (UT > 24) UT -= 24;
    
    return UT;
}

// Palo Alto, CA
local longitude = -122.1381;
local latitude = 37.4292;
local timezone = -7; // PDT
local dayofyear = date().yday+1;

local rise = riseorset(true, dayofyear, longitude, latitude, timezone);
server.log("rise "+format("%02d:%02d", math.abs(rise), (rise-math.abs(rise))*60));
    
local set = riseorset(false, dayofyear,  longitude, latitude, timezone);
server.log("set "+format("%02d:%02d", math.abs(set), (set-math.abs(set))*60));

…obviously if you want it in UTC (which you may if you are using this to do things based on the imp’s UTC clock) then pass the timezone as zero.

8

Perfect, Thanks!

9

If anyone is still interested, I figured out most of the issues in the NOAA squirrel code above and have it working. The main problem is that the trig functions expect radians and not degrees. So math.sin(latitude), should be math.sin(latitude * (PI/180)). Whenever latitude is used, substitute latitude * (PI/180) and 90.8333 should be 90.8333*(PI/180). (Make the same change for calculations using longitude.) The arc cosine also outputs radians and a following equation expects the results in degrees so math.acos(…) should be (180/PI)*math.acos(…).

The following equations should be substituted for the existing sunrise2 and sunset2, respectively. Sunset2 and Sunrise2 are in minutes so will have to be converted to hours and minutes.

local sunrise2 = 720+(4*(longitude+newha)-eqtime)

local sunset2 = 720+(4*(longitude-newha)-eqtime);

I found that I still have to add the offset from UTC time to get the correct values. As an example, add 5 hours for Eastern Standard Time where the offset from UTC Time is -5 hours. Don’t have a reason for this change yet, but it might be due how the timezone is represented in: local timeoffset = eqtime-(4longitude)+(60timezone);