Sunrise and sunset calculation

The following code implements the Wikipedia calculation for sunrise and sunset at a given latitude and longitude. I wrote it after an OP published an incorrect version and asked for help.

It can be added to to compensate for elevation, etc- see web references. Results should be within a fraction of a minute, but assume clear horizon with no elevated land on the line-of-sight towards the Sun. Please note too that there are various different definitions of sunrise!!

The photo is a soft-focus of sunset at Rhosilli- one of my favourite places. Note the left-hand sun-dog- sign of stable air at high altitude carrying ice crystals. These optical effects are fascinating- see my programming example on how rainbows form, and also Flickr photos of Brockenspectre.


This code requires Liberty BASIC ( a free version is available) to run

    '                                       ************************************
    '                                       **         JDquery.bas            **
    '                                       **   Implemented from Wikipedia   **
    '                                       ************************************

    '                       >>>>>>>>>>>>>>>>  NB calculations are at about 9sf   <<<<<<<<<<<<<<

    mainwin 80 24

    lonOmega    =     0.0           '   Longitude ( West  =+ve) of the observer in degrees.
    phi         =    51.1           '   Latitude  ( North =+ve)

    year        =  2014.0           '   Gregorian version of 14 May 2014 at 00:00:00
    month       =     5.0
    day         =    20.0
    hour        =    12.0           '   This is the time JD refers to.
    minute      =     0.0
    second      =     0.0

    month$      ="JanFebMarAprMayJunJulAugSepOctNovDec"

    a           =int( ( 14 -month) /12)
    y           =year  +4800  -a
    m           =month +12 *a -3

    print: print   " Calculation is for Gregorian date and time of "; as2digits$( hour); ":"; as2digits$( minute); ":"; as2digits$( second); " on "; day; "/"; month; "/"; year; "..."

    '   JD  is Julian Date, which is fractional.   ( JDN is Julian Day Number and is an integer)
'   >> I had 152 here...
    JD     =day +int( ( 153 *m +2) /5) +365 *y +int( y /4) -int( y /100) +int( y /400) -32045_
        +( hour -12) /24 +minute /1440 +second /86400                                               '   Current Julian Date as a fractional number
    print "   .... so Julian date ( and time) is "; using( "########.######", JD)                   '   agrees with on-line calculator
print "JD",, JD

    nStar       =JD -2451545.0009   -lonOmega /360

    n           =int( nStar +1 /2)                                                                  '   Current Julian Cycle since Jan 1st 2000
print "nstar, n",, nStar, n  '   NB 'n' wasn't integered

    Jstar       =2451545.0009    +lonOmega /360   +n                                                '   Approx. Solar Noon at longitude specified
print "Jstar",, Jstar

    M           =( 357.5291 +0.98560028 *( Jstar -2451545)) mod 360                                 '   Solar Mean Anomaly in .
print "M",, M; "."

    C           =1.9148 *sinrad( M) +0.0200 *sinrad( 2 *M) +0.0003 *sinrad( 3 *M)                   '   Equation of Centre
print "C",, C

    lambda      =( M +102.9372 +C +180) mod 360                                                     '   Ecliptic Longitude in .
print "lambda",, lambda; "."

    Jtransit    =Jstar +0.0053 *sinrad( M) -0.0069 *sinrad( 2 *lambda)                              '   Hour angle of Solar Transit ie of solar noon
print "Jtransit",, Jtransit

    sind        =sinrad( lambda) *sinrad( 23.45)
    delta       =asn( sind)         *57.29577951                                                    '   Declination of the Sun in .
'   NB wasn't converting 'delta' to degrees
print "sind", sind,
print "delta", delta; "."

    cosOmegaZero=( sinrad( -0.83) -sinrad( phi) *sind) /( cosrad( phi) *cosrad( delta))
    print "cosOmegaZero", cosOmegaZero,
    omegaZero   =acs( cosOmegaZero) *57.29577951                                                    '   Hour angle in .
'   NB wasn't converting 'omegaZero' to degrees
print "omegaZero", omegaZero; "."

    Jset        =2451545.0009 +( omegaZero +lonOmega) /360 +n +0.0053 *sinrad( M) -0.0069 *sinrad( 2 *lambda)
    Jrise       =Jtransit -( Jset -Jtransit)
print: print "Jrise, Jset", using( "#########.######", Jrise), using( "#########.######", Jset)
print , timefromJD$( Jrise), timefromJD$( Jset)

    print
    print " *** OP said should be ***"
    print "Jrise, Jset     2456797.729167              2456798.361806"

    print
    print "Wolfram, when asked for 'sunrise 51.1N 0W 20 May 20.14' gave" +chr$( 13) +_
        "    Input interpretation: sunrise 51 6 N, 0 West, May 20 2014" +chr$( 13) +_
        "    Result: 5:04:52 am BST Tuesday, May 20, 2014"

    end '   ____________________________________________________________________


    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

    function as2digits$( x)
        as2digits$ =right$( "00" +str$( x), 2)
    end function

    function sinrad( th)
        sinrad =sin( th /57.29577951)
    end function

    function cosrad( th)
        cosrad =cos( th /57.29577951)
    end function

    function formattedTime$( H, M)
        formattedTime$ =right$( "00" +str$( H), 2); ":"; right$( "00" +str$( M), 2)
    end function

Output looks like-

 Calculation is for Gregorian date and time of 12:00:00 on 20/5/2014...
   .... so Julian date ( and time) is  2456798.000000
JD                          2456798
nstar, n                    5252.9991     5253
Jstar                       2456798.0
M                           134.888258.
C                           1.33681713
lambda                      59.162275.
Jtransit                    2456798.0
sind          0.34168768    delta         19.9797298.
cosOmegaZero  -0.47512187   omegaZero     118.367285.

Jrise, Jset     2456797.669783              2456798.327379
              05/20/2014 04:04:29         05/20/2014 19:51:26