Created
June 30, 2023 15:50
-
-
Save obstruse/499a6eb1ab7b28cc7bf45815dc18a376 to your computer and use it in GitHub Desktop.
Sunrise Calculation
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
/* Copyright (GPL) 2004 Mike Chirico mchirico@comcast.net | |
Updated: Sun Nov 28 15:15:05 EST 2004 | |
Program adapted by Mike Chirico mchirico@comcast.net | |
Reference: | |
http://prdownloads.sourceforge.net/souptonuts/working_with_time.tar.gz?download | |
http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html | |
Compile: | |
gcc -o sunrise -Wall -W -O2 -s -pipe -lm sunrise.c | |
or for debug output | |
gcc -o sunrise -DDEBUG=1 -Wall -W -O2 -s -pipe -lm sunrise.c | |
This can also go in a batch job to calculate the next | |
20 days as follows: | |
#!/bin/bash | |
lat=39.95 | |
long=75.15 | |
for (( i=0; i <= 20; i++)) | |
do | |
./sunrise `date -d "+$i day" "+%Y %m %d"` $lat $long | |
done | |
*/ | |
/* gcc -DDEBUG=1 .. */ | |
#ifndef DEBUG | |
#define DEBUG 0 | |
#endif | |
#include <sys/time.h> | |
#include <time.h> | |
#include <stdlib.h> | |
#include <stdio.h> | |
#include <math.h> | |
#include <string.h> | |
#include <libgen.h> | |
double calcSunEqOfCenter(double t); | |
/* Convert degree angle to radians */ | |
double degToRad(double angleDeg) | |
{ | |
return (M_PI * angleDeg / 180.0); | |
} | |
double radToDeg(double angleRad) | |
{ | |
return (180.0 * angleRad / M_PI); | |
} | |
double calcMeanObliquityOfEcliptic(double t) | |
{ | |
double seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * (0.001813))); | |
double e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0; | |
return e0; // in degrees | |
} | |
double calcGeomMeanLongSun(double t) | |
{ | |
double L = 280.46646 + t * (36000.76983 + 0.0003032 * t); | |
while ((int)L > 360) | |
{ | |
L -= 360.0; | |
} | |
while (L < 0) | |
{ | |
L += 360.0; | |
} | |
return L; // in degrees | |
} | |
double calcObliquityCorrection(double t) | |
{ | |
double e0 = calcMeanObliquityOfEcliptic(t); | |
double omega = 125.04 - 1934.136 * t; | |
double e = e0 + 0.00256 * cos(degToRad(omega)); | |
return e; // in degrees | |
} | |
double calcEccentricityEarthOrbit(double t) | |
{ | |
double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t); | |
return e; // unitless | |
} | |
double calcGeomMeanAnomalySun(double t) | |
{ | |
double M = 357.52911 + t * (35999.05029 - 0.0001537 * t); | |
return M; // in degrees | |
} | |
double calcEquationOfTime(double t) | |
{ | |
double epsilon = calcObliquityCorrection(t); | |
double l0 = calcGeomMeanLongSun(t); | |
double e = calcEccentricityEarthOrbit(t); | |
double m = calcGeomMeanAnomalySun(t); | |
double y = tan(degToRad(epsilon) / 2.0); | |
y *= y; | |
double sin2l0 = sin(2.0 * degToRad(l0)); | |
double sinm = sin(degToRad(m)); | |
double cos2l0 = cos(2.0 * degToRad(l0)); | |
double sin4l0 = sin(4.0 * degToRad(l0)); | |
double sin2m = sin(2.0 * degToRad(m)); | |
double Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m; | |
return radToDeg(Etime) * 4.0; // in minutes of time | |
} | |
double calcTimeJulianCent(double jd) | |
{ | |
double T = (jd - 2451545.0) / 36525.0; | |
return T; | |
} | |
double calcSunTrueLong(double t) | |
{ | |
double l0 = calcGeomMeanLongSun(t); | |
double c = calcSunEqOfCenter(t); | |
double O = l0 + c; | |
return O; // in degrees | |
} | |
double calcSunApparentLong(double t) | |
{ | |
double o = calcSunTrueLong(t); | |
double omega = 125.04 - 1934.136 * t; | |
double lambda = o - 0.00569 - 0.00478 * sin(degToRad(omega)); | |
return lambda; // in degrees | |
} | |
double calcSunDeclination(double t) | |
{ | |
double e = calcObliquityCorrection(t); | |
double lambda = calcSunApparentLong(t); | |
double sint = sin(degToRad(e)) * sin(degToRad(lambda)); | |
double theta = radToDeg(asin(sint)); | |
return theta; // in degrees | |
} | |
double calcHourAngleSunrise(double lat, double solarDec) | |
{ | |
double latRad = degToRad(lat); | |
double sdRad = degToRad(solarDec); | |
double HA = (acos(cos(degToRad(90.833)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad))); | |
return HA; // in radians | |
} | |
double calcHourAngleSunset(double lat, double solarDec) | |
{ | |
double latRad = degToRad(lat); | |
double sdRad = degToRad(solarDec); | |
double HA = (acos(cos(degToRad(90.833)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad))); | |
return -HA; // in radians | |
} | |
double calcJD(int year, int month, int day) | |
{ | |
if (month <= 2) | |
{ | |
year -= 1; | |
month += 12; | |
} | |
int A = floor(year / 100); | |
int B = 2 - A + floor(A / 4); | |
double JD = floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1)) + day + B - 1524.5; | |
return JD; | |
} | |
double calcJDFromJulianCent(double t) | |
{ | |
double JD = t * 36525.0 + 2451545.0; | |
return JD; | |
} | |
double calcSunEqOfCenter(double t) | |
{ | |
double m = calcGeomMeanAnomalySun(t); | |
double mrad = degToRad(m); | |
double sinm = sin(mrad); | |
double sin2m = sin(mrad + mrad); | |
double sin3m = sin(mrad + mrad + mrad); | |
double C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289; | |
return C; // in degrees | |
} | |
double calcSunriseUTC(double JD, double latitude, double longitude) | |
{ | |
double t = calcTimeJulianCent(JD); | |
// *** First pass to approximate sunrise | |
double eqTime = calcEquationOfTime(t); | |
double solarDec = calcSunDeclination(t); | |
double hourAngle = calcHourAngleSunrise(latitude, solarDec); | |
double delta = longitude - radToDeg(hourAngle); | |
double timeDiff = 4 * delta; // in minutes of time | |
double timeUTC = 720 + timeDiff - eqTime; // in minutes | |
double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0); | |
eqTime = calcEquationOfTime(newt); | |
solarDec = calcSunDeclination(newt); | |
hourAngle = calcHourAngleSunrise(latitude, solarDec); | |
delta = longitude - radToDeg(hourAngle); | |
timeDiff = 4 * delta; | |
timeUTC = 720 + timeDiff - eqTime; // in minutes | |
return timeUTC; | |
} | |
double calcSunsetUTC(double JD, double latitude, double longitude) | |
{ | |
double t = calcTimeJulianCent(JD); | |
// *** First pass to approximate sunset | |
double eqTime = calcEquationOfTime(t); | |
double solarDec = calcSunDeclination(t); | |
double hourAngle = calcHourAngleSunset(latitude, solarDec); | |
double delta = longitude - radToDeg(hourAngle); | |
double timeDiff = 4 * delta; // in minutes of time | |
double timeUTC = 720 + timeDiff - eqTime; // in minutes | |
double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0); | |
eqTime = calcEquationOfTime(newt); | |
solarDec = calcSunDeclination(newt); | |
hourAngle = calcHourAngleSunset(latitude, solarDec); | |
delta = longitude - radToDeg(hourAngle); | |
timeDiff = 4 * delta; | |
timeUTC = 720 + timeDiff - eqTime; // in minutes | |
// printf("************ eqTime = %f \nsolarDec = %f \ntimeUTC = %f\n\n",eqTime,solarDec,timeUTC); | |
return timeUTC; | |
} | |
int main(int argc, char **argv) | |
{ | |
struct tm *ptm = NULL; | |
struct tm tm; | |
time_t seconds; // temp | |
time_t midSeconds; // local midnight | |
time_t nowSeconds; // current time | |
time_t sunSeconds; // calculated sun rise time | |
time_t setSeconds; // calculated sun set time | |
int delta = 7; // hours from here to UTC | |
int year, month, day; | |
int isdst = 1; // uses DST. Should always be "1" | |
double JD; | |
double latitude = 33.8358; // convert to just degrees. No min/sec | |
double longitude = 118.3406; | |
nowSeconds = time(0); | |
if (argc == 6) | |
{ | |
year = atoi(argv[1]); | |
month = atoi(argv[2]); | |
day = atoi(argv[3]); | |
latitude = atof(argv[4]); | |
longitude = atof(argv[5]); | |
} | |
else | |
{ | |
ptm = localtime(&nowSeconds); | |
year = ptm->tm_year + 1900; | |
month = ptm->tm_mon + 1; | |
day = ptm->tm_mday; | |
} | |
// kludge: for some reason the calcs were a day behind published times | |
// adding 1 day to JD | |
// JD = calcJD(year,month,day+1); | |
JD = calcJD(year, month, day); | |
tm.tm_year = year - 1900; | |
tm.tm_mon = month - 1; /* Jan = 0, Feb = 1,.. Dec = 11 */ | |
tm.tm_mday = day; | |
tm.tm_hour = 0; | |
tm.tm_min = 0; | |
tm.tm_sec = 0; | |
// tm.tm_isdst=-1; | |
tm.tm_isdst = isdst; // seems that isdst means that DST is used when necessary...so should be "1" | |
midSeconds = mktime(&tm); | |
sunSeconds = midSeconds + calcSunriseUTC(JD, latitude, longitude) * 60; | |
sunSeconds -= delta * 3600; | |
setSeconds = midSeconds + calcSunsetUTC(JD, latitude, longitude) * 60; | |
setSeconds -= delta * 3600; | |
if (argc == 6) | |
{ | |
ptm = localtime(&sunSeconds); | |
printf("%2d-%2d-%d Rise: %02d:%02d:%02d ", ptm->tm_mon + 1, ptm->tm_mday, ptm->tm_year + 1900, ptm->tm_hour, ptm->tm_min, ptm->tm_sec); | |
ptm = localtime(&setSeconds); | |
printf("Set: %02d:%02d:%02d ", ptm->tm_hour, ptm->tm_min, ptm->tm_sec); | |
seconds = nowSeconds - sunSeconds; | |
printf("Time: %02d:%02d:%02d\n", (int)((seconds / 60 / 60) % 24), (int)((seconds / 60) % 60), (int)(seconds % 60)); | |
} | |
else | |
{ | |
if (!strcmp(basename(argv[0]), "sunset")) | |
{ | |
ptm = localtime(&setSeconds); | |
printf("%02d%02d\n", ptm->tm_hour, ptm->tm_min); | |
} | |
else | |
{ | |
ptm = localtime(&sunSeconds); | |
printf("%02d%02d\n", ptm->tm_hour, ptm->tm_min); | |
} | |
} | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment