简介：这份MATLAB代码资源聚焦太阳高度角与方位角的高精度计算，面向气象学、建筑采光、地理信息系统及太阳能工程等领域的科研人员与工程师，用于解决特定经纬度与时间下太阳位置精确求解的问题。压缩包共3个文件，包含1个m脚本、1个pdf文档和1个txt授权说明，整体约1.28MB；其中m文件为核心计算程序，pdf提供理论背景与公式推导，txt则明确使用与分发条件。代码基于日期时间函数与地球坐标参数，结合国际天文联合会推荐公式推算太阳位置，可用于优化光伏板朝向、分析建筑日照与预测辐射变化。目前已有3714人学习下载，读者可借此掌握从理论公式到MATLAB实现的完整链路，并理解天文算法在工程场景中的落地方式。

1. 太阳高度角方位角计算：从天文算法到MATLAB工程实现

做光伏电站选址评估时，我需要在没有专业天文软件的环境下，快速算出一整年逐小时的太阳高度角和方位角。市面上的在线计算器要么精度不够，要么无法批量处理，于是我开始找一套能直接在MATLAB里跑、精度足够、还能自己改的代码。这份太阳高度角方位角高精度MATLAB计算代码就是在这个场景下进入视野的——它用天文算法实现了太阳位置的高精度求解，输入经纬度和时间，输出高度角、方位角、赤纬角、时角等参数。适合做光伏发电量估算、建筑日照分析、太阳能跟踪系统开发，以及需要批量计算太阳轨迹的科研和工程场景。代码结构清晰，参数可调，拿来就能跑，改起来也不费劲。

2. 算法核心拆解：太阳位置是怎么一步步算出来的

2.1 从儒略日到太阳几何位置

太阳位置计算的第一步是把日历时间转成儒略日（Julian Day），这是天文计算的通用时间基准。很多新手直接拿本地时间代入公式，结果偏差大到离谱——问题就出在没做时区转换和儒略日换算。

% 将公历日期转换为儒略日 % 输入：year, month, day, hour, minute, second (UTC时间) % 输出：jd - 儒略日 function jd = cal2jd(year, month, day, hour, minute, second) if month <= 2 year = year - 1; month = month + 12; end A = floor(year / 100); B = 2 - A + floor(A / 4); jd = floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1)) ... + day + B - 1524.5; jd = jd + (hour + minute/60 + second/3600) / 24; end

这段代码的关键在于：月份小于等于2时要调整年份和月份，这是儒略日算法的标准处理；B项是格里高利历的修正项，少了它会有十几天的累积误差。参数hour、minute、second必须传UTC时间，如果你手头是北京时间，记得先减8小时再传进来。

拿到儒略日后，下一步是计算儒略世纪数T，然后求太阳平黄经、平近点角、黄赤交角等中间量。这些量在代码里通常封装成一个结构体或数组返回，方便后续调用。

% 计算太阳几何位置参数 % 输入：jd - 儒略日 % 输出：struct，包含平黄经、平近点角、黄赤交角等 function sun = sun_geometry(jd) T = (jd - 2451545.0) / 36525; % 儒略世纪数 % 太阳平黄经（度） L0 = mod(280.46646 + 36000.76983 * T + 0.0003032 * T^2, 360); % 太阳平近点角（度） M = mod(357.52911 + 35999.05029 * T - 0.0001537 * T^2, 360); % 地球轨道偏心率 e = 0.016708634 - 0.000042037 * T - 0.0000001267 * T^2; % 黄赤交角（度） epsilon = 23.439291 - 0.0130042 * T - 1.64e-7 * T^2 + 5.04e-7 * T^3; sun.L0 = L0; sun.M = M; sun.e = e; sun.epsilon = epsilon; end

这里每个系数都是天文历表拟合出来的，T的幂次越高代表越精细的修正。实际使用中，T的范围一般在-1到1之间（即1900年到2100年），超出这个范围精度会下降。如果你只需要粗略估算，可以只保留T的一次项，但既然这份代码主打高精度，建议完整保留。

2.2 赤纬角和时角的求解

赤纬角是太阳直射点纬度，时角是太阳相对于当地子午线的角度。这两个量直接决定了太阳高度角和方位角的计算结果。

% 计算太阳赤纬角和时角 % 输入：sun - 几何位置结构体，jd - 儒略日，lon - 经度(东经为正) % 输出：dec - 赤纬角(度)，ha - 时角(度) function [dec, ha] = sun_dec_ha(sun, jd, lon) % 中心差（度） C = (1.914602 - 0.004817 * sun.T - 0.000014 * sun.T^2) * sind(sun.M) ... + (0.019993 - 0.000101 * sun.T) * sind(2 * sun.M) ... + 0.000289 * sind(3 * sun.M); % 太阳真黄经 true_lon = sun.L0 + C; % 太阳真近点角 true_anomaly = sun.M + C; % 日地距离（天文单位） R = 1.000001018 * (1 - sun.e^2) / (1 + sun.e * cosd(true_anomaly)); % 赤纬角 dec = asind(sind(sun.epsilon) * sind(true_lon)); % 时角：先求格林尼治时角，再转当地时角 gmst = mod(280.46061837 + 360.98564736629 * (jd - 2451545.0), 360); ha = mod(gmst + lon - true_lon, 360); if ha > 180 ha = ha - 360; end end

赤纬角计算里用到了中心差C，它修正了地球椭圆轨道导致的太阳黄经偏差。时角计算中，gmst是格林尼治平恒星时，加上经度得到当地恒星时，再减去太阳真黄经得到时角。注意时角需要归一化到-180到180度之间，否则后续三角函数计算会出问题。

2.3 高度角和方位角的最终输出

有了赤纬角和时角，高度角和方位角就是球面三角的常规计算了。

% 计算太阳高度角和方位角 % 输入：dec - 赤纬角(度)，ha - 时角(度)，lat - 纬度(度) % 输出：alt - 高度角(度)，az - 方位角(度，北为0，顺时针) function [alt, az] = sun_alt_az(dec, ha, lat) sin_alt = sind(dec) * sind(lat) + cosd(dec) * cosd(lat) * cosd(ha); alt = asind(sin_alt); cos_az = (sind(dec) - sind(alt) * sind(lat)) / (cosd(alt) * cosd(lat)); cos_az = max(-1, min(1, cos_az)); % 防止数值溢出 az = acosd(cos_az); if ha > 0 az = 360 - az; % 下午方位角翻转 end end

方位角的象限判断靠时角正负来区分：时角为负（上午）方位角取acosd结果，时角为正（下午）用360度减。这个细节很多简化代码会忽略，导致下午的方位角算出来是镜像的。另外cos_az做了一次钳位，防止浮点误差导致acosd输入超出[-1,1]范围而报NaN。

3. 工程化落地：从单点计算到批量处理

3.1 主函数封装与参数设计

把上面几个函数串起来，形成一个对外的主接口。我一般会设计成输入经纬度、时区、起止时间、时间步长，输出一个表格或结构体数组。

% 主函数：批量计算太阳高度角和方位角 % 输入：lat, lon - 纬度、经度(度) % tz - 时区(如北京为8) % t_start, t_end - 起始和结束时间(datetime) % step_min - 时间步长(分钟) % 输出：result - 包含时间、高度角、方位角的table function result = solar_position_batch(lat, lon, tz, t_start, t_end, step_min) times = t_start : minutes(step_min) : t_end; n = length(times); alt = zeros(n, 1); az = zeros(n, 1); for i = 1:n t_utc = times(i) - hours(tz); jd = cal2jd(t_utc.Year, t_utc.Month, t_utc.Day, ... t_utc.Hour, t_utc.Minute, t_utc.Second); sun = sun_geometry(jd); [dec, ha] = sun_dec_ha(sun, jd, lon); [alt(i), az(i)] = sun_alt_az(dec, ha, lat); end result = table(times', alt, az, 'VariableNames', {'Time', 'Altitude', 'Azimuth'}); end

这个批量接口的好处是一次性算出整段时间序列，直接喂给光伏模型或日照分析模块。step_min设成1就是逐分钟，设成60就是逐小时。注意时间步长越小计算量越大，但太阳位置变化是连续的，一般1分钟精度已经绰绰有余。

3.2 时区与经度符号的坑

时区和经度符号是翻车重灾区。经度东经为正、西经为负，这个在代码里必须统一。时区方面，输入时间如果是本地时间，要先减去时区偏移转成UTC再算儒略日。我见过有人把北京时间直接当UTC代入，结果太阳高度角差了十几度，排查了半天才发现是时区问题。

提示：如果你的应用场景只关心相对变化趋势，时区误差影响不大；但如果要做光伏发电量精确估算，时区必须严格对齐。

3.3 输出结果的可视化验证

算完数据后，画一张太阳轨迹图是最直观的验证方式。

% 绘制太阳轨迹图 result = solar_position_batch(39.9, 116.4, 8, datetime(2024,6,21,0,0,0), ... datetime(2024,6,21,23,59,0), 10); idx = result.Altitude > 0; % 只画白天 polarplot(deg2rad(result.Azimuth(idx)), result.Altitude(idx), 'LineWidth', 1.5); title('太阳轨迹图 (夏至日)');

极坐标图里，角度是方位角，半径是高度角。夏至日中午高度角应该接近当地纬度与赤纬角之和的余角，比如北纬40度左右，夏至日正午高度角大约73度。如果画出来明显不对，优先检查时区和经度符号。

4. 避坑与排查：那些年我踩过的太阳位置计算坑

4.1 高度角出现负值但白天时间对不上

现象：算出来的高度角在白天也有负值，或者日出日落时间和实际差了一两个小时。原因：最常见的是时区没转UTC，其次是经度符号搞反了。还有一种情况是datetime类型的小时数读取错误，比如把AM/PM搞混。解决：在cal2jd入口处打印一下传入的UTC时间，和实际UTC时间对一下。经度方面，东经116度就传116，不要传-116。

4.2 方位角在正午前后跳变

现象：方位角在正午附近从接近0度突然跳到接近360度。原因：这是正常的，正午时太阳在正南（北半球），方位角就是180度附近，上午从东往南转，下午从南往西转。但如果跳变幅度异常大，可能是时角归一化没做好。解决：检查时角是否归一化到-180到180度。如果时角用了0到360度的表示，方位角判断逻辑会出错。

4.3 批量计算速度慢

现象：逐分钟算一年数据要跑好几分钟。原因：循环里反复调用函数，MATLAB的循环开销加上函数调用开销累积起来很可观。解决：把时间序列向量化，儒略日计算和几何位置计算都可以用数组运算替代循环。或者至少把cal2jd和sun_geometry改成接受向量输入。

4.4 高纬度地区结果异常

现象：在纬度70度以上地区，某些日期算不出结果或结果明显错误。原因：极昼极夜期间，太阳可能整天在地平线以上或以下，方位角公式中的cos_az可能超出[-1,1]范围。解决：代码里已经做了钳位处理，但如果你的应用需要精确处理极昼极夜，建议增加判断逻辑：当sin_alt大于1或小于-1时直接返回±90度。

4.5 闰年和世纪年处理

现象：算2000年2月29日的数据时报错或结果偏差。原因：儒略日算法对闰年的处理依赖格里高利历修正项B，如果B算错了，2月29日附近会有1天误差。解决：检查cal2jd里的A和B计算，A = floor(year/100)，B = 2 - A + floor(A/4)。这个公式对1900年和2100年这类世纪年也适用，但要注意1900年不是闰年。

5. 二次开发：把这套代码改造成你要的样子

拿到这份代码后，最直接的改造方向是接入你自己的数据流。比如你有一个光伏电站的SCADA系统，想把太阳位置作为特征输入到发电量预测模型里，那就把solar_position_batch的输出直接拼到特征表里。我一般会加一列清晰度指数或云量修正，因为纯天文计算得到的是理论值，实际辐照度还受天气影响。

另一个常见需求是实时计算。如果你的系统需要每秒更新太阳位置，可以把儒略日计算简化——因为时间间隔很短，T的变化量极小，可以只算一次基准值然后做线性外推。这样单次计算耗时能降到微秒级。

还有一个改造点是精度分级。这份代码默认用了较高阶的修正项，如果你做的是建筑遮阳分析，精度要求没那么高，可以把T的二次项和三次项去掉，计算量能减少三分之一，精度损失在0.01度以内，完全够用。

验证改造结果的方法很简单：拿几个已知的基准点对一下。比如春分日正午，太阳高度角应该等于90度减去当地纬度；夏至日正午，高度角等于90度减纬度加23.44度。如果这几个点对上了，说明核心算法没问题。

从那以后我每次改完太阳位置代码，都会先跑一遍春分、夏至、秋分、冬至四个基准点，确认高度角误差在0.1度以内再往下做。这个习惯帮我省了很多返工时间。希望帮到你。

本文还有配套的精品资源，点击获取