Expression for length of sunrise/sunset as function of latitude and day-of-year

Question

I would like to know how the length of sunrise/sunset (and also time between civil/nautical/astronomical twilight) varies with latitude and with day of year. Ideally what I am looking for is an expression that takes latitude and day-of-year as arguments and outputs the time taken for the sun to rise or set at that latitude on that day.

Answer 1

The equations necessary for a correct answer are somewhat lengthy, for a web interface that accepts the location and prints a calendar try: https://sunrise-sunset.org/.

For a lengthy explanation of the mathematics see Wikipedia's Calculating the equation of time and Daytime webpages for information about the duration of rising and setting, along with information about the differences in solar noon and the three types of twilight.

This is a QBasic implementation by Keith Burnett for twilight times:

'********************************************************************
'   Sun rise/set /twilight
'   From Explanatory Supplement
'   definitions and functions
DEFDBL A-Z
pi = 4 * ATN(1)
tpi = 2 * pi
twopi = tpi
degs = 180 / pi
rads = pi / 180
'
'   the function below returns the true integer part,
'   even for negative numbers
'
DEF FNipart (x) = SGN(x) * INT(ABS(x))
'
'   FNday only works between 1901 to 2099 - see Meeus chapter 7
'
DEF FNday (y, m, d, h) = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
'
'   define some arc cos and arc sin functions
'
DEF FNacos (x)
    s = SQR(1 - x * x)
    FNacos = ATN(s / x)
END DEF
DEF FNasin (x)
    c = SQR(1 - x * x)
    FNasin = ATN(x / c)
END DEF
'
'   the function below returns an angle in the range
'   0 to two pi
'
DEF FNrange (x)
    b = x / tpi
    a = tpi * (b - FNipart(b))
    IF a < 0 THEN a = tpi + a
    FNrange = a
END DEF
'
'       the function below returns the time in 24 hour
'   notation - hhmm - given a decimal hours figure.
DEF FNdegmin$ (d)
c = ABS(d)
a = INT(c)
b = 60 * (c - a)
d$ = STR$(INT(a * 100 + b + .5))
d$ = RIGHT$(d$, LEN(d$) - 1)
WHILE LEN(d$) < 4
    d$ = "0" + d$
