Overview

If you compare the console sunrise/sunset with an autnoritative source you'll probably find that it's a little off. While it's close enough for most uses, there's an opportunity for improvement. Rather than live with the console's caclulation, we can easily present a more accurate result by using the resources of the host.

One problem with the console is that the clock isn't syncronized. To correct that on a Linux system we need only install a small utility to syncronize the system clock. At the same time we'll also install a shared library called libnova to do all the tedious calculations.

$ apt-get install ntp ntpdate libnova-dev

Included Header Files, Constants and Shared Data

To start, we'll need to include some header files for the library.

#include <libnova/solar.h>
#include <libnova/julian_day.h>
#include <libnova/rise_set.h>
#include <libnova/transform.h>
#include <libnova/utility.h>
#include <libnova/refraction.h>

Next we'll define some constants.

/* station lat/lon in decimal degrees */
#define LAT    ( 40.0163889) // 40 deg 0' 59"
#define LON    (-75.8183333) // 75 deg 49' 6"
/* station elevation in feet above sea level */
#define ELEV   (490.0)

Finally, some common data.

static struct ln_lnlat_posn observer;
static double JD;
static struct ln_helio_posn pos;
static struct ln_equ_posn equ;			
			

Calculations

Instead of performing the common calculations every time we need them we'll call calc_common once to store the results in the shared data structures. Since the Linux host is syncronized, the Julian Date should be pretty accurate. We'll also need a couple utility and unit conversion functions.

void calc_common() {

    observer.lat = LAT;
    observer.lng = LON;
    JD = ln_get_julian_from_sys();
    ln_get_solar_geom_coords(JD, &pos);
    ln_get_solar_equ_coords(JD, &equ);
}

float convFtoC(float tempF) {

    return (float) (5.0 * (tempF - 32) / 9.0);
}

float convinHgtomb(float pressureinHg) {

    return (float) (33.8639 * pressureinHg);
}

float calc_daylight_hours(struct ln_zonedate *r, struct ln_zonedate *s) {

    struct tm br, bs;
    time_t cr, cs;
    float h = 0;
    
    br.tm_year = bs.tm_year = s->years - 1900;
    br.tm_mon = bs.tm_mon = s->months - 1;
    br.tm_mday = bs.tm_mday = s->days;
    br.tm_hour = r->hours;
    br.tm_min = r->minutes;
    br.tm_sec = round(r->seconds);
    bs.tm_hour = s->hours;
    bs.tm_min = s->minutes;
    bs.tm_sec = round(s->seconds);
    cr = mktime(&br);
    cs = mktime(&bs);
    if (cr != -1 && cs != -1)
        h = (float) (cs - cr) / 3600.0;

    return (h);
}

Now we're ready to calculate the standard rise, set, transit and hours of daylight with better accuracy. We'll repeat the calculation for the civil, nautical and astronomical horizons.

int calc_rst(uint16_t *dest) {

    int modptr = 0;
    struct ln_rst_time rst;
    struct ln_zonedate rise, set, transit;

    if (ln_get_solar_rst_horizon(JD, &observer, LN_SOLAR_STANDART_HORIZON, &rst) != 1) {
        ln_get_local_date(rst.rise, &rise);
        ln_get_local_date(rst.transit, &transit);
        ln_get_local_date(rst.set, &set);
        dest[modptr++] = transit.hours;
        dest[modptr++] = transit.minutes;
        dest[modptr++] = (int) round(transit.seconds);
        dest[modptr++] = rise.hours;
        dest[modptr++] = rise.minutes;
        dest[modptr++] = (int) round(rise.seconds);
        dest[modptr++] = set.hours;
        dest[modptr++] = set.minutes;
        dest[modptr++] = (int) round(set.seconds);
        modptr += vp2_pack_float(calc_daylight_hours(&rise, &set), &dest[modptr]);
        }			
		
	...
	
}

We'll also calculate sun's position and an estimate of the clear sky global radiation.

int calc_rad(uint16_t *dest, struct vp2_conditions *cond) {

    int modptr = 0;
    struct ln_hrz_posn hrz;
    
    /* radius vector */
    modptr += vp2_pack_float(cond->R = pos.R, &dest[modptr]);

    /* right ascension and declination */
    ln_get_solar_equ_coords(JD, &equ);
    modptr += vp2_pack_float(equ.ra, &dest[modptr]);
    modptr += vp2_pack_float(equ.dec, &dest[modptr]);

    /* local azimuth and elevation corrected for refraction */
    ln_get_hrz_from_equ(&equ, &observer, JD, &hrz);
    hrz.alt += ln_get_refraction_adj(hrz.alt, convinHgtomb(cond->bar), convFtoC(cond->temp));
    cond->az = ln_range_degrees(hrz.az + 180.0);        /* orient to 0 degrees = North */
    cond->alt = hrz.alt;
    modptr += vp2_pack_float(cond->az, &dest[modptr]);
    modptr += vp2_pack_float(cond->alt, &dest[modptr]);
    
    /* clear sky global radiation */
    cond->csg = (0.79 - (3.75 / cond->alt)) * (Gsc / (pos.R * pos.R)) * sin(ln_deg_to_rad(cond->alt));
    if (cond->csg < 0.0) cond->csg = 0.0;
    modptr += vp2_pack_float(cond->csg, &dest[modptr]);
		
	...
	
}