Przepisałem do C#, ale gdzieś się pomyliłem, bo nie działa. Słońce praktycznie stoi w miejscu przy zmianie czasu. Szczerze to wolałbym wzory z dokładnym opisem i objaśnieniem, bo nie szukam gotowca.
using System;
namespace Simulation
{
public class SunCalc
{
const double J1970 = 2440588;
const double J2000 = 2451545;
const double Rad2Deg = 57.29578;
const double M0 = 357.5291 * Rad2Deg;
const double M1 = 0.98560028 * Rad2Deg;
const double J0 = 0.0009;
const double J1 = 0.0053;
const double J2 = -0.0069;
const double C1 = 1.9148 * Rad2Deg;
const double C2 = 0.0200 * Rad2Deg;
const double C3 = 0.0003 * Rad2Deg;
const double P = 102.9372 * Rad2Deg;
const double e = 23.45 * Rad2Deg;
const double th0 = 280.1600 * Rad2Deg;
const double th1 = 360.9856235 * Rad2Deg;
const double h0 = -0.83 * Rad2Deg; //sunset angle
const double d0 = 0.53 * Rad2Deg; //sun diameter
const double h1 = -6 * Rad2Deg; //nautical twilight angle
const double h2 = -12 * Rad2Deg; //astronomical twilight angle
const double h3 = -18 * Rad2Deg; //darkness angle
const int msInDay = 86400000;
public double DateTimeToJulianDate(DateTime date)
{
return date.ToOADate() / msInDay - 0.5 + J1970;
}
public DateTime JulianDateToDateTime(double julianDate)
{
double unixTime = (julianDate - 2440587.5) * 86400;
DateTime dtDateTime = new DateTime(1970, 1, 1, 0, 0, 0, 0, DateTimeKind.Utc);
dtDateTime = dtDateTime.AddSeconds(unixTime);
return dtDateTime;
}
public double GetJulianCycle(double J, double lw)
{
return Math.Round(J - J2000 - J0 - lw / (2 * Math.PI));
}
public double GetApproxSolarTransit(double Ht, double lw, double n)
{
return J2000 + J0 + (Ht + lw) / (2 * Math.PI) + n;
}
public double GetSolarMeanAnomaly(double Js)
{
return M0 + M1 * (Js - J2000);
}
public double GetEquationOfCenter(double M)
{
return C1 * Math.Sin(M) + C2 * Math.Sin(2 * M) + C3 * Math.Sin(3 * M);
}
public double GetEclipticLongitude(double M, double C)
{
return M + P + C + Math.PI;
}
public double GetSolarTransit(double Js, double M, double Lsun)
{
return Js + (J1 * Math.Sin(M)) + (J2 * Math.Sin(2 * Lsun));
}
public double GetSunDeclination(double Lsun)
{
return Math.Asin(Math.Sin(Lsun) * Math.Sin(e));
}
public double GetRightAscension(double Lsun)
{
return Math.Atan2(Math.Sin(Lsun) * Math.Cos(e), Math.Cos(Lsun));
}
public double GetSiderealTime(double J, double lw)
{
return th0 + th1 * (J - J2000) - lw;
}
public double GetAzimuth(double th, double a, double phi, double d)
{
var H = th - a;
return Math.Atan2(Math.Sin(H), Math.Cos(H) * Math.Sin(phi) -
Math.Tan(d) * Math.Cos(phi));
}
public double GetAltitude(double th, double a, double phi, double d)
{
var H = th - a;
return Math.Asin(Math.Sin(phi) * Math.Sin(d) +
Math.Cos(phi) * Math.Cos(d) * Math.Cos(H));
}
public double GetHourAngle(double h, double phi, double d)
{
return Math.Acos((Math.Sin(h) - Math.Sin(phi) * Math.Sin(d)) /
(Math.Cos(phi) * Math.Cos(d)));
}
public double GetSunsetJulianDate(double w0, double M, double Lsun, double lw, double n )
{
return GetSolarTransit(GetApproxSolarTransit(w0, lw, n), M, Lsun);
}
public double GetSunriseJulianDate(double Jtransit, double Jset )
{
return Jtransit - (Jset - Jtransit);
}
public (double Azimuth, double Altitude) GetSunPosition(double J, double lw, double phi )
{
var M = GetSolarMeanAnomaly(J);
var C = GetEquationOfCenter(M);
var Lsun = GetEclipticLongitude(M, C);
var d = GetSunDeclination(Lsun);
var a = GetRightAscension(Lsun);
var th = GetSiderealTime(J, lw);
return (GetAzimuth(th, a, phi, d), GetAltitude(th, a, phi, d));
}
public Info GetDayInfo(DateTime date, double lat, double lng, bool detailed)
{
var lw = -lng * Rad2Deg;
var phi = lat * Rad2Deg;
var J = DateTimeToJulianDate(date);
var n = GetJulianCycle(J, lw);
var Js = GetApproxSolarTransit(0, lw, n);
var M = GetSolarMeanAnomaly(Js);
var C = GetEquationOfCenter(M);
var Lsun = GetEclipticLongitude(M, C);
var d = GetSunDeclination(Lsun);
var Jtransit = GetSolarTransit(Js, M, Lsun);
var w0 = GetHourAngle(h0, phi, d);
var w1 = GetHourAngle(h0 + d0, phi, d);
var Jset = GetSunsetJulianDate(w0, M, Lsun, lw, n);
var Jsetstart = GetSunsetJulianDate(w1, M, Lsun, lw, n);
var Jrise = GetSunriseJulianDate(Jtransit, Jset);
var Jriseend = GetSunriseJulianDate(Jtransit, Jsetstart);
var w2 = GetHourAngle(h1, phi, d);
var Jnau = GetSunsetJulianDate(w2, M, Lsun, lw, n);
var Jciv2 = GetSunriseJulianDate(Jtransit, Jnau);
Info info = new Info()
{
Dawn = JulianDateToDateTime(Jciv2),
SunriseStart = JulianDateToDateTime(Jrise),
SunriseEnd = JulianDateToDateTime(Jriseend),
Transit = JulianDateToDateTime(Jtransit),
SunsetStart = JulianDateToDateTime(Jsetstart),
SunsetEnd = JulianDateToDateTime(Jset),
Dusk = JulianDateToDateTime(Jnau)
};
if (detailed)
{
var w3 = GetHourAngle(h2, phi, d);
var w4 = GetHourAngle(h3, phi, d);
var Jastro = GetSunsetJulianDate(w3, M, Lsun, lw, n);
var Jdark = GetSunsetJulianDate(w4, M, Lsun, lw, n);
var Jnau2 = GetSunriseJulianDate(Jtransit, Jastro);
var Jastro2 = GetSunriseJulianDate(Jtransit, Jdark);
info.MorningTwilight = new MorningTwilight() {
AstronomicalStart = JulianDateToDateTime(Jastro2),
AstronomicalEnd = JulianDateToDateTime(Jnau2),
NauticalStart = JulianDateToDateTime(Jnau2),
NauticalEnd = JulianDateToDateTime(Jciv2),
CivilStart = JulianDateToDateTime(Jciv2),
CivilEnd = JulianDateToDateTime(Jrise)
};
info.NightTwilight = new NightTwilight() {
CivilStart = JulianDateToDateTime(Jset),
CivilEnd = JulianDateToDateTime(Jnau),
NauticalStart = JulianDateToDateTime(Jnau),
NauticalEnd = JulianDateToDateTime(Jastro),
AstronomicalStart = JulianDateToDateTime(Jastro),
AstronomicalEnd = JulianDateToDateTime(Jdark)
};
}
return info;
}
public (double Azimuth, double Altitude) GetSunPosition(DateTime date, double lng, double lat)
{
return GetSunPosition(DateTimeToJulianDate(date), -lng * Rad2Deg, lat * Rad2Deg);
}
}
public class Info
{
public DateTime Dawn { get; set; }
public DateTime SunriseStart { get; set; }
public DateTime SunriseEnd { get; set; }
public DateTime Transit { get; set; }
public DateTime SunsetStart { get; set; }
public DateTime SunsetEnd { get; set; }
public DateTime Dusk { get; set; }
public MorningTwilight MorningTwilight { get; set; }
public NightTwilight NightTwilight { get; set; }
}
public class MorningTwilight
{
public DateTime AstronomicalStart { get; set; }
public DateTime AstronomicalEnd { get; set; }
public DateTime NauticalStart { get; set; }
public DateTime NauticalEnd { get; set; }
public DateTime CivilStart { get; set; }
public DateTime CivilEnd { get; set; }
}
public class NightTwilight
{
public DateTime CivilStart { get; set; }
public DateTime CivilEnd { get; set; }
public DateTime NauticalStart { get; set; }
public DateTime NauticalEnd { get; set; }
public DateTime AstronomicalStart { get; set; }
public DateTime AstronomicalEnd { get; set; }
}
}