寿星万年历---java算法实现

寿星万年历是我在网上见到的一份极高高精度的万年历，其采用先进的算法实现，其精度堪比刘安国教授为中国科学院国家授时中心制作的日梭万年历。但网络上只有javascript版本。于是自己将其翻译为java程序，并公布于此，方便大家使用。

 

寿星万年历相关信息：http://www.fjptsz.com/xxjs/xjw/rj/113.htm  设计:许剑伟

 

日梭万年历相关信息：http://www.time.ac.cn/calendar/calendar.htm 作者：刘安国

 

关于寿星万年历相关简要描述如下：

寿星万年历是一款采用现代天文算法制作的农历历算程序，含有公历与回历信息，可以很方便的进行公、农、回三历之间的转换。提供公元-4712年到公元9999年的日期查询功能。其中1500年到1940农历数据已经与陈垣的《二十史朔闰表》核对；含有从公420元(南北朝/宋武帝元年)到今的基本年号。在过去几百年中，寿星万年历的误差是非常小的，节气时刻计算及日月合朔时刻的平均误差小于1秒，太阳坐标的最大可能误差为0.2角秒，月亮坐标的最大可能误差为3角秒，平均误差为误差的1/6。万年历中含有几百个国内城市的经纬度，并且用户可根据自已的需要扩展经纬度数据。

 

 

 代码如下：

 

/**
 * @author lxslove
 * @mail moodlxs@163.com
 *
 */
public class SolarTerm {
 // ========角度变换===============
 private static final double rad = 180 * 3600 / Math.PI; // 每弧度的角秒数
 private static final double RAD = 180 / Math.PI; // 每弧度的角度数
 // ================日历计算===============
 private static final double J2000 = 2451545; // 2000年前儒略日数(2000-1-1
 // 12:00:00格林威治平时)

 // =========黄赤交角及黄赤坐标变换===========
 private static final double hcjjB[] = { 84381.448, -46.8150, -0.00059,
   0.001813 };// 黄赤交角系数表
 private static final double preceB[] = { 0, 50287.92262, 111.24406,
   0.07699, -0.23479, -0.00178, 0.00018, 0.00001 };// Date黄道上的岁差p

 private double Y = 2000;
 private double M = 1;
 private double D = 1;
 private double h = 12;
 private double m = 0;
 private double s = 0;

 private static final double[] dts = {
   // 世界时与原子时之差计算表
   -4000, 108371.7, -13036.80, 392.000, 0.0000, -500, 17201.0,
   -627.82, 16.170, -0.3413, -150, 12200.6, -346.41, 5.403, -0.1593,
   150, 9113.8, -328.13, -1.647, 0.0377, 500, 5707.5, -391.41, 0.915,
   0.3145, 900, 2203.4, -283.45, 13.034, -0.1778, 1300, 490.1, -57.35,
   2.085, -0.0072, 1600, 120.0, -9.81, -1.532, 0.1403, 1700, 10.2,
   -0.91, 0.510, -0.0370, 1800, 13.4, -0.72, 0.202, -0.0193, 1830,
   7.8, -1.81, 0.416, -0.0247, 1860, 8.3, -0.13, -0.406, 0.0292, 1880,
   -5.4, 0.32, -0.183, 0.0173, 1900, -2.3, 2.06, 0.169, -0.0135, 1920,
   21.2, 1.69, -0.304, 0.0167, 1940, 24.2, 1.22, -0.064, 0.0031, 1960,
   33.2, 0.51, 0.231, -0.0109, 1980, 51.0, 1.29, -0.026, 0.0032, 2000,
   64.7, -1.66, 5.224, -0.2905, 2150, 279.4, 732.95, 429.579, 0.0158,
   6000 };

 // 取整数部分
 public static double int2(double v) {
  v = Math.floor(v);
  if (v < 0)
   return v + 1;

  return v;
 }

 // 对超过0-2PI的角度转为0-2PI
 public static double rad2mrad(double v) {
  v = v % (2 * Math.PI);
  if (v < 0)
   return v + 2 * Math.PI;

  return v;
 }

 // 计算世界时与原子时之差,传入年
 public double deltatT(double y) {
  int i = 0;
  for (i = 0; i < 100; i += 5)
   if (y < dts[i + 5] || i == 95)
    break;

  double t1 = (y - dts[i]) / (dts[i + 5] - dts[i]) * 10;
  double t2 = t1 * t1;
  double t3 = t2 * t1;
  return dts[i + 1] + dts[i + 2] * t1 + dts[i + 3] * t2 + dts[i + 4] * t3;
 }

 // 传入儒略日(J2000起算),计算UTC与原子时的差(单位:日)
 public double deltatT2(double jd) {
  return this.deltatT(jd / 365.2425 + 2000) / 86400.0;
 }

 // 公历转儒略日,UTC=1表示原日期是UTC
 public double toJD(boolean UTC) {
  double y = this.Y; // 取出年月
  double m = this.M;
  double n = 0;

  if (m <= 2) {
   m += 12;
   y--;
  }

  if (this.Y * 372 + this.M * 31 + this.D >= 588829) {
   // 判断是否为格里高利历日1582*372+10*31+15
   n = int2(y / 100);
   n = 2 - n + int2(n / 4);// 加百年闰
  }

  n += int2(365.2500001 * (y + 4716)); // 加上年引起的偏移日数
  n += int2(30.6 * (m + 1)) + this.D; // 加上月引起的偏移日数及日偏移数
  n += ((this.s / 60 + this.m) / 60 + this.h) / 24 - 1524.5;
  if (UTC)
   return n + this.deltatT2(n - J2000);

  return n;
 }

 // 儒略日数转公历,UTC=1表示目标公历是UTC
 public void setFromJD(double jd, boolean UTC) {
  if (UTC)
   jd -= this.deltatT2(jd - J2000);

  jd += 0.5;

  // 取得日数的整数部份A及小数部分F
  double A = int2(jd);
  double F = jd - A;
  double D;

  if (A > 2299161) {
   D = int2((A - 1867216.25) / 36524.25);
   A += 1 + D - int2(D / 4);
  }
  A += 1524; // 向前移4年零2个月
  this.Y = int2((A - 122.1) / 365.25);// 年
  D = A - int2(365.25 * this.Y); // 去除整年日数后余下日数
  this.M = int2(D / 30.6001); // 月数
  this.D = D - int2(this.M * 30.6001);// 去除整月日数后余下日数
  this.Y -= 4716;
  this.M--;
  if (this.M > 12)
   this.M -= 12;
  if (this.M <= 2)
   this.Y++;
  // 日的小数转为时分秒
  F *= 24;
  this.h = int2(F);
  F -= this.h;
  F *= 60;
  this.m = int2(F);
  F -= this.m;
  F *= 60;
  this.s = F;
 }

 // 设置时间,参数例:"20000101 120000"或"20000101"
 public void setFromStr(String s) {
  this.Y = Double.parseDouble(s.substring(0, 4));
  this.M = Double.parseDouble(s.substring(4, 2));
  this.D = Double.parseDouble(s.substring(6, 2));
  this.h = Double.parseDouble(s.substring(9, 2));
  this.m = Double.parseDouble(s.substring(11, 2));
  this.s = Double.parseDouble(s.substring(13, 2)); /* 将5改为了2 */
 }

