File Coverage

blib/lib/Date/Qreki.pm

Criterion	Covered	Total	%
statement	353	406	86.9
branch	37	60	61.6
condition	15	36	41.6
subroutine	15	16	93.7
pod	2	12	16.6
total	422	530	79.6

=encoding UTF-8

=head1 NAME

Date::Qreki - convert Gregorian to Japanese "kyureki" dates.

=head1 SYNOPSIS

This module is currently undocumented, so to use it you need to read
the source code. If you are interested in helping write documentation
for the module or would like to discuss, please contact the module
author or visit
L
and vote for this project.

=head1 DESCRIPTION

=head1 FUNCTIONS

=head2 calc_kyureki

=head2 get_rokuyou

=head1 SEE ALSO

=head1 COPYRIGHT

Date::Qreki is copyright H. Takano, N. Ueno.

=head1 LICENSE

The original licence in Japanese reads as follows:

19.  配布規定について

    本当は、「一般常識の許す範囲で...」としたい所ですが、世の中には色々な人
  がいますので、ある程度の枠組みが必要なのは仕方の無いことなのでしょう。 と
  もあれ拙者が希望するのは以下のようなものです。

     i   本スクリプト・説明書を再配布する場合はオリジナルのアーカイブファイ
       ルに含まれるファイルを全て含み、オリジナルのスクリプト・説明書を改変
       しないで下さい。 また、 頒布の都合でアーカイブ形式（.arc  .zip  .zoo
       等）を変更する場合も同様に扱って下さい。

    ii   再配布する際の媒体に要するコストを除いて一切の金銭等の授受は禁止い
       たします。    （他の言語に移植したものを配布する場合にも適用いたしま
       す。）

   iii   他の言語に移植したり、  本スクリプトの一部または全部を引用して作成
       し、これを配布する場合には、オリジナルのスクリプトと本説明書を必ず同
       梱して下さい。
       その場合の著作権につきましては、引用した部分のみ著者に帰属し、その外
       は製作者に帰属します。

    iv   配布によって著作者が一切の制限を受ける可能性がないようにして下さ
       い。

  なお、 いづれの場合におきましても、承諾・転載の報告など一切不要といたしま
  す。

Translation:

19 Conditions for redistribution

I would prefer to release this under the conditions of "whatever
normal common sense allows", but since there are a variety of people
in the world, unfortunately I think some form of framework is
necessary. My wishes are as follows

i. If this script is redistributed, please also include all the
original scripts, and do not modify the original script and its
documentation. This also applies if the archive format (.arc, .zip,
.zoo, etc.) is changed.

ii. I forbid any resale of this program beyond the cost of
distribution media. (This also applies if the program is translated to
a different programming language).

iii. If this is translated to a different programming language, or if
part or all of this script is used as part of another program, please
be sure to include the original script and this explanation with
it. If part of this program is used, I retain the copyright of the
quoted material only and the remaining program remains copyrighted by
its author.

iv. Please do not allow the author of this software to be placed under
any restrictions by your redistribution.

In any case, there is absolutely no need to obtain the author's
consent or notify the author.

End of artistic translation.

The original Awk script and its documentation are in F in
the distribution.

=head1 AUTHOR

Original Awk script by H. Takano. Perl conversion by N. Ueno. This
CPAN-ification of N. Ueno's Perl script was performed by Ben Bullock.

For enquiries about this Perl module, please contact Ben Bullock
.

=cut

package Date::Qreki;
require Exporter;
@ISA = qw(Exporter);
@EXPORT_OK = qw/calc_kyureki get_rokuyou/;
use warnings;
use strict;
our $VERSION = 0.05;

#=========================================================================
# 旧暦計算サンプルプログラム  $Revision:   1.1  $
# Coded by H.Takano 1993,1994
#
# Arranged for Perl Script by N.Ueno
#
#
# オリジナルのスクリプトは高野氏のAWKです。下記より入手できます。
# http://www.vector.co.jp/soft/dos/personal/se016093.html
#
#
#========================================================================

#-----------------------------------------------------------------------
# 円周率の定義と（角度の）度からラジアンに変換する係数の定義
#-----------------------------------------------------------------------
use constant PI => 3.141592653589793238462;
use constant k => PI/180.0;

sub deg_cos
{
    my ($angle) = @_;
    return cos ($angle * k);
}

