【C语言】计算实时太阳角度（高度角、方位角），以及使用stm32单片机实时获取时间戳

整体计算方法

在编写该代码的过程中寻找了多篇博文和论文，综合所有文章且按网上的以0时的方位角的0°，且随时间累加累加至360度。我修改了博文和论文的一些角度的计算方法。得到一下代码与网站计算的方位角相互验证过，误差不超过1°
验证网站

太阳赤纬角

在这里插入图片描述
n为该天在当年的序号
另一种方法没尝试成功：

代码如下读者可以自行尝试，我没尝试成功

	int n, n0;double t, theta, delta;// 计算n和n0n = month * 30 + day;n0 = 79.6764 + 0.2422 * (year - 1985) - ((year - 1985) / 4);// 计算tt = n ;// 计算θ（日角）theta = 2 * M_PI * t / 365.2422;// 计算δ（太阳赤纬）delta = 0.3723 + 23.2567 * sin(theta) + 0.1149 * sin(2 * theta) - 0.1712 * sin(3 * theta)- 0.758 * cos(theta) + 0.3656 * cos(2 * theta) + 0.0201 * cos(3 * theta);

太阳时角

在这里插入图片描述
时角这里的话中午时时角为0下午时角为负，可以简化为时角=15*(t-12)
不过这里计算的是时区的时间，不是当地时间所以存在误差，我怀疑代码计算的误差大部分都是出自这里
当地时角位于东7区的换算方法: (120度 – 当地经度) * 4，1个小时15度，平均下来，1度4分钟。
既时角=15*(t-12)+(120-经度)

太阳高度角

在这里插入图片描述
通过上述两个方程，即可求出高度角的各个参数。

太阳方位角

在这里插入图片描述
h为太阳高度角
Φ为纬度数据
σ为太阳的赤纬角

代码部分

C语言三角函数运算

因为角度计算涉及大量三角函数，所以C语言的三角函数计算不可避免需要使用。

C语言的三角函数计算函数

需要包含头文件#include <math.h>

  double __cdecl sin(double _X);double __cdecl cos(double _X);double __cdecl tan(double _X);double __cdecl asin(double _X);double __cdecl acos(double _X);double __cdecl atan(double _X);double __cdecl atan2(double _Y,double _X);

其中sin、cos、tan 函数的参数是一个用弧度表示的角度，分别返回这个角度的正弦、余弦、正切值。度和弧度都是衡量角的大小的单位。度用 ° 来表示，弧度用rad表示。

1度＝π/180 弧度 ( ≈0.017453弧度 )

asin()、acos、atan 函数分别返回它们的参数的反正弦、反余弦、反正切值。功能刚好和sin、cos、tan 函数相反。

赤纬角计算

通过公式：
在这里插入图片描述

/********************************************************************************************************************************************
*@*名称 ： Declination_Angle_Calculation
*@*功能 ： 计算太阳赤纬角
*@*形参 ： YearData年数据	DayData日数据
*@*返回值 ： DeclinationAngle太阳赤纬角 
*********************************************************************************************************************************************/
static double Declination_Angle_Calculation(double Daydata,short YearData)
{double DeclinationAngle;DeclinationAngle = 23.45 * sin((Daydata-80)/370.0 * 360*PI/180);		//赤纬角公式 return DeclinationAngle;												//返回赤纬角 
}

时角计算

通过公式：
时角=15*(t-12)

/********************************************************************************************************************************************
*@*名称 ： Hour_Angle_Calculation
*@*功能 ： 计算太阳时角
*@*形参 ： HourData小时数据 
*@*返回值 ： HourAngle太阳时角 
*********************************************************************************************************************************************/
static double Hour_Angle_Calculation(double HourData)
{double HourAngle;HourAngle = 15*(HourData-12);return HourAngle;}

太阳高度角计算

通过公式
在这里插入图片描述

