通过圭表计算影子长度的模型——python实现

学校的地理老师发了一个表，需要统计一年的影子长度，作业要等一年的时间再交.......不知道一年后都忘到哪里去了，索性现在就通过编程来搞定它，一起来！

本程序可实现的核心功能：
日期输入：接受MM-DD格式的日期输入
地点固定：北京（纬度39.9042°，经度116.4074°）
圭表高度：固定 10cm
影长计算：基于太阳高度角的几何计算

import math
from datetime import datetime, timedelta
import numpy as np

class GnomonShadowCalculator:
    def __init__(self):
        # 北京的地理坐标（天安门广场）
        self.latitude = 39.9042  # 纬度，单位：度
        self.longitude = 116.4074  # 经度，单位：度
        self.timezone = 8  # 时区：东八区

        # 圭表参数
        self.gnomon_height = 10.0  # 圭表高度，单位：厘米

        # 天文常数
        self.obliquity = 23.4397  # 黄赤交角，单位：度

    def date_to_julian(self, month, day, year=2024):
        """将日期转换为儒略日"""
        # 简化计算，假设为2024年（闰年）
        month_days = [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

        # 计算年内的第几天
        day_of_year = sum(month_days[:month-1]) + day

        # 计算儒略日（简化版，从春分3月20日开始计数）
        # 实际计算中，春分通常在3月20日或21日
        spring_equinox = 80  # 3月20日是年内的第80天（2024年）
        julian_day = day_of_year - spring_equinox

        return julian_day

    def calculate_solar_declination(self, julian_day):
        """计算太阳赤纬"""
        # 使用简化公式：δ = ε * sin(360/365.24 * n)
        # 其中ε是黄赤交角，n是春分后的天数
        n = julian_day
        # 将角度转换为弧度
        angle_rad = math.radians(360 * n / 365.24)

        # 计算赤纬（单位：度）
        declination = self.obliquity * math.sin(angle_rad)

        return declination

    def calculate_equation_of_time(self, julian_day):
        """计算时差（太阳时与平太阳时的差值）"""
        # 简化计算时差
        n = julian_day
        b = math.radians(360 * (n - 81) / 365)

        # 简化时差公式（单位：分钟）
        eot = 9.87 * math.sin(2*b) - 7.53 * math.cos(b) - 1.5 * math.sin(b)

        return eot / 60.0  # 转换为小时

    def calculate_solar_noon(self, date_str, eot):
        """计算正午时间（太阳时）"""
        # 解析日期
        month, day = map(int, date_str.split('-'))

        # 平太阳时正午（地方时12:00）
        solar_noon_local = 12.0

        # 经度修正（北京经度116.4°，东八区中心经度120°）
        longitude_correction = (self.longitude - 120) * 4 / 60  # 每度4分钟

        # 计算真太阳时正午
        true_solar_noon = solar_noon_local - longitude_correction - eot

        # 确保时间在0-24之间
        if true_solar_noon < 0:
            true_solar_noon += 24
        elif true_solar_noon >= 24:
            true_solar_noon -= 24

        return true_solar_noon

    def calculate_solar_altitude(self, date_str, hour=12):
        """计算太阳高度角"""
        # 解析日期
        month, day = map(int, date_str.split('-'))
        julian_day = self.date_to_julian(month, day)

        # 计算太阳赤纬
        declination = self.calculate_solar_declination(julian_day)

        # 计算时角（基于正午）
        # 计算时差
        eot = self.calculate_equation_of_time(julian_day)

        # 计算真太阳时
        true_solar_time = hour - eot

        # 时角计算（每小时15度，正午为0）
        hour_angle = 15 * (true_solar_time - 12)

        # 将角度转换为弧度
        lat_rad = math.radians(self.latitude)
        dec_rad = math.radians(declination)
        ha_rad = math.radians(hour_angle)

        # 计算太阳高度角公式：
        # sin(alt) = sin(lat)*sin(dec) + cos(lat)*cos(dec)*cos(ha)
        sin_alt = (math.sin(lat_rad) * math.sin(dec_rad) + 
                   math.cos(lat_rad) * math.cos(dec_rad) * math.cos(ha_rad))

        # 确保sin_alt在有效范围内
        sin_alt = max(min(sin_alt, 1.0), -1.0)

        # 计算高度角（单位：度）
        altitude = math.degrees(math.asin(sin_alt))

        return altitude

    def calculate_shadow_length(self, date_str, hour=12):
        """计算影长"""
        # 计算太阳高度角
        altitude = self.calculate_solar_altitude(date_str, hour)

        # 如果太阳在地平线以下，返回None
        if altitude <= 0:
            return None, altitude

        # 将高度角转换为弧度
        alt_rad = math.radians(altitude)

        # 计算影长：影长 = 圭表高度 / tan(太阳高度角)
        shadow_length = self.gnomon_height / math.tan(alt_rad)

        return shadow_length, altitude

    def calculate_daily_shadow_curve(self, date_str):
        """计算一天的影长变化曲线"""
        results = []

        for hour in range(6, 19):  # 从6点到18点
            for minute in [0, 30]:  # 每小时两个点
                time_decimal = hour + minute/60.0
                shadow_length, altitude = self.calculate_shadow_length(date_str, time_decimal)

                if shadow_length is not None:
                    time_str = f"{hour:02d}:{minute:02d}"
                    results.append({
                        'time': time_str,
                        'solar_altitude': round(altitude, 2),
                        'shadow_length': round(shadow_length, 2)
                    })

        return results

    def find_shortest_shadow(self, date_str):
        """查找最短影长（正午时刻）"""
        # 计算正午影长
        shadow_length, altitude = self.calculate_shadow_length(date_str, 12)

        if shadow_length is None:
            return None

        # 计算时差以确定真太阳时正午
        month, day = map(int, date_str.split('-'))
        julian_day = self.date_to_julian(month, day)
        eot = self.calculate_equation_of_time(julian_day)
        true_noon = self.calculate_solar_noon(date_str, eot)

        return {
            'date': date_str,
            'true_solar_noon': f"{int(true_noon):02d}:{int((true_noon % 1) * 60):02d}",
            'solar_altitude': round(altitude, 2),
            'shadow_length': round(shadow_length, 2),
            'gnomon_height': self.gnomon_height
        }

def main():
    """主程序"""
    print("=" * 50)
    print("圭表影长测量程序")
    print("地点：北京 (纬度: 39.9042°N, 经度: 116.4074°E)")
    print(f"圭表高度: 10 cm")
    print("=" * 50)

    # 创建计算器实例
    calculator = GnomonShadowCalculator()

    while True:
        print("\n请选择操作:")
        print("1. 计算特定日期的影长")
        print("2. 查看一天的影长变化")
        print("3. 退出程序")

        choice = input("请输入选择 (1-3): ").strip()

        if choice == '1':
            # 输入日期
            date_input = input("请输入日期 (格式: MM-DD, 例如: 06-21): ").strip()

            try:
                # 验证日期格式
                datetime.strptime(date_input, "%m-%d")

                # 计算最短影长（正午）
                result = calculator.find_shortest_shadow(date_input)

                if result:
                    print("\n" + "=" * 50)
                    print(f"日期: {result['date']}")
                    print(f"真太阳时正午: {result['true_solar_noon']}")
                    print(f"太阳高度角: {result['solar_altitude']}°")
                    print(f"圭表高度: {result['gnomon_height']} cm")
                    print(f"正午影长: {result['shadow_length']} cm")
                    print("=" * 50)

                    # 计算其他时间的影长
                    print("\n其他时间的影长:")
                    print("-" * 40)
                    print(f"{'时间':<8} {'太阳高度角(°)':<15} {'影长(cm)':<12}")
                    print("-" * 40)

                    for hour in [9, 10, 11, 12, 13, 14, 15]:
                        shadow, alt = calculator.calculate_shadow_length(date_input, hour)
                        if shadow:
                            print(f"{hour:02d}:00   {alt:<15.2f} {shadow:<12.2f}")
                else:
                    print("无法计算影长（太阳可能在地平线以下）")

            except ValueError:
                print("日期格式错误，请使用 MM-DD 格式")

        elif choice == '2':
            # 输入日期查看影长变化
            date_input = input("请输入日期 (格式: MM-DD, 例如: 06-21): ").strip()

            try:
                # 验证日期格式
                datetime.strptime(date_input, "%m-%d")

                # 计算一天的影长变化
                daily_data = calculator.calculate_daily_shadow_curve(date_input)

                if daily_data:
                    print(f"\n{date_input} 影长变化曲线:")
                    print("=" * 60)
                    print(f"{'时间':<8} {'太阳高度角(°)':<15} {'影长(cm)':<12} {'相对长度':<20}")
                    print("=" * 60)

                    # 找到最大影长用于归一化显示
                    max_shadow = max(item['shadow_length'] for item in daily_data)

                    for item in daily_data:
                        # 创建简单的条形图
                        bar_length = int(20 * item['shadow_length'] / max_shadow) if max_shadow > 0 else 0
                        bar = "█" * bar_length

                        print(f"{item['time']:<8} {item['solar_altitude']:<15.2f} "
                              f"{item['shadow_length']:<12.2f} {bar}")
                else:
                    print("无法计算影长数据")

            except ValueError:
                print("日期格式错误，请使用 MM-DD 格式")

        elif choice == '3':
            print("感谢使用圭表影长测量程序！")
            break

        else:
            print("无效选择，请重新输入")

if __name__ == "__main__":
    # 示例：计算夏至日（6月21日）的影长
    calculator = GnomonShadowCalculator()

    # 测试几个重要日期
    test_dates = ["03-21", "06-21", "09-23", "12-22"]

    print("重要节气日期的正午影长（北京）:")
    print("-" * 60)
    print(f"{'节气':<10} {'日期':<10} {'太阳高度角(°)':<15} {'影长(cm)':<12}")
    print("-" * 60)

    for date in test_dates:
        result = calculator.find_shortest_shadow(date)
        if result:
            # 给节气命名
            if date == "03-21":
                solar_term = "春分"
            elif date == "06-21":
                solar_term = "夏至"
            elif date == "09-23":
                solar_term = "秋分"
            elif date == "12-22":
                solar_term = "冬至"
            else:
                solar_term = "未知"

            print(f"{solar_term:<10} {date:<10} {result['solar_altitude']:<15.2f} "
                  f"{result['shadow_length']:<12.2f}")

    print("\n")

    # 运行主程序
    main()