在我的C应用程序中，我想计算给定日期，纬度和经度的日出/日落时间。我一直在网上搜索，但我找不到工作样本。
我试图实现这个样本：http://souptonuts.sourceforge.net/code/sunrise.c.html
但是这个样本没有正常工作。
是否有一个简单的C源代码或方法，我可以在我的应用程序中轻松实现？
编辑： 我在这个link上实现了代码，但它给了我错误的日落/日出值。我也尝试了Saul的链接here，但它也给了我错误的结果。 我有41N，28E的位置。当我尝试这些代码时，两个样本都说日出值大约是10:13而日落是23:24。但正确的值是06：06,20：13。 我无法理解这个问题。
示例代码似乎在VC ++ 2010中有一些小的更改：
#include <sys/time.h>线。#define _USE_MATH_DEFINES，以便定义M_PI。%T调用中的两个strftime()更改为%X。既然您有一个工作样本，您可以调试工作版本和您的版本，以查看计算开始不同的地方并缩小问题范围。无论是单步执行程序还是自由使用临时printf()调用就像样本一样。
如果您需要特定帮助，则必须发布您的代码（指向整个文件的链接或您需要帮助的特定代码段）。
这是适应性的，非常准确。您将获得所有组件，然后您需要计算的是Zenith到Horizon角度的反余弦。是的，有更简单的方法，但你不是基本上想跟踪太阳？它有朝一日会派上用场。 http://www.nrel.gov/midc/spa/
//a practical approximation for a microcontroller using lookup table
//seems to only use a small fraction of the resources required by the trigonometric functions
//ignores daylight saving: idea is for the device to approximate and trigger actual sunrise, sunset event as opposed to actual (politically correct local hhmm)
//using globals gps.Lat gps.LatDir(0:S 1:N) gps.Day (day of month) gps.Mth
//needs solar noon offset,latitude,dayOfMonth,Month
//LocalNoon-SolarNoon (may be off +/- up to ?? 2 hours) depending on local timezone and longitude (? typical discrepancy about 20min)
void SunriseSunsetCalculations(u8 SolarNoonOffsetInDeciHours,u8 SolarNoonOffsetDirection){
if(SolarNoonOffsetInDeciHours>20){SolarNoonOffsetInDeciHours=0;}//limit offest to 2 hours: discard offset on error
//--------references:
//https://www.orchidculture.com/COD/daylength.html
//http://aa.usno.navy.mil/data/docs/Dur_OneYear.php
//SolarTime is up two two hours off in either direction depending on timezone:
//http://blog.poormansmath.net/how-much-is-time-wrong-around-the-world/
//------------------
//lookUpTable Of daylength in decihours(6min)(fits well in u8) - 15day And 5 DegLat resolution
//SolarNoonOffsetDirection: 0: negative(b4 local noon), 1: +ve (solar noon after local noon)
u8 a[13][12]={//length of day in decihours
// Dec Nov Oct Sep Aug Jul
//Jan Feb Mar Apr May June
{120,120,120,120,120,120,120,120,120,120,120,120}, //lat 0____0
{126,125,124,123,122,121,119,118,116,115,115,114}, //lat 10____1
{129,128,127,125,123,121,119,117,115,113,112,111}, //lat 15____2
{132,131,129,127,124,121,118,115,113,110,109,108}, //lat 20____3
{135,134,131,128,125,122,118,114,111,108,106,105}, //lat 25____4
{139,137,134,130,127,122,117,113,108,105,102,101}, //lat 30____5
{143,141,137,133,128,123,117,111,106,102, 98, 97}, //lat 35____6
{148,145,141,135,130,123,116,109,103, 98, 94, 92}, //lat 40____7
{154,150,145,138,132,124,115,107,100, 93, 88, 86}, //lat 45____8
{161,157,150,142,134,124,114,105, 96, 88, 82, 79}, //lat 50____9
{170,165,156,146,137,125,113,102, 91, 81, 73, 69}, //lat 55___10
{183,176,165,152,140,126,112, 98, 84, 72, 61, 56}, //lat 60___11
{200,185,171,152,134,121,101, 84, 65, 53, 40, 33} //lat 65___12
};
u8 b[]={6,12,17,22,27,32,37,42,47,52,57,62,90}; // latitude limit cutoffs to get index of lookUpTable
u8 lat=gps.Lat/10000000; //lat stored in u32 to 7 decimals resolution (32bit unsigned integer)
u8 i=0; while(b[i]<lat){i++;} //get row index for daylength table
u8 k,ix; //k: 15 day offset; ix: column index for daylength table
k=gps.Day/15;if(k){k=1;}//which half of the month (avoid k=2 eg:31/15)
if(!gps.LatDir){ //0:southern latitudes
if(gps.Mth<7){
ix=(gps.Mth-1)*2+k; //2 fields per month (k to select)
}else{ //beyond june, use row in reverse
ix=11-(gps.Mth-7)*2-k; //2 fields per month (k to select)
}
}else{ //1:northern latitudes
if(gps.Mth<7){ //with first six month read rows in reverse
ix=11-(gps.Mth-1)*2-k; //2 fields per month (k to select)
}else{ //beyond june, use same row forward
ix=(gps.Mth-7)*2+k; //2 fields per month (k to select)
}
}
//debug only:...dcI("\r\ni", i ,Blue,Red,1);dcI("ix", ix ,Blue,Red,1);dcI("a[i][ix]", a[i][ix] ,Blue,Red,1);
u8 h=a[i][ix]/2; //HalfTheDayLightLength in deciHours
u8 j[]={-1,1}; //multiplier: 0:(-) 1:(+)
u8 sn=120+SolarNoonOffsetInDeciHours*j[SolarNoonOffsetDirection]; //Solar noon
u8 srdh=sn-h; //sunrise in deciHours = solarNoon minus HalfTheDayLightLength
u8 ssdh=sn+h; //sunset in deciHours = solarNoon plus HalfTheDayLightLength
//debug only:...dcI("\r\nSunRiseDeciHours", srdh ,Blue,Red,1);dcI("SunSetDeciHours", ssdh ,Blue,Red,1);
gps.HmSunRise=deciHourTohhmm(srdh);
gps.HmSunSet =deciHourTohhmm(ssdh);
}
u16 deciHourTohhmm(u8 dh){ //return unsigned integer from 0 to 2400
u16 h=(dh/10)*100; //hours: hh00
u8 r= dh%10; //fraction hour remainder to be converted to minutes
u8 m= 6*r; //60*r/10
return(h+m);
}
/*
Example Output: (!!! solarNoonOffset kept at 0(ignored) for the below example)
:_(08474300)___:_(A)___:_(381234567)___:_(S)___:_(1431234567)___:_(E)___
:_(GPS OK)___:_(Sat)___:_(12)___:_(Aug)___:_(2017)___hhmm:_(847)___
//...........
i:_(7)___ix:_(9)___a[i][ix]:_(98)___
SunRiseDeciHours:_(71)___SunSetDeciHours:_(169)___
HmSunRise:_(706)___
HmSunSet:_(1654)___
//............same as above but change LatDir to 1 (northern hemisphere)
i:_(7)___ix:_(2)___a[i][ix]:_(141)___
SunRiseDeciHours:_(50)___SunSetDeciHours:_(190)___
HmSunRise:_(500)___
HmSunSet:_(1900)___
..........ignore dcI(...) it's just a custom function printing through the serial port
*/
给出日期，纬度和经度，计算sunrise / sunset time的十个简单步骤
5A。计算太阳的正确提升
RA = atan(0.91764 * tan(L))
NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
5B。右提升值需要与L在同一象限
Lquadrant = (floor( L/90)) * 90
RAquadrant = (floor(RA/90)) * 90
RA = RA + (Lquadrant - RAquadrant)
5C。正确的提升值需要转换成小时
RA = RA / 15
7A。计算太阳的当地小时角
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。完成计算H并转换为小时
if if rising time is desired:
H = 360 - acos(cosH)
if setting time is desired:
H = acos(cosH)
H = H / 15
有一个c解决方案，其中包含日出/设置的目标c包装：https://github.com/berkley/ObjectiveCUtil
使用这个guide（由@BenjaminMonate首次发布，我假设是@Geetha使用的那个），我构建了这个C函数，它似乎正常工作。
#include <math.h>
#define PI 3.1415926
#define ZENITH -.83
float calculateSunrise(int year,int month,int day,float lat, float lng,int localOffset, int daylightSavings) {
/*
localOffset will be <0 for western hemisphere and >0 for eastern hemisphere
daylightSavings should be 1 if it is in effect during the summer otherwise it should be 0
*/
//1. first calculate the day of the year
float N1 = floor(275 * month / 9);
float N2 = floor((month + 9) / 12);
float N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
float N = N1 - (N2 * N3) + day - 30;
//2. convert the longitude to hour value and calculate an approximate time
float lngHour = lng / 15.0;
float t = N + ((6 - lngHour) / 24); //if rising time is desired:
//float t = N + ((18 - lngHour) / 24) //if setting time is desired:
//3. calculate the Sun's mean anomaly
float M = (0.9856 * t) - 3.289;
//4. calculate the Sun's true longitude
float L = fmod(M + (1.916 * sin((PI/180)*M)) + (0.020 * sin(2 *(PI/180) * M)) + 282.634,360.0);
//5a. calculate the Sun's right ascension
float RA = fmod(180/PI*atan(0.91764 * tan((PI/180)*L)),360.0);
//5b. right ascension value needs to be in the same quadrant as L
float Lquadrant = floor( L/90) * 90;
float 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
float sinDec = 0.39782 * sin((PI/180)*L);
float cosDec = cos(asin(sinDec));
//7a. calculate the Sun's local hour angle
float cosH = (sin((PI/180)*ZENITH) - (sinDec * sin((PI/180)*lat))) / (cosDec * cos((PI/180)*lat));
/*
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
float H = 360 - (180/PI)*acos(cosH); // if if rising time is desired:
//float H = acos(cosH) // if setting time is desired:
H = H / 15;
//8. calculate local mean time of rising/setting
float T = H + RA - (0.06571 * t) - 6.622;
//9. adjust back to UTC
float UT = fmod(T - lngHour,24.0);
//10. convert UT value to local time zone of latitude/longitude
return UT + localOffset + daylightSavings;
}
void printSunrise() {
float localT = calculateSunrise(/*args*/);
double hours;
float minutes = modf(localT,&hours)*60;
printf("%.0f:%.0f",hours,minutes);
}
scottmrogowski提供的代码很有用，但有两个问题
干杯马雷克
也许试试这段代码。它经过测试和运作。希望你喜欢...
#include "stdafx.h"
#include <iostream>
#include <math.h>
#include <time.h>
using namespace std;
//STANDARD CONSTANTS
double pi = 3.1415926535; // Pi
double solarConst = 1367; // solar constant W.m-2
// Function to convert radian to hours
double RadToHours (double tmp)
{
//double pi = 3.1415926535; // Pi
return (tmp * 12 / pi);
}
// Function to convert hours to radians
double HoursToRads (double tmp)
{
//double pi = 3.1415926535; // Pi
return (tmp * pi / 12);
}
// Function to calculate the angle of the day
double AngleOfDay (int day, // number of the day
int month, // number of the month
int year // year
)
{ // local vars
int i, leap;
int numOfDays = 0; // number of Day 13 Nov=317
int numOfDaysofMonths[12] = {0,31,28,31,30,31,30,31,31,30,31,30}; // Number of days per month
int AllYearDays; // Total number of days in a year 365 or 366
double DayAngle; // angle of the day (radian)
//double pi = 3.1415926535; // Pi
// leap year ??
leap = 0;
if ((year % 400)==0)
{ AllYearDays = 366;
leap = 1;
}
else if ((year % 100)==0) AllYearDays = 365;
else if ((year % 4)==0)
{ AllYearDays = 366;
leap = 1;
}
else AllYearDays = 365;
// calculate number of day
for (i=0;i<month;i++) numOfDays += numOfDaysofMonths[i];
if ( (month > 2) && leap) numOfDays++;
numOfDays += day;
// calculate angle of day
DayAngle = (2*pi*(numOfDays-1)) / AllYearDays;
return DayAngle;
}
// Function to calculate declination - in radian
double Declination (double DayAngle // angle day in radian
)
{
double SolarDeclination;
// Solar declination (radian)
SolarDeclination = 0.006918
- 0.399912 * cos (DayAngle)
+ 0.070257 * sin (DayAngle)
- 0.006758 * cos (2*DayAngle)
+ 0.000907 * sin (2*DayAngle)
- 0.002697 * cos (3*DayAngle)
+ 0.00148 * sin (3*DayAngle);
return SolarDeclination;
}
// Function to calculate Equation of time ( et = TSV - TU )
double EqOfTime (double DayAngle // angle day (radian)
)
{
double et;
// Equation of time (radian)
et = 0.000075
+ 0.001868 * cos (DayAngle)
- 0.032077 * sin (DayAngle)
- 0.014615 * cos (2*DayAngle)
- 0.04089 * sin (2*DayAngle);
// Equation of time in hours
et = RadToHours(et);
return et;
}
// Calculation of the duration of the day in radian
double DayDurationRadian (double _declination, // _declination in radian
double lat // latitude in radian
)
{
double dayDurationj;
dayDurationj = 2 * acos( -tan(lat) * tan(_declination) );
return dayDurationj;
}
// Function to calculate Day duration in Hours
double DayDuratInHours (double _declination // _declination in radian
, double lat // latitude in radian
)
{
double dayDurationj;
dayDurationj = DayDurationRadian(_declination, lat);
dayDurationj = RadToHours(dayDurationj);
return dayDurationj;
}
// Function to calculate the times TSV-UTC
double Tsv_Tu (double rlong // longitude en radian positive a l est.
,double eqOfTime // Equation of times en heure
)
{
double diffUTC_TSV; double pi = 3.1415926535; // Pi
// diffUTC_TSV Solar time as a function of longitude and the eqation of time
diffUTC_TSV = rlong * (12 / pi) + eqOfTime;
// difference with local time
return diffUTC_TSV;
}
// Calculations of the orbital excentricity
double Excentricity(int day,
int month,
int year)
{
double dayAngleRad, E0;
// calculate the angle of day in radian
dayAngleRad = AngleOfDay(day, month, year);
// calculate the excentricity
E0 = 1.000110 + 0.034221 * cos(dayAngleRad)
+ 0.001280 * sin(dayAngleRad)
+0.000719 * cos(2*dayAngleRad)
+0.000077 * sin(2*dayAngleRad);
return E0;
}
// Calculate the theoretical energy flux for the day radiation
double TheoreticRadiation(int day, int month, int year,
double lat // Latitude in radian !
)
{
double RGth; // Theoretical radiation
double decli; // Declination
double E0;
double sunriseHourAngle; // Hour angle of sunset
// Calculation of the declination in radian
decli = Declination (AngleOfDay(day, month, year));
// Calcuate excentricity
E0 = Excentricity(day, month, year);
// Calculate hour angle in radian
sunriseHourAngle = DayDurationRadian(decli, lat) / 2;
// Calculate Theoretical radiation en W.m-2
RGth = solarConst * E0 * (cos(decli)*cos(lat)*sin(sunriseHourAngle)/sunriseHourAngle + sin(decli)*sin(lat));
return RGth;
}
// Function to calculate decimal hour of sunrise: result in local hour
double CalclulateSunriseLocalTime(int day,
int month,
int year,
double rlong,
double rlat)
{
// local variables
int h1, h2;
time_t hour_machine;
struct tm *local_hour, *gmt_hour;
double result;
// Calculate the angle of the day
double DayAngle = AngleOfDay(day, month, year);
// Declination
double SolarDeclination = Declination(DayAngle);
// Equation of times
double eth = EqOfTime(DayAngle);
// True solar time
double diffUTC_TSV = Tsv_Tu(rlong,eth);
// Day duration
double dayDurationj = DayDuratInHours(SolarDeclination,rlat);
// local time adjust
time( &hour_machine ); // Get time as long integer.
gmt_hour = gmtime( &hour_machine );
h1 = gmt_hour->tm_hour;
local_hour = localtime( &hour_machine ); // local time.
h2 = local_hour->tm_hour;
// final result
result = 12 - fabs(dayDurationj / 2) - diffUTC_TSV + h2-h1;
return result;
}
// Function to calculate decimal hour of sunset: result in local hour
double CalculateSunsetLocalTime(int day,
int month,
int year,
double rlong,
double rlat)
{
// local variables
int h1, h2;
time_t hour_machine;
struct tm *local_hour, *gmt_hour;
double result;
// Calculate the angle of the day
double DayAngle = AngleOfDay(day, month, year);
// Declination
double SolarDeclination = Declination(DayAngle);
// Equation of times
double eth = EqOfTime(DayAngle);
// True solar time
double diffUTC_TSV = Tsv_Tu(rlong,eth);
// Day duration
double dayDurationj = DayDuratInHours(SolarDeclination,rlat);
// local time adjust
time( &hour_machine ); // Get time as long integer.
gmt_hour = gmtime( &hour_machine );
h1 = gmt_hour->tm_hour;
local_hour = localtime( &hour_machine ); // local time.
h2 = local_hour->tm_hour;
// resultat
result = 12 + fabs(dayDurationj / 2) - diffUTC_TSV + h2-h1;
return result;
}
// Function to calculate decimal hour of sunrise: result universal time
double CalculateSunriseUniversalTime(int day,
int month,
int year,
double rlong,
double rlat)
{
double result;
// Calculate the angle of the day
double DayAngle = AngleOfDay(day, month, year);
// Declination
double SolarDeclination = Declination(DayAngle);
// Equation of times
double eth = EqOfTime(DayAngle);
// True solar time
double diffUTC_TSV = Tsv_Tu(rlong,eth);
// Day duration
double dayDurationj = DayDuratInHours(SolarDeclination,rlat);
// resultat
result = 12 - fabs(dayDurationj / 2) - diffUTC_TSV;
return result;
}
// Function to calculate decimal hour of sunset: result in universal time
double CalculateSunsetUniversalTime(int day,
int month,
int year,
double rlong,
double rlat)
{
double result;
// Calculate the angle of the day
double DayAngle = AngleOfDay(day, month, year);
// Declination
double SolarDeclination = Declination(DayAngle);
// Equation of times
double eth = EqOfTime(DayAngle);
// True solar time
double diffUTC_TSV = Tsv_Tu(rlong,eth);
// Day duration
double dayDurationj = DayDuratInHours(SolarDeclination,rlat);
// resultat
result = 12 + fabs(dayDurationj / 2) - diffUTC_TSV;
return result;
}
// Function to calculate the height of the sun in radians the day to day j and hour TU
double SolarHeight (int tu, // universal times (0,1,2,.....,23)
int day,
int month,
int year,
double lat, // latitude in radian
double rlong // longitude in radian
)
{
// local variables
double pi = 3.1415926535; // Pi
double result, tsvh;
// angle of the day
double DayAngle = AngleOfDay(day, month, year);
// _declination
double decli = Declination(DayAngle);
// eq of time
double eq = EqOfTime(DayAngle);
// calculate the tsvh with rlong positiv for the east and negative for the west
tsvh = tu + rlong*180/(15*pi) + eq;
// hour angle per hour
double ah = acos( -cos((pi/12)*tsvh) );
// final result
result = asin( sin(lat)*sin(decli) + cos(lat)*cos(decli)*cos(ah) );
return result;
}
///////////EXTRA FUNCTIONS/////////////////////////////
//Julian day conversion for days calculations
//Explanation for this sick formula...for the curious guys...
//http://www.cs.utsa.edu/~cs1063/projects/Spring2011/Project1/jdn-explanation.html
int julian(int year, int month, int day) {
int a = (14 - month) / 12;
int y = year + 4800 - a;
int m = month + 12 * a - 3;
if (year > 1582 || (year == 1582 && month > 10) || (year == 1582 && month == 10 && day >= 15))
return day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
else
return day + (153 * m + 2) / 5 + 365 * y + y / 4 - 32083;
}
int _tmain(int argc, _TCHAR* argv[])
{
int day = 14;
int month = 11;
int year = 2013;
double lat = 39.38;
double lon = 22.75;
double rlat = 39.38 * pi/180;
double rlong = 22.75 * pi/180;
double _AngleOfDay = AngleOfDay ( day , month , year );
cout << "Angle of day: " << _AngleOfDay << "\n";
double _Declinaison = Declination (_AngleOfDay);
cout << "Declination (Delta): " << _Declinaison << "\n";
double _EqOfTime = EqOfTime (_AngleOfDay);
cout << "Declination (Delta): " << _EqOfTime << "\n";
double _DayDuratInHours = DayDuratInHours (_Declinaison, rlat);
cout << "Day duration: " << _DayDuratInHours << "\n";
double _Excentricity = Excentricity(day, month, year);
cout << "Excentricity: " << _Excentricity << "\n";
double _TheoreticRadiation = TheoreticRadiation(day, month, year, rlat);
cout << "Theoretical radiation: " << _TheoreticRadiation << "\n";
double _CalclulateSunriseLocalTime = CalclulateSunriseLocalTime
(day, month, year, rlong, rlat);
cout << "Sunrise Local Time: " << _CalclulateSunriseLocalTime << "\n";
double _CalculateSunsetLocalTime = CalculateSunsetLocalTime
(day, month, year, rlong, rlat);
cout << "Sunrise Local Time: " << _CalculateSunsetLocalTime << "\n";
return 0;
}
纯粹的C实施显然是稀缺的，但如果你愿意从C++或C#移植，那么有几个选择：
我将Tomoyose的C代码移植到C＃。虽然需要手动补偿，但由于与在线资源的差异大约为4分钟，因此工作正常。这可能是由于日落或上升的构成存在差异，但正式情况下应该是当太阳从它下方触及地平线时。我希望细节隐藏在某个地方的常数参数中 - 超出我的范围，尝试正确的修复。如果其他人管理它我想知道;-)
它使用NodaTime，这是Jon Skeet对Java的JodaTime的实现。 NodaTime在nuGet中可用。
using System;
using NodaTime;
namespace YourNamespaceHere
{
public static class MySunset
{
// Based on Tomoyose's
// http://stackoverflow.com/questions/7064531/sunrise-sunset-times-in-c
private static DateTimeZone mUtcZone = DateTimeZoneProviders.Tzdb["Etc/UTC"];
private static int mSecondsInDay = 24 * 60 * 60;
// Convert radian to hours
private static double RadiansToHours(double radians)
{
return (radians * 12.0 / Math.PI);
}
// Convert hours to radians
private static double HoursToRadians(double hours)
{
return (hours * Math.PI / 12.0);
}
// Calculate the angle of the day
private static double AngleOfDay(int year, int month, int day)
{
DateTime date = new DateTime(year, month, day);
int daysInYear = DateTime.IsLeapYear(year) ? 366 : 365;
// Return angle of day in radians
return (2.0 * Math.PI * ((double)date.DayOfYear - 1.0)) / (double)daysInYear;
}
// Calculate declination in radians
private static double SolarDeclination(double angleOfDayRadians)
{
// Return solar declination in radians
return 0.006918
- 0.399912 * Math.Cos(angleOfDayRadians)
+ 0.070257 * Math.Sin(angleOfDayRadians)
- 0.006758 * Math.Cos(2.0 * angleOfDayRadians)
+ 0.000907 * Math.Sin(2.0 * angleOfDayRadians)
- 0.002697 * Math.Cos(3.0 * angleOfDayRadians)
+ 0.00148 * Math.Sin(3.0 * angleOfDayRadians);
}
// Calculate Equation of time ( eot = TSV - TU )
private static double EquationOfTime(double angleOfDayRadians)
{
// Equation of time (radians)
double et = 0.000075
+ 0.001868 * Math.Cos(angleOfDayRadians)
- 0.032077 * Math.Sin(angleOfDayRadians)
- 0.014615 * Math.Cos(2.0 * angleOfDayRadians)
- 0.04089 * Math.Sin(2.0 * angleOfDayRadians);
// Return equation-of-time in hours
return RadiansToHours(et);
}
// Calculate the duration of the day in radians
private static double DayDurationRadians(double declinationRadians, double latitudeRadians)
{
return 2.0 * Math.Acos(-Math.Tan(latitudeRadians) * Math.Tan(declinationRadians));
}
// Calculate day duration in hours
private static double DayDurationHours(double declinationRadians, double latitudeRadians)
{
return RadiansToHours(
DayDurationRadians(declinationRadians, latitudeRadians)
);
}
// Calculate the times TSV-UTC
private static double Tsv_Tu(double longitudeRadians, double equationOfTime)
{
// Solar time as a function of longitude and the equation of time
return longitudeRadians * (12.0 / Math.PI) + equationOfTime;
}
private static void GetDayParameters(int year, int month, int day, double latitude, double longitude,
out double dayDuration, out double diffUTC_TSV)
{
double latitudeRadians = latitude * Math.PI / 180.0;
double longitudeRadians = longitude * Math.PI / 180.0;
// Calculate the angle of the day
double dayAngle = AngleOfDay(year, month, day);
// Declination
double solarDeclination = SolarDeclination(dayAngle);
// Equation of times
double equationOfTime = EquationOfTime(dayAngle);
// True solar time
diffUTC_TSV = Tsv_Tu(longitudeRadians, equationOfTime);
// Day duration
dayDuration = DayDurationHours(solarDeclination, latitudeRadians);
}
// Calculate decimal UTC hour of sunrise.
private static double CalculateSunriseUTC(int year, int month, int day, double latitude, double longitude)
{
double dayDuration;
double diffUTC_TSV;
GetDayParameters(year, month, day, latitude, longitude, out dayDuration, out diffUTC_TSV);
return 12.0 - Math.Abs(dayDuration / 2.0) - diffUTC_TSV;
}
// Calculate decimal UTC hour of sunset.
private static double CalculateSunsetUTC(int year, int month, int day, double latitude, double longitude)
{
double dayDuration;
double diffUTC_TSV;
GetDayParameters(year, month, day, latitude, longitude, out dayDuration, out diffUTC_TSV);
return 12.0 + Math.Abs(dayDuration / 2.0) - diffUTC_TSV;
}
// Public methods to return sun rise and set times.
public static Tuple<ZonedDateTime, ZonedDateTime> GetSunRiseSet(LocalDate dateAtLocation, DateTimeZone locationZone, double latitude, double longitude)
{
// latitude-longitude must lie within zone of locationZone
// Get UTC rise and set
double dayDuration;
double diffUTC_TSV;
GetDayParameters(dateAtLocation.Year, dateAtLocation.Month, dateAtLocation.Day, latitude, longitude, out dayDuration, out diffUTC_TSV);
double sunriseUtcDecimal = 12.0 - Math.Abs(dayDuration / 2.0) - diffUTC_TSV;
double sunsetUtcDecimal = 12.0 + Math.Abs(dayDuration / 2.0) - diffUTC_TSV;
// Convert decimal UTC to UTC dates
// If a UTC time is negative then it means the date before in the UTC timezone.
// So if negative need to minus 1 day from date and set time as 24 - decimal_time.
LocalDateTime utcRiseLocal;
LocalDateTime utcSetLocal;
if (sunriseUtcDecimal < 0)
{
LocalDate utcDateAdjusted = dateAtLocation.PlusDays(-1);
// Normalize() is important here; otherwise only have access to seconds.
Period utcTimeAdjusted = Period.FromSeconds((long)((24.0 + sunriseUtcDecimal) / 24.0 * mSecondsInDay)).Normalize(); // + a negative
utcRiseLocal = new LocalDateTime(utcDateAdjusted.Year, utcDateAdjusted.Month, utcDateAdjusted.Day,
(int)utcTimeAdjusted.Hours, (int)utcTimeAdjusted.Minutes, (int)utcTimeAdjusted.Seconds);
}
else
{
Period utcTime = Period.FromSeconds((long)(sunriseUtcDecimal / 24.0 * mSecondsInDay)).Normalize();
utcRiseLocal = new LocalDateTime(dateAtLocation.Year, dateAtLocation.Month, dateAtLocation.Day,
(int)utcTime.Hours, (int)utcTime.Minutes, (int)utcTime.Seconds);
}
if (sunsetUtcDecimal < 0) // Maybe not possible?
{
LocalDate utcDateAdjusted = dateAtLocation.PlusDays(-1);
Period utcTimeAdjusted = Period.FromSeconds((long)((24.0 + sunsetUtcDecimal) / 24.0 * mSecondsInDay)).Normalize();
utcSetLocal = new LocalDateTime(utcDateAdjusted.Year, utcDateAdjusted.Month, utcDateAdjusted.Day,
(int)utcTimeAdjusted.Hours, (int)utcTimeAdjusted.Minutes, (int)utcTimeAdjusted.Seconds);
}
else
{
Period utcTime = Period.FromSeconds((long)(sunsetUtcDecimal / 24.0 * mSecondsInDay)).Normalize();
utcSetLocal = new LocalDateTime(dateAtLocation.Year, dateAtLocation.Month, dateAtLocation.Day,
(int)utcTime.Hours, (int)utcTime.Minutes, (int)utcTime.Seconds);
}
// FIX: always about 4 minutes later/earlier than other sources
utcRiseLocal = utcRiseLocal.PlusMinutes(-4);
utcSetLocal = utcSetLocal.PlusMinutes(4);
// Get zoned datetime from UTC local datetimes
ZonedDateTime utcRiseZoned = new LocalDateTime(utcRiseLocal.Year, utcRiseLocal.Month, utcRiseLocal.Day, utcRiseLocal.Hour, utcRiseLocal.Minute, utcRiseLocal.Second).InZoneLeniently(mUtcZone);
ZonedDateTime utcSetZoned = new LocalDateTime(utcSetLocal.Year, utcSetLocal.Month, utcSetLocal.Day, utcSetLocal.Hour, utcSetLocal.Minute, utcSetLocal.Second).InZoneLeniently(mUtcZone);
// Return zoned UTC to zoned local
return new Tuple<ZonedDateTime, ZonedDateTime>
(
new ZonedDateTime(utcRiseZoned.ToInstant(), locationZone),
new ZonedDateTime(utcSetZoned.ToInstant(), locationZone)
);
}
public static Tuple<ZonedDateTime, ZonedDateTime> GetSunRiseSet(int year, int month, int day, string locationZone, double latitude, double longitude)
{
return GetSunRiseSet(new LocalDate(year, month, day), DateTimeZoneProviders.Tzdb[locationZone], latitude, longitude);
}
public static Tuple<ZonedDateTime, ZonedDateTime> GetSunRiseSet(ZonedDateTime zonedDateAtLocation, double latitude, double longitude)
{
return GetSunRiseSet(zonedDateAtLocation.LocalDateTime.Date, zonedDateAtLocation.Zone, latitude, longitude);
}
}
}