#=========================================================================
# 六曜算出関数
#
# 引数：新暦年月日
# 戻値：0:大安 1:赤口 2:先勝 3:友引 4:先負 5:仏滅
#
#=========================================================================
sub get_rokuyou
{

	my ($year,$mon,$day) = @_;
	my ($tm0,$q_year,$q_mon,$q_day,$uruu,$q_yaer);

	($q_yaer,$uruu,$q_mon,$q_day) = calc_kyureki($year,$mon,$day);

	return(($q_mon + $q_day) % 6);
}

#=========================================================================
# 新暦に対応する、旧暦を求める。
#
# 呼び出し時にセットする変数
# 引　数　year : 計算する日付
#         mon
#         day
#
# 戻り値　kyureki : 答えの格納先（配列に答えをかえす）
#         　　  kyureki[0] : 旧暦年
#         　　  kyureki[1] : 平月／閏月 flag .... 平月:0 閏月:1
#         　　  kyureki[2] : 旧暦月
#         　　  kyureki[3] : 旧暦日
#
#=========================================================================
sub calc_kyureki
{
	my ($year,$mon,$day) = @_;
	my (@kyureki,$tm,@saku,$lap,@a,$i,@m);

	my $tm0 = YMDT2JD($year,$mon,$day,0,0,0);

#-----------------------------------------------------------------------
# 計算対象の直前にあたる二分二至の時刻を求める
# chu[0,0]:二分二至の時刻  chu[0,1]:その時の太陽黄経
#-----------------------------------------------------------------------
        my @chu;
	($chu[0][0],$chu[0][1]) = before_nibun($tm0);

#-----------------------------------------------------------------------
# 中気の時刻を計算（４回計算する）
# chu[i,0]:中気の時刻  chu[i,1]:太陽黄経
#-----------------------------------------------------------------------
	for($i=1;$i<4;$i++){
		($chu[$i][0],$chu[$i][1]) = calc_chu($chu[$i-1][0]+32.0);
	}

#-----------------------------------------------------------------------
#  計算対象の直前にあたる二分二至の直前の朔の時刻を求める
#-----------------------------------------------------------------------
	$saku[0] = calc_saku($chu[0][0]);


#-----------------------------------------------------------------------
# 朔の時刻を求める
#-----------------------------------------------------------------------
	for($i=1;$i<5;$i++){
		$tm=$saku[$i-1];
		$tm += 30.0;
		$saku[$i]=calc_saku($tm);

# 前と同じ時刻を計算した場合（両者の差が26日以内）には、初期値を
# +33日にして再実行させる。
		if( abs( int($saku[$i-1])-int($saku[$i]) ) <= 26.0 ){
			$saku[$i]=calc_saku($saku[$i-1]+35.0);
		}
	}

#-----------------------------------------------------------------------
# saku[1]が二分二至の時刻以前になってしまった場合には、朔をさかのぼり過ぎ
# たと考えて、朔の時刻を繰り下げて修正する。
# その際、計算もれ（saku[4]）になっている部分を補うため、朔の時刻を計算
# する。（近日点通過の近辺で朔があると起こる事があるようだ...？）
#-----------------------------------------------------------------------
	if( int($saku[1]) <= int($chu[0][0]) ){
		for($i=0;$i<5;$i++){
			$saku[$i]=$saku[$i+1];
		}
		$saku[4] = calc_saku($saku[3]+35.0);
	}

#-----------------------------------------------------------------------
# saku[0]が二分二至の時刻以後になってしまった場合には、朔をさかのぼり足
# りないと見て、朔の時刻を繰り上げて修正する。
# その際、計算もれ（saku[0]）になっている部分を補うため、朔の時刻を計算
# する。（春分点の近辺で朔があると起こる事があるようだ...？）
#-----------------------------------------------------------------------
	elsif( int($saku[0]) > int($chu[0][0]) ){
		for($i=4;$i>0;$i--){
			$saku[$i] = $saku[$i-1];
		}
		$saku[0] = calc_saku($saku[0]-27.0);
	}

#-----------------------------------------------------------------------
# 閏月検索Ｆｌａｇセット
# （節月で４ヶ月の間に朔が５回あると、閏月がある可能性がある。）
# lap=0:平月  lap=1:閏月
#-----------------------------------------------------------------------
	if(int($saku[4]) <= int($chu[3][0]) ){
		$lap=1;
	}else{
		$lap=0;
	}

#-----------------------------------------------------------------------
# 朔日行列の作成
# m[i,0] ... 月名（1:正月 2:２月 3:３月 ....）
# m[i,1] ... 閏フラグ（0:平月 1:閏月）
# m[i,2] ... 朔日のjd
#-----------------------------------------------------------------------
	$m[0][0]=int($chu[0][1]/30.0) + 2;
	if(defined $m[0][1] && $m[0][1] > 12 ){
		$m[0][0]-=12;
	}
	$m[0][2]=int($saku[0]);
	$m[0][1]=0;

	for($i=1;$i<5;$i++){
		if($lap == 1 && $i !=1 ){
			if( int($chu[$i-1][0]) <= int($saku[$i-1]) || int($chu[$i-1][0]) >= int($saku[$i]) ){
				$m[$i-1][0] = $m[$i-2][0];
				$m[$i-1][1] = 1;
				$m[$i-1][2] = int($saku[$i-1]);
				$lap=0;
			}
		}
		$m[$i][0] = $m[$i-1][0]+1;
		if( $m[$i][0] > 12 ){
			$m[$i][0]-=12;
		}
		$m[$i][2]=int($saku[$i]);
		$m[$i][1]=0;
	}

#-----------------------------------------------------------------------
# 朔日行列から旧暦を求める。
#-----------------------------------------------------------------------
	my $state=0;
	for($i=0;$i<5;$i++){
		if(int($tm0) < int($m[$i][2])){
			$state=1;
			last;
		}elsif(int($tm0) == int($m[$i][2])){
			$state=2;
			last;
		}
	}
	if($state==0||$state==1){
		$i--;
	}

	$kyureki[1]=$m[$i][1];
	$kyureki[2]=$m[$i][0];
	$kyureki[3]=int($tm0)-int($m[$i][2])+1;

#-----------------------------------------------------------------------
# 旧暦年の計算
# （旧暦月が10以上でかつ新暦月より大きい場合には、
#   まだ年を越していないはず...）
#-----------------------------------------------------------------------

	@a = JD2YMDT($tm0);
	$kyureki[0] = $a[0];
	if($kyureki[2] > 9 && $kyureki[2] > $a[1]){
		$kyureki[0]--;
	}

	return($kyureki[0],$kyureki[1],$kyureki[2],$kyureki[3]);

}

