简介：这份资源提供一套基于NOAA算法的太阳方位角与高度角计算C代码，面向从事太阳能追日、测地定位、气象观测等需要精确太阳位置解算的开发者与工程技术人员。代码输入经纬度、日期、时间与时区即可输出太阳方位角与高度角，方位角以正北为零度参考，南北东西半球统一适用，计算精度较高。压缩包内共1个docx文件，约20KB，内容以C源码为主，涵盖角度与弧度转换、积日计算、日期转儒略日、儒略世纪数换算、太阳几何平黄经等核心函数，并附有NOAA官方算法参考链接，便于对照验证。目前已有1973人学习下载，适合需要将太阳位置计算模块集成到嵌入式或上位机项目中的读者参考，可直接移植或按需裁剪，快速完成追日控制、光照分析等功能的开发。

1. 太阳方位角高度角计算C代码：从天文算法到嵌入式落地的完整路径

做光伏追光、无人机航迹规划或者建筑日照分析时，绕不开一个基础问题：给定时间、经纬度，太阳此刻在天上的哪个位置？太阳方位角高度角计算C代码就是解决这件事的。它把天文历法里的太阳位置算法用C语言实现，输入年月日时分和经纬度，输出高度角（仰角）和方位角。听起来简单，但真正落地时会遇到一堆问题：时间用世界时还是本地时？大气折射要不要修正？精度要求到角秒还是度？嵌入式平台浮点性能够不够？我做过几个光伏追踪和户外设备姿态校准的项目，这套算法从PC验证到MCU部署踩了不少坑，下面把完整路径拆开讲清楚，包括算法选型、代码实现、参数配置和常见翻车点。

2. 太阳位置算法选型：为什么SPA和NOAA简化算法各有适用场景

2.1 三种主流算法的精度与开销对比

太阳位置计算不是只有一种方法。从最简单的低精度近似到天文台级别的SPA（Solar Position Algorithm），精度和计算量差了好几个数量级。选错了要么精度不够导致追光偏几度，要么在MCU上跑不动。

NOAA简化算法是我在工程中用得最多的。它基于Jean Meeus的《Astronomical Algorithms》，把太阳赤纬、时差、高度角、方位角的计算压缩到几十行C代码，在Cortex-M3上跑一次不到1ms。SPA精度高但涉及大量三角函数和多项式展开，在资源受限的MCU上不太现实。

选型建议很直接：光伏追光、气象站、户外设备姿态校准用NOAA简化算法足够；如果做天文望远镜指向或者需要角秒级精度，才考虑SPA或者查表插值。

2.2 时间系统：UTC、本地时和儒略日的转换关系

这是最容易翻车的地方。太阳位置计算的所有公式都基于世界时（UTC），但你的输入大概率是本地时间。如果不做时区转换，算出来的方位角能偏几十度。

核心转换链路是：本地时间 → UTC → 儒略日（Julian Day）→ 儒略世纪数T。儒略日是自公元前4713年1月1日12时起算的连续天数，是所有天文计算的时间基准。

// 本地时间转UTC再转儒略日 // year, month, day, hour, minute, second 为本地时间 // timezone 为时区偏移，东八区传 8.0 double local_to_julian_day(int year, int month, int day, int hour, int minute, double second, double timezone) { // 先转UTC double ut_hour = hour + minute / 60.0 + second / 3600.0 - timezone; // 处理跨日：如果UTC小时为负或超过24，需要调整日期 int day_offset = 0; if (ut_hour < 0.0) { ut_hour += 24.0; day_offset = -1; } else if (ut_hour >= 24.0) { ut_hour -= 24.0; day_offset = 1; } // 调整日期（简化处理，实际需考虑月末年末） day += day_offset; // 如果月份是1或2月，按Meeus算法视为上一年的13或14月 int y = year; int m = month; if (m <= 2) { y -= 1; m += 12; } // 格里高利历儒略日公式 int A = y / 100; int B = 2 - A + A / 4; double jd = (int)(365.25 * (y + 4716)) + (int)(30.6001 * (m + 1)) + day + B - 1524.5 + ut_hour / 24.0; return jd; }

