Sunrise calculation

sunrise.cpp

/* source: http://aa.quae.nl/en/reken/zonpositie.html */

// source: http://en.wikipedia.org/wiki/Julian_day
int julian_jdn(int year, int month, int day) { // convert Gregorian calendar date to Julian Day Number
  int a = (14 - month) / 12;
  int y = year + 4800 - a;
  int m = month + 12 * a - 3;
  return day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
}

double to_rad(double angle) {
  return PI * angle / 180;
}

double to_deg(double angle) {
  return 180 * angle / PI;
}

double full_circle(double angle) { // draw a full circle (370' to (370'- 360') = 10')
  return angle - (int)(angle / 360) * 360;
}

double earth_mean_anomaly(double jdn) {
  return full_circle(357.5291 + 0.98560028 * (jdn - 2451545));
}

double earth_equation_of_center(double mean_anomaly) {
  return 1.9148 * sin(to_rad(mean_anomaly)) +
         0.0200 * sin(2 * to_rad(mean_anomaly)) +
         0.0003 * sin(3 * to_rad(mean_anomaly));
}

double earth_true_anomaly(double jdn) {
  return earth_mean_anomaly(jdn) + earth_equation_of_center(earth_mean_anomaly(jdn));
}

double sun_ecliptic_longitude(double true_anomaly) {
  return full_circle(true_anomaly + 102.9372 + 180);
}

double earth_declination(double ecliptic_longitude) {
  return to_deg(asin(
           sin(to_rad(ecliptic_longitude)) * sin(to_rad(23.45))
         ));
}

double sun_sunrise(double ln, double declination, double jdn, double le, double mean_anomaly) {
  double h = -full_circle(to_deg(acos(
               (sin(to_rad(-0.83)) - sin(to_rad(ln)) * sin(to_rad(declination))) /
               (cos(to_rad(ln)) * cos(to_rad(declination)))
             )));

  return 2451545 + 0.0009 + (h - le) / 360 +
       (int)((jdn - 2451545 - 0.0009) / 1 - (h - le) / 360) +
       0.0053 * sin(to_rad(mean_anomaly)) -
       0.0069 * sin(to_rad(2 * full_circle(mean_anomaly + 102.9372 + 180)));
}

void calc_sunrise(double * &sunrise, int &sun_s, double la = 47, double lo = 19, int start_year = 2019, int start_month = 4, int start_day = 15, int end_year = 2019, int end_month = 8, int end_day = 31) {
  double ta, l, delta;
  int j = julian_jdn(start_year, start_month, start_day);
  int ej = julian_jdn(end_year, end_month, end_day);
  sun_s = ej - j + 1;
  sunrise = new double[sun_s];
  for (int i = 0; j+i <= ej; i++) {
    ta = earth_true_anomaly(j+i);
    l = sun_ecliptic_longitude(ta);
    delta = earth_declination(l);
    sunrise[i] = sun_sunrise(la, delta, j+i, lo, ta);
  }
}

sunrise.py

from math import sin, cos, asin, acos

PI = 3.141592653589793238462643383279
defaults = [47, 49, 2019, 4, 15, 2019, 8, 31]

def julian_jdn(year, month, day): # convert Gregorian calendar date to Julian Day Number
  a = int((14 - month) / 12)
  y = int(year + 4800 - a)
  m = int(month + 12 * a - 3)
  return int(day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045)

def to_rad(angle):
  global PI
  return PI * angle / 180

def to_deg(angle):
  global PI
  return 180 * angle / PI

def full_circle(angle): # draw a full circle (370' to (370'- 360') = 10')
  return angle - int(angle / 360) * 360

def earth_mean_anomaly(jdn):
  return full_circle(357.5291 + 0.98560028 * (jdn - 2451545))

def earth_equation_of_center(mean_anomaly):
  return 1.9148 * sin(to_rad(mean_anomaly)) +\
         0.0200 * sin(2 * to_rad(mean_anomaly)) +\
         0.0003 * sin(3 * to_rad(mean_anomaly))

def earth_true_anomaly(jdn):
  return earth_mean_anomaly(jdn) + earth_equation_of_center(earth_mean_anomaly(jdn))

def sun_ecliptic_longitude(true_anomaly):
  return full_circle(true_anomaly + 102.9372 + 180)

def earth_declination(ecliptic_longitude):
  return to_deg(asin(sin(to_rad(ecliptic_longitude)) * sin(to_rad(23.45))))

