万年历 - java算法实现(5)

14286.1503796870,-6.0860358E-05, 6.2714140E-08,-1.9984990E-10}; private static final double M1n[] = { 3.81034392032, 8.39968473021E+03,-3.31919929753E-05, //月球平黄经系数

3.20170955005E-08,-1.53637455544E-10};

// ==================日位置计算=================== private double EnnT = 0; // 调用Enn前先设置EnnT时间变量

// 计算E10,E11,E20等,即:某一组周期项或泊松项算出,计算前先设置EnnT时间 public double Enn(double[] F) { double v = 0;

for (int i = 0; i < F.length; i += 3)

v += F[i] * Math.cos(F[i + 1] + EnnT * F[i + 2]); return v; }

// 返回地球位置,日心Date黄道分点坐标 public double[] earCal(double jd) { EnnT = jd / 365250; double llr[] = new double[3];

double t1 = EnnT, t2 = t1 * t1, t3 = t2 * t1, t4 = t3 * t1, t5 = t4 * t1;

llr[0] = Enn(E10) + Enn(E11) * t1 + Enn(E12) * t2 + Enn(E13) * t3 + Enn(E14) * t4 + Enn(E15) * t5; llr[1] = Enn(E20) + Enn(E21) * t1;

llr[2] = Enn(E30) + Enn(E31) * t1 + Enn(E32) * t2 + Enn(E33) * t3; llr[0] = rad2mrad(llr[0]); return llr; }

// 传回jd时刻太阳的地心视黄经及黄纬 public double[] sunCal2(double jd) { double[] sun = earCal(jd); sun[0] += Math.PI;

sun[1] = -sun[1]; // 计算太阳真位置 ZD d = nutation(jd);

sun[0] = rad2mrad(sun[0] + d.Lon); // 补章动 addGxc(jd, sun); // 补周年黄经光行差 return sun; // 返回太阳视位置 }

// ==================月位置计算=================== private double MnnT = 0; // 调用Mnn前先设置MnnT时间变量 // 计算M10,M11,M20等,计算前先设置MnnT时间 public double Mnn(double[] F) {

double v = 0, t1 = MnnT, t2 = t1 * t1, t3 = t2 * t1, t4 = t3 * t1; for (int i = 0; i < F.length; i += 6) v += F[i]

* Math.sin(F[i + 1] + t1 * F[i + 2] + t2 * F[i + 3] + t3 * F[i + 4] + t4 * F[i + 5]); return v; }

// 返回月球位置,返回地心Date黄道坐标 public double[] moonCal(double jd) { MnnT = jd / 36525;

double t1 = MnnT, t2 = t1 * t1, t3 = t2 * t1, t4 = t3 * t1; double[] llr = new double[3];

llr[0] = (Mnn(M10) + Mnn(M11) * t1 + Mnn(M12) * t2) / rad; llr[1] = (Mnn(M20) + Mnn(M21) * t1) / rad;

llr[2] = (Mnn(M30) + Mnn(M31) * t1) * 0.999999949827; llr[0] = llr[0] + M1n[0] + M1n[1] * t1 + M1n[2] * t2 + M1n[3] * t3 + M1n[4] * t4;

llr[0] = rad2mrad(llr[0]); // 地心Date黄道原点坐标(不含岁差) addPrece(jd, llr); // 补岁差 return llr; }

// 传回月球的地心视黄经及视黄纬 public double[] moonCal2(double jd) {

double[] moon = moonCal(jd); ZD d = nutation(jd);

moon[0] = rad2mrad(moon[0] + d.Lon); // 补章动 return moon; }

// 传回月球的地心视赤经及视赤纬 public double[] moonCal3(double jd) { double[] moon = moonCal(jd); HCconv(moon, hcjj1(jd));

nutationRaDec(jd, moon); // 补赤经及赤纬章动

// 如果黄赤转换前补了黄经章动及交章动,就不能再补赤经赤纬章动 return moon; }

// ==================地心坐标中的日月位置计算=================== public double jiaoCai(int lx, double t, double jiao) {