这段代码有几个关键点。第一，时区偏移用小数表示，东八区是8.0，印度是5.5。第二，跨日处理不能省，否则晚上算出来的儒略日会差一天。第三，1月和2月要当作上一年的13月和14月处理，这是Meeus算法的约定。第四，B值只在格里高利历（1582年10月15日之后）有效，如果做历史数据回溯需要判断历法。

注意：儒略日公式里的day包含了日期调整，如果跨月跨年需要额外处理。工程中我一般先把本地时间转成Unix时间戳，再从时间戳算儒略日，这样跨日跨月跨年都不用自己处理。

2.3 赤纬角和时差：两个决定精度的中间量

太阳赤纬角（Declination）是太阳直射点相对于赤道的纬度，范围在±23.44°之间。时差（Equation of Time）是真太阳时和平太阳时的差值，范围约±16分钟。这两个量直接决定后续高度角和方位角的精度。

// 计算太阳赤纬角和时差（单位：度） // jd 为儒略日 void calc_declination_eot(double jd, double *declination, double *eot) { // 儒略世纪数 double T = (jd - 2451545.0) / 36525.0; // 太阳平黄经（度） double L0 = fmod(280.46646 + T * (36000.76983 + T * 0.0003032), 360.0); // 太阳平近点角（度） double M = 357.52911 + T * (35999.05029 - 0.0001537 * T); double M_rad = M * M_PI / 180.0; // 地球轨道偏心率修正 double C = sin(M_rad) * (1.914602 - T * (0.004817 + 0.000014 * T)) + sin(2 * M_rad) * (0.019993 - 0.000101 * T) + sin(3 * M_rad) * 0.000289; // 太阳真黄经 double true_longitude = L0 + C; // 黄赤交角（度） double epsilon = 23.439291 - 0.0130042 * T; // 赤纬角 double lambda_rad = true_longitude * M_PI / 180.0; double eps_rad = epsilon * M_PI / 180.0; *declination = asin(sin(eps_rad) * sin(lambda_rad)) * 180.0 / M_PI; // 时差计算 double y = tan(eps_rad / 2.0) * tan(eps_rad / 2.0); double L0_rad = L0 * M_PI / 180.0; double e = 0.016708634 - T * (0.000042037 + 0.0000001267 * T); double Etime = y * sin(2 * L0_rad) - 2 * e * sin(M_rad) + 4 * e * y * sin(M_rad) * cos(2 * L0_rad) - 0.5 * y * y * sin(4 * L0_rad) - 1.25 * e * e * sin(2 * M_rad); *eot = Etime * 4.0 * 180.0 / M_PI; // 转成分钟对应的度数 }

赤纬角计算里，L0是太阳平黄经，M是平近点角，C是中心差。这三个量决定了太阳在黄道上的真实位置。时差计算用的是Meeus的简化公式，精度在±0.5分钟以内，对绝大多数工程应用足够。

参数说明：T是儒略世纪数，从J2000.0历元起算。L0的系数36000.76983是每世纪太阳平黄经的变化率。M的系数35999.05029是每世纪平近点角的变化率。这些常数不要改，改了精度就崩了。

3. 高度角和方位角的C代码实现：从公式到可运行函数

3.1 完整计算流程与函数接口设计

把前面的模块串起来，整个计算流程是：本地时间 → UTC → 儒略日 → 儒略世纪数 → 赤纬角和时差 → 时角 → 高度角 → 方位角。

接口设计上我习惯用一个结构体打包输入输出，这样在嵌入式里传参方便，也便于后续扩展。

// 太阳位置计算结果结构体 typedef struct { double altitude; // 高度角，单位：度，地平线以上为正 double azimuth; // 方位角，单位：度，正北为0，顺时针增加 double declination; // 赤纬角，单位：度 double hour_angle; // 时角，单位：度 } SunPosition; // 输入参数结构体 typedef struct { int year, month, day; int hour, minute; double second; double latitude; // 纬度，北纬为正，单位：度 double longitude; // 经度，东经为正，单位：度 double timezone; // 时区偏移，东八区为8.0 } SunInput; // 主计算函数 int calc_sun_position(const SunInput *in, SunPosition *out) { // 参数合法性检查 if (in->latitude < -90.0 || in->latitude > 90.0) return -1; if (in->longitude < -180.0 || in->longitude > 180.0) return -2; if (in->month < 1 || in->month > 12) return -3; // 第一步：本地时间转儒略日 double jd = local_to_julian_day(in->year, in->month, in->day, in->hour, in->minute, in->second, in->timezone); // 第二步：计算赤纬角和时差 double declination, eot; calc_declination_eot(jd, &declination, &eot); // 第三步：计算真太阳时和时角 double T = (jd - 2451545.0) / 36525.0; double L0 = fmod(280.46646 + T * (36000.76983 + T * 0.0003032), 360.0); // 真太阳时（分钟） double true_solar_time = in->hour * 60.0 + in->minute + in->second / 60.0 + eot + 4.0 * in->longitude - 60.0 * in->timezone; // 归一化到0-1440分钟 true_solar_time = fmod(true_solar_time, 1440.0); if (true_solar_time < 0) true_solar_time += 1440.0; // 时角（度），正午为0，下午为正 double hour_angle = true_solar_time / 4.0 - 180.0; if (hour_angle < -180.0) hour_angle += 360.0; // 第四步：计算高度角 double lat_rad = in->latitude * M_PI / 180.0; double dec_rad = declination * M_PI / 180.0; double ha_rad = hour_angle * M_PI / 180.0; double sin_alt = sin(lat_rad) * sin(dec_rad) + cos(lat_rad) * cos(dec_rad) * cos(ha_rad); double altitude = asin(sin_alt) * 180.0 / M_PI; // 第五步：计算方位角 double cos_alt = cos(altitude * M_PI / 180.0); double azimuth; if (fabs(cos_alt) < 1e-10) { // 太阳在天顶或天底，方位角无定义 azimuth = 0.0; } else { double sin_az = -sin(ha_rad) * cos(dec_rad) / cos_alt; double cos_az = (sin(dec_rad) - sin(lat_rad) * sin_alt) / (cos(lat_rad) * cos_alt); azimuth = atan2(sin_az, cos_az) * 180.0 / M_PI; // 归一化到0-360 azimuth = fmod(azimuth + 360.0, 360.0); } // 输出 out->altitude = altitude; out->azimuth = azimuth; out->declination = declination; out->hour_angle = hour_angle; return 0; }

这个函数是整个方案的核心。逻辑说明：先做参数检查，然后依次算儒略日、赤纬角时差、真太阳时、时角，最后用球面三角公式算高度角和方位角。方位角用atan2而不是atan，因为atan只能返回±90°，atan2能覆盖全象限。

参数说明：latitude和longitude用十进制度，北纬东经为正。timezone东八区传8.0，西五区传-5.0。second用double是为了支持亚秒级精度，如果不需要可以改成int。

3.2 大气折射修正：什么时候必须加，什么时候可以省

太阳高度角低于85°时，大气折射会让太阳看起来比实际位置高。地平线附近折射量约0.57°，天顶附近接近0。做光伏追光时如果不修正，早晨和傍晚的追踪角度会偏。

// 大气折射修正（单位：度） // 输入为几何高度角，输出为视高度角 double atmospheric_refraction(double altitude) { // 低于-1度时不修正（太阳已在地平线以下） if (altitude < -1.0) return 0.0; double refraction; if (altitude > 85.0) { refraction = 0.0; } else { // 使用Saemundsson公式 double alt_rad = altitude * M_PI / 180.0; double te = tan(alt_rad); refraction = 1.02 / te; // 单位：角分 // 更精确的表达式 refraction = 1.02 / tan((altitude + 10.3 / (altitude + 5.11)) * M_PI / 180.0); } return refraction / 60.0; // 角分转度 }

Saemundsson公式在-1°到85°范围内精度约0.1角分。如果只做光伏追光，修正到0.1°就够了，用简化公式refraction = 1.02 / tan((altitude + 10.3/(altitude+5.11)) * M_PI/180.0)即可。如果做天文观测，需要用更复杂的折射模型，还要考虑温度气压。

注意：折射修正只影响高度角，不影响方位角。修正后的高度角是视高度角，追光系统应该用视高度角来对准。

3.3 方位角象限判断：atan2的坑和归一化处理

