《Numerical Recipes》给的第一个程序是计算月相、儒略日和普通日之间的互相转换这样的。但貌似以前都居然不知道有儒略日这个东西，这样太伤我心了。查一下：

儒略日是啥呢？

儒略日一种不用年和月却只用日的纪年法。它的起点，也就是（儒略日=0）那天是B.C.4713年1月1日12:00到B.C.4713年1月2日12:00（格林威治时间）

对于计算机的好处就是它可以直接用int储存。

书中的第一个例程是void flmoon(const int n, const int phase, int * const julian_day, double * const frac);

它的作用是传递给它n和phase作为输入，其中phase是你想要计算的月相，n是让它计算出（从1900.01开始的）第n个phase月相的儒略日。

julian_day作为儒略日的输出，frac是那天的一个校准值，可以得到什么时候出现的那种月相，也就是几点几分这样的。

比较悲剧的是，它并没有告诉我解释公式是怎么回事，也不告诉我它是怎么做到的。大家先来看一下：

/**
*@file calendar.h
*/
#ifndef CALENDAR_H_INCLUDED
#define CALENDAR_H_INCLUDED
/**
*@brief	This routine calculates the phases of the moon.
*@param	n	Input for n-th such phase to be calculated
*@param	phase	Input for the phase to be calculated (phase = 0 for new moon, 1 for first quarter, 2 for full, 3 for last quarter)
*@param	julian_day	Output the Julian Day Number of that day
*@param	frac	Output the fractional part of that day
*Introduction:
*Given an integer n and a code phase for the phase desired
*(phase = 0 for new moon, 1 for first quarter, 2 for full, 3 for last quarter),
*the routine returns the Julian Day Number julian_day, and the fractional part of a day frac to be
*added to it, of the n-th such phase since January, 1900. Greenwich Mean Time is assumed.
*/
void flmoon(const int n, const int phase, int * const julian_day, double * const frac){
	const double RAD = 3.141592653589793238/180.0;
	int i;
	double am, as, c, t, t2, xtra;
	c = n + phase/4.0;
	t = c/1236.85;
	t2 = t*t;
	as = 359.2242 + 29.105356 * c;
	am = 306.0253 + 385.816918 * c + 0.010730 * t2;
	*julian_day = 2415020 + 28 * n + 7 * phase;
	xtra = 0.75933 + 1.53058868 * c + ((1.178e-4)-(1.55e-7) * t) * t2;
	if(phase == 0 || phase == 2)
		xtra += (0.1734 - 3.93e-4 * t) * sin(RAD * as) - 0.4068 * sin(RAD * am);
	else if(phase == 1 || phase == 3)
		xtra += (0.1721 - 4.0e-4 * t) * sin(RAD * as) - 0.6280 * sin(RAD * am);
	else printf("phase is unknown in flmoon\n");
	i = (int)(xtra >= 0.0 ? floor(xtra) : ceil(xtra - 1.0));
	*julian_day += i;
	*frac = xtra - i;
}

第二个例程是int julday(const int month, const int day, const int year);

它的作用是把普通的日期转换到儒略日去。

值得一提的是，数字365.25是四年一闰年的平均每年有的天数，但30.6001却让我难以理解、缺少指明性。根据第二个参考链接，可能与月份的计数方式与蔡勒公式相同所致。

/**
*@brief	Convert calendar date into Julian Day Number
*@param	month	the month of the calendar date
*@param	day	the day of the calendar date
*@param	year	the year of the calendar date
*@return	the corresponding Julian Day
*In this routine julday returns the Julian Day Number that begins
*at noon of the calendar date specified by month, day, and year,
*all integer variables. Positive year signifies A.D.; negative, B.C.
*@attention Remember that the year after 1 B.C. was 1 A.D.
*/
int julday(const int month, const int day, const int year){
	const int IGREG = 15+31*(10+12*1582); //Gregorian Calendat adopted Oct. 15, 1582.
	int ja, jul, jy = year, jm;
	if(jy == 0)
		printf("julday: there is no year zero.");
	if(jy < 0)
		++jy;
	if(month > 2){
		jm = month+1;
	} else{
		--jy;
		jm = month+13;
	}
	jul = (int)(floor(365.25*jy) + floor(30.6001*jm) + day + 1720995);
	if(day + 31*(month+12*year) >= IGREG){
		//Test whether to change to Gregorian Calendar
		ja = (int)(0.01*jy);
		jul += 2 - ja + (int)(0.25*ja);
	}
	return jul;
}

第三个例程是第二个例程的反函数。将儒略日转换回普通日期。

/**
*@brief	Convert Julian Day into Calendar Date
*@param	julian	Input for the Julian Day
*@param	month	Output the corresponding month
*@param	day	Output the corresponding day
*@param	year	Output the corresponding year
*Inverse of the function julday given above. Here julian is input as a Julian Day number,
*and the routine outputs month, day, and year as the month, day, and year on which the
*specified Julian Day stared at noon.
*/
void caldat(const int julian, int * const month, int * const day, int * const year){
	const int IGREG = 2299161;
	int ja, jalpha, jb, jc, julian_day, je;

	if(julian >= IGREG){
		//Cross-over to Gregorian Calendar produces this correction.
		jalpha = (int)((((double)(julian - 1867216)) - 0.25)/36524.25);
		ja = julian + 1 + jalpha - (int)(0.25*jalpha);
	} else if (julian < 0){
		/*
		*Make day number positive by adding integer number of Julian centuries,
		*then subtract them off at the end.
		*/
		ja = julian + 36525 * (1 - julian/36525);
	} else {
		ja = julian;
	}
	jb = ja + 1524;
	jc = (int)(6680.0 + (((double)(jb-2439870))-122.1)/365.25);
	julian_day = (int)(365 * jc + (0.25*jc));
	je = (int)((jb - julian_day)/30.6001);
	*day = jb - julian_day - (int)(30.6001*je);
	*month = je - 1;
	if(*month > 12)
		*month -= 12;
	*year = jc - 4715;
	if(*month > 2)
		--(*year);
	if(julian < 0)
		*year -= 100 * (1 - julian / 36525);
}

#endif // CALENDAR_H_INCLUDED

原书上的代码需要包含nr3.h，这样的话就必须要拿到他的nr3.h的头文件才行，也就是得买他的光盘了（有点小舍不得），所以呢，就将其中会用到的它自己实现的特性还原回了C的代码。这样的话，至少大家不会因为调试不出而直接失去对它的兴趣。

参考地址：

http://zh.wikipedia.org/wiki/%E5%84%92%E7%95%A5%E6%97%A5

http://www.fjptsz.com/xxjs/xjw/rj/117/07.htm

• Relative Posts
• Global Relative Posts