/********************************************************************************************************************************************
*@*名称 ： Height_Angle_Calculation
*@*功能 ： 计算太阳高度角
*@*形参 ： LatitudeData纬度数据		HourData小时数据 	DeclinationAngle赤纬角数据 
*@*返回值 ： ElevatingAngle太阳高度角 
*********************************************************************************************************************************************/
static double Height_Angle_Calculation(double LatitudeData,double HourAngle,double DeclinationAngle)
{double ElevatingAngle;													//高度角double ElevatingAngle_sin;												//高度角正弦值double DeclinationAngle_sin = sin(DeclinationAngle * PI / 180);			//赤纬角的正弦值 double DeclinationAngle_cos = cos(DeclinationAngle * PI / 180);			//赤纬角余弦值 double LatitudeData_sin = sin(LatitudeData * PI / 180);					//纬度正弦值 double LatitudeData_cos = cos(LatitudeData * PI / 180);					//纬度余弦值 double HourAngle_cos = cos(HourAngle * PI/180);ElevatingAngle_sin = (DeclinationAngle_sin*LatitudeData_sin) + (DeclinationAngle_cos*LatitudeData_cos*HourAngle_cos);ElevatingAngle = asin(ElevatingAngle_sin)*180/PI;return ElevatingAngle;
}

太阳方位角计算

通过公式
在这里插入图片描述

