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.
' ************************************
' ** 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
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
