「VCS/AstroSim7」の編集履歴（バックアップ）一覧はこちら

VCS/AstroSim7 - (2013/01/11 (金) 22:06:54) の１つ前との変更点

追加された行は緑色になります。

削除された行は赤色になります。

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

参考
-[[国立天文台 暦計算室>http://eco.mtk.nao.ac.jp/koyomi/]]
-[[均時差 (equation of time)>http://mail2.nara-edu.ac.jp/~asait/kuiper_belt/time/equation_of_time.htm]]

Program.cs
#highlight(c#){{
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秒間に回転する角度

    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 da;
    double dx;

    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 delta = Math.Sin((0.533 / 2 + (35 * 60 + 8) / 3600.0) * Math.PI / 180);
        da = Math.Sin(latRad) * delta;
        dx = Math.Cos(latRad) * delta;
    }

    // 日出・日没・南中時の計算
    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 += 0.5; // 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 eclipticLon = (2 * Math.PI) * ((MJD / solarYear - 0.3399541) % 1.0);
        if (Math.PI < eclipticLon) eclipticLon -= (2 * Math.PI);

        // 平均近点角（概算）近日点から次の近日点までの角度
        double M = (2 * Math.PI) * ((MJD / anomalisticYear - 0.1242853) % 1.0);
        if (Math.PI < M) M -= (2 * Math.PI);

        // 楕円効果と傾斜効果
        int ellipseEffect = (int)Math.Round(-2 * e * Math.Sin(M) / earthRotationPerSec);
        double t = Math.Tan(obliquity * Math.PI / 360);
        int obliquityEffect = (int)Math.Round(t * t * Math.Sin(2 * eclipticLon) / earthRotationPerSec);
        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 + da) * -Math.Tan(latRad) - dx) / 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)));
    }
}
}}

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

参考
-[[国立天文台 暦計算室>http://eco.mtk.nao.ac.jp/koyomi/]]
-[[均時差 (equation of time)>http://mail2.nara-edu.ac.jp/~asait/kuiper_belt/time/equation_of_time.htm]]

Program.cs
#highlight(c#){{
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;
    }
}
}}