#=========================================================================
# 中気の時刻を求める
#
# 呼び出し時にセットする変数
# tm ........ 計算対象となる時刻（ユリウス日）
# chu ....... 戻り値を格納する配列のポインター
# i ......... 戻り値を格納する配列の要素番号
# 戻り値 .... 中気の時刻、その時の黄経を配列で渡す
#
#=========================================================================
sub calc_chu
{
	my ($tm) = @_;
	my ($tm1,$tm2,$t,$rm_sun0,$rm_sun,$delta_t1,$delta_t2,$delta_rm);
	my (@temp);
#-----------------------------------------------------------------------
#時刻引数を分解する
#-----------------------------------------------------------------------
	$tm1 = int( $tm );
	$tm2 = $tm - $tm1;

#-----------------------------------------------------------------------
# JST ==> DT （補正時刻=0.0sec と仮定して計算）
#-----------------------------------------------------------------------
	$tm2-=9.0/24.0;

#-----------------------------------------------------------------------
# 中気の黄経 λsun0 を求める
#-----------------------------------------------------------------------
	$t=($tm2+0.5) / 36525.0;
	$t=$t + ($tm1-2451545.0) / 36525.0;
	$rm_sun = LONGITUDE_SUN( $t );

	$rm_sun0 = 30.0*int($rm_sun/30.0);

#-----------------------------------------------------------------------
# 繰り返し計算によって中気の時刻を計算する
# （誤差が±1.0 sec以内になったら打ち切る。）
#-----------------------------------------------------------------------
        $delta_t1 = 0;
	for( $delta_t2 = 1.0 ; abs( $delta_t1 + $delta_t2 ) > ( 1.0 / 86400.0 ) ; ){

#-----------------------------------------------------------------------
# λsun を計算
#-----------------------------------------------------------------------
		$t =($tm2+0.5) / 36525.0;
		$t =$t + ($tm1-2451545.0) / 36525.0;
		$rm_sun=LONGITUDE_SUN( $t );

#-----------------------------------------------------------------------
# 黄経差 Δλ＝λsun −λsun0
#-----------------------------------------------------------------------
		$delta_rm = $rm_sun - $rm_sun0 ;

#-----------------------------------------------------------------------
# Δλの引き込み範囲（±180°）を逸脱した場合には、補正を行う
#-----------------------------------------------------------------------
		if( $delta_rm > 180.0 ){
			$delta_rm-=360.0;
		}elsif( $delta_rm < -180.0 ){
			$delta_rm+=360.0;
		}

#-----------------------------------------------------------------------
# 時刻引数の補正値 Δt
# delta_t = delta_rm * 365.2 / 360.0;
#-----------------------------------------------------------------------
		$delta_t1 = int($delta_rm * 365.2 / 360.0);
		$delta_t2 = $delta_rm * 365.2 / 360.0;
		$delta_t2 -= $delta_t1;

#-----------------------------------------------------------------------
# 時刻引数の補正
# tm -= delta_t;
#-----------------------------------------------------------------------
		$tm1 = $tm1 - $delta_t1;
		$tm2 = $tm2 - $delta_t2;
		if($tm2 < 0){
			$tm2+=1.0;$tm1-=1.0;
		}
	}

#-----------------------------------------------------------------------
# 戻り値の作成
# chu[i,0]:時刻引数を合成するのと、DT ==> JST 変換を行い、戻り値とする
# （補正時刻=0.0sec と仮定して計算）
# chu[i,1]:黄経
#-----------------------------------------------------------------------
	$temp[0] = $tm2+9.0/24.0;
	$temp[0] += $tm1;
	$temp[1] = $rm_sun0;

	return(@temp);
}

#=========================================================================
# 直前の二分二至の時刻を求める
#
# 呼び出し時にセットする変数
# tm ........ 計算対象となる時刻（ユリウス日）
# nibun ..... 戻り値を格納する配列のポインター
# 戻り値 .... 二分二至の時刻、その時の黄経を配列で渡す
# （戻り値の渡し方がちょっと気にくわないがまぁいいや。）
#=========================================================================
sub before_nibun
{
	my ($tm) = @_;
	my (@nibun,$tm1,$tm2,$t,$rm_sun0,$rm_sun,$delta_t1,$delta_t2,$delta_rm);

#-----------------------------------------------------------------------
#時刻引数を分解する
#-----------------------------------------------------------------------
	$tm1 = int( $tm );
	$tm2 = $tm - $tm1;

#-----------------------------------------------------------------------
# JST ==> DT （補正時刻=0.0sec と仮定して計算）
#-----------------------------------------------------------------------
	$tm2-=9.0/24.0;

#-----------------------------------------------------------------------
# 直前の二分二至の黄経 λsun0 を求める
#-----------------------------------------------------------------------
	$t=($tm2+0.5) / 36525.0;
	$t=$t + ($tm1-2451545.0) / 36525.0;
	$rm_sun=LONGITUDE_SUN( $t );
	$rm_sun0=90*int($rm_sun/90.0);

#-----------------------------------------------------------------------
# 繰り返し計算によって直前の二分二至の時刻を計算する
# （誤差が±1.0 sec以内になったら打ち切る。）
#-----------------------------------------------------------------------
        $delta_t1 = 0;
	for( $delta_t2 = 1.0 ; abs( $delta_t1+$delta_t2 ) > ( 1.0 / 86400.0 ) ; ){

#-----------------------------------------------------------------------
# λsun を計算
#-----------------------------------------------------------------------
		$t=($tm2+0.5) / 36525.0;
		$t=$t + ($tm1-2451545.0) / 36525.0;
		$rm_sun=LONGITUDE_SUN( $t );

#-----------------------------------------------------------------------
# 黄経差 Δλ＝λsun −λsun0
#-----------------------------------------------------------------------
		$delta_rm = $rm_sun - $rm_sun0 ;

#-----------------------------------------------------------------------
# Δλの引き込み範囲（±180°）を逸脱した場合には、補正を行う
#-----------------------------------------------------------------------
		if( $delta_rm > 180.0 ){
			$delta_rm-=360.0;
		}elsif( $delta_rm < -180.0){
			$delta_rm+=360.0;
		}

#-----------------------------------------------------------------------
# 時刻引数の補正値 Δt
# delta_t = delta_rm * 365.2 / 360.0;
#-----------------------------------------------------------------------
		$delta_t1 = int($delta_rm * 365.2 / 360.0);
		$delta_t2 = $delta_rm * 365.2 / 360.0;
		$delta_t2 -= $delta_t1;

#-----------------------------------------------------------------------
# 時刻引数の補正
# tm -= delta_t;
#-----------------------------------------------------------------------
		$tm1 = $tm1 - $delta_t1;
		$tm2 = $tm2 - $delta_t2;
		if($tm2 < 0){
			$tm2+=1.0;$tm1-=1.0;
		}

	}

#-----------------------------------------------------------------------
# 戻り値の作成
# nibun[0,0]:時刻引数を合成するのと、DT ==> JST 変換を行い、戻り値とする
# （補正時刻=0.0sec と仮定して計算）
# nibun[0,1]:黄経
#-----------------------------------------------------------------------
	$nibun[0] = $tm2+9.0/24.0;
	$nibun[0] += $tm1;
	$nibun[1] = $rm_sun0;

	return(@nibun);

}

#=========================================================================
# 朔の計算
# 与えられた時刻の直近の朔の時刻（JST）を求める
#
# 呼び出し時にセットする変数
# tm ........ 計算対象となる時刻（ユリウス日）
# 戻り値 .... 朔の時刻
#
# ※ 引数、戻り値ともユリウス日で表し、時分秒は日の小数で表す。
#
#=========================================================================
sub calc_saku
{
	my ($tm) = @_;
	my ($lc,$t,$tm1,$tm2,$rm_sun,$rm_moon,$delta_rm,$delta_t1,$delta_t2);

#-----------------------------------------------------------------------
# ループカウンタのセット
#-----------------------------------------------------------------------
	$lc=1;

#-----------------------------------------------------------------------
#時刻引数を分解する
#-----------------------------------------------------------------------
	$tm1 = int( $tm );
	$tm2 = $tm - $tm1;

#-----------------------------------------------------------------------
# JST ==> DT （補正時刻=0.0sec と仮定して計算）
#-----------------------------------------------------------------------
	$tm2-=9.0/24.0;

#-----------------------------------------------------------------------
# 繰り返し計算によって朔の時刻を計算する
# （誤差が±1.0 sec以内になったら打ち切る。）
#-----------------------------------------------------------------------
        $delta_t1 = 0;
	for( $delta_t2 = 1.0 ; abs( $delta_t1+$delta_t2 ) > ( 1.0 / 86400.0 ) ; $lc++){

#-----------------------------------------------------------------------
# 太陽の黄経λsun ,月の黄経λmoon を計算
# t = (tm - 2451548.0 + 0.5)/36525.0;
#-----------------------------------------------------------------------
		$t=($tm2+0.5) / 36525.0;
		$t=$t + ($tm1-2451545.0) / 36525.0;
		$rm_sun=LONGITUDE_SUN( $t );
		$rm_moon=LONGITUDE_MOON( $t );

#-----------------------------------------------------------------------
# 月と太陽の黄経差Δλ
# Δλ＝λmoon−λsun
#-----------------------------------------------------------------------
		$delta_rm = $rm_moon - $rm_sun ;

#-----------------------------------------------------------------------
# ループの１回目（lc=1）で delta_rm < 0.0 の場合には引き込み範囲に
# 入るように補正する
#-----------------------------------------------------------------------
		if( $lc==1 && $delta_rm < 0.0 ){
			$delta_rm = NORMALIZATION_ANGLE( $delta_rm );
		}
#-----------------------------------------------------------------------
#   春分の近くで朔がある場合（0 ≦λsun≦ 20）で、月の黄経λmoon≧300 の
#   場合には、Δλ＝ 360.0 − Δλ と計算して補正する
#-----------------------------------------------------------------------
		elsif( $rm_sun >= 0 && $rm_sun <= 20 && $rm_moon >= 300 ){
			$delta_rm = NORMALIZATION_ANGLE( $delta_rm );
			$delta_rm = 360.0 - $delta_rm;
		}
#-----------------------------------------------------------------------
# Δλの引き込み範囲（±40°）を逸脱した場合には、補正を行う
#-----------------------------------------------------------------------
		elsif( abs( $delta_rm ) > 40.0 ) {
			$delta_rm = NORMALIZATION_ANGLE( $delta_rm );
		}

#-----------------------------------------------------------------------
# 時刻引数の補正値 Δt
# delta_t = delta_rm * 29.530589 / 360.0;
#-----------------------------------------------------------------------
		$delta_t1 = int($delta_rm * 29.530589 / 360.0);
		$delta_t2 = $delta_rm * 29.530589 / 360.0;
		$delta_t2 -= $delta_t1;

#-----------------------------------------------------------------------
# 時刻引数の補正
# tm -= delta_t;
#-----------------------------------------------------------------------
		$tm1 = $tm1 - $delta_t1;
		$tm2 = $tm2 - $delta_t2;
		if($tm2 < 0.0){
			$tm2+=1.0;$tm1-=1.0;
		}

#-----------------------------------------------------------------------
# ループ回数が15回になったら、初期値 tm を tm-26 とする。
#-----------------------------------------------------------------------
		if($lc == 15 && abs( $delta_t1+$delta_t2 ) > ( 1.0 / 86400.0 ) ){
			$tm1 = int( $tm-26 );
			$tm2 = 0;
		}

#-----------------------------------------------------------------------
# 初期値を補正したにも関わらず、振動を続ける場合には初期値を答えとして
# 返して強制的にループを抜け出して異常終了させる。
#-----------------------------------------------------------------------
		elsif( $lc > 30 && abs( $delta_t1+$delta_t2 ) > ( 1.0 / 86400.0 ) ){
			$tm1=$tm;$tm2=0;
			last;
		}
	}

#-----------------------------------------------------------------------
# 時刻引数を合成するのと、DT ==> JST 変換を行い、戻り値とする
# （補正時刻=0.0sec と仮定して計算）
#-----------------------------------------------------------------------

	return($tm2+$tm1+9.0/24.0);
}

#=========================================================================
#  角度の正規化を行う。すなわち引数の範囲を ０≦θ＜３６０ にする。
#=========================================================================
sub NORMALIZATION_ANGLE
{
	my ($angle) = @_;
	my ($angle1,$angle2);

	if( $angle < 0.0 ){
		$angle1 = -$angle;
		$angle2 = int( $angle1 / 360.0 );
		$angle1 -= 360.0 * $angle2;
		$angle1 = 360.0 - $angle1;
	}else{
		$angle1 = int( $angle / 360.0 );
		$angle1 = $angle - 360.0 * $angle1;
	}

	return($angle1);
}

#=========================================================================
# 太陽の黄経 λsun を計算する
#=========================================================================
sub LONGITUDE_SUN
{
	my ($t) = @_;
	my ($th,$ang);

#-----------------------------------------------------------------------
# 摂動項の計算
#-----------------------------------------------------------------------
	$ang = NORMALIZATION_ANGLE(  31557.0 * $t + 161.0 );
	$th =       .0004 * deg_cos( $ang );
	$ang = NORMALIZATION_ANGLE(  29930.0 * $t +  48.0 );
	$th = $th +  .0004 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(   2281.0 * $t + 221.0 );
	$th = $th +  .0005 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(    155.0 * $t + 118.0 );
	$th = $th +  .0005 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  33718.0 * $t + 316.0 );
	$th = $th +  .0006 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(   9038.0 * $t +  64.0 );
	$th = $th +  .0007 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(   3035.0 * $t + 110.0 );
	$th = $th +  .0007 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  65929.0 * $t +  45.0 );
	$th = $th +  .0007 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  22519.0 * $t + 352.0 );
	$th = $th +  .0013 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  45038.0 * $t + 254.0 );
	$th = $th +  .0015 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE( 445267.0 * $t + 208.0 );
	$th = $th +  .0018 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(     19.0 * $t + 159.0 );
	$th = $th +  .0018 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  32964.0 * $t + 158.0 );
	$th = $th +  .0020 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  71998.1 * $t + 265.1 );
	$th = $th +  .0200 * deg_cos ($ang );
	$ang = NORMALIZATION_ANGLE(  35999.05 * $t + 267.52 );
	$th = $th - 0.0048 * $t * deg_cos ($ang ) ;
	$th = $th + 1.9147     * deg_cos ($ang ) ;

#-----------------------------------------------------------------------
# 比例項の計算
#-----------------------------------------------------------------------
	$ang = NORMALIZATION_ANGLE( 36000.7695 * $t );
	$ang = NORMALIZATION_ANGLE( $ang + 280.4659 );
	$th  = NORMALIZATION_ANGLE( $th + $ang );

	return($th);
}

