satoh

知っているようで知らないデジタル図形処理 ~任意形状図形の選択から交点・接点・接線計算まで、方程式では解けないノウハウを紹介

   

日の出時間と日没時間の計算方法~プログラム編

日の出、日没、そして南中時刻を具体的に計算してみる

以前に日の出日没時刻計算の概要を示したので、今回は具体的な計算を示してみましょう。
既知なのは、日付と緯度・経度です。これらから日の出・日没・南中の時刻を求めます。
まず、以前の記事では省略してしまった、

    Et    : 均時差[h]
    δ     : 太陽赤緯[°]

これらを計算しなければなりません。
最初に通算日(1月1日を1とした1~365)から係数 χ を決めます。
通算日を day としましょう。2月1日ならば32です。

    day : 通算日
    π    : 円周率
    χ    : 係数

とすれば、係数χは、

    χ = (day – 1) * 2 * π / 365

です。これにより均時差 Et、太陽赤緯 δ は

    Et = (0.0172 + 0.4281 * cos(χ) – 7.3515 * sin(χ)
                          – 3.3495 * cos(2 * χ) – 9.3619 * sin(2 * χ)) / 60.0
    δ = 0.006918 – 0.399912 * cos(χ) + 0.070257 * sin(χ)
                           – 0.006758 * cos(2 * χ) + 0.000907 * sin(2 * χ)
                           – 0.002697 * cos(3 * χ) + 0.001480 * sin(3 * χ)

となります。
これで各時刻を計算する準備ができました。
緯度・経度は、°で与えるようにします。14度23分15秒であれば、14 + 23 / 60 + 15 / 3600ですね。

    φ  : 緯度[°]
    L   : 経度[°]

得られる時刻を、

    hour   : 南中時刻
    hourS : 日の出時刻
    hourE : 日没時刻

としましょう。
まず簡単な南中時刻を計算しましょう。太陽が真南に来る時刻ですね。
(南中時刻は日によって変化します。一年を通して一定ではありません。詳しくは国立天文台のサイトで確認してください)

    hour  = ((9.0 – Et + 12.0) * 15.0 – L) / 15.0

南中時刻は経度だけで決まるのでこれだけです。
次に日の出・日没時刻を求めます。太陽高度が0°の時がそれに対応しますが、以前の記事で示したのは太陽の視半径と屈折理知を考慮したほうがより現実に近い、といして、-50’/60とそたと思います。
今回は、より現実に近く、[-大気差-眼高差-視半径+地心視差] ≒ 0.899 °を使ってみましょう。

    ω   : 太陽の時角

は、

    ω = acos((sin(-0.899) – sin(φ) * sin(δ)) / (cos(φ) * cos(δ)))

となります。
時角がわかれば、

    hourS = (-ω + (9.0 – Et + 12.0) * 15.0 – L) / 15
    hourE = ( ω + (9.0 – Et + 12.0) * 15.0 – L) / 15

で完了です。
おまけですが、南中時の太陽高度も求めてみましょう。
このとき、時角 ω は南中なので、0°ですから。cos(ω) = 1 、つまり

    h = asin(sin(φ) * sin(δ) + cos(φ) * cos(δ) * cos(ω))
       =  90° – φ + δ

となります。

プログラム

では早速プログラムにしてみます。

#define		M_PI        3.14159265358979323846
#define		DEG(a)		((a) * 180 / M_PI)
#define		RAD(a)		((a) * M_PI / 180)

//----------------------------------------------------------------------------------------------
//		南中時の太陽位置と時刻計算
//----------------------------------------------------------------------------------------------
int CalcSun(long	day,		//---------- I 通算日
			double	phi,		//---------- I 緯度
			double	L,			//---------- I 経度
			double	&hourS,		//---------- O 日出時刻
			double	&hourE,		//---------- O 日没時刻
			double	&hour)		//---------- O 南中時刻
{
	double		khi ;
	double		omegaS ;
	double		Et ;				//------ 均時差[h]
	double		delta ;				//------ 太陽赤緯[°]
	double		omega ;				//------ 太陽の時角[°]
	double		alpha ;				//------ 太陽の方位角[°]
	double		h ;					//------ 太陽高度[°](=90°-θz)

//----------------- °で渡されるのでラジアンに変換する
	phi		= RAD(phi) ;
	L		= RAD(L) ;

//----------------- δ,χの計算
	khi		= (day - 1) * 2.0 * M_PI / 365.0 ;
	delta	= 0.006918 - 0.399912 * cos(khi) + 0.070257 * sin(khi)
					   - 0.006758 * cos(2 * khi) + 0.000907 * sin(2 * khi)
					   - 0.002697 * cos(3 * khi) + 0.001480 * sin(3 * khi) ;

//----------------- Et, ωの計算
	Et		= (0.0172 + 0.4281 * cos(khi) - 7.3515 * sin(khi)
					  - 3.3495 * cos(2 * khi) - 9.3619 * sin(2 * khi)) / 60.0 ;

//----------------- 南中時刻
	hour	= ((9.0 - Et + 12.0) * 15.0 - DEG(L)) / 15.0 ;

//----------------- 南中時の太陽高度(南中時は ω= 0)
//	h		= asin(sin(phi) * sin(delta) + cos(phi) * cos(delta) * cos(omega)) ;
	h		= RAD(90) - phi + delta ;		//------- cos(omega = 0) == 1

//----------------- 日出・日没時刻 : -大気差-眼高差-視半径+地心視差 ≒ -0.899
	omegaS	= acos((sin(RAD(-0.899)) - sin(phi) * sin(delta)) / (cos(phi) * cos(delta))) ;

	hourS = (DEG(-omegaS) + (9.0 - Et + 12.0) * 15.0 - DEG(L)) / 15.0 ;
	hourE = (DEG( omegaS) + (9.0 - Et + 12.0) * 15.0 - DEG(L)) / 15.0 ;

	return 0 ;
}

 - C, ソーラー, 計算