C语言中的日出日落时间

Question

C语言中的日出日落时间

6

在我的C应用程序中，我想计算给定日期、纬度和经度的日出/日落时间。我已经在网上搜索过，但找不到可行的示例。

我尝试实现这个示例：http://souptonuts.sourceforge.net/code/sunrise.c.html 但是这个示例没有正确工作。

是否有一个简单的C源代码或方法，我可以轻松地在我的应用程序中实现？

编辑：我在这个链接link上实现了代码，但是它给了我错误的日落/日出值。我也尝试了Saul的链接here，但是它也给了我错误的结果。我位于41N，28E位置。当我尝试这些代码时，两个示例都说日出时间约为10:13，日落时间为23:24。但正确的值是06:06，20:13。我无法理解问题所在。

- Fer

4

"didn't work correctly" 的意思是指某事物没有按照预期的方式正常运作。 - Oliver Charlesworth

2

“Did not work correctly”是一个糟糕的问题描述。如果您想要帮助让代码正常工作，请提供具体哪些部分出现了问题（假设与代码有关）。否则，本网站不是搜索引擎。 - Mat

13个回答

8

这似乎很容易实现：
http://edwilliams.org/sunrise_sunset_algorithm.htm

- Benjamin Monate

6

使用这个指南（最初由@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

4

scottmrogowski提供的代码非常有用，但是有两个问题

//float H = (180/PI)*acos(cosH) //如果需要设置时间：
float localT=fmod(24 + calculateSunrise(/* args */),24.0); //在printSunrise函数中

谢谢 Marek

- user2018761

2

也许可以尝试这段代码。它已经经过测试并且可行。希望您喜欢它...

#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;
}

- Tomoyose

本地时间变量在午夜时有一个小错误。首先按以下方式计算dif，而不是h2-h1： int dif = h2-h1; if (dif<-12) dif += 24; - Adion

2

这里有一个关于日出/日落的 C 语言解决方案，还带有 Objective-C 封装，链接在这里：https://github.com/berkley/ObjectiveCUtil。

- cberkley

2

我将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
        // https://dev59.com/NlrUa4cB1Zd3GeqPiUYB

        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);
        }
    }
}

- user1350972

2

我最近写了一个库来完成这个任务。它同时提供了C和Python的API接口 - https://github.com/adonmo/daylight

以下是日出时间计算的摘录：

Original Answer翻译成"最初的回答"

time_t Sunclock::sunrise(time_t date) {
  date = date + tz_offset * 60 * 60;
  struct tm *t = gmtime(&date);
  double _time_of_day = time_of_day(date);
  double _julian_day = julian_day(t, _time_of_day, tz_offset);
  double _julian_century = julian_century(_julian_day);
  double _mean_obliq_ecliptic = mean_obliq_ecliptic(_julian_century);
  double _mean_long_sun = mean_long_sun(_julian_century);
  double _mean_anom_sun = mean_anom_sun(_julian_century);
  double _sun_eq_of_centre = sun_eq_of_centre(_mean_anom_sun, _julian_century);
  double _sun_true_long = sun_true_long(_mean_long_sun, _sun_eq_of_centre);
  double _obliq_corr = obliq_corr(_mean_obliq_ecliptic, _julian_century);
  double _sun_app_long = sun_app_long(_sun_true_long, _julian_century);
  double _eccent_earth_orbit = eccent_earth_orbit(_julian_century);
  double _var_y = var_y(_obliq_corr);
  double _eq_of_time =
      eq_of_time(_var_y, _mean_long_sun, _eccent_earth_orbit, _mean_anom_sun);
  double _declination = declination(_obliq_corr, _sun_app_long);
  double _hour_angle_sunrise = hour_angle_sunrise(_declination);

  double noon_decimal_day =
      (720 - 4 * longitude - _eq_of_time + tz_offset * 60) / 1440;
  double decimal_day = noon_decimal_day - _hour_angle_sunrise * 4 / 1440;
  return time_from_decimal_day(date, decimal_day) - tz_offset * 60 * 60;
}

上面的代码应该给出了一个大致的想法，但其他功能的代码以及完整的代码实际上可以从这里阅读：https://github.com/adonmo/daylight/blob/master/source/Sunclock.cpp。最初的回答。

- chaitan94

1

//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
*/

- Denes Borsos

1

纯 C 实现似乎很少，但如果您愿意从 C++ 或 C# 进行移植，那么有一些选择：

C++ 日出/日落时间计算
C# 根据纬度和经度获取日出和日落时间

- Saul

我在第一个链接上运行了你的C++示例，但是项目显示错误的结果。请阅读我上面的编辑。 - Fer

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Geetha · Accepted Answer

按照以下十个简单步骤计算给定日期、纬度和经度的日出/日落时间

首先计算一年中的天数

N1 = floor(275 * month / 9) N2 = floor((month + 9) / 12) N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3)) N = N1 - (N2 * N3) + day - 30
将经度转换为小时值并计算近似时间

lngHour = longitude / 15

如果需要日出时间： t = N + ((6 - lngHour) / 24) 如果需要日落时间： t = N + ((18 - lngHour) / 24)
计算太阳的平近点角

M = (0.9856 * t) - 3.289
计算太阳的真实黄经

L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634 注意：L可能需要通过加/减360来调整到[0,360)范围内

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

计算太阳的赤纬

sinDec = 0.39782 * sin(L) cosDec = cos(asin(sinDec))

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

计算日出/日落的当地平均时间

T = H + RA - (0.06571 * t) - 6.622
调整回协调世界时（UTC）

UT = T - lngHour 注意：UT可能需要通过加/减24来调整到[0,24)范围内
将UT值转换为纬度/经度的本地时区时间

localT = UT + localOffset