Last active
March 21, 2023 20:24
-
-
Save AlessandroMinali/f5813341618aa71a429a6f2d6df82056 to your computer and use it in GitHub Desktop.
Calculate Sun Rise and Set times based. Takes two args: <lat> <lon>
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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; | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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