# 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; "°."

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

delta       =asn( sind)         *57.29577951                                                    '   Declination of the Sun in °.
'   NB wasn't converting 'delta' to degrees
print "sind", sind,
print "delta", 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

end function

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
``` 