Skip to content

Instantly share code, notes, and snippets.

@AlessandroMinali
Last active March 21, 2023 20:24
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save AlessandroMinali/f5813341618aa71a429a6f2d6df82056 to your computer and use it in GitHub Desktop.
Save AlessandroMinali/f5813341618aa71a429a6f2d6df82056 to your computer and use it in GitHub Desktop.
Calculate Sun Rise and Set times based. Takes two args: <lat> <lon>
// reference: https://en.wikipedia.org/wiki/Sunrise_equation#Complete_calculation_on_Earth
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
double rad(double deg) {
// PI / 180
return deg * 0.017453292519943295;
}
double deg(double rad) {
// 180 / PI
return rad * 57.29577951308232;
}
int main(int argc, char *argv[]) {
if (argc != 3) { printf("Provide <lat> and <lon>.\n"); exit(1); }
time_t rawtime;
struct tm *info;
char buffer[80];
double lat = strtod(argv[1], NULL);
double lon = strtod(argv[2], NULL);
// Calculate current Julian day
int n = ceil(time(&rawtime) / 86400.0 - 10957.4992);
// Mean solar noon
double J_star = n - (lon / 360.0);
// Solar mean anomaly
double M = fmod(357.5291 + 0.98560028 * J_star, 360);
// Equation of the center
double C = (1.9148 * sinl(rad(M))) + (0.0200 * sin(rad(2 * M))) +(0.0003 * sin(rad(3 * M)));
// Ecliptic longitude
double lam = fmod(M + C + 282.9372, 360);
// Solar transit
double J_transit = 2451545.0 + J_star + (0.0053 * sin(rad(M))) - (0.0069 * sin(rad(2 * lam)));
// Declination of the Sun
// rad(23.44)
double delta = asin(sin(rad(lam)) * sin(0.40910517666747087));
// Hour angle
// rad(-0.83)
double t = sin(-0.014486232791552934) - sin(rad(lat)) * sin(delta);
double b = cos(rad(lat)) * cos(delta);
double omega = deg(acos(t / b));
// Calculate sunrise and sunset
double rise = J_transit - (omega / 360);
double set = J_transit + (omega / 360);
info = localtime(&rawtime);
strftime(buffer, 80, "%Y-%m-%d", info);
printf("Date:\t\t%s\n", buffer);
rawtime = (rise - 2440587.5) * 86400;
info = localtime(&rawtime);
strftime(buffer, 80, "%H:%M", info);
printf("Sunrise:\t%s\n", buffer);
rawtime = (set - 2440587.5) * 86400;
info = localtime(&rawtime);
strftime(buffer, 80, "%H:%M", info);
printf("Sunset:\t\t%s\n", buffer);
return 0;
}
# reference: https://en.wikipedia.org/wiki/Sunrise_equation#Complete_calculation_on_Earth
require 'date'
def rad(deg)
deg * Math::PI / 180.0
end
def deg(rad)
rad * 180.0 / Math::PI
end
lat = ARGV[0] || 43.720480
lon = ARGV[1] || -79.713230
# Calculate current Julian day
n = Date.today.jd - 2_451_545.0 + 0.0008
# Mean solar noon
J_star = n - (lon / 360.0)
# Solar mean anomaly
M = (357.5291 + 0.98560028 * J_star) % 360
# Equation of the center
C = (1.9148 * Math.sin(rad(M))) + (0.0200 * Math.sin(rad(2 * M))) +( 0.0003 * Math.sin(rad(3 * M)))
# Ecliptic longitude
lam = (M + C + 180 + 102.9372) % 360
# Solar transit
J_transit = 2_451_545.0 + J_star + (0.0053 * Math.sin(rad(M))) - (0.0069 * Math.sin(rad(2 * lam)))
# Declination of the Sun
delta = Math.asin(Math.sin(rad(lam)) * Math.sin(rad(23.44)))
# Hour angle
t = Math.sin(rad(-0.83)) - Math.sin(rad(lat)) * Math.sin(delta)
b = Math.cos(rad(lat)) * Math.cos(delta)
omega = deg(Math.acos(t / b))
# Calculate sunrise and sunset
rise = J_transit - (omega / 360)
set = J_transit + (omega / 360)
offset = ARGV[2] || '-4'
puts "Date: #{Date.today}"
puts "Sunrise: #{DateTime.jd(rise).new_offset(offset).yield_self {|i| i + 0.5}.strftime '%H:%M'}"
puts "Sunset: #{DateTime.jd(set).new_offset(offset).yield_self {|i| i + 0.5}.strftime '%H:%M'}"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment