Skip to main content
Here is an implementation of the Sunrise/Sunset Algorithm in VBA.

Official times of sunrise og sunset are calculated with variables depending on location given in latitude and longitude and corrected for difference between UCT/GMT and local time. Tested in Excel 2007.

Source of pseudocode algorithm: Almanac for Computers, 1990; Nautical Almanac Office; United States Naval Observatory; Washington, DC 20392

Inputs: day (integer), month (integer), year (double), latitude (double), longitude (double), localOffset (difference between local time and UCT/GMT), status (return time of official sun rise=0; return time of official sun set=1)

day, month and year are return values of standard Excel functions DAY, MONTH and YEAR applied to a cell formatted as an date.

My own location: 57.07 north and 9.94 east (eastern coordinates are positive, western negative)
This year the sun will rise at 5:49:04 AM and set at 9:00:03 PM at my home address on my birthday.

Option Explicit

Public Const cPI As Double = 3.14159265358979

Public Function IntervalCorrection(a, b)
    Dim Temp As Double
    Temp = a
    If Temp < 0 Then
    Temp = Temp + b
    End If
    If Temp > b Then
    Temp = Temp - b
    End If
    IntervalCorrection = Temp
End Function

'Implementation of trigonometric functions taking degrees as argument instead of radians
Public Function cosDeg(degree As Double)
cosDeg = Cos(degree * cPI / 180#)
End Function
Public Function sinDeg(degree As Double)
sinDeg = Sin(degree * cPI / 180#)
End Function
Public Function tanDeg(degree As Double)
tanDeg = Tan(degree * cPI / 180#)
End Function
'Implementation of trigonometric functions returning degrees instead of radians
Public Function acosDeg(arg As Double)
acosDeg = 180# * Application.WorksheetFunction.Acos(arg) / cPI
End Function
Public Function asinDeg(arg As Double)
asinDeg = 180# * Application.WorksheetFunction.Asin(arg) / cPI
End Function
Public Function atanDeg(arg As Double)
atanDeg = 180# * Atn(arg) / cPI
End Function

Public Function suntimes(day As Integer, month As Integer, year As Integer, latitude As Double, longitude As Double, localOffset As Integer, state As Double)
'The day of the year and associated intermediate variables
Dim N As Double
Dim N1 As Double
Dim N2 As Double
Dim N2A As Double
Dim N3 As Double
'Hour value and approximate time
Dim lngHour As Double
Dim t As Double
'Suns mean anomaly
Dim M As Double
'Suns true longitude and intermediate variables
Dim L As Double
Dim L1 As Double
Dim L2 As Double
'Suns right ascension
Dim RA As Double
'Intermediate variables in calculation of ascension value
Dim Lquadrant As Double
Dim RAquadrant As Double
'Suns declination
Dim sinDec As Double
Dim cosDec As Double
'Suns local hour value and intermediate values
Dim H As Double
Dim cosH As Double
Dim cosH1 As Double
Dim cosH2 As Double
Dim cosH3 As Double
'Local mean time of rising
Dim LT As Double
'Coordinated Universal Time of rising
Dim UT As Double
'Definition of zenith - official 90.0 degrees, civil 96 50' degrees, nautical 102 degrees, astronomical 108 degrees
Dim zenith As Double
zenith = 90#

'Calculate the day of the year
N1 = WorksheetFunction.Floor(275# * month / 9#, 1)
N2 = WorksheetFunction.Floor((month + 9#) / 12#, 1)
N2A = WorksheetFunction.Floor(year / 4#, 1)
N3 = 1 + WorksheetFunction.Floor((year - 4# * N2A + 2#) / 3#, 1)
N = N1 - (N2 * N3) + day - 30#
Debug.Print "The value of variable N is: " & N
'Convert the longitude to hour value and calculate an approximate time
lngHour = longitude / 15#
t = N + ((6# - lngHour) / 24#)
If (state > 0) Then
t = N + ((18# - lngHour) / 24#)
End If

'Calculate the Sun's mean anomaly
M = (0.9856 * t) - 3.289
Debug.Print "The value of variable M is: " & M
'Calculate the Suns true longitude
L1 = 1.916 * sinDeg(M)
L2 = 0.02 * sinDeg(2# * M)
L = M + L1 + L2 + 282.634
L = IntervalCorrection(L, 360)

'Calculate the Suns right ascension
RA = atanDeg(0.91764 * tanDeg(L))
RA = IntervalCorrection(RA, 360)

'Right ascension needs to be in the same quadrant as L
Lquadrant = (Application.WorksheetFunction.Floor(L / 90#, 1)) * 90#
RAquadrant = (Application.WorksheetFunction.Floor(RA / 90#, 1)) * 90#
RA = RA + (Lquadrant - RAquadrant)

'Right ascension value needs to be converted into hours
RA = RA / 15#

'Calculate the Suns declination
sinDec = 0.39782 * sinDeg(L)
cosDec = cosDeg(asinDeg(sinDec))

'Calculate the Suns local hour angle
cosH1 = cosDeg(zenith)
cosH2 = sinDec * sinDeg(latitude)
cosH3 = cosDec * cosDeg(latitude)
cosH = (cosH1 - cosH2) / cosH3

If (cosH > 1) Then
    MsgBox ("The sun do not rise")
    Exit Function
 End If

 If (cosH < -1) Then
    MsgBox ("The sun do not set")
    Exit Function
 End If

'Finish calculating H and convert into hours
H = (360 - acosDeg(cosH)) / 15#
If (state > 0) Then
H = acosDeg(cosH) / 15#
End If

'Calculate local mean time of rising
LT = H + RA - (0.06571 * t) - 6.622
'Adjust back to UTC
UT = LT - lngHour
UT = IntervalCorrection(UT, 24)
Debug.Print "The value of variable cPI is: " & cPI
Debug.Print "The value of variable UT is: " & UT
Debug.Print "The value of variable localOffset is: " & localOffset
'Return value
suntimes = UT + localOffset
End Function


Popular posts from this blog

HackRF on Windows 8

This technical note is based on an extract from thread. I have made several changes and added recommendations. I have experienced lot of latency using GnuRadio and HackRF on Pentoo Linux, so I wanted to try out GnuRadio on Windows.

HackRF One is a transceiver, so besides SDR capabilities, it can also transmit signals, inkluding sweeping a given range, uniform and Gaussian signals. Pentoo Linux provides the most direct access to HackRF and toolboxes. Install Pentoo Linux on a separate drive, then you can use osmocom_siggen from a terminal to transmit signals such as near-field GSM bursts, which will only be detectable within a meter.

Installation of MGWin and cmake: Download and install the following packages:
- MinGW Setup (Go to the Installer directory and download setup file)
- CMake (I am using CMake 3.2.2 and I installed it in C:\CMake, this path is important in the commands we must send in the MinGW shell)
Download and extract the packages respectively in the path C:\MinGW\msys\…

Example: Beeswarm plot in R


data <- read.dta("C:/Users/hellmund/Documents/MyStataDataFile.dta")





png(file="C:/Users/hellmund/Documents/il6.png", bg="transparent")

beeswarm(data$il6~data$group,data=data, method=c("swarm"),pch=16,pwcol=data$Gender,xlab='',ylab='il6',ylim=c(0,20))


boxplot(data$il6~data$group, data=data, add = T, names = c("","",""), col="#0000ff22")

Example: Business cards typeset with LaTeX

So you enjoy the quality of a professional typesetting system? You got Avery labels, a working MikTeX and the ticket package installed...
You might find some assistance from a half criminal paranoid zealot system administrator, willing to guide you through a dinosaur kingdom of TeX ... but that kind of assistance might also just leave you with nothing.

It was easy to get the layout of the labels with the option zw32010, but how about page margins? I tried to set things straight with the layouts package (\usepackage{layouts}\currentpage \pagedesign), but then there was still some unwanted white space and margins...

To make things less complicated I decided to make a single card. The solution is a hack because it needs customization (with voffset and hoffset as you see n the TeX code below) but the adjustment is more straightforward, especially if you use the boxed option with ticket.

The card was converted to png with Ghostscript and I could easily print the business cards with Averys …