日の出計算 その１ - その２があるかはわかりませんが.
2012/10/21
java
計算
天文
日の出
Javaです。
とりあえずそれっぽい感じできたので。
コードはかなり意味不明です。
あとから直すということで。
インターネットでいろいろあさるといろんな人が日の出計算をおこなっているようですね。
もとネタはこの本ですが、正直理解していないので、ぼくのコードはかなりあやういです。
日の出・日の入りの計算―天体の出没時刻の求め方
長沢 工
import static java.lang.Math.*;
public class A20121021 {
static class YMDH {
int Y, M, D, H;
public YMDH(int yyyy, int mm, int dd, int hh) {
this.Y = yyyy;
this.M = mm;
this.D = dd;
this.H = hh;
}
public double d() {
return H / 24.;
}
public String toString() {
StringBuilder builder = new StringBuilder();
builder.append("y:[").append(Y).append("] ").append("m:[")
.append(M).append("] ").append("d:[").append(D)
.append("] ").append("h:[").append(H).append("] ")
.append("d[").append(d()).append("] ");
return builder.toString();
}
}
static class Degree {
double deg;
public Degree(double argdeg) {
this.deg = argdeg;
}
public double toR() {
return toRadians(deg);
}
public String toString() {
StringBuilder builder = new StringBuilder();
builder.append("deg:[").append(deg).append("] ");
return builder.toString();
}
}
static class LongitudeAndLatitude {
/** 緯度 */
Degree latitude;
/** 経度 */
Degree longitude;
public LongitudeAndLatitude(double a, double b) {
this.latitude = new Degree(a);
this.longitude = new Degree(b);
}
public String toString() {
StringBuilder builder = new StringBuilder();
builder.append("lat:[").append(latitude.deg).append("] ")
.append("lng[").append(longitude.deg).append("]");
return builder.toString();
}
}
static class Sun {
/** a:赤経 */
Degree a;
/** δ:赤緯 デルタ */
Degree δ;
/** r:距離 */
double r;
/** th:恒星時 θ シータ */
Degree th;
public String toString() {
StringBuilder builder = new StringBuilder();
builder.append("a:[").append(a).append("] ").append("δ:[")
.append(δ).append("] ").append("r:[").append(r)
.append("] ").append("th:[").append(th).append("] ");
return builder.toString();
}
}
public static void main(String[] args) {
// for(int i=1;i<=12;i++){
// for(int j=1;j<=31;j++){
// tokyoSunrise(2012, i, j);
// }
// }
sunrise("", 2012, 7, 9, 31.36, 130.33);
sunrise("", 2012, 7, 9, 32.48, 130.43);
sunrise("", 2012, 7, 9, 31.54, 131.25);
sunrise("", 2012, 7, 9, 33.14, 131.37);
sunrise("根室", 2012, 7, 9, 43.20, 145.35);
}
static void sunrise(String location, int Y, int M, int D, double longitude,
double latitude) {
YMDH ymdh = new YMDH(Y, M, D, 6);
LongitudeAndLatitude longitudeAndLatitude = new LongitudeAndLatitude(
longitude, latitude);
double d = b(ymdh, longitudeAndLatitude);
// System.out.println(d);
d = d * 24;
double d2 = d % 1;
// System.out.println(d);
// System.out.println(d2*60);
System.out.println(location + "\t" + ymdh.Y + "年" + ymdh.M + "月"
+ ymdh.D + "日" + (int) floor(d) + "時" + (int) floor(d2 * 60)
+ "分 \t" + (int) floor(d) + " " + (int) floor(d2 * 60));
}
static double b(YMDH ymdh, LongitudeAndLatitude longitudeAndLatitude) {
Sun sun = a2(ymdh.Y, ymdh.M, ymdh.D, ymdh.H,
longitudeAndLatitude.longitude.deg);
// (3) 出没高度k
double k = k(sun.r);
// (4)
double tk = tk(sun.δ.deg, longitudeAndLatitude.latitude.deg, k);
tk = toDegrees(acos(tk));
double t = t(sun.a.deg, sun.th.deg);
if (0 < tk)
tk *= -1;
// 補正値Δd
double dd = deltaD(tk, t);
// System.out.println(dd);
double d = ymdh.d() + dd;// 補正値
// System.out.println("補正値:" + d);
return d;
}
static Sun a2(int Y, int M, int D, int H, double λ2) {
{
// T 時刻変数
double T = 時刻変数T(Y, M, D);
double λs = λs(T);
double q = q(T);
double r = r(q);
// double a = toDegrees(λs);
double e = toRadians(23.439291) - toRadians(0.000130042) * T;
double tana = tan(λs) * cos(e);
double sinδ = sin(λs) * sin(e);
double _λs = toDegrees(λs) % 360.;
tana = atan(tana);
if (0 <= _λs && _λs < 90) {
// System.out.print((int)_λs + " bingo1 ***\t");
tana += toRadians(360);
} else if (90 <= _λs && _λs < 180) {
// System.out.print((int)_λs + " bingo2 ***\t");
tana += toRadians(180);
} else if (180 <= _λs && _λs < 270) {
tana += toRadians(180);
// System.out.print((int)_λs + " bingo3 ***\t");
} else if (270 <= _λs && _λs < 360) {
tana += toRadians(360);
// System.out.print((int)_λs + " bingo4 ***\t");
}
// System.out.print((int)toDegrees(tana)+ "\t");
sinδ = asin(sinδ);
// d:時刻変数
double d = H / 24.;
// 恒星時を求める
double th = siderealTimeθ(T, d, λ2);
Sun sun = new Sun();
{
// System.out.println("tana:" + tana);
// System.out.println("sinδ:" + sinδ);
sun.a = new Degree(toDegrees(tana));
sun.δ = new Degree(toDegrees(sinδ));
sun.r = r;
sun.th = new Degree(th);
}
return sun;
}
}
/**
* 仮定時刻に対する補正値Δdの計算
*
* @param r
* @return
*/
static double deltaD(double tk, double t) {
double deltaD = 0;
deltaD = (tk - t) / 360.;
return deltaD;
}
/**
* 太陽の出没高度kを計算する
*
* @param r
* @return
*/
static double k(double r) {
// 太陽の視半径
final double S = 0.266994 / r;
// 大気差
final double R = 0.585556;
// 視差
final double Π = 0.0024428 / r;
final double E0 = 0;
double k = -S - E0 - R + Π;
return k;
}
/**
* 時角を求める φ
*
* @return
*/
static double tk(double delta, double phi, double k) {
double tk = 0;
double sink = sin(toRadians(k));
double cosDelta = cos(toRadians(delta));
double sinDelta = sin(toRadians(delta));
double cosPhi = cos(toRadians(phi));
double sinPhi = sin(toRadians(phi));
tk = (sink - sinDelta * sinPhi) / (cosDelta * cosPhi);
return tk;
}
static double tk(Degree delta, Degree phi, Degree k) {
double tk = 0;
double sink = sin(k.toR());
double cosDelta = cos(delta.toR());
double sinDelta = sin(delta.toR());
double cosPhi = cos(phi.toR());
double sinPhi = sin(phi.toR());
tk = (sink - sinDelta * sinPhi) / (cosDelta * cosPhi);
return tk;
}
/**
* 恒星時から太陽の時角tを計算する
*
* @param delta
* @param phi
* @param k
* @return
*/
static double t(double a, double th) {
double t = 0;
t = th - a;
return t;
}
static double b(double b, double bb) {
double a = 60;
double aa = 60 * 60;
return b / a + bb / aa;
}
static double r(double q) {
return pow(10., q);
}
/**
*
* @param T
* @param d
* @param λ
* @return
*/
static double siderealTimeθ(double T, double d, double λ) {
double th = 0;
th = 325.4606 + 360.007700536 * T + 0.00000003879 * pow(T, 2) + 360.
* d + λ;
double th2 = th;
if (th2 < 0) {
while (th2 < 0) {
th2 += 360;
}
} else if (360 <= th2) {
while (360 <= th2) {
th2 -= 360;
}
}
return th2;
}
static double q(double T) {
double[] L = new double[6];
L[0] = (0.007256 - 0.0000002 * T)
* sin(toRadians(267.54) + toRadians(359.991) * T);
L[1] = (0.000091) * sin(toRadians(265.1) + toRadians(719.98) * T);
L[2] = (0.000030) * sin(toRadians(90.0));
L[3] = (0.000013) * sin(toRadians(27.8) + toRadians(4452.67) * T);
L[4] = (0.000007) * sin(toRadians(254) + toRadians(450.4) * T);
L[5] = (0.000007) * sin(toRadians(156) + toRadians(329.6) * T);
double q = 0;
for (double d : L) {
q += d;
}
return q;
}
static double λs(double T) {
double[] L = new double[20];
L[0] = toRadians(280.4603) + toRadians(360.00769) * T;
L[1] = (toRadians(1.9146) - toRadians(0.00005))
* sin(toRadians(357.538) + toRadians(359.991) * T);
L[2] = (toRadians(0.0200))
* sin(toRadians(355.05) + toRadians(719.981) * T);
L[3] = (toRadians(0.0048))
* sin(toRadians(234.95) + toRadians(19.341) * T);
L[4] = (toRadians(0.0020))
* sin(toRadians(247.1) + toRadians(329.64) * T);
L[5] = (toRadians(0.0018))
* sin(toRadians(297.8) + toRadians(4452.67) * T);
L[6] = (toRadians(0.0018))
* sin(toRadians(251.3) + toRadians(0.20) * T);
L[7] = (toRadians(0.0015))
* sin(toRadians(343.2) + toRadians(450.37) * T);
L[8] = (toRadians(0.0013))
* sin(toRadians(81.4) + toRadians(225.18) * T);
L[9] = (toRadians(0.0008))
* sin(toRadians(132.5) + toRadians(659.29) * T);
L[10] = (toRadians(0.0007))
* sin(toRadians(153.3) + toRadians(90.38) * T);
L[11] = (toRadians(0.0007))
* sin(toRadians(206.8) + toRadians(30.35) * T);
L[12] = (toRadians(0.0006))
* sin(toRadians(29.8) + toRadians(337.18) * T);
L[13] = (toRadians(0.0005))
* sin(toRadians(207.4) + toRadians(1.50) * T);
L[14] = (toRadians(0.0005))
* sin(toRadians(291.2) + toRadians(22.81) * T);
L[15] = (toRadians(0.0004))
* sin(toRadians(234.9) + toRadians(315.56) * T);
L[16] = (toRadians(0.0004))
* sin(toRadians(157.3) + toRadians(299.30) * T);
L[17] = (toRadians(0.0004))
* sin(toRadians(21.1) + toRadians(720.02) * T);
L[18] = (toRadians(0.0003))
* sin(toRadians(352.5) + toRadians(1079.97) * T);
L[19] = (toRadians(0.0003))
* sin(toRadians(329.7) + toRadians(44.43) * T);
double λs = 0;
for (double d : L) {
λs += d;
}
return λs;
}
static double 時刻変数T(int Y, int M, int D) {
double deltaT = deltaT(Y, M);
double T = 0.;
double K = 経過日数K(Y, M, D, deltaT);
T = K / 365.25;
return T;
}
static double 経過日数K(int Y, int M, int D, double deltaT) {
// 1月と２月は前年の１３月、１４月として考える
if (M == 1 || M == 2) {
M += 12;
Y -= 1;
}
// J2000.0からの経過日数を求める
// (2000 + Y)年M月D日として計算する
Y = Y - 2000;
return 経過日数K2(Y, M, D, deltaT);
}
static double 経過日数K2(int Y, int M, int D, double deltaT) {
double K = 0;
double KK = 0;
KK = 365 * Y + 30 * M + D - 33.875 + Math.floor(3 * (M + 1) / 5.)
+ Math.floor(Y / 4.);
K = KK + 6 / 24. + deltaT / 86400;
return K;
}
static double deltaT(int yyyy, int mm) {
// y = year + (month - 0.5)/12
double y = yyyy + (mm - 0.5) / 12.;
if (1961 <= y && y <= 1986) {
// Between years 1961 and 1986, calculate:
// ΔT = 45.45 + 1.067*t - t^2/260 - t^3 / 718
// where: t = y - 1975
double t = y - 1975;
double dT = 45.45 + 1.067 * t - Math.pow(t, 2) / 260.
- Math.pow(t, 3) / 718.;
return dT;
}
if (1986 <= y && y <= 2005) {
// Between years 1986 and 2005, calculate:
// ΔT = 63.86 + 0.3345 * t - 0.060374 * t^2 + 0.0017275 * t^3 +
// 0.000651814 * t^4
// + 0.00002373599 * t^5
// where: t = y - 2000
double t = y - 2000;
double dT = 63.86 + 0.3345 * t - 0.060374 * Math.pow(t, 2)
+ 0.0017275 * Math.pow(t, 3) + 0.000651814 * Math.pow(t, 4)
+ 0.00002373599 * Math.pow(t, 5);
return dT;
}
if (2005 <= y && y <= 2050) {
// Between years 2005 and 2050, calculate:
// ΔT = 62.92 + 0.32217 * t + 0.005589 * t^2
// where: t = y - 2000
double t = y - 2000;
double dT = 62.92 + 0.32217 * t + 0.005589 * Math.pow(t, 2);
return dT;
}
return 0;
}
static double deltaT(int yyyy) {
if (1961 <= yyyy && yyyy <= 1986) {
// Between years 1961 and 1986, calculate:
// ΔT = 45.45 + 1.067*t - t^2/260 - t^3 / 718
// where: t = y - 1975
double t = yyyy - 1975;
double dT = 45.45 + 1.067 * t - Math.pow(t, 2) / 260
- Math.pow(t, 3) / 718;
return dT;
}
if (1986 <= yyyy && yyyy <= 2005) {
// Between years 1986 and 2005, calculate:
// ΔT = 63.86 + 0.3345 * t - 0.060374 * t^2 + 0.0017275 * t^3 +
// 0.000651814 * t^4
// + 0.00002373599 * t^5
// where: t = y - 2000
double t = yyyy - 2000;
double dT = 63.86 + 0.3345 * t - 0.060374 * Math.pow(t, 2)
+ 0.0017275 * Math.pow(t, 3) + 0.000651814 * Math.pow(t, 4)
+ 0.00002373599 * Math.pow(t, 5);
return dT;
}
if (2005 <= yyyy && yyyy <= 2050) {
// Between years 2005 and 2050, calculate:
// ΔT = 62.92 + 0.32217 * t + 0.005589 * t^2
// where: t = y - 2000
double t = yyyy - 2000;
double dT = 62.92 + 0.32217 * t + 0.005589 * Math.pow(t, 2);
return dT;
}
return 0;
}
}
: