（转）根据经纬度计算日出、日落、中天、天亮、天黑和昼长时间

Jean Meeus的《天文算法》（Astronomical Algorithms，2nd Edition）第二版中第7章第60页内有详细介绍计算儒略日的方法：

设Y为给定年份，M为月份，D为该月日期（带小数，把时：分：秒折算成日的形式）。运算符INT表示为取所给数的整数部分，也即小数点前的部分。

1.若M > 2，Y和M不变。

若 M =1或2，以Y–1代Y，以M+12代M。

换句话说，如果日期在1月或2月，则被看作是在前一年的13月或14月。

2.对格里高利历(即1582年10月15日以后)，有

A = INT(Y/100),

B = 2 - A + INT(A/4).

另外，对于儒略历(即1582年10月15日之前)，取B=0。

3.所求的儒略日即为：

JD=INT(365.25(Y+4716))+INT(30.60001(M+1))+D+B-1524.5

import java.io.IOException;
import java.util.Date;
import java.util.HashMap;
import java.util.Map;
import java.util.Properties;

public class DayTime {

    private static Properties propt;

    private static double RAD = 180.0 * 3600 / Math.PI;

    private static double midDayTime;

    private static double dawnTime;

    private static String codeStr = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";

    static {
        propt = new Properties();
        try {
            propt.load(DayTime.class.getClassLoader().getResourceAsStream("jwd.properties"));
        } catch (IOException e) {
            e.printStackTrace();
        }
    }

    @SuppressWarnings("deprecation")
    public static Map<DayTimeType, String> dailytime(String area) {
        Date date = new Date();
        return dailytime(area, date.getYear(), date.getMonth(), date.getDay(), 0, 0, 0,8.0);
    }

    public static Map<DayTimeType, String> dailytime(String area, int year, int month, int day, int hour, int min, int sec,double tz) {
        Map<DayTimeType, String> map = new HashMap<DayTimeType, String>();
        String jwd = decodeJWD(area);
        System.out.println(jwd);
        Double wd = (Double.parseDouble(jwd.substring(0, 2)) + Double.parseDouble(jwd.substring(2, 4)) / 60) / 180 * Math.PI;
        Double jd = -(Double.parseDouble(jwd.substring(4, 7)) + Double.parseDouble(jwd.substring(7)) / 60) / 180 * Math.PI;
        double richu = timeToDouble(year, month, day, hour, min, sec) - 2451544.5;
        for (int i = 0; i < 10; i++)
            richu = sunRiseTime(richu, jd, wd, tz / 24);// 逐步逼近法算10次
        // 日出
        map.put(DayTimeType.SUNRISE, doubleToStr(richu));
        // 日落
        map.put(DayTimeType.SUNSET, doubleToStr(midDayTime + midDayTime - richu));
        // 中天
        map.put(DayTimeType.MIDTIME, doubleToStr(midDayTime));
        // 天亮
        map.put(DayTimeType.DAWNTIME, doubleToStr(dawnTime));
        // 天黑
        map.put(DayTimeType.NIGHTTIME, doubleToStr(midDayTime + midDayTime - dawnTime));
        // 昼长
        map.put(DayTimeType.DAYTIME, doubleToStr((midDayTime - dawnTime) * 2 - 0.5));
        return map;
    }

    // 加密
    public static String decodeJWD(String encode) {
        StringBuilder jwd = new StringBuilder();
        for (int i = 0; i < 4; i++)
            if (2 == i)
                jwd.append(String.format("%03d", codeStr.indexOf(encode.charAt(i)) + 73));
            else
                jwd.append(String.format("%02d", codeStr.indexOf(encode.charAt(i))));
        return jwd.toString();
    }

    // 解密
    public static String encodeJWD(Integer decode) {
        StringBuilder jwd = new StringBuilder();
        int i = 230811316;
        int ge = i % 100;
        int shi = i % 100000 - ge;
        int bai = i % 10000000 - shi;
        int qian = i % 1000000000 - bai;
        shi = shi / 100 - 73;
        bai = bai / 100000;
        qian = qian / 10000000;
        jwd.append(codeStr.charAt(qian)).append(codeStr.charAt(bai)).append(codeStr.charAt(shi)).append(codeStr.charAt(ge));
        return jwd.toString();
    }
    /**
     * 
     * @param date
     *            儒略日平午
     * @param lo
     *            地理经度
     * @param la
     *            地理纬度
     * @param tz
     *            时区
     * @return 太阳升起时间
     */
    public static Double sunRiseTime(double date, double lo, double la, double tz) {
        date = date - tz;
        // 太阳黄经以及它的正余弦值
        double t = date / 36525;
        double j = sunHJ(t);
        // 太阳黄经以及它的正余弦值
        double sinJ = Math.sin(j);
        double cosJ = Math.cos(j);
        // 其中2*PI*(0.7790572732640 + 1.00273781191135448*jd)恒星时（子午圈位置）
        double gst = 2 * Math.PI * (0.779057273264 + 1.00273781191135 * date) + (0.014506 + 4612.15739966 * t + 1.39667721 * t * t) / RAD;
        double E = (84381.406 - 46.836769 * t) / RAD; // 黄赤交角
        double a = Math.atan2(sinJ * Math.cos(E), cosJ);// '太阳赤经
        double D = Math.asin(Math.sin(E) * sinJ); // 太阳赤纬
        double cosH0 = (Math.sin(-50 * 60 / RAD) - Math.sin(la) * Math.sin(D)) / (Math.cos(la) * Math.cos(D)); // 日出的时角计算，地平线下50分
        double cosH1 = (Math.sin(-6 * 3600 / RAD) - Math.sin(la) * Math.sin(D)) / (Math.cos(la) * Math.cos(D)); // 天亮的时角计算，地平线下6度，若为航海请改为地平线下12度
        // 严格应当区分极昼极夜，本程序不算
        if (cosH0 >= 1 || cosH0 <= -1)
            return -0.5;// 极昼
        double H0 = -Math.acos(cosH0); // 升点时角（日出）若去掉负号 就是降点时角，也可以利用中天和升点计算
        double H1 = -Math.acos(cosH1);
        double H = gst - lo - a; // 太阳时角
        midDayTime = date - degree(H) / Math.PI / 2 + tz; // 中天时间
        dawnTime = date - degree(H - H1) / Math.PI / 2 + tz;// 天亮时间
        return date - degree(H - H0) / Math.PI / 2 + tz; // 日出时间，函数返回值
    }

