/*
 * calendar.c -- 日付関係
 * Copyright (C) 1989, 1990, 2006 by Candy
 */
/*
 * $Id: calendar.c,v 2.1 1992/06/05 15:47:10 candy Exp candy $
 */
#include "calendar.h"
#ifndef __GNUC__
#define inline
#endif
/*
 * 西暦 year 年が、閏年なら non 0 を返す。そうでなければ 0 を返す。
 * year が 400 の倍数であるか
 * year が 4 の倍数であるが 100 の倍数でない
 * ならば non 0 を返す。
 */
extern int
isleap(int year)
{
	return (year % 4 == 0 && year % 100 != 0) || year % 400 == 0;
}/* isleap */

/*
 * 西暦 year 年 month 月の日数を返す。閏年の 2 月ならば 29 を返す。
 */
extern int
day_of_manth(int year, int month)
{
	static int dom[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
	int days;
	if (month == 2 && isleap(year))
		days = dom[2 - 1] + 1;
	else
		days = dom[month - 1];
	return days;
}/* day_of_manth */

/*
 * 日付と基準日からの日数との変換を行います。
 * 基準日からの日数は、西暦 org_year 年 1 月 1 日を 0 日とします。
 * 結果が負になるような場合も正しく返します。
 *
 * Note:
 *
 * org_year について
 * org_year は 400 の倍数でなければなりません。
 *
 * 内部処理の単位
 * Y 年 3 月〜 12 月、Y + 1 年 1 月〜 2 月を
 * 内部では Y 年' 0 月'〜 11 月' として処理します。
 * d 日は内部では d - 1 日' とします。
 * すると閏日は 11 月' 28 日' になります。
 * 日数も内部では 0 年' 0 月' 0 日' を 0 とします。
 *
 * sod[i] について
 * sod[i] は i 月' 0 日' の 0 月' 0 日' からの日数です。
 *
 *     3   4   5   6   7   8   9  10  11  12   1   2 [月]
 * dom[12] = {31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31, 28}; とすると、
 *
 * sod[i] = { 0 (i = 0 の時)
 *            sod[i - 1] + dom[i - 1] (i > 0 の時) }
 *
 * と定義されます。
 *
 * 計算
 * days = "0年'0月'0日'から yy 年' mm 月' dd 日'までの日数"
 * uu = "0年'0月'0日'から yy 年'0月'0日'迄に閏日があった回数"
 * とすると、
 * days = yy * 365 + sod[mm] + dd + uu
 * となります。
 *
 * オフセットについて
 * ・org_year 年1月1日=-1年'10月'0日'
 * ・0年'0月'0日'から -1年'10月'0日'迄には-1回の閏日があった。
 * だか
 * = -1 * 365 + sod[10] +  0 + (-1)
 * = -1 * 365 + 306 + 0 - 1
 * = -60
 * 即ち、-60がその値です。
 *
 */
#define OFFSET (-60)
static int org_year = 1600;
static int sod[12] = { 0, 31, 61, 92,122,153,184,214,245,275,306,337};

/* day 値 - 修正ユリウス日(MJD) */
#define DAY_MINUS_MJD (94553)
#define AJD_MINUS_MJD (2400000.5)

/*
 * AJD = MJD + 2494553.5
 *
 * フリーゲルの公式
 * MJD = floor(365.25 * Y) + floor(Y / 400) - floor(Y / 100)
 *	+ floor(30.59 * (m - 2)) + d - 678912
 * 1 月 2 月は前年の 13 月 14 月とする。d = 1..31
 */

/*
 * 除算
 * dst / src を越えない最大の整数を返す。
 * 1990/03/25 の gnulib.a の ___modsi3 がバグってて、
 * dst < 0 の時バグります。あ〜あ。
 */
static inline long
my_divl(long dst, long src)
{
	long q, r;
	q = dst / src;
	r = dst % src;
	if (r < 0) {
		r = r + src;
		q = q - 1;
	}
	return q;
}/* my_divl */

/*
 * 基準日から西暦 year 年 month 月 day 日までの日数を返します。
 * month = (1 .. 12);
 * day = (1 .. 31);
 */
extern long
ymd_to_mjd(int year, int month, int day)
{
	int yy, mm, dd, uu;
	long mjd;
	if (month >= 3) {
		yy = year - org_year;
		mm = month - 3;
	}
	else {
		yy = year - org_year - 1;
		mm = month - 3 + 12;
	}
	dd = day - 1;
	if (yy < 0)
		uu = my_divl(yy, 4) - my_divl(yy, 100) + my_divl(yy, 400);
	else
		uu = yy / 4 - yy / 100 + yy / 400;
	mjd = (long)yy * 365 + sod[mm] + dd + uu - OFFSET - DAY_MINUS_MJD;
	return mjd;
}/* ymd_to_mjd */

/*
 * ymd_to_mjd() で求めた基準日からの日数 date を西暦年月日に変換し、
 * それぞれ *yearp, *monthp, *dayp に返します。
 *
 *Note:
 * 年' = y400 * 400 + y100 * 100 + y4 * 4 + y1
 * の形で表します。
 *
 * 年'  * y100 y4     y1
 * 0..2 00..23 0..2  365 * 3 * 24 * 3 = 78840
 * 0..2 00..23    3  366     * 24 * 3 = 26352
 * 0..2     24 0..2  365 * 3      * 3 =  3285
 * 0..2     24    3  365          * 3 =  1095
 *    3 00..23 0..2  365 * 3 * 24     = 26280
 *    3 00..23    3  366     * 24     =  8784
 *    3     99 0..2  365 * 3          =  1095
 *    3     99    3  366              =   366
 *
 * D1   =                   365   0..0 年'の日数
 * D4   = D1 * 4 + 1   =   1461   0..3 年'の日数
 * D100 = D4 * 25 - 1  =  36524  0..99 年'の日数
 * D400 = D100 * 4 + 1 = 146097 0..399 年'の日数
 *
 */
#define D1 365
#define D4 1461
#define D100 36524
#define D400 146097

extern void
mjd_to_ymd(long mjd, int *yearp, int *monthp, int *dayp)
{
	int y400, y100, y4, y1;
	long r400, r100, r4, r1;
	int yy, mm, dd;
	mjd += OFFSET + DAY_MINUS_MJD;
	if (mjd < 0)
		y400 = my_divl(mjd, D400);
	else
		y400 = mjd / D400;
	r400 = mjd - y400 * D400; /* r400 = (0 .. D400 - 2, D400 - 1) */
	y100 = r400 / D100; /* y100 = (0 .. 3, 4) */
	r100 = r400 % D100; /* r100 = (0 .. D100 - 1) */
	y4 = r100 / D4; /* y4 = (0 .. 24) */
	r4 = r100 % D4; /* r4 = (0 .. D4 - 2, D4 - 1) */
	y1 = r4 / D1; /* y1 = (0 .. 3, 4) */
	r1 = r4 % D1; /* r1 = (0 .. D1 - 1) */
	yy = ((y400 * 4 + y100) * 25 + y4) * 4 + y1;
	mm = 1;
	while (mm < 12 && sod[mm] <= r1)
		mm++;
	mm--;
	dd = r1 - sod[mm];
	/*
	 * r400 == D400 - 1 の時(400年に一度の閏日)
	 * r4 == D4 - 1 の時(4年に一度の閏日)
	 * の場合の修	 */
	if (y100 == 4 || y1 == 4) {
		yy--;
		mm = 11;
		dd = 28;
	}
	if (mm < 10) {
		*yearp = yy + org_year;
		*monthp = mm + 3;
	}
	else {
		*yearp = yy + org_year + 1;
		*monthp = mm + 3 - 12;
	}
	*dayp = dd + 1;
}/* mjd_to_ymd */

/*
 * ymd_to_mjd() で求めた基準日からの日数 date の曜日を返す。
 * 日曜日なら0、土曜日なら6を返す。
 * 400年間の日数は 146097 日で7の倍数ですから、ありがたいこと
 *
 */
extern int
week_of_mjd(long mjd)
{
	long q;
	int r;
	static char wod[7] = {6, 0, 1, 2, 3, 4, 5};
	mjd += DAY_MINUS_MJD;
	if (mjd < 0)
		q = my_divl(mjd, 7L);
	else
		q = mjd / 7;
	r = mjd - q * 7;
	return wod[r];
}/* week_of_mjd */

extern long
mjd_to_ajd(long mjd)
{
	return mjd + (AJD_MINUS_MJD + 0.5);
}/* mjd_to_ajd */

extern long
ajd_to_mjd(long ajd)
{
	return ajd - (AJD_MINUS_MJD + 0.5);
}/* ajd_to_mjd */