Thursday, February 14, 2019

Solar Position Functions, Google Sheets Edition

Using the sun as a navigation aid is both an obvious idea and one fraught with difficulty. Making sense of the sun's position throughout the day involves understanding your latitude, longitude, timezone offset and a bunch of non-trivial (to me!) math. Fortunately, I've been able to untangle this math. The next step was to make leveraging these calclations a breeze. For this, I turned to Google Sheets.

I popped into Google Sheets, opened up the Script Editor and started translating scheme code to JavaScript. When I was done, I had the following custom Sheets functions implemented:

# Sun's Position
elevation(lat,lng,ts,tzOffset=null) - sun's position in degrees above the horizon
azimuth(lat,lng,ts,tzOffset=null) - heading of the sun's position in degrees 

# Notable times when the sun is 0° above the horizon
sunrise(lat,lng,ts,tzOffset=null)
sunset(lat,lng,ts,tzOffset=null)

You can find the code that implements these functions here. The ts (timestamp) parameter passed in above uses the internal Google Sheets timestamp format. This allows the above code to gracefully integrate with other Sheets functions. The tzOffset parameter is optional; if unset it falls back to the user's current timezone preference. This auto behavior is handy for seamlessly dealing with daylight saving changes throughout the year.

Consider these examples:

These demonstrate printing a full day's worth of elevations and azimuths. Using these charts I can trivially derive my current direction by using the sun and my watch, or derive the time by using the sun and a compass.

Stay tuned for more ideas as to how these functions can be used.

Update: apparently my attempt to link to the Google Sheet's App Script was a bust. Here's the code that powers the spreadsheet:

/*
* See:
* http://www.blogbyben.com/2019/01/go-towards-light-direction-and-time.html
* https://github.com/benjisimon/code/blob/master/programming-praxis/sun-position.scm
* for an explanation as to what this is all about.
*/

eval(UrlFetchApp.fetch('https://cdnjs.cloudflare.com/ajax/libs/moment.js/2.24.0/moment.js').getContentText());

function solarPosition(lat, lng, ts, tzOffset) {
  var doy = moment(ts).dayOfYear();
  var hod = moment(ts).hour() + (moment(ts).minute() / 60);
  var tz  = tzOffset === undefined ? (moment(ts).utcOffset() / 60) : tzOffset;
  
  var usingDegrees = function(fn) {
    return function(dd) { return fn(dd * 0.0174533); };
  }
  var toDegrees = function(fn) {
    return function(v) { return fn(v) * 57.2958; };
  }
  
  var cos = usingDegrees(Math.cos);
  var sin = usingDegrees(Math.sin);
  var tan = usingDegrees(Math.tan);
  var acos = toDegrees(Math.acos);
  var asin = toDegrees(Math.asin);
  var atan = toDegrees(Math.atan);
  
  var lstm = tz * 15;
  
  var eot = (function() {
    var b = (360 / 365) * (doy - 81);
    return (9.87 * sin(2*b)) -
           (7.53 * cos(b)) -
           (1.5 * sin(b));
  })();
  
  var tcf = (4 * (lng - lstm)) + eot;
  var lst = hod + (tcf / 60);
  var hra = 15 * (lst - 12);
    
  var d = 23.45 * sin((360/365) * (doy - 81));
  
  var e = asin(((sin(d)) * sin(lat)) +
               ((cos(d)) * cos(lat) * cos(hra)));
  
  var a = (function() {
    var a = acos(((sin(d) * cos(lat)) - (cos(d) * sin(lat) * cos(hra))) /
                 (cos(e)));
    return hra < 0 ? a : (360 - a);
  })();
  
  var sr = 12 - (((1/15) * acos(-1 * tan(lat) * tan(d)))) - (tcf / 60);
  var ss = 12 + (((1/15) * acos(-1 * tan(lat) * tan(d)))) - (tcf / 60);
  
  return {
    elevation: e,
    azimuth: a,
    sunrise: sr / 24,
    sunset: ss / 24,
  };
}


/**
 * Elevation of the sun
 *
 * @param lat latitude of interest
 * @param lng longitude of interest
 * @param ts timestamp
 * @param tzOffset optional timezone offset from UTC. 
 * @return elevation in degrees
 * @customfunction
 */
function elevation(lat, lng, ts, tzOffset) {
  return solarPosition(lat, lng, ts, tzOffset).elevation;
}

/**
 * Azimuth (compass direction) of the sun
 *
 * @param lat latitude of interest
 * @param lng longitude of interest
 * @param ts timestamp
 * @param tzOffset optional timezone offset from UTC. 
 * @return azimuth in degrees
 * @customfunction
 */
function azimuth(lat, lng, ts, tzOffset) {
  return solarPosition(lat, lng, ts, tzOffset).azimuth;
}

/**
 * Sunrise
 *
 * @param lat latitude of interest
 * @param lng longitude of interest
 * @param ts timestamp
 * @param tzOffset optional timezone offset from UTC. 
 * @return sunrise as Google Sheets time
 * @customfunction
 */
function sunrise(lat, lng, ts, tzOffset) {
  return solarPosition(lat, lng, ts, tzOffset).sunrise;
}

/**
 * Sunset
 *
 * @param lat latitude of interest
 * @param lng longitude of interest
 * @param ts timestamp
 * @param tzOffset optional timezone offset from UTC. 
 * @return sunset as Google Sheets time
 * @customfunction
 */
function sunset(lat, lng, ts, tzOffset) {
  return solarPosition(lat, lng, ts, tzOffset).sunset;
}

1 comment:

  1. Nice job. I wonder if, for the elevation, you've considered atmospheric refraction.

    ReplyDelete