def sun_sunrise(ln, declination, jdn, le, mean_anomaly):
  h = -full_circle(to_deg(acos(\
        (sin(to_rad(-0.83)) - sin(to_rad(ln)) * sin(to_rad(declination))) /\
        (cos(to_rad(ln)) * cos(to_rad(declination)))\
      )))

  return 2451545 + 0.0009 + (h - le) / 360 +\
       int((jdn - 2451545 - 0.0009) / 1 - (h - le) / 360) +\
       0.0053 * sin(to_rad(mean_anomaly)) -\
       0.0069 * sin(to_rad(2 * full_circle(mean_anomaly + 102.9372 + 180)))

def calc_sunrise(la=defaults[0], lo=defaults[1],\
                 start_year=defaults[2], start_month=defaults[3], start_day=defaults[4],\
                 end_year=defaults[5], end_month=defaults[6], end_day=defaults[7]):
  j = julian_jdn(start_year, start_month, start_day)
  ej = julian_jdn(end_year, end_month, end_day)

  def f(i):
    ta = earth_true_anomaly(j+i)
    delta = earth_declination(sun_ecliptic_longitude(ta))
    return sun_sunrise(la, delta, j+i, lo, ta)

  return map(f, range(ej-j))

if __name__ == '__main__':
  import sys
  print('\n'.join(str(e) for e in calc_sunrise(*list(map(float, sys.argv[1:])))))

sunrise_console.cpp

/* source: http://aa.quae.nl/en/reken/zonpositie.html */

// g++ -O3 sunrise.cpp -o sunrise
// sunrise > napkelte_2019_apr_15_aug_31_JD_UTC.txt

#include <stdio.h>
#include <math.h>
#define PI 3.141592653589793238462643383279

// source: http://en.wikipedia.org/wiki/Julian_day
int julian_jdn(int day, int month, int year) { // convert Gregorian calendar date to Julian Day Number
  int a = (14 - month) / 12;
  int y = year + 4800 - a;
  int m = month + 12 * a - 3;
  return day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
}

double to_rad(double angle) {
  return PI * angle / 180;
}

double to_deg(double angle) {
  return 180 * angle / PI;
}

double full_circle(double angle) { // draw a full circle (370' to (370'- 360') = 10')
  return angle - (int)(angle / 360) * 360;
}

double earth_mean_anomaly(double jdn) {
  return full_circle(357.5291 + 0.98560028 * (jdn - 2451545));
}

double earth_equation_of_center(double mean_anomaly) {
  return 1.9148 * sin(to_rad(mean_anomaly)) +
         0.0200 * sin(2 * to_rad(mean_anomaly)) +
         0.0003 * sin(3 * to_rad(mean_anomaly));
}

double earth_true_anomaly(double jdn) {
  return earth_mean_anomaly(jdn) + earth_equation_of_center(earth_mean_anomaly(jdn));
}

double sun_ecliptic_longitude(double true_anomaly) {
  return full_circle(true_anomaly + 102.9372 + 180);
}

double earth_declination(double ecliptic_longitude) {
  return to_deg(asin(
           sin(to_rad(ecliptic_longitude)) * sin(to_rad(23.45))
         ));
}

double sun_sunrise(double ln, double declination, double jdn, double le, double mean_anomaly) {
  double h = -full_circle(to_deg(acos(
               (sin(to_rad(-0.83)) - sin(to_rad(ln)) * sin(to_rad(declination))) /
               (cos(to_rad(ln)) * cos(to_rad(declination)))
             )));

  return 2451545 + 0.0009 + (h - le) / 360 +
       (int)((jdn - 2451545 - 0.0009) / 1 - (h - le) / 360) +
       0.0053 * sin(to_rad(mean_anomaly)) -
       0.0069 * sin(to_rad(2 * full_circle(mean_anomaly + 102.9372 + 180)));
}

int main() {
  double ln = 47;  /* north latitude */
  double le = 19;  /* east longitude */
  int year = 2019, month = 4, je = julian_jdn(31, month + 4, year);
  double ta, delta;

  for (int i = julian_jdn(15, month, year); i <= je; i++) {
    ta = earth_true_anomaly(i);
    delta = earth_declination(sun_ecliptic_longitude(ta));
    printf("%.9f\n", sun_sunrise(ln, delta, i, le, ta));
    /*t  = (sunrise - (int)sunrise) > 0.5 ? (sunrise - (int)sunrise) - 0.5 : (sunrise - (int)sunrise) + 0.5;
    printf("%2d:%02d\n", (int)(t * 24) + (gmt + dst), (int)(t * 24 * 60) % 60);*/
  }

  return 0;
}

Source:
https://blog.sleeplessbeastie.eu/2013/05/21/how-to-determine-the-sunrise-and-sunset-times/

To get sunset, just change -full_circle to full_circle in sun_sunrise, and add 1 + after return.