方位角计算最容易出错的地方是象限判断。用atan算出来的值范围是-90°到90°，无法区分东南西北。必须用atan2，而且要注意参数顺序。

atan2(y, x)返回的是点(x, y)与x轴正方向的夹角。在太阳方位角计算中，我们以正北为0°，顺时针增加。但数学坐标系里正北是y轴正方向，正东是x轴正方向，所以需要做坐标变换。

// 方位角归一化到0-360度，正北为0，顺时针 double normalize_azimuth(double az) { az = fmod(az, 360.0); if (az < 0.0) az += 360.0; return az; } // 另一种方位角计算方式，用球面三角余弦公式 double calc_azimuth_cosine(double latitude, double declination, double altitude, double hour_angle) { double lat_rad = latitude * M_PI / 180.0; double dec_rad = declination * M_PI / 180.0; double alt_rad = altitude * M_PI / 180.0; double cos_az = (sin(dec_rad) - sin(lat_rad) * sin(alt_rad)) / (cos(lat_rad) * cos(alt_rad)); // 数值裁剪，防止浮点误差导致acos参数越界 if (cos_az > 1.0) cos_az = 1.0; if (cos_az < -1.0) cos_az = -1.0; double az = acos(cos_az) * 180.0 / M_PI; // 根据时角判断东西方向 if (hour_angle > 0) { az = 360.0 - az; // 下午，太阳在西边 } return normalize_azimuth(az); }

余弦公式法在高度角接近90°时会出现除零问题，因为cos(altitude)趋近于0。实际工程中我一般用atan2法，稳定性更好。余弦公式法可以作为交叉验证。

4. 嵌入式部署与精度验证：从PC到MCU的移植要点

4.1 浮点运算优化：定点化还是保留float

在Cortex-M0/M3这类没有FPU的MCU上，double运算靠软件模拟，一次太阳位置计算可能要几十毫秒。如果追光系统要求每秒更新，就必须优化。

常见做法有三种：一是改用float，精度损失约0.01°，对追光足够；二是定点化，把角度放大10000倍用int32存储；三是查表加插值，把赤纬角和时差按天做成表，运行时插值。

// 定点化版本：角度放大10000倍，用int32表示 #define ANGLE_SCALE 10000 #define DEG_TO_FIXED(d) ((int32_t)((d) * ANGLE_SCALE)) #define FIXED_TO_DEG(f) ((double)(f) / ANGLE_SCALE) // 定点化正弦查表，0-90度，步长1度 static const int32_t sin_table[91] = { 0, 175, 349, 523, 698, 872, 1045, 1219, 1392, 1564, // ... 省略中间数据，实际使用时需要补全 10000 }; int32_t fixed_sin(int32_t angle_fixed) { // 归一化到0-360度 int32_t deg = (angle_fixed / ANGLE_SCALE) % 360; if (deg < 0) deg += 360; // 利用对称性查表 if (deg <= 90) return sin_table[deg]; if (deg <= 180) return sin_table[180 - deg]; if (deg <= 270) return -sin_table[deg - 180]; return -sin_table[360 - deg]; }

定点化能把计算时间压到float的1/3左右，但代码复杂度上升。我的经验是：如果MCU有FPU（如Cortex-M4F），直接用float，别折腾定点；如果没有FPU且计算频率低于1Hz，用float也能接受；只有高频更新且无FPU时才考虑定点化。

4.2 用已知数据验证：三个必须通过的测试用例

代码写完必须验证。我一般用三个测试用例：春分正午、夏至正午、冬至正午。春分时太阳直射赤道，北纬某地正午高度角等于90°减去纬度。夏至和冬至分别加减23.44°。

// 测试用例：北纬39.9度，东经116.4度，东八区 void run_test_cases(void) { SunInput in = {0}; SunPosition out = {0}; // 测试1：春分正午 in.year = 2024; in.month = 3; in.day = 20; in.hour = 12; in.minute = 0; in.second = 0; in.latitude = 39.9; in.longitude = 116.4; in.timezone = 8.0; calc_sun_position(&in, &out); // 预期高度角约50.1度（90 - 39.9） printf("春分正午: 高度角=%.2f, 方位角=%.2f\n", out.altitude, out.azimuth); // 测试2：夏至正午 in.month = 6; in.day = 21; calc_sun_position(&in, &out); // 预期高度角约73.5度（90 - 39.9 + 23.44） printf("夏至正午: 高度角=%.2f, 方位角=%.2f\n", out.altitude, out.azimuth); // 测试3：冬至正午 in.month = 12; in.day = 21; calc_sun_position(&in, &out); // 预期高度角约26.7度（90 - 39.9 - 23.44） printf("冬至正午: 高度角=%.2f, 方位角=%.2f\n", out.altitude, out.azimuth); }

