Nachylenie Ziemi względem Słońca

0

Cześć chciałbym w przestrzeni 3D obliczyć na podstawie podanej daty i czasu aktualne kąty nachylenia Ziemi względem Słońca by odzwierciedlić przybliżony kierunek padania światła na Ziemię zakładając, że np. Słońce znajduje się w punkcie (0,0,0) a Ziemia np (100,100,100), chciałbym ją obrócić tak by padające światło na model Ziemi odzwierciedlało miejsca w których obecnie jest jasno.

Ktoś mógłby mnie nakierować na jakieś algorytmy? Dzięki.

1

http://suncalc.net/
kod jest w javascript na pierwszy rzut oka tutaj:
http://suncalc.net/scripts/suncalc.js

0

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

0

Ok to działa https://gist.github.com/paulhayes/54a7aa2ee3cccad4d37bb65977eb19e2 Trzeba tylko podać odpowiednie kąty i obraca się wokół globu jak należy.

1 użytkowników online, w tym zalogowanych: 0, gości: 1