VCS/AstroSim7 - Sample is best!

開発環境	Microsoft Visual C# 2010 Express (SP1)
実行環境	Microsoft Windows XP Home Edition (SP3)
プロジェクトの種類	空のプロジェクト
プロジェクト名	AstroSim4

参考

• 国立天文台 暦計算室
• 均時差 (equation of time)

Program.cs

using System;
 
class Program
{
    const double solarYear = 365.24219; // 太陽年(day) 365d5h48m45s
    const double anomalisticYear = 365.259643; // 近点年(day) 365d6h13m53.1552s
    const double e = 0.01671022; // 離心率(Orbital eccentricity)
    const double earthRotationPerSec = (2 * Math.PI) / 86400; // 地球が1秒間に回転する角度
    readonly double K = Math.Sqrt((1 + e) / (1 - e));
    const double epsilon = 1.0e-14;
 
    const double lat = 35.4500; // 緯度
    const double lon = 139.6500; // 経度
    const double ST = 9; // 標準時
    const double latRad = lat * Math.PI / 180;
    readonly int equationOfTimeRev = (int)Math.Round(240 * lon - 3600 * ST); // 南中時補正
 
    double ra;
    double rx;
 
    static void Main()
    {
        Program program = new Program();
        Console.WriteLine(string.Format("lat:{0:f4} lon:{1:f4} ST:{2}", lat, lon, ST));
        DateTime dateStart = new DateTime(2013, 1, 1);
        DateTime dateEnd = dateStart.AddYears(1);
        for (DateTime dt = dateStart; dt < dateEnd; dt = dt.AddDays(1))
        {
            program.Calc(dt);
        }
        Console.ReadLine();
    }
 
    Program()
    {
        // 太陽の視角(0.533deg)と大気差(35m8s)による昼の長さの補正
        double revision = Math.Sin((0.533 / 2 + (35 * 60 + 8) / 3600.0) * Math.PI / 180);
        ra = Math.Sin(latRad) * revision;
        rx = Math.Cos(latRad) * revision;
    }
 
    // 日出・日没・南中時の計算
    void Calc(DateTime dt)
    {
        // 修正ユリウス日
        int y = dt.Year;
        int m = dt.Month;
        int d = dt.Day;
        if (m < 3)
        {
            y--;
            m += 12;
        }
        double MJD = (int)(365.25 * y) + (y / 400) - (y / 100) + (int)(30.59 * (m - 2)) + d - 678912;
        MJD += (12 - ST) / 24; // 標準時における12:00
 
        // 黄道傾斜角
        double T = (MJD - 51554.5) / 36525; // 2000/1/1 12:00(UT)からのユリウス世紀(36525日)
        double obliquity = (84381.406 - 46.836769 * T - 0.00059 * T * T + 0.001813 * T * T * T) / 3600;
 
        // 平均近点角（概算）近日点から次の近日点までの角度
        double Ma = ModAngle((2 * Math.PI) * ((MJD / anomalisticYear - 0.1242853) % 1.0));
        double E;
        double Ta; // 真近点角
        KeplersEquation(Ma, out E, out Ta);
 
        // 春分点(vernal equinox)の真近点角
        double MJDv = ((int)(MJD / solarYear - 0.3399541) + 0.3399541) * solarYear;
        double Mv = (2 * Math.PI) * ((MJDv / anomalisticYear - 0.1242853) % 1.0);
        double Tv;
        KeplersEquation(Mv, out E, out Tv);
 
        // 黄径（概算）春分点から次の春分点までの角度
        //double eclipticLon = ModAngle((2 * Math.PI) * ((MJD / solarYear - 0.3399541) % 1.0));
        double eclipticLon = ModAngle(Ta - Tv);
 
        // 楕円効果と傾斜効果
        //int ellipseEffect2 = (int)Math.Round(-2 * e * Math.Sin(M) / earthRotationPerSec);
        //double t = Math.Tan(obliquity * Math.PI / 360);
        //int obliquityEffect2 = (int)Math.Round(t * t * Math.Sin(2 * eclipticLon) / earthRotationPerSec);
        int ellipseEffect = (int)Math.Round(Ma / earthRotationPerSec - Ta / earthRotationPerSec);
        int obliquityEffect = CalcObliquityEffect(obliquity, eclipticLon);
        int equationOfTime = ellipseEffect + obliquityEffect; // 均時差
        int transit = 43200 - equationOfTime - equationOfTimeRev; // 南中時
 
        // 天球上の太陽軌道の高さと半径
        double solarDecl = Math.Sin(eclipticLon) * obliquity * Math.PI / 180; // 太陽の赤緯
        double solarAlt = Math.Sin(solarDecl);
        double solarRad = Math.Cos(solarDecl);
 
        // 天球上の太陽軌道と地平面の交点＝日出・日没
        double x = ((solarAlt + ra) * -Math.Tan(latRad) - rx) / solarRad;
        int halfDaytime = (int)Math.Round(43200 * Math.Acos(x) / Math.PI);
        int rising = transit - halfDaytime;
        int setting = transit + halfDaytime;
 
        Console.WriteLine(string.Format(@"{0} {1:hh\:mm} {2:hh\:mm} {3:hh\:mm}",
            dt.ToShortDateString(), new TimeSpan(0, 0, rising + 30),
            new TimeSpan(0, 0, transit + 30), new TimeSpan(0, 0, setting + 30)));
    }
 
    double ModAngle(double angle)
    {
        while (angle <= -Math.PI) angle += (2 * Math.PI);
        while (Math.PI < angle) angle -= (2 * Math.PI);
        return angle;
    }
 
    // 傾斜効果の計算
    int CalcObliquityEffect(double obliquity, double eclipticLon)
    {
        // 地球を基準とした太陽の公転
        double x = Math.Cos(eclipticLon);
        double r = Math.Sin(eclipticLon);
        double y = Math.Cos(obliquity * Math.PI / 180) * r;
        double celestialEquator = Math.Atan2(y, x); // 天の赤道上の角度
        return (int)Math.Round(eclipticLon / earthRotationPerSec - celestialEquator / earthRotationPerSec);
    }
 
    // 漸化式によりケプラー方程式を解く
    // M 平均近点角(mean anomaly)
    // E 離心近点角(Eccentric anomaly)
    // T 真近点角(true anomaly)
    void KeplersEquation(double M, out double E, out double T)
    {
        double E0 = M; // 初項
        for (int i = 0; ; )
        {
            i++;
            E = M + e * Math.Sin(E0);
            if ((E0 - epsilon < E) && (E < E0 + epsilon))
            {
                break;
            }
            if (10 <= i)
            {
                Console.WriteLine(string.Format("計算打ち切り M={0} E={1}", M, E));
                break;
            }
            E0 = E;
        }
        T = Math.Atan(K * Math.Tan(E / 2)) * 2;
    }
}