#=========================================================================
# 月の黄経 λmoon を計算する
#=========================================================================
sub LONGITUDE_MOON
{
	my ($t) = @_;
	my ($th,$ang);

#-----------------------------------------------------------------------
# 摂動項の計算
#-----------------------------------------------------------------------
  $ang = NORMALIZATION_ANGLE( 2322131.0  * $t + 191.0  );
   $th =      .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(    4067.0  * $t +  70.0  );
   $th = $th + .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  549197.0  * $t + 220.0  );
   $th = $th + .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1808933.0  * $t +  58.0  );
   $th = $th + .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  349472.0  * $t + 337.0  );
   $th = $th + .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  381404.0  * $t + 354.0  );
   $th = $th + .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  958465.0  * $t + 340.0  );
   $th = $th + .0003 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   12006.0  * $t + 187.0  );
   $th = $th + .0004 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   39871.0  * $t + 223.0  );
   $th = $th + .0004 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  509131.0  * $t + 242.0  );
   $th = $th + .0005 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1745069.0  * $t +  24.0  );
   $th = $th + .0005 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1908795.0  * $t +  90.0  );
   $th = $th + .0005 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 2258267.0  * $t + 156.0  );
   $th = $th + .0006 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  111869.0  * $t +  38.0  );
   $th = $th + .0006 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   27864.0  * $t + 127.0  );
   $th = $th + .0007 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  485333.0  * $t + 186.0  );
   $th = $th + .0007 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  405201.0  * $t +  50.0  );
   $th = $th + .0007 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  790672.0  * $t + 114.0  );
   $th = $th + .0007 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1403732.0  * $t +  98.0  );
   $th = $th + .0008 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  858602.0  * $t + 129.0  );
   $th = $th + .0009 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1920802.0  * $t + 186.0  );
   $th = $th + .0011 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1267871.0  * $t + 249.0  );
   $th = $th + .0012 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1856938.0  * $t + 152.0  );
   $th = $th + .0016 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  401329.0  * $t + 274.0  );
   $th = $th + .0018 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  341337.0  * $t +  16.0  );
   $th = $th + .0021 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   71998.0  * $t +  85.0  );
   $th = $th + .0021 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  990397.0  * $t + 357.0  );
   $th = $th + .0021 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  818536.0  * $t + 151.0  );
   $th = $th + .0022 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  922466.0  * $t + 163.0  );
   $th = $th + .0023 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   99863.0  * $t + 122.0  );
   $th = $th + .0024 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1379739.0  * $t +  17.0  );
   $th = $th + .0026 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  918399.0  * $t + 182.0  );
   $th = $th + .0027 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(    1934.0  * $t + 145.0  );
   $th = $th + .0028 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  541062.0  * $t + 259.0  );
   $th = $th + .0037 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1781068.0  * $t +  21.0  );
   $th = $th + .0038 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(     133.0  * $t +  29.0  );
   $th = $th + .0040 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1844932.0  * $t +  56.0  );
   $th = $th + .0040 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1331734.0  * $t + 283.0  );
   $th = $th + .0040 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  481266.0  * $t + 205.0  );
   $th = $th + .0050 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   31932.0  * $t + 107.0  );
   $th = $th + .0052 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  926533.0  * $t + 323.0  );
   $th = $th + .0068 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  449334.0  * $t + 188.0  );
   $th = $th + .0079 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  826671.0  * $t + 111.0  );
   $th = $th + .0085 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1431597.0  * $t + 315.0  );
   $th = $th + .0100 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1303870.0  * $t + 246.0  );
   $th = $th + .0107 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  489205.0  * $t + 142.0  );
   $th = $th + .0110 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1443603.0  * $t +  52.0  );
   $th = $th + .0125 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   75870.0  * $t +  41.0  );
   $th = $th + .0154 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  513197.9  * $t + 222.5  );
   $th = $th + .0304 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  445267.1  * $t +  27.9  );
   $th = $th + .0347 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  441199.8  * $t +  47.4  );
   $th = $th + .0409 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  854535.2  * $t + 148.2  );
   $th = $th + .0458 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 1367733.1  * $t + 280.7  );
   $th = $th + .0533 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  377336.3  * $t +  13.2  );
   $th = $th + .0571 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   63863.5  * $t + 124.2  );
   $th = $th + .0588 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  966404.0  * $t + 276.5  );
   $th = $th + .1144 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(   35999.05 * $t +  87.53 );
   $th = $th + .1851 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  954397.74 * $t + 179.93 );
   $th = $th + .2136 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  890534.22 * $t + 145.7  );
   $th = $th + .6583 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE(  413335.35 * $t +  10.74 );
   $th = $th + 1.2740 * deg_cos ($ang );
  $ang = NORMALIZATION_ANGLE( 477198.868 * $t + 44.963 );
   $th = $th + 6.2888 * deg_cos ($ang );

#-----------------------------------------------------------------------
# 比例項の計算
#-----------------------------------------------------------------------
  $ang = NORMALIZATION_ANGLE(  481267.8809 * $t );
  $ang = NORMALIZATION_ANGLE(  $ang + 218.3162 );
  $th  = NORMALIZATION_ANGLE(  $th  +  $ang );

  return($th);
}

