日の出計算 その２
2012/10/24
java
天文
Javaです。ここんとこずーと日の出の計算をおこなっていました。
計算結果を国立天文台の値とくらべてみます。
結果は以下です。１分違いの不一致が188個でした。どうでしょうか。使えない値でしょうか。左の値が、国立天文台の時間です。
年月日 出 時分 比較 2012/1/1 6:50 6:50 ○ 2012/1/2 6:51 6:50 × 2012/1/3 6:51 6:50 × 2012/1/4 6:51 6:50 × 2012/1/5 6:51 6:51 ○ 2012/1/6 6:51 6:51 ○ 2012/1/7 6:51 6:51 ○ 2012/1/8 6:51 6:51 ○ 2012/1/9 6:51 6:51 ○ 2012/1/10 6:51 6:51 ○ 2012/1/11 6:51 6:50 × 2012/1/12 6:51 6:50 × 2012/1/13 6:51 6:50 × 2012/1/14 6:50 6:50 ○ 2012/1/15 6:50 6:50 ○ 2012/1/16 6:50 6:50 ○ 2012/1/17 6:50 6:49 × 2012/1/18 6:49 6:49 ○ 2012/1/19 6:49 6:49 ○ 2012/1/20 6:49 6:48 × 2012/1/21 6:48 6:48 ○ 2012/1/22 6:48 6:47 × 2012/1/23 6:47 6:47 ○ 2012/1/24 6:47 6:46 × 2012/1/25 6:46 6:46 ○ 2012/1/26 6:46 6:45 × 2012/1/27 6:45 6:45 ○ 2012/1/28 6:45 6:44 × 2012/1/29 6:44 6:43 × 2012/1/30 6:43 6:43 ○ 2012/1/31 6:42 6:42 ○ 2012/2/1 6:42 6:41 × 2012/2/2 6:41 6:41 ○ 2012/2/3 6:40 6:40 ○ 2012/2/4 6:39 6:39 ○ 2012/2/5 6:39 6:38 × 2012/2/6 6:38 6:37 × 2012/2/7 6:37 6:36 × 2012/2/8 6:36 6:35 × 2012/2/9 6:35 6:35 ○ 2012/2/10 6:34 6:34 ○ 2012/2/11 6:33 6:33 ○ 2012/2/12 6:32 6:32 ○ 2012/2/13 6:31 6:31 ○ 2012/2/14 6:30 6:30 ○ 2012/2/15 6:29 6:28 × 2012/2/16 6:28 6:27 × 2012/2/17 6:27 6:26 × 2012/2/18 6:26 6:25 × 2012/2/19 6:24 6:24 ○ 2012/2/20 6:23 6:23 ○ 2012/2/21 6:22 6:22 ○ 2012/2/22 6:21 6:21 ○ 2012/2/23 6:20 6:19 × 2012/2/24 6:19 6:18 × 2012/2/25 6:17 6:17 ○ 2012/2/26 6:16 6:16 ○ 2012/2/27 6:15 6:14 × 2012/2/28 6:14 6:13 × 2012/2/29 6:12 6:12 ○ 2012/3/1 6:11 6:11 ○ 2012/3/2 6:10 6:09 × 2012/3/3 6:08 6:08 ○ 2012/3/4 6:07 6:07 ○ 2012/3/5 6:06 6:05 × 2012/3/6 6:04 6:04 ○ 2012/3/7 6:03 6:03 ○ 2012/3/8 6:02 6:01 × 2012/3/9 6:00 6:00 ○ 2012/3/10 5:59 5:58 × 2012/3/11 5:57 5:57 ○ 2012/3/12 5:56 5:56 ○ 2012/3/13 5:55 5:54 × 2012/3/14 5:53 5:53 ○ 2012/3/15 5:52 5:51 × 2012/3/16 5:50 5:50 ○ 2012/3/17 5:49 5:49 ○ 2012/3/18 5:48 5:47 × 2012/3/19 5:46 5:46 ○ 2012/3/20 5:45 5:44 × 2012/3/21 5:43 5:43 ○ 2012/3/22 5:42 5:41 × 2012/3/23 5:40 5:40 ○ 2012/3/24 5:39 5:39 ○ 2012/3/25 5:38 5:37 × 2012/3/26 5:36 5:36 ○ 2012/3/27 5:35 5:34 × 2012/3/28 5:33 5:33 ○ 2012/3/29 5:32 5:31 × 2012/3/30 5:30 5:30 ○ 2012/3/31 5:29 5:29 ○ 2012/4/1 5:28 5:27 × 2012/4/2 5:26 5:26 ○ 2012/4/3 5:25 5:24 × 2012/4/4 5:23 5:23 ○ 2012/4/5 5:22 5:21 × 2012/4/6 5:21 5:20 × 2012/4/7 5:19 5:19 ○ 2012/4/8 5:18 5:17 × 2012/4/9 5:16 5:16 ○ 2012/4/10 5:15 5:15 ○ 2012/4/11 5:14 5:13 × 2012/4/12 5:12 5:12 ○ 2012/4/13 5:11 5:11 ○ 2012/4/14 5:10 5:09 × 2012/4/15 5:08 5:08 ○ 2012/4/16 5:07 5:07 ○ 2012/4/17 5:06 5:05 × 2012/4/18 5:05 5:04 × 2012/4/19 5:03 5:03 ○ 2012/4/20 5:02 5:01 × 2012/4/21 5:01 5:00 × 2012/4/22 5:00 4:59 × 2012/4/23 4:58 4:58 ○ 2012/4/24 4:57 4:57 ○ 2012/4/25 4:56 4:55 × 2012/4/26 4:55 4:54 × 2012/4/27 4:54 4:53 × 2012/4/28 4:52 4:52 ○ 2012/4/29 4:51 4:51 ○ 2012/4/30 4:50 4:50 ○ 2012/5/1 4:49 4:49 ○ 2012/5/2 4:48 4:48 ○ 2012/5/3 4:47 4:46 × 2012/5/4 4:46 4:45 × 2012/5/5 4:45 4:44 × 2012/5/6 4:44 4:43 × 2012/5/7 4:43 4:42 × 2012/5/8 4:42 4:42 ○ 2012/5/9 4:41 4:41 ○ 2012/5/10 4:40 4:40 ○ 2012/5/11 4:39 4:39 ○ 2012/5/12 4:38 4:38 ○ 2012/5/13 4:38 4:37 × 2012/5/14 4:37 4:36 × 2012/5/15 4:36 4:35 × 2012/5/16 4:35 4:35 ○ 2012/5/17 4:35 4:34 × 2012/5/18 4:34 4:33 × 2012/5/19 4:33 4:33 ○ 2012/5/20 4:32 4:32 ○ 2012/5/21 4:32 4:31 × 2012/5/22 4:31 4:31 ○ 2012/5/23 4:31 4:30 × 2012/5/24 4:30 4:30 ○ 2012/5/25 4:30 4:29 × 2012/5/26 4:29 4:28 × 2012/5/27 4:29 4:28 × 2012/5/28 4:28 4:28 ○ 2012/5/29 4:28 4:27 × 2012/5/30 4:27 4:27 ○ 2012/5/31 4:27 4:26 × 2012/6/1 4:27 4:26 × 2012/6/2 4:26 4:26 ○ 2012/6/3 4:26 4:25 × 2012/6/4 4:26 4:25 × 2012/6/5 4:25 4:25 ○ 2012/6/6 4:25 4:25 ○ 2012/6/7 4:25 4:25 ○ 2012/6/8 4:25 4:24 × 2012/6/9 4:25 4:24 × 2012/6/10 4:25 4:24 × 2012/6/11 4:25 4:24 × 2012/6/12 4:25 4:24 × 2012/6/13 4:25 4:24 × 2012/6/14 4:25 4:24 × 2012/6/15 4:25 4:24 × 2012/6/16 4:25 4:24 × 2012/6/17 4:25 4:24 × 2012/6/18 4:25 4:24 × 2012/6/19 4:25 4:25 ○ 2012/6/20 4:25 4:25 ○ 2012/6/21 4:25 4:25 ○ 2012/6/22 4:26 4:25 × 2012/6/23 4:26 4:25 × 2012/6/24 4:26 4:26 ○ 2012/6/25 4:27 4:26 × 2012/6/26 4:27 4:26 × 2012/6/27 4:27 4:27 ○ 2012/6/28 4:28 4:27 × 2012/6/29 4:28 4:27 × 2012/6/30 4:28 4:28 ○ 2012/7/1 4:29 4:28 × 2012/7/2 4:29 4:29 ○ 2012/7/3 4:30 4:29 × 2012/7/4 4:30 4:30 ○ 2012/7/5 4:31 4:30 × 2012/7/6 4:31 4:31 ○ 2012/7/7 4:32 4:31 × 2012/7/8 4:32 4:32 ○ 2012/7/9 4:33 4:32 × 2012/7/10 4:33 4:33 ○ 2012/7/11 4:34 4:34 ○ 2012/7/12 4:35 4:34 × 2012/7/13 4:35 4:35 ○ 2012/7/14 4:36 4:35 × 2012/7/15 4:37 4:36 × 2012/7/16 4:37 4:37 ○ 2012/7/17 4:38 4:37 × 2012/7/18 4:39 4:38 × 2012/7/19 4:39 4:39 ○ 2012/7/20 4:40 4:40 ○ 2012/7/21 4:41 4:40 × 2012/7/22 4:41 4:41 ○ 2012/7/23 4:42 4:42 ○ 2012/7/24 4:43 4:42 × 2012/7/25 4:44 4:43 × 2012/7/26 4:44 4:44 ○ 2012/7/27 4:45 4:45 ○ 2012/7/28 4:46 4:45 × 2012/7/29 4:47 4:46 × 2012/7/30 4:47 4:47 ○ 2012/7/31 4:48 4:48 ○ 2012/8/1 4:49 4:48 × 2012/8/2 4:50 4:49 × 2012/8/3 4:50 4:50 ○ 2012/8/4 4:51 4:51 ○ 2012/8/5 4:52 4:52 ○ 2012/8/6 4:53 4:52 × 2012/8/7 4:54 4:53 × 2012/8/8 4:54 4:54 ○ 2012/8/9 4:55 4:55 ○ 2012/8/10 4:56 4:55 × 2012/8/11 4:57 4:56 × 2012/8/12 4:57 4:57 ○ 2012/8/13 4:58 4:58 ○ 2012/8/14 4:59 4:59 ○ 2012/8/15 5:00 4:59 × 2012/8/16 5:01 5:00 × 2012/8/17 5:01 5:01 ○ 2012/8/18 5:02 5:02 ○ 2012/8/19 5:03 5:02 × 2012/8/20 5:04 5:03 × 2012/8/21 5:05 5:04 × 2012/8/22 5:05 5:05 ○ 2012/8/23 5:06 5:06 ○ 2012/8/24 5:07 5:06 × 2012/8/25 5:08 5:07 × 2012/8/26 5:08 5:08 ○ 2012/8/27 5:09 5:09 ○ 2012/8/28 5:10 5:09 × 2012/8/29 5:11 5:10 × 2012/8/30 5:11 5:11 ○ 2012/8/31 5:12 5:12 ○ 2012/9/1 5:13 5:12 × 2012/9/2 5:14 5:13 × 2012/9/3 5:14 5:14 ○ 2012/9/4 5:15 5:15 ○ 2012/9/5 5:16 5:15 × 2012/9/6 5:17 5:16 × 2012/9/7 5:17 5:17 ○ 2012/9/8 5:18 5:18 ○ 2012/9/9 5:19 5:18 × 2012/9/10 5:20 5:19 × 2012/9/11 5:20 5:20 ○ 2012/9/12 5:21 5:21 ○ 2012/9/13 5:22 5:21 × 2012/9/14 5:23 5:22 × 2012/9/15 5:23 5:23 ○ 2012/9/16 5:24 5:24 ○ 2012/9/17 5:25 5:25 ○ 2012/9/18 5:26 5:25 × 2012/9/19 5:27 5:26 × 2012/9/20 5:27 5:27 ○ 2012/9/21 5:28 5:28 ○ 2012/9/22 5:29 5:28 × 2012/9/23 5:30 5:29 × 2012/9/24 5:30 5:30 ○ 2012/9/25 5:31 5:31 ○ 2012/9/26 5:32 5:31 × 2012/9/27 5:33 5:32 × 2012/9/28 5:33 5:33 ○ 2012/9/29 5:34 5:34 ○ 2012/9/30 5:35 5:35 ○ 2012/10/1 5:36 5:35 × 2012/10/2 5:37 5:36 × 2012/10/3 5:37 5:37 ○ 2012/10/4 5:38 5:38 ○ 2012/10/5 5:39 5:38 × 2012/10/6 5:40 5:39 × 2012/10/7 5:41 5:40 × 2012/10/8 5:41 5:41 ○ 2012/10/9 5:42 5:42 ○ 2012/10/10 5:43 5:43 ○ 2012/10/11 5:44 5:43 × 2012/10/12 5:45 5:44 × 2012/10/13 5:46 5:45 × 2012/10/14 5:46 5:46 ○ 2012/10/15 5:47 5:47 ○ 2012/10/16 5:48 5:48 ○ 2012/10/17 5:49 5:49 ○ 2012/10/18 5:50 5:49 × 2012/10/19 5:51 5:50 × 2012/10/20 5:52 5:51 × 2012/10/21 5:53 5:52 × 2012/10/22 5:54 5:53 × 2012/10/23 5:54 5:54 ○ 2012/10/24 5:55 5:55 ○ 2012/10/25 5:56 5:56 ○ 2012/10/26 5:57 5:57 ○ 2012/10/27 5:58 5:58 ○ 2012/10/28 5:59 5:59 ○ 2012/10/29 6:00 5:59 × 2012/10/30 6:01 6:00 × 2012/10/31 6:02 6:01 × 2012/11/1 6:03 6:02 × 2012/11/2 6:04 6:03 × 2012/11/3 6:05 6:04 × 2012/11/4 6:06 6:05 × 2012/11/5 6:07 6:06 × 2012/11/6 6:08 6:07 × 2012/11/7 6:09 6:08 × 2012/11/8 6:10 6:09 × 2012/11/9 6:11 6:10 × 2012/11/10 6:12 6:11 × 2012/11/11 6:13 6:12 × 2012/11/12 6:14 6:13 × 2012/11/13 6:15 6:14 × 2012/11/14 6:16 6:15 × 2012/11/15 6:17 6:16 × 2012/11/16 6:18 6:17 × 2012/11/17 6:19 6:18 × 2012/11/18 6:20 6:19 × 2012/11/19 6:21 6:20 × 2012/11/20 6:22 6:21 × 2012/11/21 6:23 6:22 × 2012/11/22 6:24 6:23 × 2012/11/23 6:25 6:24 × 2012/11/24 6:25 6:25 ○ 2012/11/25 6:26 6:26 ○ 2012/11/26 6:27 6:27 ○ 2012/11/27 6:28 6:28 ○ 2012/11/28 6:29 6:29 ○ 2012/11/29 6:30 6:30 ○ 2012/11/30 6:31 6:31 ○ 2012/12/1 6:32 6:32 ○ 2012/12/2 6:33 6:32 × 2012/12/3 6:34 6:33 × 2012/12/4 6:35 6:34 × 2012/12/5 6:36 6:35 × 2012/12/6 6:36 6:36 ○ 2012/12/7 6:37 6:37 ○ 2012/12/8 6:38 6:38 ○ 2012/12/9 6:39 6:38 × 2012/12/10 6:40 6:39 × 2012/12/11 6:40 6:40 ○ 2012/12/12 6:41 6:41 ○ 2012/12/13 6:42 6:41 × 2012/12/14 6:42 6:42 ○ 2012/12/15 6:43 6:43 ○ 2012/12/16 6:44 6:43 × 2012/12/17 6:44 6:44 ○ 2012/12/18 6:45 6:45 ○ 2012/12/19 6:46 6:45 × 2012/12/20 6:46 6:46 ○ 2012/12/21 6:47 6:46 × 2012/12/22 6:47 6:47 ○ 2012/12/23 6:48 6:47 × 2012/12/24 6:48 6:48 ○ 2012/12/25 6:49 6:48 × 2012/12/26 6:49 6:48 × 2012/12/27 6:49 6:49 ○ 2012/12/28 6:50 6:49 × 2012/12/29 6:50 6:49 × 2012/12/30 6:50 6:50 ○ 2012/12/31 6:50 6:50 ○ 不一致:188
ソースコード
import static java.lang.Math.*;
import java.text.SimpleDateFormat;
import java.util.Calendar;
import java.util.Date;
/**
*
* @author nakawakashigeto
*
*/
public class A20121021 {
public static void main(String[] args) {
Calendar calendar = Calendar.getInstance();
calendar.set(2012, 1 - 1, 1, 0, 0, 0);
// 検証のために1年分の出力
while(true){
int yyyy = calendar.get(Calendar.YEAR);
int m = calendar.get(Calendar.MONTH) + 1;
int d = calendar.get(Calendar.DAY_OF_MONTH);
sunrise("東京", yyyy, m, d, 35.6581 , 139.7414);
// 一日進める
calendar.add(Calendar.DAY_OF_YEAR, 1);
if(2013 <= calendar.get(Calendar.YEAR)) break;
}
}
/**
* 日付表現クラス
*
*/
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 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("] ");
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 {
/** ra:赤経 α right ascension */
Degree ra;
/** δ:赤緯 デルタ declination */
Degree dec;
/** r:距離 */
double r;
/** th:恒星時 θ シータ sidereal time */
Degree th;
public String toString() {
StringBuilder builder = new StringBuilder();
builder.append("a:[").append(ra).append("] ").append("δ:[")
.append(dec).append("] ").append("r:[").append(r)
.append("] ").append("th:[").append(th).append("] ");
return builder.toString();
}
}
/**
* 日の出計算を行う
*
* @param location
* @param Y
* @param M
* @param D
* @param longitude
* @param latitude
*/
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 = calcSunrise(ymdh, longitudeAndLatitude);
// 補正なぜこれが必要になるのか？
// 書籍には少数点の値のみでよいと書いてある
if(d >= 1)
d -= floor(d);
// System.out.println(d);
d = d * 24.;
double d2 = d % 1;
String s = location + "\t" + ymdh.Y + "年" + ymdh.M + "月"
+ ymdh.D + "日" + (int) floor(d) + "時" + (int) floor(d2 * 60)
+ "分 \t";
String s2 = String.format("%d:%02d", (int) floor(d), (int) floor(d2 * 60));
System.out.println(s + s2);
}
/**
* 仮定時刻に対する日の出計算を行い、補正値を返す。
*
* @param ymdh
* @param longitudeAndLatitude
* @return
*/
static double calcSunrise(YMDH ymdh, LongitudeAndLatitude longitudeAndLatitude) {
// (1) 日の出時刻を仮定し、日単位の時刻変数dに換算する
// d:時刻変数
double d = ymdh.H / 24.;
// (2) 太陽の赤経、赤緯、恒星時を求める
Sun sun = a2(ymdh.Y, ymdh.M, ymdh.D, ymdh.H, d,
longitudeAndLatitude.longitude.deg);
// (3) 出没高度k
double k = k(sun.r);
// (4) 時角tk
double tk = tk(sun.dec.deg, longitudeAndLatitude.latitude.deg, k);
tk = toDegrees(acos(tk));
// (5) 恒星時から太陽の時角tを計算する
double t = t(sun.ra.deg, sun.th.deg);
//
if (0 < tk)
tk *= -1;
// (6) 仮定時刻に対する補正値Δdの計算をする
double dd = deltaD(tk, t);
// System.out.print("ra:" + (int)sun.ra.deg +"\tt:"+ t + "\ttk:" + tk + "\td:" + d + " dd:" + dd + " ");
// 補正値を加えて新しく仮定時刻を決定する
double newd = d + dd;// 補正値
return newd;
}
/**
*
* @param Y
* @param M
* @param D
* @param H
* @param λ2
* @return
*/
static Sun a2(int Y, int M, int D, int H, double d, double λ2) {
// T 時刻変数
double T = 時刻変数T(Y, M, D);
double λs = λs(T);
double q = q(T);
double r = r(q);
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) {
tana += toRadians(360);
} else if (90 <= _λs && _λs < 180) {
tana += toRadians(180);
} else if (180 <= _λs && _λs < 270) {
tana += toRadians(180);
} else if (270 <= _λs && _λs < 360) {
tana += toRadians(360);
}
// System.out.print((int)toDegrees(tana)+"\t");
sinδ = asin(sinδ);
// 恒星時を求める
double th = siderealTimeθ(T, d, λ2);
Sun sun = new Sun();
{
sun.ra = new Degree(toDegrees(tana));
sun.dec = new Degree(toDegrees(sinδ));
sun.r = r;
sun.th = new Degree(th);
}
// System.out.print(sun+"\t");
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;
}
/**
*
* @param delta
* @param phi
* @param k
* @return
*/
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;
}
/**
* 距離 r
* @param q
* @return
*/
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;
}
}
: