日出日落时间和年均光照时长计算 java

最新推荐文章于 2023-07-18 15:10:30 发布

红尘不归尘

最新推荐文章于 2023-07-18 15:10:30 发布

阅读量3.6k

点赞数 3

文章标签： java 几何学开发语言

本文链接： https://blog.csdn.net/Empty_mountain/article/details/121911623

版权

至入行多年仍是新手的我

文章目录

前言
一、天文知识
二、计算方法
总结

前言

太阳系行星运行的规律早就被前人研究的彻彻底底，区域气象、楼宇采光等方面对太阳相对地球的运转规律应用还是比较多的。最近对有建筑物遮挡情况下区域光照时长，基于地域观测值对太阳辐射量的计算，做了一些研究，记录一下这段时间学习的过程。

不同专业间的第一道壁垒基础定义

一、天文知识

这段主要是写定义

2、太阳方位角

自观察者所在地朝正南方向的水平线与太阳光线在地平面上的投影之间的夹角，通常规定，上午的太阳方位角为正，下午的太阳方位角为负。

4、时角

地球自转-周为一天，即24h，不同的时间有不同的时角，以Ω表示之。地球转一周为360°、因而每小时的时角为15°，即以，Ω= 15t，t表示时数。

定义很多，不在赘述，较为重要的是上面三个。

太阳高度角： $h_{s}$ :
$\sin(h_{s})=\sin(\varphi)\sin(\delta)+\cos(\varphi)\cos(\delta)\cos(\Omega)$
太阳方位角 $A_{s}$ :
$\cos(A_{s})=\frac{\sin(h_{s})\sin(\varphi)-sin(\delta)}{\cos(h_s)\cos(\varphi)}$
地理纬度： $\varphi$
赤纬角： $\delta$
$\delta=23.45sin(2\pi*\frac{284+n}{365})$ (单位：°)
n表示天数，1月1日取为1，12月31日取为365。

时角： $\Omega$
角度的单位要统一下面的是按照word体系插入的公式

2、日出日落时间

当太阳高度角为0时，对应时角即为日出日落时刻的时角。
对应 $\cos(\Omega)=-\tan(\delta)\tan(\varphi)$

计算得到时角后，需要转化为角度制，除以15，得到半天的日照时长；用12+或者减去这个值就是日出日落时间。

代码如下（示例）：

public class SolarUtil {
    /**
     * 计算赤纬角 δ
     * @param day 一年里面的第几天 如1月1日 1; 12月31 365；
     * @return 角度制 赤纬角
     */
    public static double declination(int day){
        double args = 23.45;
        double season = 2* PI * (284 + day) / 365;
        return args * sin(season);
    }
    
    public static double solarTime2(double latitude,int day){
        // 判断极昼或者极夜
        double flag = latitude * SolarUtil.declination(day);
        double t = Math.acos( Math.tan(latitude * Math.PI / 180) * Math.tan(SolarUtil.declination(day) * Math.PI / 180));
        // 出现极昼或者是极夜
        if (Double.toString(t).equals("NaN")){
            if (flag < 0)
                return 0.0;
            else return 180.0;
        }
        double toDegrees = Math.toDegrees(t);
        return toDegrees;
    }
}

3、年均光照时长

目标是计算建筑物顶部每单元位置处的光照时间，不考虑建筑物的复杂结构，按照简单的长方体进行处理。
包结构

1.建筑物实体类

用左上角和右下角的横纵坐标定义建筑物。

package com.dwind.pojo;

import java.util.ArrayList;

/**
 * @author lcdwind
 * @date 2021/10/29 16:54
 */

public class Building {
    private double height;
    private double RDPointX;
    private double RDPointY;
    private double LUPointX;
    private double LUPointY;

    public Building() {
    }

    public Building(double height, double RDPointX, double RDPointY, double LUPointX, double LUPointY) {
        this.height = height;
        this.RDPointX = RDPointX;
        this.RDPointY = RDPointY;
        this.LUPointX = LUPointX;
        this.LUPointY = LUPointY;
    }

    public double getHeight() {
        return height;
    }

    public void setHeight(double height) {
        this.height = height;
    }

    public double getRDPointX() {
        return RDPointX;
    }

    public void setRDPointX(double RDPointX) {
        this.RDPointX = RDPointX;
    }

    public double getRDPointY() {
        return RDPointY;
    }

    public void setRDPointY(double RDPointY) {
        this.RDPointY = RDPointY;
    }

    public double getLUPointX() {
        return LUPointX;
    }

    public void setLUPointX(double LUPointX) {
        this.LUPointX = LUPointX;
    }

    public double getLUPointY() {
        return LUPointY;
    }

    public void setLUPointY(double LUPointY) {
        this.LUPointY = LUPointY;
    }

    @Override
    public String toString() {
        return "Building{" +
                "height=" + height +
                ", RDPointX=" + RDPointX +
                ", RDPointY=" + RDPointY +
                ", LUPointX=" + LUPointX +
                ", LUPointY=" + LUPointY +
                '}';
    }
}

3.工具类

本来觉得会经常用到，实际上没怎么用，反倒是多写了许多没有用的。

package com.dwind.util;

import static java.lang.Math.*;

/**
 * @author lcdwind
 * @date 2021/10/29 15:45
 */
public class SolarUtil {

    // 计算赤纬角 δ

    /**
     * 计算赤纬角 δ
     * @param day 一年里面的第几天 如1月1日 1; 12月31 365；
     * @return 角度制 赤纬角
     */
    public static double declination(int day){
        double args = 23.45;
        double season = 2* PI * (284 + day) / 365;
        return args * sin(season);
    }

    /**
     * 将当地时间转换为对应时刻得地球自转时角
     * @param tureSoloarTime 真太阳时 例如 12.5 表示12：30
     * @return 角度制 时角
     */
    public static double hourAngle(double tureSoloarTime){
        return  15.0 * (tureSoloarTime - 12);
    }

    /**
     * 返回对应经度下，北京时间相应的当地时间，也就是太阳时。
     * @param longitude 计算位置的经度
     * @param CST 北京时间
     * @param day 一年当中第几天
     * @return 当地太阳时
     */
    public static double trueSoloarTime(double longitude,double CST,double day){
        double B = 2* PI*(day-81) / 364;
        return  CST +(9.87 * sin(2*  B) - 7.53*cos(B) - 1.5*sin(B))/60+ (longitude - 120) / 15.0;
    }

    /**
     *
     * @param longitude 计算位置的经度
     * @param solarTime 当地太阳时
     * @param day 一年中第几天 用于修正地球自转时进动和变动产生的差异
     * @return 北京时间
     */
    public static double CSTTime(double longitude,double solarTime,double day){
        double B = 2* PI*(day-81) / 364;
        return solarTime + (120 - longitude) / 15.0 - (9.87 * sin(2*B) - 7.53*cos(B) - 1.5*sin(B))/60;
    }

    /**
     * 返回太阳高度角
     * @param latitude 纬度 弧度制
     * @param declination 赤纬角 弧度制
     * @param hourAngle 时角 弧度制
     * @return 太阳高度角 弧度制
     */
    public static double solarHight(double latitude, double declination, double hourAngle){
        return asin(sin(latitude) * sin(declination) + cos(latitude) * cos(declination) * cos(hourAngle));
    }

    /**
     * 返回太阳方位角
     * @param solarHight 太阳高度角
     * @param latidude 纬度
     * @param declination 赤纬角
     * @param time 当地太阳时，用于修正太阳方位角的正负值（与12作比较）
     * @return 太阳方位角
     */
    public static double solarAzimuth(double solarHight,double latidude, double declination,double time){

        if (time == 0){
            if(latidude>declination){
                return 0;
            } else {
                return PI;
            }
        }
        double a = (sin(solarHight) * sin(latidude) - sin(declination)) / cos(solarHight) / cos(latidude);
        if (time > 0){
            // time>12 表示下午，方位角为正。
            return acos(a);
        }else {
            // time < 12 上午，方位角为负。
            return -acos(a);
        }

    }


    /**
     * 根据太阳高度角计算太阳时角
     * @param xDiff 用于计算太阳高度角
     * @param yDiff 用于计算太阳高度角
     * @param highDiff 用于计算太阳高度角
     * @param latitude 此地的纬度
     * @param day 一年中的第几天
     * @return 角度制太阳时角
     */
    public static double solarTime(double xDiff,double yDiff,double highDiff,double latitude,int day){
        double sin_Hs = highDiff / Math.sqrt(xDiff * xDiff + yDiff * yDiff + highDiff * highDiff);


        double flag = Math.sin(latitude * Math.PI / 180)* Math.sin(SolarUtil.declination(day) * Math.PI /180);
        double t = Math.acos(sin_Hs - flag /( Math.cos(latitude * Math.PI / 180) * Math.cos(SolarUtil.declination(day) * Math.PI / 180)));
        if (Double.toString(t).equals("NaN")){
            if (flag < 0)
                return 0.0;
            else return 180.0;
        }
        double toDegrees = Math.toDegrees(t);
        return toDegrees;
    }
    public static double solarTime2(double latitude,int day){

        // 判断极昼或者极夜
        double flag = latitude * SolarUtil.declination(day);
        double t = Math.acos( Math.tan(latitude * Math.PI / 180) * Math.tan(SolarUtil.declination(day) * Math.PI / 180));
        // 出现极昼或者是极夜
        if (Double.toString(t).equals("NaN")){
            if (flag < 0)
                return 0.0;
            else return 180.0;
        }
        double toDegrees = Math.toDegrees(t);
        return toDegrees;
    }

    public static double sunDownTime(double latitude, double declination){
        return acos(-tan(latitude)*tan(declination));
    }

    /**
     * 所判断点在两点组成的直线的哪一侧，采用向量法计算，向量方向由点1指向点2
     * @param x1 点1横坐标
     * @param y1 点2纵坐标
     * @param x2 点2横坐标
     * @param y2 点2纵坐标
     * @param tempX 所判断点横坐标
     * @param tempY 所判断点纵坐标
     * @return 返回true表示线上或者右侧，false表示左侧
     */
    public static boolean lineSide(double x1,double y1, double x2,double y2, double tempX,double tempY){

//        if (x2 == x1){
//            return tempX >= x1;
//        }
        if(((tempY-y1)*(x2-x1) - (y2-y1)*(tempX-x1))>0) {
            // 左侧
            return false;
        }else {
            // 线上或者右侧
            return true;
        }

    }

    /**
     * 返回判断点是否在多边形内部
     * @param points 多边形的顶点数组
     * @param tempX 判断点横坐标
     * @param tempY 判断点纵坐标
     * @return 在多边形内或者上面 true，外面false
     */
    public static boolean isCovered(double[][] points,double tempX,double tempY){
        boolean flag = true;
        for (int i =0;i<points.length-1;i ++){
            boolean b = lineSide(points[i][0], points[i][1], points[i + 1][0], points[i + 1][1], tempX, tempY);
            if (!b){
                flag = false;
                break;
            }
        }
        return flag;
    }
}

4.计算部分

package com.dwind.building;

import com.dwind.pojo.Building;
import com.dwind.util.SolarUtil;

import java.util.ArrayList;
import java.util.HashMap;

/**
 * @author lcdwind
 * @date 2021/10/29 16:35
 */
public class Sunshine {
    public HashMap<Building, Double> sunTime(Building target, ArrayList<Building> buildings, double latitude, double longitude) {
        HashMap<String, Double> map = new HashMap<>();
        // 用于记录每个小单元的年均光照时间
        HashMap<Building, Double> cubeMap = new HashMap<>();
        // 切分的单位长度
        double cube = 1.0;
        // target切分成一个个小的单元 单位cube；
        ArrayList<Building> buildingCube = new ArrayList<>();
        for (double xMark = target.getLUPointX(), yMark = target.getRDPointY(); yMark < target.getLUPointY(); yMark += cube) {
            for (; xMark < target.getRDPointX(); xMark += cube) {
                Building building = new Building(target.getHeight(), xMark + cube, yMark, xMark, yMark + cube);
                buildingCube.add(building);
            }
            xMark = target.getLUPointX();
        }


        // 循环计算每一小块儿的年均光照
        for (Building building : buildingCube) {
            map.put("sunHour", 0.0);
            // 循环计算每天的光照情况
            for (int day = 1; day <= 365; day++) {

//遮挡物在建筑物的东边还是西边有区别没有？
                // 这里只判断高度角得出时角是不够精确的，并且在目标建筑物左右两侧的遮挡物并没有做区分处理，两者会被最高最近的结果所覆盖。
                // 有一种思路是根据建筑物的光影图棒长图计算出当日的日照影子时间变化，然后按照3种情况做分析，最后取中心值的点做为标记点也是不错的方法

                // sunUpTime 日出时刻的时角
                double sunUpTime = SolarUtil.solarTime(1, 1, 0, latitude, day);
                // sunUP,sunDown 日出日落时刻
                double sunUp = 12 - sunUpTime / 15;
                double sunDown = 12 + sunUpTime / 15;
                // 一天中没有被遮挡的时间。
                double noCoverTime = 0;
                // 每间隔 0.01*60 分钟采样一次，判断目标位置是否被覆盖
                for (double time = sunUp; time <= sunDown; time += 0.01) {
                    boolean total = true;
                    // 逐个建筑物判断遮挡情况
                    for (Building value : buildings) {
                        // 高度差，太阳高度角，太阳方位角，影子多边形数组
                        double highDiff = value.getHeight() - building.getHeight();
                        double sunHeight = SolarUtil.solarHight(latitude * Math.PI / 180, SolarUtil.declination(day) * Math.PI / 180, Math.abs(time - 12) * 15 * Math.PI / 180);
                        double sunAzimuth = SolarUtil.solarAzimuth(sunHeight, latitude * Math.PI / 180, SolarUtil.declination(day) * Math.PI / 180, time - 12);
                        double[][] polygon = null;
                        double shadowX, shadowY;
                        double shadowOfY = Math.abs(highDiff / Math.tan(sunHeight) * Math.cos(sunAzimuth));
                        double shadowOfX = highDiff / Math.tan(sunHeight) * Math.sin(sunAzimuth);
                        if (time > 12) {
                            // 上午的情况
                            shadowX = value.getLUPointX() - shadowOfX;
                            shadowY = value.getRDPointY() - shadowOfY;
                            polygon = new double[][]{
                                    {value.getRDPointX(), value.getRDPointY()},
                                    {shadowX + (value.getRDPointX() - value.getLUPointX()), shadowY},
                                    {shadowX, shadowY},
                                    {shadowX, shadowY + (value.getLUPointY() - value.getRDPointY())},
                                    {value.getLUPointX(), value.getRDPointY()},
                                    {value.getRDPointX(), value.getRDPointY()}
                            };
                        } else {
                            // 下午的情况
                            shadowX = value.getRDPointX() - shadowOfX;
                            shadowY = value.getRDPointY() - shadowOfY;
                            polygon = new double[][]{
                                    {value.getRDPointX(), value.getLUPointY()},
                                    {shadowX, shadowY + (value.getLUPointY() - value.getRDPointY())},
                                    {shadowX, shadowY},
                                    {shadowX - (value.getRDPointX() - value.getLUPointX()), shadowY},
                                    {value.getLUPointX(), value.getRDPointY()},
                                    {value.getRDPointX(), value.getLUPointY()}
                            };

                        }
                        // 有没有被覆盖
                        boolean covered = SolarUtil.isCovered(polygon, (building.getRDPointX() + building.getLUPointX()) / 2, (building.getLUPointY() + building.getRDPointY()) / 2);

                        if (time <= sunUp || time >= sunDown) {
                            covered = true;
                        }
                        if (covered) {
                            // 任意建筑物覆盖即为覆盖
                            total = false;
                            break;
                        }
                    }
                    if (total) {
                        noCoverTime += 0.01;
                    }
                }
                map.put("sunHour", map.get("sunHour") + noCoverTime);

            }
            cubeMap.put(building, map.get("sunHour"));

        }
        return cubeMap;
    }

}

该处主要是先将周围建筑物的阴影区计算出来，边界点存入数组中，构成一个多边形，判断目标位置是否在多边形区域内，是就意味着在影子范围内，被遮挡。
循环判断每个周边建筑物遮挡情况，循环计算一年中每天的光照时长，最后循环计算每个目标位置的年光照时长。循环非常多，计算非常耗费算力，我的测试用例大概需要23s的运行时长。

6、测试

package com.dwind.building;

import com.dwind.pojo.Building;
import com.dwind.util.SolarUtil;
import org.junit.jupiter.api.Test;

import java.util.ArrayList;
import java.util.HashMap;

import static org.junit.jupiter.api.Assertions.*;

/**
 * @author lcdwind
 * @date 2021/11/1 9:31
 */
class SunshineTest {

    @Test
    void sunTime() {
        Building target = new Building(10, 5, 0, 0, 5);

        Building building = new Building(20, 2, 5, -5, 10);
        Building building1 = new Building(20, 5, 6, 3, 10);

        ArrayList<Building> buildings = new ArrayList<>();
        buildings.add(building);
        buildings.add(building1);

        double latitude = 34.5;
        double longitude = 113.5;
//        70.9 6108.000000000003 79.9 7524.0 85.9 8267.999999999996

        long start = System.currentTimeMillis();
        Sunshine sunshine = new Sunshine();
        HashMap<Building, Double> map = sunshine.sunTime(target, buildings, latitude,longitude);
        System.out.println(map);
        long end = System.currentTimeMillis();
        System.out.println(end - start);
    }
}

测试结果：

{Building{height=10.0, RDPointX=3.0, RDPointY=2.0, LUPointX=2.0, LUPointY=3.0}=2723.4999999999563, Building{height=10.0, RDPointX=4.0, RDPointY=4.0, LUPointX=3.0, LUPointY=5.0}=1665.539999999981, Building{height=10.0, RDPointX=4.0, RDPointY=3.0, LUPointX=3.0, LUPointY=4.0}=2292.6099999999683, Building{height=10.0, RDPointX=5.0, RDPointY=2.0, LUPointX=4.0, LUPointY=3.0}=2907.0799999999567, Building{height=10.0, RDPointX=5.0, RDPointY=3.0, LUPointX=4.0, LUPointY=4.0}=2533.6099999999633, Building{height=10.0, RDPointX=4.0, RDPointY=1.0, LUPointX=3.0, LUPointY=2.0}=3120.6399999999494, Building{height=10.0, RDPointX=5.0, RDPointY=0.0, LUPointX=4.0, LUPointY=1.0}=3385.0299999999424, Building{height=10.0, RDPointX=5.0, RDPointY=4.0, LUPointX=4.0, LUPointY=5.0}=2101.9799999999723, Building{height=10.0, RDPointX=2.0, RDPointY=2.0, LUPointX=1.0, LUPointY=3.0}=2588.4099999999594, Building{height=10.0, RDPointX=2.0, RDPointY=4.0, LUPointX=1.0, LUPointY=5.0}=1172.2899999999863, Building{height=10.0, RDPointX=1.0, RDPointY=1.0, LUPointX=0.0, LUPointY=2.0}=2923.4099999999553, Building{height=10.0, RDPointX=1.0, RDPointY=3.0, LUPointX=0.0, LUPointY=4.0}=1865.1199999999715, Building{height=10.0, RDPointX=2.0, RDPointY=3.0, LUPointX=1.0, LUPointY=4.0}=1956.239999999974, Building{height=10.0, RDPointX=2.0, RDPointY=0.0, LUPointX=1.0, LUPointY=1.0}=3285.299999999948, Building{height=10.0, RDPointX=1.0, RDPointY=2.0, LUPointX=0.0, LUPointY=3.0}=2502.2399999999625, Building{height=10.0, RDPointX=3.0, RDPointY=1.0, LUPointX=2.0, LUPointY=2.0}=3080.9399999999505, Building{height=10.0, RDPointX=5.0, RDPointY=1.0, LUPointX=4.0, LUPointY=2.0}=3174.8799999999483, Building{height=10.0, RDPointX=4.0, RDPointY=2.0, LUPointX=3.0, LUPointY=3.0}=2791.96999999996, Building{height=10.0, RDPointX=3.0, RDPointY=4.0, LUPointX=2.0, LUPointY=5.0}=1539.9799999999827, Building{height=10.0, RDPointX=1.0, RDPointY=0.0, LUPointX=0.0, LUPointY=1.0}=3225.899999999951, Building{height=10.0, RDPointX=3.0, RDPointY=3.0, LUPointX=2.0, LUPointY=4.0}=2182.5799999999695, Building{height=10.0, RDPointX=2.0, RDPointY=1.0, LUPointX=1.0, LUPointY=2.0}=2997.8999999999523, Building{height=10.0, RDPointX=1.0, RDPointY=4.0, LUPointX=0.0, LUPointY=5.0}=1096.6699999999862, Building{height=10.0, RDPointX=3.0, RDPointY=0.0, LUPointX=2.0, LUPointY=1.0}=3339.789999999945, Building{height=10.0, RDPointX=4.0, RDPointY=0.0, LUPointX=3.0, LUPointY=1.0}=3362.059999999947}
20191