简介太阳高度角与方位角的高精度计算代码以MATLAB实现面向从事气象分析、建筑节能、光伏系统设计或GIS开发的研究人员与工程师解决追光、日照模拟及辐射估算中的太阳位置求解问题。代码以核心的M脚本为载体基于日期、经纬度、海拔等输入参数结合IAU推荐的天文算法计算太阳位置可直接运行并输出可靠结果适用于太阳能电池板朝向优化、建筑采光模拟和光伏电站选址评估等场景。压缩包共3个文件包含1个M脚本、1个PDF理论说明文档和1个TXT许可证文件整体仅1.28MB附件精简、便于快速查验和二次开发。PDF文档对Haversine、Vincenty公式及MATLAB调用方式做了展开帮助读者理解求算原理与误差控制从而迁移到其他地理或能源项目中。已有3712人学习使用适合需要快速获得高精度太阳位置计算工具并掌握底层方法的中高级MATLAB用户。1. 太阳高度角方位角高精度计算先算对你头顶那颗太阳的“表”再谈 MATLAB 代码做光伏逐时仿真、建筑日照分析或遥感影像校正时太阳高度角和方位角是所有计算的第一级输入。很多开发者从网上下载简化公式一算就发现上午十点的方位角偏出二十度正午前后的高度角还不对称——问题不在公式本身而在公式背后的时间系统没有对齐。高精度计算的核心不是把某个级数展开背得更熟而是把 UTC、平太阳时、真太阳时、赤纬角、时角这条链路每个环节都校正好。这篇笔记面向需要在 MATLAB 里落地逐时、逐分钟太阳位置的工程人员给出可直接复用的算法、参数边界和验证方法。2. 高精度计算背后的天文链条真太阳时、赤纬角与方位角判向2.1 为什么简化公式会差出十几分钟平太阳时、真太阳时与时差方程太阳位置计算里最容易翻车的是拿“钟表中午”当成“太阳中天”。钟表走的是平太阳时太阳实际过中天的时刻叫真太阳时两者之差就是时差方程Equation of TimeEoT。地球公转轨道是椭圆的且黄道面与赤道面之间存在约 23.44° 的交角两个因素叠加使真太阳时与平太阳时全年最多相差约 16 分钟。光伏功率曲线在正午附近斜率很陡16 分钟的错位足以让某几个时刻的仿真功率偏离实测值 10% 以上。平太阳时的当地版本由经度决定某地的平太阳时等于 UTC 加上 经度/15 小时东经为正。比如东经 120° 的地方正午平太阳时对应 UTC 时间 04:00。真太阳时再叠加 EoT即真太阳时 UTC lon/15 EoT/60EoT 的工程近似式常用EoT ≈ 9.87 * sin(2B) - 7.53 * cos(B) - 1.5 * sin(B)其中 B 360 * (N - 81) / 364N 为一年中的日序号。这个式子在全行业都常见优点是手算快、适合查表缺点是误差大约在几十秒的量级对应太阳角度误差约零点零几度。如果只是白天/黑夜判断完全够用但要做到“高精度”三个字就必须回到天文算法定义由太阳视赤经与平黄经之差反推 EoT。这个差别我会在第 3 章代码里体现。这里有一个关键认知时差方程不是玄学它由开普勒方程主导的轨道项和黄赤交角主导的投影项构成。理解了这一点你就知道为什么任何“全天 24 点统一加 5 分钟”的土办法都不可靠——EoT 的变化率不是常数春秋季和冬夏季的日变化形态完全不同。2.2 赤纬角的高精度取法从单正弦到黄道级数展开赤纬角 δ 是太阳直射点的地理纬度全年在 ±23.44° 之间摆动。网上流传最广的 Cooper 公式是δ 23.45 * sin(360 * (284 N) / 365)它只有一个正弦项把一个椭圆轨道里的太阳视运动强行拟合成匀速正弦误差最大能到 1° 左右。1° 的赤纬误差在正午高度角公式里就直接带来 1° 的高度误差对光伏阵列倾角优化这种要求 0.1° 量级精度的场景显然不合格。高精度做法的链条是先算儒略日世纪数 T再求太阳几何平黄经 L0 和平近点角 M由 M 通过中心差级数修正得到太阳真黄经 λ最后用黄道坐标到赤道坐标的转换求赤纬δ arcsin( sin(ε) * sin(λ) )其中 ε 是黄赤交角λ 是太阳真黄经。中心差 C 是级数展开前三项为C (1.914602 - 0.004817T) * sin(M) (0.019993 - 0.000101T) * sin(2M) 0.000289 * sin(3M)第一项来自开普勒方程的偏心修正第二、第三项是二阶和三阶周期项。对 2020 年前后的日期这一组级数可以把太阳黄经做到 0.01° 量级。如果还要更高可以在 λ 里追加光行差修正 -0.00569° 和章动主项 -0.00478 * sin(Ω)其中 Ω 是月球升交点黄经。这两项对结果影响约 0.01°普通工程可以忽略但代码里我会留出可选项。2.3 高度角与方位角公式用 atan2 做掉最麻烦的象限分支太阳高度角 h 的标准公式是h arcsin( sin(φ) * sin(δ) cos(φ) * cos(δ) * cos(ω) )φ 是地理纬度北纬为正ω 是时角正午为 0下午为正。这个公式本身非常稳定难点在方位角。围绕方位角有很多版本必须注意零度起算方向。下面这个版本以正北为 0°顺时针增加A atan2( -sin(ω), tan(δ) * cos(φ) - sin(φ) * cos(ω) )使用 atan2 而不是 acos 版本是有原因的acos 版本先算 cos(A)反余弦的结果天然落在 0° 到 180°同一个值对应上午和下午两个可能的方位你必须再拿时角符号去人为修正。而 atan2 把分子分母同时交给函数时角 ω 的符号直接决定方位角落在哪个象限少一道手工分支就少一类边界错误。还有一个很容易踩的定义坑不少经典文献和海外软件里方位角从正南起算向东为正、向西为负。同一组角度数据从“北为 0°”换成“南为 0°”数值上有的要加 180°有的要取负转换错一个符号结果就从东南跑到西南。所以拿到别人的计算函数第一件事不是看公式而是确认它的零度定义和旋转方向。3. 用 MATLAB 实现高精度太阳位置函数完整代码与逐段参数说明3.1 函数接口与时间转换输入本地时刻统一到 UTC 儒略日我习惯把整个计算拆成两层一层负责把“当地钟表时间”换算成 UTC 儒略日另一层负责纯天文计算。这样做的血泪经验是一旦把时区换算混进天文公式出错时你根本分不清是 JD 不对还是赤纬不对。function JD local_time_to_jd(y, mo, d, hh, mm, ss, tz_offset) % 本地钟表时间 - UTC儒略日 % 输入y年, mo月, d日, hh时, mm分, ss秒, tz_offset时区东八区传8西五区传-5 % 输出JD为UTC儒略日标准天文儒略日正午起点格里高利历 t_local datetime(y, mo, d, hh, mm, ss); % 无时区的墙钟时间 t_utc t_local - hours(tz_offset); % 换算到UTC墙钟 JD juliandate(t_utc); % datetime - 标准JD end这段代码的逻辑很直白先把年月日时分秒构造成一个不带时区的 datetime再用hours(tz_offset)把本地时间减回去得到 UTC 的墙钟时间。juliandate会把无时区的 datetime 按字面值当作 UTC 解释所以t_utc的值必须已经是 UTC。这里最容易错的点就是时区方向东八区的本地时间比 UTC 早 8 小时所以要减 8而不是加 8。参数表如下参数含义单位/约定lat_deg地理纬度度北纬为正南纬为负lon_deg地理经度度东经为正西经为负tz_offset时区偏移小时UTC8 传 8UTC-5 传 -5JD儒略日连续天数UTC 基准正午起点h_deg太阳高度角度地平以上为正A_deg太阳方位角度正北0顺时针3.2 核心计算代码太阳赤纬、时差、高度角与方位角下面这个函数接收 UTC 儒略日和经纬度输出高度角、方位角和赤纬角。它只依赖 JD不碰任何时区问题因此可以放心批量调用。function [h_deg, A_deg, delta_deg] solar_angles_jd(JD, lat_deg, lon_deg) % 高精度太阳高度角/方位角计算基于经典天文算法 % 输入JD为UTC儒略日可接受向量lat_deg纬度lon_deg经度 % 输出h_deg高度角度A_deg方位角度正北0顺时针delta_deg赤纬角度 % 1. 儒略日世纪数 T (JD - 2451545.0) / 36525.0; % 2. 太阳几何平黄经和平近点角度归一化到0~360 L0 mod(280.46646 36000.76983 * T 0.0003032 * T^2, 360); M mod(357.52911 35999.05029 * T - 0.0001537 * T^2 4.0e-8 * T^3, 360); % 3. 太阳中心差与真黄经 Mrad deg2rad(M); C (1.914602 - 0.004817 * T - 0.000014 * T^2) * sin(Mrad) ... (0.019993 - 0.000101 * T) * sin(2 * Mrad) ... 0.000289 * sin(3 * Mrad); lambda mod(L0 C, 360); % 太阳真黄经 % 可选更高精度的视黄经修正光行差章动主项 % Omega 125.04 - 1934.136 * T; % lambda lambda - 0.00569 - 0.00478 * sind(Omega); % lambda mod(lambda, 360); % 4. 黄赤交角、赤纬、视赤经 eps 23.439291 - 0.0130042 * T; % 黄赤交角 delta_deg asind(sind(eps) * sind(lambda)); % 赤纬 alpha atan2d(sind(lambda) * cosd(eps), cosd(lambda)); % 视赤经 alpha mod(alpha 360, 360); % 5. 时差方程用视赤经反推单位度再转分钟 EoT_deg L0 - 0.0057183 - alpha; EoT_deg mod(EoT_deg 180, 360) - 180; % 归一到[-180, 180) EoT_min EoT_deg * 4; % 6. 真太阳时与时角 utc_h mod(JD 0.5, 1) * 24; % 当日UTC小时0~24 true_solar_h mod(utc_h lon_deg / 15 EoT_min / 60, 24); omega 15 * (true_solar_h - 12); % 时角度下午为正 % 7. 高度角 sinh sind(lat_deg) * sind(delta_deg) cosd(lat_deg) * cosd(delta_deg) * cosd(omega); h_deg asind(sinh); % 8. 方位角正北0度顺时针 x tand(delta_deg) * cosd(lat_deg) - sind(lat_deg) * cosd(omega); y -sind(omega); A_deg mod(atan2d(y, x), 360); end代码逻辑按天文算法标准链路展开JD 先转世纪数 TT 决定了所有周期性参数L0 和 M 是描述太阳轨道的基本量C 修正地球公转椭圆带来的不均匀性得到太阳真黄经 λ再由黄道坐标转到赤道坐标得到赤纬 δ 和赤经 α。EoT 由 L0、固定光行差常数和 α 共同决定而不是用日序号近似这是“高精度”和我们平时抄来的简易公式最本质的区别。几个参数要点utc_h的mod(JD 0.5, 1)是天文 JD 转小时的老写法原因是标准 JD 以正午为起点加 0.5 才是从午夜算起的日分数改错这个偏移会让所有时刻差 12 小时EoT 归一化到 [-180, 180) 是为了避免跨 0° 时数值跳变方位角最后mod(A_deg 360, 360)或mod(atan2d(y,x), 360)都能把 -180 到 180 转到 0 到 360我建议直接写后者少一步加法。3.3 调用示例与输出验证单点检查与逐时序列批量计算单点计算用来验证数学批量计算才真正用于工程。先看单点计算某个北方城市北纬 39.9°东经 116.4°在 2024 年夏至日中午 12 点东八区的太阳位置。JD local_time_to_jd(2024, 6, 21, 12, 0, 0, 8); [h, A] solar_angles_jd(JD, 39.9, 116.4); fprintf(高度角: %.3f 度, 方位角: %.3f 度\n, h, A);注意这里“中午 12 点”是钟表时间不是真太阳时正午因此计算出的高度角并不是当天最大值方位角也不会正好是 180°。这是初学者最容易误读的输出钟表 12 点和太阳中天之间差着一截经度修正和 EoT。想找真正的太阳中天得在 12 点前后扫描高度角取最大值后面第 5 章会用到这个技巧。批量场景是光伏逐时仿真最常见的形式算某个地点全年 8760 小时的太阳位置。得益于 MATLAB 的数组运算JD 直接传向量函数内部所有mod和三角函数都逐元素工作一行代码就能出全年结果。lat 31.2; lon 121.5; tz 8; % 某东部沿海城市 t_all datetime(2024, 1, 1, 0, 0, 0) hours(0:24*366-1); % 2024闰年逐小时 JD_all juliandate(t_all - hours(tz)); % 本地转UTC [h_all, A_all] solar_angles_jd(JD_all, lat, lon); % 按天reshape画出第1、7、14...天的高度角日曲线 h_mat reshape(h_all, 24, []); plot(0:23, h_mat(:, 1:7:end)); xlabel(钟表小时东八区); ylabel(太阳高度角度);这段代码里有一个容易被忽略的事实reshape(h_all, 24, [])的 24 行对应的是 UTC 小时还是本地小时因为t_all是从东八区 0 点开始的每个 24 点块天然对齐到当地日期而JD_all只用于计算不参与 reshape所以画出来的日曲线是东八区钟表时间的视角。这样处理既正确也符合使用习惯但如果你把t_all换成 UTC 起点日曲线就会整体错位这是时区问题的另一个隐蔽形态。4. 太阳位置计算排查实录五个高频翻车现场与避坑方法4.1 跨年结果跳变儒略日里的月份调整现象某工程师用自写的 JD 公式算太阳位置发现 12 月 31 日和次年 1 月 1 日的高度角曲线在年界处明显不连续差了有半度到一度。原因经典儒略日公式要求月份小于等于 2 时把年份减 1、月份加 12当作上一年的 13、14 月处理。很多人抄公式时漏掉这条前置分支直接用 1 月、2 月的原始数字套公式导致 1 月的 JD 出现系统性偏差。解决优先使用datetime加juliandate它内部已经处理了历法边界如果一定要维护手写 JD把月份调整这个 if 分支放在所有计算之前并专门写一个跨年用例回归测试。我这个函数里没有手写 JD就是想把这个最容易出错的环节交给 MATLAB 本身。4.2 方位角象限错乱同一个高度角有两个方位解现象某日照分析项目里早上 8 点算出的方位角落在西北方向200 多度下午反而跑到东北方向整个日轨迹图像是被镜像过。原因用了反余弦版本的方位角公式。反余弦值域只有 0 到 180 度同一高度角对应的上午解和下午解全靠人工判断只要时角符号的转换规则写错象限就翻车。更隐蔽的是有些资料里的方位角以正南为 0 度早上为正、下午为负照抄数值再和以正北为 0 度且顺时针的行业标准比较结果自然全错。解决用第 3 章的 atan2 版本且明确在代码注释里写下“正北0顺时针”。如果你要和国外开源库对算先统一零度定义再谈数值差异。4.3 正午峰值不对称时区和经度符号混用现象某光伏项目的逐时功率仿真正午峰值出现在下午 1 点 30 分且上午和下午曲线不对称换一个地点偏移量还不一样。原因直接用“本地钟表时间”代替了“当地平太阳时”。比如项目地在东经 116°却用东八区基准 120° 的钟表经度差 4° 造成 16 分钟固定偏移再叠加时差方程正午就可能偏移到 13 点以后。解决内部所有计算统一走 UTC真太阳时由UTC lon/15 EoT/60单独求出不要在 UTC 里混入时区再减回来。还要记住时区是政治定义经度才是物理量新疆、西藏地区用东八区钟表日出日落和太阳位置能差出一小时以上这种地方做日照分析绝不能靠时区换太阳时。4.4 日出日落判断不一致要不要加大气折射修正现象按高度角大于 0 度判断“太阳已升起”但实测日出时间总是比这个判据早几分钟而且清晨傍晚的方位角与观测值对不上。原因高度角 0 度是太阳中心的几何位置地面观测者看到太阳上边缘跃出地平线时太阳中心实际在地平线下约 0.83 度其中约 0.57 度来自大气折射约 0.27 度来自太阳半径。解决做光伏发电仿真这类工程计算用几何位置即可但要在判断白天时统一采用高度角大于 0 度的定义做观星、日照时长校核时把临界值改为 -0.83 度再判断。两种口径没有对错关键是整条计算链一致不要几何位置和视位置混用。4.5 MATLAB datetime 的隐性时区偏移现象同一批 UTC 时刻数据有人用datetime直接转换后得到的时间比预期快或慢了 8 小时回头查 JD 又查不出问题。原因datetime(y, mo, d, ...)创建的是不带时区的墙钟时间juliandate会把它按字面值当 UTC 解释但用户心里想的是本地时间直接喂进去于是 JD 混入时区偏移赤纬和时角全部错位。解决在函数入口处约定“本函数只认 UTC”并在调用前显式写t_utc t_local - hours(tz_offset)。我第 3 章的local_time_to_jd专门把这个换算写成一行的原因就在这里每次调用都过一遍时区减法比把时区藏在函数内部更不容易出错。5. 用三个基准把精度验到心里有数再固化成一测就用的脚本5.1 基准一正午高度角和日出日落时刻找真太阳时正午不是一个简单的时间点正确做法是在钟表正午前后扫描高度角取最大值作为当天太阳中天高度。以某北方城市夏至日为例理论正午高度角约等于 90° 减去纬度与赤纬之差的绝对值即 90 - |39.9 - 23.44| ≈ 73.5°。scan datetime(2024, 6, 21, 10, 0, 0) minutes(0:240); % 东八区10:00-14:00 JD_scan juliandate(scan - hours(8)); [h_scan, ~] solar_angles_jd(JD_scan, 39.9, 116.4); [h_max, i_max] max(h_scan); fprintf(太阳中天高度: %.3f°, 时刻: %s\n, h_max, string(scan(i_max)));这段代码本身就是一个好习惯不拿钟表 12 点冒充正午而是用扫描的方式把中天时刻找出来。如果输出的最高高度与理论值偏差超过 0.1°说明 JD 或赤纬链路有问题中天时刻偏离 12 点越远说明经度或时差环节越可能出错。5.2 基准二二分二至日的赤纬角核对赤纬角是最独立的校验对象。春分秋分时太阳赤纬应接近 0°夏至约 23.44°冬至约 -23.44°。这个数据不需要经纬度参与直接调用函数取第三个输出即可。只要这四个点的赤纬误差在 0.01° 以内黄经链路基本可信。如果你要做得更细可以把当地正午附近的太阳轨迹画成一张全年逐时分布图看日出方位角从冬至的东南方向一步步移到夏至的东北方向再移回去。方向连贯、没有突变说明方位角象限处理是稳的如果某一天突然跳 180°基本可以断定又是零度定义或归一化出了问题。5.3 我把验证写进单元测试的习惯我会把这些基准固化成test_solar_position.m脚本每次改了代码或换了机器第一件事就是跑一遍回归。断言写成assert(abs(h_max - 73.5) 0.1, ...)这样的形式让失败信息直接指出哪个环节越界。这套组合拳我用了很久挡住过多次自以为是的优化——要知道太阳位置计算里一个“看起来更精细”的修正项有时反而会引入更大的误差。先验基准再动代码是我在这个方向上最实用的教训。希望帮到你。本文还有配套的精品资源点击获取