// lx=1时计算t时刻日月角距与jiao的差, lx=0计算t时刻太阳黄经与jiao的差 double[] sun = earCal(t); // 计算太阳真位置(先算出日心坐标中地球的位置) sun[0] += Math.PI;

sun[1] = -sun[1]; // 转为地心坐标 addGxc(t, sun); // 补周年光行差 if (lx == 0) { ZD d = nutation(t);

sun[0] += d.Lon; // 补黄经章动 return rad2mrad(jiao - sun[0]); }

double[] moon = moonCal(t); // 日月角差与章动无关 return rad2mrad(jiao - (moon[0] - sun[0])); }

// ==================已知位置反求时间=================== public double jiaoCal(double t1, double jiao, int lx) { // t1是J2000起算儒略日数 // 已知角度(jiao)求时间(t)

// lx=0是太阳黄经达某角度的时刻计算(用于节气计算)

// lx=1是日月角距达某角度的时刻计算(用于定朔望等) // 传入的t1是指定角度对应真时刻t的前一些天

// 对于节气计算,应满足t在t1到t1+360天之间,对于Y年第n个节气(n=0是春分),t1可取值Y*365.2422+n*15.2

// 对于朔望计算,应满足t在t1到t1+25天之间,在此范围之外,求右边的根 double t2 = t1, t = 0, v; if (lx == 0)

t2 += 360; // 在t1到t2范围内求解(范气360天范围),结果置于t else t2 += 25;

jiao *= Math.PI / 180; // 待搜索目标角 // 利用截弦法计算

double v1 = jiaoCai(lx, t1, jiao); // v1,v2为t1,t2时对应的黄经 double v2 = jiaoCai(lx, t2, jiao); if (v1 < v2)

v2 -= 2 * Math.PI; // 减2pi作用是将周期性角度转为连续角度 double k = 1, k2; // k是截弦的斜率

for (int i = 0; i < 10; i++) { // 快速截弦求根,通常截弦三四次就已达所需精度 k2 = (v2 - v1) / (t2 - t1); // 算出斜率 if (Math.abs(k2) > 1e-15) k = k2; // 差商可能为零,应排除 t = t1 - v1 / k;

v = jiaoCai(lx, t, jiao);// 直线逼近法求根(直线方程的根) if (v > 1)

v -= 2 * Math.PI; // 一次逼近后,v1就已接近0,如果很大,则应减1周 if (Math.abs(v) < 1e-8) break; // 已达精度 t1 = t2; v1 = v2; t2 = t;

v2 = v; // 下一次截弦 } return t; }

//==================节气计算=================== public static final String jqB[] = { //节气表

\春分\清明\谷雨\立夏\小满\芒种\夏至\小暑\大暑\立秋\处暑\白露\ \秋分\寒露\霜降\立冬\小雪\大雪\冬至\小寒\大寒\立春\雨水\惊蛰\ public void JQtest(int y) { // 节气使计算范例,y是年分,这是个测试函数 double jd = 365.2422 * (y - 2000), q; String s1, s2;

for (int i = 0; i < 24; i++) { q = jiaoCal(jd + i * 15.2, i * 15, 0);

q = q + J2000 + (double)8 / 24; // 计算第i个节气(i=0是春风),结果转为北京时 setFromJD(q, true);

s1 = toStr(); // 将儒略日转成世界时 setFromJD(q, false);

s2 = toStr(); // 将儒略日转成日期格式(输出日期形式的力学时) System.out.println(jqB[i] + \显示 } }

// =================定朔弦望计算======================== public void dingSuo(int y, double arc) { // 这是个测试函数 double jd = 365.2422 * (y - 2000), q; String s1, s2;

System.out.println(\月份:世界时 原子时\ for (int i = 0; i < 12; i++) {

q = jiaoCal(jd + 29.5 * i, arc, 1) + J2000 + 8 / 24; // 计算第i个节气(i=0是春风),结果转为北京时

setFromJD(q, true);

s1 = toStr(); // 将儒略日转成世界时 setFromJD(q, false);

s2 = toStr(); // 将儒略日转成日期格式(输出日期形式的力学时) System.out.println((i + 1) + \月 : \显示 } }

//=================农历计算======================== /*****