 // 日期转为串
 public String toStr() {
  String Y = "     " + (int)this.Y;
  String M = "0" + (int)this.M;
  String D = "0" + (int)this.D;

  double h = this.h, m = this.m, s = Math.floor(this.s + .5);
  if (s >= 60) {
   s -= 60;
   m++;
  }
  if (m >= 60) {
   m -= 60;
   h++;
  }
  String sh = "0" + (int)h, sm = "0" + (int)m, ss = "0" + (int)s;
  Y = Y.substring(Y.length() - 5, Y.length());
  M = M.substring(M.length() - 2, M.length());
  D = D.substring(D.length() - 2, D.length());
  sh = sh.substring(sh.length() - 2, sh.length());
  sm = sm.substring(sm.length() - 2, sm.length());
  ss = ss.substring(ss.length() - 2, ss.length());
  return Y + "-" + M + "-" + D + " " + sh + ":" + sm + ":" + ss;
 }

 // 算出:jd转到当地UTC后,UTC日数的整数部分或小数部分
 // 基于J2000力学时jd的起算点是12:00:00时,所以跳日时刻发生在12:00:00,这与日历计算发生矛盾
 // 把jd改正为00:00:00起算,这样儒略日的跳日动作就与日期的跳日同步
 // 改正方法为jd=jd+0.5-deltatT+shiqu/24
 // 把儒略日的起点移动-0.5(即前移12小时)
 // 式中shiqu是时区,北京的起算点是-8小时,shiqu取8
 public double Dint_dec(double jd, int shiqu, boolean dec) {
  double u = jd + 0.5 - this.deltatT2(jd) + shiqu / 24;
  if (dec)
   return Math.floor(u); // 返回整数部分
  else
   return u - Math.floor(u); // 返回小数部分
 }

 // 计算两个日期的相差的天数,输入字串格式日期,如:"20080101"
 double d1_d2(String d1, String d2) {
  double Y = this.Y, M = this.M, D = this.D, h = this.h, m = this.m, s = this.s; // 备份原来的数据

  this.setFromStr(d1.substring(0, 8) + " 120000");
  double jd1 = this.toJD(false);
  this.setFromStr(d2.substring(0, 8) + " 120000");
  double jd2 = this.toJD(false);

  this.Y = Y;
  this.M = M;
  this.D = D;
  this.h = h;
  this.m = m;
  this.s = s; // 还原
  if (jd1 > jd2)
   return Math.floor(jd1 - jd2 + .0001);
  else
   return -Math.floor(jd2 - jd1 + .0001);
 }

 // 返回黄赤交角(常规精度),短期精度很高
 public static double hcjj1(double t) {
  double t1 = t / 36525;
  double t2 = t1 * t1;
  double t3 = t2 * t1;
  return (hcjjB[0] + hcjjB[1] * t1 + hcjjB[2] * t2 + hcjjB[3] * t3) / rad;
 }

 // 黄赤转换(黄赤坐标旋转)
 public static void HCconv(double[] JW, double E) {
  // 黄道赤道坐标变换,赤到黄E取负
  double HJ = rad2mrad(JW[0]), HW = JW[1];
  double sinE = Math.sin(E), cosE = Math.cos(E);
  double sinW = cosE * Math.sin(HW) + sinE * Math.cos(HW) * Math.sin(HJ);
  double J = Math.atan2(Math.sin(HJ) * cosE - Math.tan(HW) * sinE, Math
    .cos(HJ));
  JW[0] = rad2mrad(J);
  JW[1] = Math.asin(sinW);
 }

 // 补岁差
 public static void addPrece(double jd, double[] zb) {
  int i;
  double t = 1, v = 0, t1 = jd / 365250;
  for (i = 1; i < 8; i++) {
   t *= t1;
   v += preceB[i] * t;
  }
  zb[0] = rad2mrad(zb[0] + (v + 2.9965 * t1) / rad);
 }

 // ===============光行差==================
 private static final double GXC_e[] = { 0.016708634, -0.000042037,
   -0.0000001267 }; // 离心率
 private static final double GXC_p[] = { 102.93735 / RAD, 1.71946 / RAD,
   0.00046 / RAD }; // 近点
 private static final double GXC_l[] = { 280.4664567 / RAD,
   36000.76982779 / R