#=========================================================================
# 年月日、時分秒（世界時）からユリウス日（JD）を計算する
#
# ※ この関数では、グレゴリオ暦法による年月日から求めるものである。
#    （ユリウス暦法による年月日から求める場合には使用できない。）
#=========================================================================
sub YMDT2JD
{
	my ($year,$month,$day,$hour,$min,$sec) = @_;
	my ($jd,$t);

	if( $month < 3.0 ){
		$year -= 1.0;
		$month += 12.0;
	}

	$jd  = int( 365.25 * $year );
	$jd += int( $year / 400.0 );
	$jd -= int( $year / 100.0 );
	$jd += int( 30.59 * ( $month-2.0 ) );
	$jd += 1721088;
	$jd += $day;

	$t  = $sec / 3600.0;
	$t += $min /60.0;
	$t += $hour;
	$t  = $t / 24.0;

	$jd += $t;

	return( $jd );

}

#=========================================================================
# ユリウス日（JD）から年月日、時分秒（世界時）を計算する
#
# 戻り値の配列TIME[]の内訳
# TIME[0] ... 年  TIME[1] ... 月  TIME[2] ... 日
# TIME[3] ... 時  TIME[4] ... 分  TIME[5] ... 秒
#
# ※ この関数で求めた年月日は、グレゴリオ暦法によって表されている。
#
#=========================================================================
sub JD2YMDT
{

	my ($JD) = @_;
	my (@TIME,$x0,$x1,$x2,$x3,$x4,$x5,$x6,$tm);

	$x0 = int( $JD+68570.0);
	$x1 = int( $x0/36524.25 );
	$x2 = $x0 - int( 36524.25*$x1 + 0.75 );
	$x3 = int( ( $x2+1 )/365.2425 );
	$x4 = $x2 - int( 365.25*$x3 )+31.0;
	$x5 = int( int($x4) / 30.59 );
	$x6 = int( int($x5) / 11.0 );

	$TIME[2] = $x4 - int( 30.59*$x5 );
	$TIME[1] = $x5 - 12*$x6 + 2;
	$TIME[0] = 100*( $x1-49 ) + $x3 + $x6;

# 2月30日の補正
	if($TIME[1]==2 && $TIME[2] > 28){
		if($TIME[0] % 100 == 0 && $TIME[0] % 400 == 0){
			$TIME[2]=29;
		}elsif($TIME[0] % 4 ==0){
			$TIME[2]=29;
		}else{
			$TIME[2]=28;
		}
	}

	$tm=86400.0*( $JD - int( $JD ) );
	$TIME[3] = int( $tm/3600.0 );
	$TIME[4] = int( ($tm - 3600.0*$TIME[3])/60.0 );
	$TIME[5] = int( $tm - 3600.0*$TIME[3] - 60*$TIME[4] );

	return(@TIME);
}

#=========================================================================
# 今日が２４節気かどうか調べる
#
# 引数　 .... 計算対象となる年月日　$year $mon $day
#
# 戻り値 .... ２４節気の名称
#
#=========================================================================
sub check_24sekki
{
	my ($year,$mon,$day) = @_;
	my ($tm1,$tm2,$t,$rm_sun_today,$rm_sun_today0,$rm_sun_tommorow,$rm_sun_tommorow0);

#-----------------------------------------------------------------------
# ２４節気の定義
#-----------------------------------------------------------------------
	my (@sekki24) = ("春分","清明","穀雨","立夏","小満","芒種","夏至","小暑","大暑","立秋","処暑","白露",
	                   "秋分","寒露","霜降","立冬","小雪","大雪","冬至","小寒","大寒","立春","雨水","啓蟄");

	my $tm = YMDT2JD($year,$mon,$day,0,0,0);

#-----------------------------------------------------------------------
#時刻引数を分解する
#-----------------------------------------------------------------------
	$tm1 = int( $tm );
	$tm2 = $tm - $tm1;
	$tm2-=9.0/24.0;
	$t=($tm2+0.5) / 36525.0;
	$t=$t + ($tm1-2451545.0) / 36525.0;

	#今日の太陽の黄経
	$rm_sun_today = LONGITUDE_SUN( $t );

	$tm++;
	$tm1 = int($tm);
	$tm2 = $tm - $tm1;
	$tm2-=9.0/24.0;
	$t=($tm2+0.5) / 36525.0;
	$t=$t + ($tm1-2451545.0) / 36525.0;

	#明日の太陽の黄経
	$rm_sun_tommorow = LONGITUDE_SUN($t);

	#
	$rm_sun_today0   = 15.0 * int($rm_sun_today / 15.0);
	$rm_sun_tommorow0 = 15.0 * int($rm_sun_tommorow / 15.0);

	if($rm_sun_today0 != $rm_sun_tommorow0){
		return($sekki24[$rm_sun_tommorow0 / 15]);
	}else{
		return('');
	}
}

1;