这三个用例能覆盖赤纬角的正负和零值。如果春分正午高度角偏差超过0.5°，说明时差或者经度修正有问题。如果夏至冬至偏差大，说明赤纬角计算有问题。方位角正午应该接近180°（北半球太阳在正南），如果差很多，检查时角符号和方位角归一化。

4.3 与在线天文台数据对比的自动化脚本

手工验证只能看几个点，要全面验证得跟权威数据对比。我一般写个Python脚本，调用C代码编译成的动态库，批量对比。

import ctypes import datetime import ephem # 天文计算库，用于生成参考数据 # 加载C动态库 lib = ctypes.CDLL('./libsunpos.so') class SunInput(ctypes.Structure): _fields_ = [ ('year', ctypes.c_int), ('month', ctypes.c_int), ('day', ctypes.c_int), ('hour', ctypes.c_int), ('minute', ctypes.c_int), ('second', ctypes.c_double), ('latitude', ctypes.c_double), ('longitude', ctypes.c_double), ('timezone', ctypes.c_double) ] class SunPosition(ctypes.Structure): _fields_ = [ ('altitude', ctypes.c_double), ('azimuth', ctypes.c_double), ('declination', ctypes.c_double), ('hour_angle', ctypes.c_double) ] def c_calc(year, month, day, hour, minute, lat, lon, tz): inp = SunInput(year, month, day, hour, minute, 0.0, lat, lon, tz) out = SunPosition() lib.calc_sun_position(ctypes.byref(inp), ctypes.byref(out)) return out.altitude, out.azimuth def ref_calc(year, month, day, hour, minute, lat, lon, tz): obs = ephem.Observer() obs.lat = str(lat) obs.lon = str(lon) obs.date = datetime.datetime(year, month, day, hour, minute) sun = ephem.Sun() sun.compute(obs) alt = float(sun.alt) * 180.0 / 3.141592653589793 az = float(sun.az) * 180.0 / 3.141592653589793 return alt, az # 批量对比 max_alt_err = 0.0 max_az_err = 0.0 for day in range(1, 366, 10): for hour in [6, 9, 12, 15, 18]: d = datetime.datetime(2024, 1, 1) + datetime.timedelta(days=day-1) c_alt, c_az = c_calc(d.year, d.month, d.day, hour, 0, 39.9, 116.4, 8.0) r_alt, r_az = ref_calc(d.year, d.month, d.day, hour, 0, 39.9, 116.4, 8.0) alt_err = abs(c_alt - r_alt) az_err = abs(c_az - r_az) if alt_err > max_alt_err: max_alt_err = alt_err if az_err > max_az_err: max_az_err = az_err print(f"最大高度角误差: {max_alt_err:.4f}度") print(f"最大方位角误差: {max_az_err:.4f}度")

这个脚本能跑出全年不同时段的误差分布。NOAA简化算法的典型误差在0.1°以内，如果超过0.5°，说明某个中间量算错了。注意ephem库本身也有误差，对比时要把它的精度考虑进去。

5. 避坑与排查：太阳位置计算中最容易翻车的五个点

5.1 时区没转UTC导致方位角整体偏移

现象：算出来的方位角在正午时不是180°而是偏几十度，而且偏差随时间规律变化。

原因：输入的是本地时间，但公式里直接当UTC用了。东八区差8小时，地球每小时转15°，8小时就是120°。

解决：在计算儒略日之前先减时区偏移。如果跨日了要调整日期。我一般先把本地时间转成Unix时间戳，再从时间戳算儒略日，这样跨日跨月跨年都不用自己处理。

5.2 经度符号搞反导致时差修正方向错误

现象：高度角基本正确，但方位角在上午偏西、下午偏东，跟实际相反。

原因：经度符号约定不统一。有的代码用西经为正，有的用东经为正。真太阳时公式里4.0 * longitude，如果经度符号反了，时差修正方向就反了。

解决：统一约定东经为正、西经为负。在函数入口做检查，如果经度绝对值超过180就报错。测试时用东经和西经各测一次，看正午方位角是否都接近180°。

5.3 浮点精度不足导致小角度计算出现跳变

现象：太阳接近天顶时方位角突然跳变，或者高度角在0°附近抖动。

原因：用了float而不是double，或者acos/asin的参数因为浮点误差略微超过±1.0，导致返回NaN。

解决：关键计算用double。在acos/asin之前做数值裁剪，把参数限制在[-1.0, 1.0]。方位角计算用atan2代替acos，稳定性更好。

5.4 闰年和月末日期处理遗漏

现象：每年2月底到3月初算出来的结果偏差突然变大，过了3月又恢复正常。

原因：儒略日计算里对闰年和月末的处理有遗漏。2月29日、12月31日这些边界日期容易出错。

解决：用Unix时间戳做中转，时间戳库会自动处理闰年和月末。如果必须自己算，把每个月的天数做成表，闰年2月特殊处理。

5.5 大气折射修正过度导致夜间计算出负高度角

现象：太阳已经落山了，但修正后的高度角还是正的，追光系统还在追。

原因：折射修正公式在高度角低于-1°时应该返回0，但有些实现没做这个判断，导致夜间算出虚假的太阳位置。

解决：在折射修正函数入口加判断，高度角低于-1°直接返回0。追光系统还应该加一个高度角阈值，低于5°就停止追踪，避免无效动作。

6. 进阶技巧：用查表加插值把计算速度再提一个量级

前面讲的都是实时计算。如果MCU性能实在不够，或者需要每秒计算上百次，可以用查表加插值。核心思路是：赤纬角和时差的变化很慢，一天之内变化不到0.5°，可以按天做成表，运行时线性插值。

// 赤纬角和时差查找表，按年积日存储 typedef struct { int16_t declination; // 赤纬角，放大100倍，单位0.01度 int16_t eot; // 时差，放大100倍，单位0.01分钟 } SunTableEntry; // 生成一年的查找表（闰年366天，平年365天） void generate_sun_table(int year, SunTableEntry *table, int *days) { int is_leap = (year % 4 == 0 && year % 100 != 0) || (year % 400 == 0); *days = is_leap ? 366 : 365; for (int d = 0; d < *days; d++) { // 计算该年积日对应的儒略日 double jd = local_to_julian_day(year, 1, 1, 12, 0, 0.0, 0.0) + d; double dec, eot; calc_declination_eot(jd, &dec, &eot); table[d].declination = (int16_t)(dec * 100.0); table[d].eot = (int16_t)(eot * 100.0); } } // 运行时插值 void interpolate_sun_table(const SunTableEntry *table, int days, double day_of_year, double *declination, double *eot) { int idx = (int)day_of_year; if (idx >= days - 1) idx = days - 2; double frac = day_of_year - idx; *declination = (table[idx].declination * (1.0 - frac) + table[idx + 1].declination * frac) / 100.0; *eot = (table[idx].eot * (1.0 - frac) + table[idx + 1].eot * frac) / 100.0; }

这个方案把赤纬角和时差的计算从每次几十次三角函数降到了两次乘加。表的大小是365×4字节≈1.5KB，对大多数MCU都能接受。精度损失约0.01°，对追光完全够用。

插值时注意边界处理：day_of_year接近年末时idx会越界，要限制在days-2以内。frac是小数部分，用于线性插值。

我一般会在PC上生成好表，导出成C数组直接编译进固件，而不是在MCU上实时生成。这样启动快，也不占运行时间。生成表的脚本用Python写，输出格式化成C数组，直接贴到工程里。

最后一个习惯：每次改完代码，我都会用第4章的自动化脚本跑一遍全年对比。有一次改时区处理逻辑，以为只是小改动，结果脚本跑出来最大方位角误差2.3°，查了半天发现是跨日判断里把>=写成了>，导致整点时刻的日期调整漏了。这种错误肉眼很难发现，但脚本一跑就露馅。希望帮到你。

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