<p>The equations necessary for a correct answer are somewhat lengthy, for a web interface that accepts the location and prints a calendar try: <a href="https://sunrise-sunset.org/" rel="nofollow noreferrer">https://sunrise-sunset.org/</a>.</p>

<p>For a lengthy explanation of the mathematics see Wikipedia's <a href="https://en.wikipedia.org/wiki/Equation_of_time#Calculating_the_equation_of_time" rel="nofollow noreferrer">Calculating the equation of time</a> and <a href="https://en.wikipedia.org/wiki/Daytime" rel="nofollow noreferrer">Daytime</a> webpages for information about the <strong>duration</strong> of rising and setting, along with information about the differences in <a href="https://en.wikipedia.org/wiki/Daytime#Variations_in_solar_noon" rel="nofollow noreferrer">solar noon</a> and the <a href="https://en.wikipedia.org/wiki/Twilight" rel="nofollow noreferrer">three types of twilight</a>.</p>

<p><a href="https://i.sstatic.net/1X8Dy.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/1X8Dy.png" alt="Twilight Types"></a></p>

<p>This is a <a href="http://www.stargazing.net/kepler/sunrise.html" rel="nofollow noreferrer">QBasic implementation by Keith Burnett</a> for twilight times:</p>

<pre><code>'********************************************************************
'   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 &lt; 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$) &lt; 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) &gt; .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 &gt; 1
    correction = 0
  CASE cosc &lt; -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
'******************************************************************
</code></pre>

<p>This is an implementation by <a href="https://astronomy.stackexchange.com/a/13053/18054">Barry Carter</a> of rising/setting <strong>duration</strong> (see <a href="https://github.com/barrycarter/bcapps/tree/master/ASTRO" rel="nofollow noreferrer">his GitHub</a> for associated files):</p>

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

#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;string.h&gt;
#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, &amp;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&lt;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) &lt; .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)&lt;1 || abs(rstart-stime)&lt;1) {continue;}

    // TODO: this really inefficient (and not even always correct?)
    if (bc_sky_elev(6, lat, lon, 0., rstart, "Sun", 0.) &gt;
    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;
}
</code></pre>

<hr>

<p>For an explanation see <a href="https://en.wikipedia.org/wiki/Sunrise_equation" rel="nofollow noreferrer">Wikipedia's Sunrise Equation</a> webpage.</p>

<p>This is the algorithm for rising/setting <strong>times</strong>:</p>

<pre><code>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 &gt;  1) 
      the sun never rises on this location (on the specified date)
    if (cosH &lt; -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
</code></pre>

<p>Here is an implementation of rising/setting <strong>times</strong> in Basic:</p>

<pre><code>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)&lt;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&gt;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&lt;0 THEN 800
620 D=SQR(D)
630 IF V0&lt;0 AND V2&gt;0 THEN PRINT R$;
640 IF V0&lt;0 AND V2&gt;0 THEN M8=1
650 IF V0&gt;0 AND V2&lt;0 THEN PRINT S$;
660 IF V0&gt;0 AND V2&lt;0 THEN W8=1
670 E=(-B+D)/(2*A)
680 IF E&gt;1 OR E&lt;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&lt;0 THEN AZ=AZ+180
770 IF AZ&lt;0 THEN AZ=AZ+360
780 IF AZ&gt;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&lt;0 THEN PRINT M3$
880 IF V2&gt;0 THEN PRINT M4$
890 RETURN
900 '
910 '   Fundamental arguments
920 '     (Van Flandern &amp;
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 --&gt; JD
1180 INPUT "Year, Month, Day";Y,M,D
1190 G=1: IF Y&lt;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&gt;=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 &amp; Telescope for August 1994, page 84.
</code></pre>

CC BY-SA 3.0
2018-01-06T05:09:09.270Z
@barrycarter - Thanks for your input. While I had satisfied Alby's question I've now added info about duration, solar noon and **three** definitions of twilight (who knew?!) and a plug for your GitHub along with a short QBasic program.
Rob
CC BY-SA 3.0
2018-01-02T17:14:46.567Z
2018-01-06T05:03:51.540Z
2018-01-06T05:03:51.540Z
Rob
<p>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.</p>

<the-sun><earth>
Expression for length of sunrise/sunset as function of latitude and day-of-year
Alby
# Expression for length of sunrise/sunset as function of latitude and day-of-year

(Asked by: username_0 on 2018-01-02)

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.

**Tags:** the-sun, earth


---

## 1 answers

### Answer by username_1 on 2018-01-02. Score: 3 | Accepted answer

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](https://en.wikipedia.org/wiki/Equation_of_time#Calculating_the_equation_of_time) and [Daytime](https://en.wikipedia.org/wiki/Daytime) webpages for information about the **duration** of rising and setting, along with information about the differences in [solar noon](https://en.wikipedia.org/wiki/Daytime#Variations_in_solar_noon) and the [three types of twilight](https://en.wikipedia.org/wiki/Twilight).

[![Twilight Types](https://i.sstatic.net/1X8Dy.png)](https://i.sstatic.net/1X8Dy.png)

This is a [QBasic implementation by Keith Burnett](http://www.stargazing.net/kepler/sunrise.html) 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](https://astronomy.stackexchange.com/a/13053/18054) of rising/setting **duration** (see [his GitHub](https://github.com/barrycarter/bcapps/tree/master/ASTRO) 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](https://en.wikipedia.org/wiki/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.

```

> **Comments:**
> **username_1**: @barrycarter - Thanks for your input. While I had satisfied username_0's question I've now added info about duration, solar noon and **three** definitions of twilight (who knew?!) and a plug for your GitHub along with a short QBasic program.  

---

