' Sunrise and Sunset calculation from http://users.electromagnetic.net/bu/astro/sunrise-set.php
' Liberty BASIC (or LBB) version by Richard Russell, http://www.rtrussell.co.uk/  24th May 2014

' We will find the sunrise and sunset for this date:

    day = 20
    month = 5
    year = 2014

' At this location:

    lw = 0.0  ' Longitude West, degrees
    ln = 51.1 ' Latitude North, degrees

' Find the Julian Day for the date in question:

    JD = 2415386 + date$(month;"/";day;"/";year)

' Calculate the Julian Cycle number:

    nStar = (JD - 2451545.0009) - (lw / 360)
    n = int(nStar + 0.5) 

' Approximate the Julian date of solar noon:

    JStar = 2451545.0009 + (lw / 360) + n 

' Calculate the mean solar anomaly:

    M = (357.5291 + 0.98560028 * (JStar - 2451545)) mod 360 

' Calculate the equation of center:

    C = (1.9148 * sinrad(M)) + (0.0200 * sinrad(2*M)) + (0.0003 * sinrad(3*M)) 

' Calculate the ecliptical longitude of the sun:

    lambda = (M + 102.9372 + C + 180) mod 360 

' Calculate an accurate Julian date for solar noon:

    Jtransit = JStar + (0.0053 * sinrad(M)) - (0.0069 * sinrad(2 * lambda)) 

' Calculate the declination of the sun:

    delta = degasn(sinrad(lambda) * sinrad(23.45)) 

' Calculate the hour angle:

    H = degacs((sinrad(-0.83) - sinrad(ln) * sinrad(delta)) / (cosrad(ln) * cosrad(delta))) 

' Go back through the approximation again:

    JStarStar = 2451545.0009 + ((H + lw) / 360) + n 

' Calculate sunset time:

    Jset = JStarStar + (0.0053 * sinrad(M)) - (0.0069 * sinrad(2 * lambda)) 

' Assume that solar noon is half-way between sunrise and sunset (valid for latitudes < 60):

    Jrise = 2 * Jtransit - Jset

' Print results

print "Sunrise: JD = "; using("#######.########", Jrise); ", UTC = "; timefromJD$(Jrise)
print "Sunset:  JD = "; using("#######.########", Jset);  ", UTC = "; timefromJD$(Jset)

print "According to Wolfram Alpha sunrise is at 04:04 UTC and sunset at 19:52 UTC"
end

function sinrad(a)
    sinrad =sin(a / 57.29577951308232088)
end function

function cosrad(a)
    cosrad =cos(a / 57.29577951308232088)
end function

function degasn(n)
    degasn = 57.29577951308232088 * asn(n)
end function

function degacs(n)
    degacs = 57.29577951308232088 * acs(n)
end function

function timefromJD$(jd)
    jd = jd + 0.5
    secs = int((jd - int(jd)) * 24 * 60 * 60 + 0.5)
    mins = int(secs / 60)
    hour = int(mins / 60)
    timefromJD$ = date$(int(jd) - 2415386);" "; _
    dig2$(hour);":";dig2$(mins mod 60);":";dig2$(secs mod 60)
end function

function dig2$(n)
    dig2$ = right$("0";n, 2)
end function