/********************************************************************************************************************************************
*@*名称 ： Azimuth_Angle_Calculation
*@*功能 ： 计算太阳方位角
*@*形参 ： ElevatingAngle太阳高度角		HourAngle太阳时角	DeclinationAngle赤纬角数据	HourData小时数据	LatitudeData纬度数据
*@*返回值 ： 太阳方位角（0点开始360度） 
*********************************************************************************************************************************************/
static double  Azimuth_Angle_Calculation(double ElevatingAngle,double HourAngle,double DeclinationAngle,double HourData,double LatitudeData)
{double AzimuthAngle,AzimuthAngle_cos;double ElevatingAngle_cos = cos(ElevatingAngle * PI / 180);				//高度角余弦值 double ElevatingAngle_sin = sin(ElevatingAngle * PI / 180);				//高度角正弦值 double LatitudeData_sin = sin(LatitudeData * PI / 180);					//纬度正弦值 double LatitudeData_cos = cos(LatitudeData * PI / 180);					//纬度余弦值  	double HourAngle_sin = sin(HourAngle * PI / 180); 						//时角正弦值 double DeclinationAngle_cos = cos(DeclinationAngle * PI / 180);			//赤纬角余弦值 double DeclinationAngle_sin = sin(DeclinationAngle * PI / 180);			//赤纬角正弦值AzimuthAngle_cos = (ElevatingAngle_sin*LatitudeData_sin-DeclinationAngle_sin)/(ElevatingAngle_cos*LatitudeData_cos);		//初步计算太阳方位角弧度值的cos值 AzimuthAngle = acos(AzimuthAngle_cos)*180/PI;																				//计算方位角的角度 if (HourData > 12) 														//不判断的话12点为180度最高，然后两边的角度都是180度递减，所以需要一个360度的角度反馈需要将方位角判断一下。 {return AzimuthAngle + 180;}else {return -AzimuthAngle + 180;}}

其中函数尾部的返回时180度加减计算得出的方位角，是为了对齐网页计算得出的方位角，以0时的方位角的0°，且随时间累加累加至360度。

角度计算到此就结束了

时间戳转换

下文内容多数引用自：博客

因为我打算用stm32f103mcu进行太阳方位角以及高度角的计算，所以需要用stm32单片机获取时间戳信息

我采用的方法是使用esp8266-01s模块来通过联网获取时间戳的方式得到时间。

配置 ESP8266 入网

AT+CWMODE=1                     // station模式，设置成3也行
AT+RST                          // 重启模块
AT+CWJAP="ssid","password"      // 连接AP

连接授时服务器

AT+CIPMUX=0                                // 设置WiFi模块为单路链接
AT+CIPSTART="UDP","1.cn.pool.ntp.org",123  
// 链接到 NTP 服务器，时间服务器默认的端口号是 123，协议为UDP.AT+CIPSEND=48                              
// 发送48byte的数据若ESP8266返回'>'则继续注意，严格遵照以下数据且使用16进制发送
//切勿发送新行E3 0 6 EC 0 0 0 0 0 0 0 0 31 4E 31 34 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
//用16进制接收数据
E3 0 6 EC 0 0 0 0 0 0 0 0 31 4E 31 34 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

不知道为什么我好像需要发送两遍下边的16进制数据，然后接收两帧数据

接收数据

执行完上面两步，此时如果顺利的话会受到一大堆16进制数据
诸如：41 54 2B 43 49 50 53 45 4E 44 3D 34 38 0D 0D 0A 0D 0A 4F 4B 0D 0A 3E 20 0D 0A 52 65 63 76 20 34 38 20 62 79 74 65 73 0D 0A 0D 0A 53 45 4E 44 20 4F 4B 0D 0A 0D 0A 2B 49 50 44 2C 34 38 3A 24 02 06 E9 00 00 00 52 00 00 08 F2 0A 45 02 22 E0 4F 57 31 79 AF 20 82 00 00 00 00 00 00 00 00 E0 4F 59 9C 2F 4A 40 07 E0 4F 59 9C 2F 4C 22 B9
此时需要找到数据中的关键帧，也就是服务器返回的48位时间数据，它们通常具有以下格式：（如果找不到的话请重试2中的3、4步）
- ① 它前面的数据为"2B 49 50 44 2C 34 38 3A"(+IPD,48:）,他后面没有东西;
- ② 它通常以 24 开头.3.对数据提取后编号，其中[40][41][42][43]为我们需要的时间数据，在本例中为"E0 4F 5C 00".

处理数据

①将”E04F5C00“转为十进制为3763297280，此即为1900距今的秒数。
②把这个时间减去2208988800（70年），得到标准的Unix时间戳。[3763297280-2208988800=1554308480]
③[1554308480/31556736+1970=2019]，得到当下年份。
④计算月份和日期需要考虑平年闰年和每月天数，太麻烦我也没搞懂。
⑤[(1554308480/3600)%24+8=24]，得到北京时间-时(别忘了满24归零)
⑥[(1554308480/60)%60=21]，得到北京时间-分
⑦[1554308480%60=20]，得到北京时间-秒，既当下时间是 24：21：20
⑧将时间戳(1554308480)输入时间戳生成网页与计算结果进行比较

时间戳解码并加入代码

因为我们计算太阳角度的代码，只需要获取年份，当天在一年中的序号，以及小时分钟信息。
所以有以下计算

/********************************************************************************************************************************************
*@*名称 ： Timestamp_Calculation
*@*功能 ： 计算时间戳
*@*形参 ： Timestamp时间戳数据 	Year_data年数据存放的地址	DayData日数据存放的地址		 Hour_data小时数据存放的地址 
*@*返回值 ： 无
*********************************************************************************************************************************************/
static void  Timestamp_Calculation(unsigned int Timestamp,short *Year_data,short *DayData,double *Hour_data)
{Timestamp  = Timestamp-0x83AA7E80;						//将网络获取的时间戳减去70年得到标准时间戳 *Year_data = Timestamp/31556736+1970;					//用时间戳获取年份 *DayData = (short)((Timestamp%31556736)/86400)+1;		//获取日期 *Hour_data = (Timestamp/3600)%24+8;						//获取小时数 if(*Hour_data==24){*Hour_data = 0;										//当时间为24点时归零 } *Hour_data += ((Timestamp/60)%60)/60.0;					//获取分钟信息并且转为小数形式加入小时数中 
}