简介这份资源提供一套基于NOAA算法的太阳方位角与高度角计算C代码面向从事太阳能追日、测地定位、气象观测等需要精确太阳位置解算的开发者与工程技术人员。代码输入经纬度、日期、时间与时区即可输出太阳方位角与高度角方位角以正北为零度参考南北东西半球统一适用计算精度较高。压缩包内共1个docx文件约20KB内容以C源码为主涵盖角度与弧度转换、积日计算、日期转儒略日、儒略世纪数换算、太阳几何平黄经等核心函数并附有NOAA官方算法参考链接便于对照验证。目前已有1973人学习下载适合需要将太阳位置计算模块集成到嵌入式或上位机项目中的读者参考可直接移植或按需裁剪快速完成追日控制、光照分析等功能的开发。1. 太阳方位角高度角计算C代码从天文算法到嵌入式落地的完整路径做光伏追光、无人机航迹规划或者建筑日照分析时绕不开一个基础问题给定时间、经纬度太阳此刻在天上的哪个位置太阳方位角高度角计算C代码就是解决这件事的。它把天文历法里的太阳位置算法用C语言实现输入年月日时分和经纬度输出高度角仰角和方位角。听起来简单但真正落地时会遇到一堆问题时间用世界时还是本地时大气折射要不要修正精度要求到角秒还是度嵌入式平台浮点性能够不够我做过几个光伏追踪和户外设备姿态校准的项目这套算法从PC验证到MCU部署踩了不少坑下面把完整路径拆开讲清楚包括算法选型、代码实现、参数配置和常见翻车点。2. 太阳位置算法选型为什么SPA和NOAA简化算法各有适用场景2.1 三种主流算法的精度与开销对比太阳位置计算不是只有一种方法。从最简单的低精度近似到天文台级别的SPASolar Position Algorithm精度和计算量差了好几个数量级。选错了要么精度不够导致追光偏几度要么在MCU上跑不动。算法典型精度计算量适用场景简易正弦近似±1°~3°极低对精度无要求的粗略判断NOAA简化算法±0.1°~0.5°低光伏追光、日照分析SPA完整算法±0.0003°高天文观测、精密定日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/(altitude5.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(daysday-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°查了半天发现是跨日判断里把写成了导致整点时刻的日期调整漏了。这种错误肉眼很难发现但脚本一跑就露馅。希望帮到你。本文还有配套的精品资源点击获取