WEND
FNdegmin$ = d$
END DEF
CLS
'
'    get the date and geographical position from user
'
INPUT "     lat : ", glat
INPUT "    long : ", glong
INPUT "    zone : ", zone
INPUT "rise/set : ", riset
INPUT " Sun alt : ", altitude
INPUT "   year  : ", y
INPUT "   month : ", m
INPUT "   day   : ", d
day = FNday(y, m, d, 0)
PRINT "  day no : "; day
utold = pi
utnew = 0
sinalt = SIN(altitude * rads) 'go for the sunrise/sunset altitude first
sinphi = SIN(glat * rads)
cosphi = COS(glat * rads)
glong = glong * rads
DO WHILE ABS(utold - utnew) > .001
  utold = utnew
  days = day + utold / tpi
  t = days / 36525
  '     get arguments of Sun's orbit
  L = FNrange(4.8949504201433# + 628.331969753199# * t)
  G = FNrange(6.2400408# + 628.3019501# * t)
  ec = .033423 * SIN(G) + .00034907# * SIN(2 * G)
  lambda = L + ec
  E = -1 * ec + .0430398# * SIN(2 * lambda) - .00092502# * SIN(4 * lambda)
  obl = .409093 - .0002269# * t
  delta = FNasin(SIN(obl) * SIN(lambda))
  GHA = utold - pi + E
  cosc = (sinalt - sinphi * SIN(delta)) / (cosphi * COS(delta))
  SELECT CASE cosc
  CASE cosc > 1
    correction = 0
  CASE cosc < -1
    correction = pi
  CASE ELSE
    correction = FNacos(cosc)
  END SELECT
  utnew = FNrange(utold - (GHA + glong + riset * correction))
LOOP
PRINT "    UT   : "; FNdegmin$(utnew * degs / 15)
PRINT "   zone  : "; FNdegmin$(utnew * degs / 15 + zone)
END
'******************************************************************

This is an implementation by Barry Carter of rising/setting duration (see his GitHub for associated files):

// http://astronomy.stackexchange.com/questions/12824/how-long-does-a-sunrise-or-sunset-take

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "SpiceUsr.h"
#include "SpiceZfc.h"
#include "/home/barrycarter/BCGIT/ASTRO/bclib.h"

// since the sun travels no more than 362 degrees per day, and has an
// angular diameter of no less than 31 minutes, the fastest possible
// sunrise is 31 minutes/362 degrees or 2 minutes; thus, 120 below in
// bc_between below is valid

int main(int argc, char **argv) {

  SpiceDouble radii[3];
  SpiceInt n;
  char direction[10];

  furnsh_c("/home/barrycarter/BCGIT/ASTRO/standard.tm");

  // the year 2015
  double stime = 1419984000, etime = 1451692800;

  // test
  // double stime =  1388559600, etime = 1419984000;

  // lat/lon from argv
  double lat = atof(argv[1]), lon = atof(argv[2]);

  // the sun's radii
  bodvrd_c("SUN", "RADII", 3, &n, radii);


  double *results = bc_between(9, lat*rpd_c(), lon*rpd_c(), 0., stime, etime,
                   "Sun", -34/60.*rpd_c(), 120., radii[0]);

  // NOTE: can NOT ignore result 0 and 1, but must check if they are borders
  for (int i=0; i<1000; i++) {

    double rstart = results[2*i], rend = results[2*i+1];

    // if we start seeing 0s, we are out of true answers
    // TODO: this may break near 1970
    if (abs(rstart) < .001) {break;}

    // if the end result is too close to etime, result is inaccurate
    // same if start result is too close to stime
    if (abs(rend-etime)<1 || abs(rstart-stime)<1) {continue;}

    // TODO: this really inefficient (and not even always correct?)
    if (bc_sky_elev(6, lat, lon, 0., rstart, "Sun", 0.) >
    bc_sky_elev(6, lat, lon, 0., rend, "Sun", 0.)) {
      strcpy(direction,"SET");
    } else {
      strcpy(direction,"RISE");
    }

    // the "day" and length of sunrise/sunset
    printf("%f %f %s %f %f %f\n", lat, lon, direction, rstart, rend, rend-rstart);
  }
  return 0;
}

---

For an explanation see Wikipedia's Sunrise Equation webpage.

This is the algorithm for rising/setting times:

Sunrise/Sunset Algorithm

Source:
    Almanac for Computers, 1990
    published by Nautical Almanac Office
    United States Naval Observatory
    Washington, DC 20392

Inputs:
    day, month, year:      date of sunrise/sunset
    latitude, longitude:   location for sunrise/sunset
    zenith:                Sun's zenith for sunrise/sunset
      offical      = 90 degrees 50'
      civil        = 96 degrees
      nautical     = 102 degrees
      astronomical = 108 degrees

    NOTE: longitude is positive for East and negative for West
        NOTE: the algorithm assumes the use of a calculator with the
        trig functions in "degree" (rather than "radian") mode. Most
        programming languages assume radian arguments, requiring back
        and forth convertions. The factor is 180/pi. So, for instance,
        the equation RA = atan(0.91764 * tan(L)) would be coded as RA
        = (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree
        answer with a degree input for L.


1. first calculate the day of the year

    N1 = floor(275 * month / 9)
    N2 = floor((month + 9) / 12)
    N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
    N = N1 - (N2 * N3) + day - 30

2. convert the longitude to hour value and calculate an approximate time

    lngHour = longitude / 15

    if rising time is desired:
      t = N + ((6 - lngHour) / 24)
    if setting time is desired:
      t = N + ((18 - lngHour) / 24)

3. calculate the Sun's mean anomaly

    M = (0.9856 * t) - 3.289

4. calculate the Sun's true longitude

    L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
    NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

5a. calculate the Sun's right ascension

    RA = atan(0.91764 * tan(L))
    NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

5b. right ascension value needs to be in the same quadrant as L

    Lquadrant  = (floor( L/90)) * 90
    RAquadrant = (floor(RA/90)) * 90
    RA = RA + (Lquadrant - RAquadrant)

5c. right ascension value needs to be converted into hours

    RA = RA / 15

6. calculate the Sun's declination

    sinDec = 0.39782 * sin(L)
    cosDec = cos(asin(sinDec))

7a. calculate the Sun's local hour angle

    cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))

    if (cosH >  1) 
      the sun never rises on this location (on the specified date)
    if (cosH < -1)
      the sun never sets on this location (on the specified date)

7b. finish calculating H and convert into hours

    if if rising time is desired:
      H = 360 - acos(cosH)
    if setting time is desired:
      H = acos(cosH)

    H = H / 15

8. calculate local mean time of rising/setting

    T = H + RA - (0.06571 * t) - 6.622

9. adjust back to UTC

    UT = T - lngHour
    NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24

10. convert UT value to local time zone of latitude/longitude

    localT = UT + localOffset

Here is an implementation of rising/setting times in Basic:

10 '         Sunrise-Sunset
20 GOSUB 300
30 INPUT "Lat, Long (deg)";B5,L5
40 INPUT "Time zone (hrs)";H
50 L5=L5/360: Z0=H/24
60 GOSUB 1170: T=(J-2451545)+F
70 TT=T/36525+1: ' TT = centuries
80 '               from 1900.0
90 GOSUB 410: T=T+Z0
100 '
110 '       Get Sun's Position
120 GOSUB 910: A(1)=A5: D(1)=D5
130 T=T+1
140 GOSUB 910: A(2)=A5: D(2)=D5
150 IF A(2)<A(1) THEN A(2)=A(2)+P2
160 Z1=DR*90.833: ' Zenith dist.
170 S=SIN(B5*DR): C=COS(B5*DR)
180 Z=COS(Z1): M8=0: W8=0: PRINT
190 A0=A(1): D0=D(1)
200 DA=A(2)-A(1): DD=D(2)-D(1)
210 FOR C0=0 TO 23
220 P=(C0+1)/24
230 A2=A(1)+P*DA: D2=D(1)+P*DD
240 GOSUB 490
250 A0=A2: D0=D2: V0=V2
260 NEXT
270 GOSUB 820: '  Special msg?
280 END
290 '
300 '        Constants
310 DIM A(2),D(2)
320 P1=3.14159265: P2=2*P1
330 DR=P1/180: K1=15*DR*1.0027379
340 S$="Sunset at  "
350 R$="Sunrise at "
360 M1$="No sunrise this date"
370 M2$="No sunset this date"
380 M3$="Sun down all day"
390 M4$="Sun up all day"
400 RETURN
410 '     LST at 0h zone time
420 T0=T/36525
430 S=24110.5+8640184.813*T0
440 S=S+86636.6*Z0+86400*L5
450 S=S/86400: S=S-INT(S)
460 T0=S*360*DR
470 RETURN
480 '
490 '  Test an hour for an event
500 L0=T0+C0*K1: L2=L0+K1
510 H0=L0-A0: H2=L2-A2
520 H1=(H2+H0)/2: '  Hour angle,
530 D1=(D2+D0)/2: '  declination,
540 '                at half hour
550 IF C0>0 THEN 570
560 V0=S*SIN(D0)+C*COS(D0)*COS(H0)-Z
570 V2=S*SIN(D2)+C*COS(D2)*COS(H2)-Z
580 IF SGN(V0)=SGN(V2) THEN 800 
590 V1=S*SIN(D1)+C*COS(D1)*COS(H1)-Z
600 A=2*V2-4*V1+2*V0: B=4*V1-3*V0-V2
610 D=B*B-4*A*V0: IF D<0 THEN 800
620 D=SQR(D)
630 IF V0<0 AND V2>0 THEN PRINT R$;
640 IF V0<0 AND V2>0 THEN M8=1
650 IF V0>0 AND V2<0 THEN PRINT S$;
660 IF V0>0 AND V2<0 THEN W8=1
670 E=(-B+D)/(2*A)
680 IF E>1 OR E<0 THEN E=(-B-D)/(2*A)
690 T3=C0+E+1/120: ' Round off
700 H3=INT(T3): M3=INT((T3-H3)*60)
710 PRINT USING "##:##";H3;M3;
720 H7=H0+E*(H2-H0)
730 N7=-COS(D1)*SIN(H7)
740 D7=C*SIN(D1)-S*COS(D1)*COS(H7)
750 AZ=ATN(N7/D7)/DR
760 IF D7<0 THEN AZ=AZ+180
770 IF AZ<0 THEN AZ=AZ+360
780 IF AZ>360 THEN AZ=AZ-360
790 PRINT USING ",  azimuth ###.#";AZ
800 RETURN
810 '
820 '   Special-message routine
830 IF M8=0 AND W8=0 THEN 870
840 IF M8=0 THEN PRINT M1$
850 IF W8=0 THEN PRINT M2$
860 GOTO 890
870 IF V2<0 THEN PRINT M3$
880 IF V2>0 THEN PRINT M4$
890 RETURN
900 '
910 '   Fundamental arguments
920 '     (Van Flandern &
930 '     Pulkkinen, 1979)
940 L=.779072+.00273790931*T
950 G=.993126+.0027377785*T
960 L=L-INT(L): G=G-INT(G)
970 L=L*P2: G=G*P2
980 V=.39785*SIN(L)
990 V=V-.01000*SIN(L-G)
1000 V=V+.00333*SIN(L+G)
1010 V=V-.00021*TT*SIN(L)
1020 U=1-.03349*COS(G)
1030 U=U-.00014*COS(2*L)
1040 U=U+.00008*COS(L)
1050 W=-.00010-.04129*SIN(2*L)
1060 W=W+.03211*SIN(G)
1070 W=W+.00104*SIN(2*L-G)
1080 W=W-.00035*SIN(2*L+G)
1090 W=W-.00008*TT*SIN(G)
1100 '
1110 '    Compute Sun's RA and Dec
1120 S=W/SQR(U-V*V)
1130 A5=L+ATN(S/SQR(1-S*S))
1140 S=V/SQR(U):D5=ATN(S/SQR(1-S*S))
1150 R5=1.00021*SQR(U)
1160 RETURN
1165 '
1170 '     Calendar --> JD
1180 INPUT "Year, Month, Day";Y,M,D
1190 G=1: IF Y<1583 THEN G=0
1200 D1=INT(D): F=D-D1-.5
1210 J=-INT(7*(INT((M+9)/12)+Y)/4)
1220 IF G=0 THEN 1260
1230 S=SGN(M-9): A=ABS(M-9)
1240 J3=INT(Y+S*INT(A/7))
1250 J3=-INT((INT(J3/100)+1)*3/4)
1260 J=J+INT(275*M/9)+D1+G*J3
1270 J=J+1721027+2*G+367*Y
1280 IF F>=0 THEN 1300
1290 F=F+1: J=J-1
1300 RETURN
1310 '
1320 '   This program by Roger W. Sinnott calculates the times of sunrise
1330 '   and sunset on any date, accurate to the minute within several
1340 '   centuries of the present.  It correctly describes what happens in the 
1350 '   arctic and antarctic regions, where the Sun may not rise or set on
1360 '   a given date.  Enter north latitudes positive, west longitudes
1370 '   negative.  For the time zone, enter the number of hours west of
1380 '   Greenwich (e.g., 5 for EST, 4 for EDT).  The calculation is
1390 '   discussed in Sky & Telescope for August 1994, page 84.