    /**
     * 保证角度∈(-π,π)
     * 
     * @param ag
     * @return ag
     */
    public static Double degree(double ag) {
        ag = mod(ag, 2 * Math.PI);
        return ag <= -Math.PI ? ag + 2 * Math.PI : ag > Math.PI ? ag - 2 * Math.PI : ag;
    }

    public static Double mod(double num1, double num2) {
        num2 = Math.abs(num2);
        // 只是取决于Num1的符号
        return num1 >= 0 ? num1 - ((int) (num1 / num2)) * num2 : ((int) (Math.abs(num1) / num2)) * num2 - Math.abs(num1);
    }
    /**
     * @param t
     *            儒略世纪数
     * @return 太阳黄经
     */
    public static Double sunHJ(double t) {
        t = t + (32.0 * (t + 1.8) * (t + 1.8) - 20) / 86400.0 / 36525.0;
        // 儒略世纪年数,力学时
        double j = 48950621.66 + 6283319653.318 * t + 53 * t * t - 994 + 334166 * Math.cos(4.669257 + 628.307585 * t) + 3489 * Math.cos(4.6261 + 1256.61517 * t) + 2060.6 * Math.cos(2.67823 + 628.307585 * t) * t;
        return j / 10000000;
    }

    /**
     * 儒略日的计算
     * 
     * @param y
     *            年
     * @param M
     *            月
     * @param d
     *            日
     * @param h
     *            小时
     * @param m
     *            分
     * @param s
     *            秒
     * @return int
     */
    public static int timeToDouble(int y, int M, int d, int h, int m, int s) {
        double time = 0;
        if (M <= 2) {
            M += 12;
            y -= 1;
        }
        if (y * 372 + M * 31 + d >= 588829) {
            time = (int) (y / 100);
            time = 2 - time + (int) (time / 4);
        }
        time += (int) Math.round(365.25 * (y + 4716) + 0.01) + (int) (30.60001 * (M + 1)) + d + (h * 3600 + m * 60 + s) / (24 * 3600) - 1524.5;
        return (int) Math.round(time);
    }

    public static String doubleToStr(double time) {
        double t = time + 0.5;
        t = (t - (int) t) * 24;
        int h = (int) t;
        t = (t - h) * 60;
        int m = (int) t;
        t = (t - m) * 60;
        int s = (int) t;
        return h + ":" + m + ":" + s;
    }

    public static void main(String[] args) {
        // System.out.println(decodeJWD("N3dS肇庆"));
        // String jwd = decodeJWD("N3dS肇庆");
        // Double wd = (Double.parseDouble(jwd.substring(0, 2)) +
        // Double.parseDouble(jwd.substring(2, 4)) / 60) / 180 * Math.PI;
        // Double jd = -(Double.parseDouble(jwd.substring(4, 7)) +
        // Double.parseDouble(jwd.substring(7)) / 60) / 180 * Math.PI;
        // double richu = timeToDouble(2012, 6, 25, 0, 0, 0) - 2451544.5;
        // System.out.println(richu + " " + timeToDouble(2012, 6, 25, 0, 0, 0));
        // richu = sunRiseTime(richu, jd, wd, 8.0 / 24);
        // System.out.println(richu);
        // richu = sunRiseTime(richu, jd, wd, 8.0 / 24);// 逐步逼近法算两次
        // System.out.println(doubleToStr(richu));
        String t = "03";
        System.out.println(Integer.parseInt(t));
        System.out.println(propt.getProperty("HI"));
        String[] strs = propt.getProperty("HI").split(" ");
        for (int i = 1; i < strs.length; i++) {
            System.out.println(strs[i] + " 日出 " + dailytime(strs[i]).get(DayTimeType.SUNRISE));
            System.out.println(strs[i] + " 日落" + dailytime(strs[i]).get(DayTimeType.SUNSET));
            System.out.println(strs[i] + " 中天" + dailytime(strs[i]).get(DayTimeType.MIDTIME));
            System.out.println(strs[i] + " 天亮" + dailytime(strs[i]).get(DayTimeType.DAWNTIME));
            System.out.println(strs[i] + " 天黑" + dailytime(strs[i]).get(DayTimeType.NIGHTTIME));
            System.out.println(strs[i] + " 昼长" + dailytime(strs[i]).get(DayTimeType.DAYTIME));
        }
        System.out.println(encodeJWD(230811316));
    }
}

枚举类型：

public enum DayTimeType {
    SUNRISE, SUNSET, MIDTIME, DAWNTIME, NIGHTTIME, DAYTIME
}