太阳高度角与方位角计算程序-stm32G070(上)

经验体会: python验证算法,然后“翻译”  C语言

变量名长些好,尽可能表达物理含义。

参考高手的程序

参考文献:

太阳时的计算_邱国全

太阳高度角和方位角的计算算法_王慧

sunPosition.h

复制代码
/*
#纬度   latitude   北半球正          -90  ~  90度
#经度   altitude   东半球正          -180 ~ 180度
#太阳真时    apparent_solar_time
#太阳时角    solar_hour_angle        -180~180度
#计算赤纬角  solar_declination_angle -23.45 ~ 23.45度
#计算高度角  solar_altitude_angle    -90~90度
#太阳方位角  solar_azimuth_angle     -180~ 180度
#太阳天顶角  solar_zenith_angle       0 -180度
#太阳升起时角solar_hour_angle_sunrise  -180~0度
*/
#ifndef __solar_position_algorithm_header
#define __solar_position_algorithm_header
#include "stm32g0xx_hal.h"
typedef struct
{
    //----------------------INPUT VALUES------------------------

    double  year;           //输入: 4-digit year,    valid range: -2000 to 6000, error code: 1
    double  month;          //输入: 2-digit month,         valid range: 1 to 12, error code: 2
    double day;             //输入: 2-digit day,           valid range: 1 to 31, error code: 3
    double hour;            //输入: Observer local hour,   valid range: 0 to 24, error code: 4
    double minute;          //输入: Observer local minute, valid range: 0 to 59, error code: 5
    double second;          //输入: Observer local second, valid range: 0 to 59, error code: 6
    double timezone;     //输入: Observer time zone (negative west of Greenwich)
                         //       valid range: -18   to   18 hours,   error code: 7

    double longitude;    // 输入:经度Observer longitude (negative west of Greenwich)
                         //           valid range: -180  to  180 degrees, error code: 8

    double latitude;     // 输入:纬度Observer latitude (negative south of equator)
                         //           valid range: -90   to   90 degrees, error code: 9
    //---------------------OUTPUT VALUES------------------------
      double  ndays;                   // 从元旦(含)计,全年第几天
      double  mean_solar_time;         //输出:平太阳时
      double  error_time;              //时差: 真太阳时与平太阳时之差
    double     apparent_solar_time;     //输出:太阳真时    0-24
      double  solar_hour_angle;        //输出:太阳时角    -180~180度
      double  solar_declination_angle; //输出:赤纬角      -23.45 ~ 23.45度
      double  solar_altitude_angle;    //输出:高度角      -90~90度
        double  solar_azimuth_angle;     //输出:太阳方位角  -180~ 180度
        double  solar_hour_angle_sunrise;//输出:太阳升起时角-180~0度
} spa_data;

//Calculate SPA output values (in structure) based on input values passed in structure
uint8_t sunpos_cal(spa_data *spa);
#endif
复制代码

sunPosition.c

复制代码
#include <math.h>
#include <sunPosition.h>

#define PI         3.1415926535897932384626433832795028841971

// 计算 平太阳时
// INPUT:  
//         Beijing_time       北京标准时间,    单位:分钟
//                 longitude_beijing  北京标准时间经度,单位:度
//                 longitude_location 当地时间经度,东半球正,单位:度
//                 latitude_locaiton  当地时间维度,北半球正,单位:度
// OUTPUT: 时差,单位:分钟
double mean_solar_time_cal(double Beijing_time,float longitude_beijing,float  longitude_location,float latitude_locaiton)
{/*计算设备所在地的平太阳时*/
    double mean_solar_time;
    if(latitude_locaiton>=0)
    {
     mean_solar_time = Beijing_time - 4*(longitude_beijing-longitude_location);
    }
    else
    {
     mean_solar_time = Beijing_time + 4*(longitude_beijing-longitude_location);
    }
    return mean_solar_time;
}

// 计算 平太阳时-真太阳时的差
// INPUT: 当年第几天,从元旦计,1~365
// OUTPUT: 时差,单位:分钟
double time_error(double nday)
{/*计算设备所在地的平太阳时*/
    double gama_big;
    double ET;
    gama_big = 2*PI*(nday-1)/365;
    ET = 0.01719+0.428145*cos(gama_big);
    ET=ET-7.3520484*sin(gama_big);
    ET=ET-3.349758*cos(gama_big*2);
    ET=ET-9.371988*sin(gama_big*2);
    return ET;
}
// FUN:    计算真太阳时 单位:分钟
// INPUT: 
//         mean_solar_time 单位:分钟
//         error_time      单位:分钟
// OUTPUT: apparent_solar_time 单位:分钟
double apparent_solar_time_cal(const double mean_solar_time,const double error_time)
{/*计算设备所在地的真太阳时*/
    double apparent_solar_time;
    apparent_solar_time = mean_solar_time + error_time;
    return apparent_solar_time;
}

// FUN:    计算真太阳时 单位:分钟
// INPUT: 
//         mean_solar_time 单位:分钟
//         error_time      单位:分钟
// OUTPUT: apparent_solar_time 单位:分钟

double solar_hour_angle_cal(double apparent_solar_time)
{/*计算设备所在地的时间角度*/
    double temp;
    temp = (apparent_solar_time-12*60)*15.00/60;
    return temp;
}

// FUN:   计算当天赤纬角度
// INPUT: nday  当年第几天,从元旦计,1~365
// OUTPUT: solar_declination_angle,单位:度
double solar_declination_angle_cal(double nday)
{/*计算当天赤纬角度*/
    double solar_declination_angle;
    solar_declination_angle = 23.45 *sin(2*PI*(284+nday)/365);
    return solar_declination_angle;
}

// FUN:   计算当天日升高度角度
// INPUT: 
//       latitude                纬度     度
//       solar_declination_angle 赤纬角   度
// OUTPUT: solar_declination_angle,单位:度
double solar_hour_angle_sunrise_cal(double latitude, double solar_declination_angle)
{/*计算计算当天日升高度角度*/
    double solar_hour_angle_sunrise;
    /*计算余弦值*/
    solar_hour_angle_sunrise = sin(PI*latitude/180)*sin(PI*solar_declination_angle/180);
    solar_hour_angle_sunrise = solar_hour_angle_sunrise/( cos(PI*latitude/180)*cos(PI*solar_declination_angle/180) );
    if(solar_hour_angle_sunrise>1)//#截断误差处理
            solar_hour_angle_sunrise=1;
    if(solar_hour_angle_sunrise<-1)//#截断误差处理
            solar_hour_angle_sunrise=-1;    
    /*计算角度值*/    
    solar_hour_angle_sunrise = (-acos(-solar_hour_angle_sunrise)*180/PI);
    return solar_hour_angle_sunrise;
}

// FUN:   计算太阳高度角度
// INPUT: 
//        latitude                 当地纬度,   单位:度
//        solar_declination_angle  当天赤纬角, 单位:度
//        solar_hour_angle         当地太阳时角, 单位:度
// OUTPUT: solar_declination_angle,单位:度
double solar_altitude_angle_cal(float latitude,float solar_declination_angle,float solar_hour_angle)
{/*计算太阳高度角度*/
    float item1,item2;
    float total;
    float solar_altitude_angle;
    item1=sin(PI*latitude/180)*sin(PI*solar_declination_angle/180);
  item2=cos(PI*latitude/180)*cos(PI*solar_declination_angle/180)*cos(PI*solar_hour_angle/180);    
    total = item1+item2;
    if(total<-1)//#截断误差处理
            total=-1;
    if(total>1)//#截断误差处理
            total=1;
    solar_altitude_angle=asin(total)*180/PI;
    return solar_altitude_angle;
}

// FUN:   计算太阳方位角度
// INPUT: 
//        latitude                 当地纬度,   单位:度
//        solar_declination_angle  当天赤纬角, 单位:度
//        solar_hour_angle         当地太阳时角,   单位:度
//    solar_altitude_angle     当地太阳高度角, 单位:度
// OUTPUT: solar_declination_angle,单位:度
double solar_azimuth_angle_cal(float latitude,float solar_declination_angle,float solar_hour_angle,float solar_altitude_angle)
{/*计算太阳方位角度*/
  float sin_azimuth;
    float cos_azimuth;
    float azimuth;
    sin_azimuth = cos(solar_declination_angle*PI/180);
  sin_azimuth = sin_azimuth*sin(solar_hour_angle*PI/180);
    sin_azimuth = sin_azimuth/cos(solar_altitude_angle*PI/180);
    
    cos_azimuth =  sin(solar_altitude_angle*PI/180);
  cos_azimuth =  cos_azimuth*sin(latitude*PI/180);
  cos_azimuth =  cos_azimuth-sin(solar_declination_angle*PI/180);    
    cos_azimuth =  cos_azimuth /cos(solar_altitude_angle*PI/180);
    cos_azimuth =  cos_azimuth/cos(latitude*PI/180);

    if(cos_azimuth>1)
            cos_azimuth=1;
    if(cos_azimuth<-1)
            cos_azimuth=-1;
    if(sin_azimuth<=0)
            azimuth=-acos(cos_azimuth )*180/PI;
    if(sin_azimuth>0)
            azimuth= acos(cos_azimuth )*180/PI;     
    return azimuth;
}
// 输入合法性检查
uint8_t validate_inputs(spa_data *spa)
{
    if ((spa->year        < -2000) || (spa->year        > 6000)) return 1;
    if ((spa->month       < 1    ) || (spa->month       > 12  )) return 2;
    if ((spa->day         < 1    ) || (spa->day         > 31  )) return 3;
    if ((spa->hour        < 0    ) || (spa->hour        > 24  )) return 4;
    if ((spa->minute      < 0    ) || (spa->minute      > 59  )) return 5;
    if ((spa->second      < 0    ) || (spa->second      > 59  )) return 6;
    if ((spa->hour        == 24  ) && (spa->minute      > 0   )) return 5;
    if ((spa->hour        == 24  ) && (spa->second      > 0   )) return 6;

    if (fabs(spa->timezone)      > 18      ) return 8;
    if (fabs(spa->longitude)     > 180     ) return 9;
    if (fabs(spa->latitude)      > 90      ) return 10;
    
    return 0;
}
// 计算Sun Position
uint8_t sunpos_cal(spa_data *spa)
{
     const uint8_t month_days[]={31,28,31,30,31,30,31,30,31,30,31,30,31};
     
     uint8_t ret;     
   uint16_t ndays;
     float    totalminutesBeijing;
     uint16_t index;
     ret =0;
     ret= validate_inputs(spa);
     if(ret!=0) return ret;
     /*计算从元旦计,天数*/
     ndays =0;
   for(index=1;index< spa->month;index++)
     {
       ndays= ndays + month_days[index-1];
     }
     spa->ndays = ndays + spa->day;
     /*计算北京时间,分钟*/     
     totalminutesBeijing = spa->hour*60 +spa->minute +spa->second/60;
     /*计算当地当天平太阳时*/         
     spa->mean_solar_time = mean_solar_time_cal(totalminutesBeijing,120,spa->longitude,spa->latitude);
     /*计算当天时差(真太阳时-平太阳时)*/
   spa->error_time =     time_error(spa->ndays); 
     /*计算当天真太阳时*/    
   spa->apparent_solar_time =apparent_solar_time_cal(spa->mean_solar_time,spa->error_time);
     /*计算当天太阳时角*/
   spa->solar_hour_angle =     solar_hour_angle_cal(spa->apparent_solar_time);
     /*计算当天赤纬角*/    
     spa->solar_declination_angle=solar_declination_angle_cal(spa->ndays);
     /*计算当时高度角*/         
   spa->solar_altitude_angle = solar_altitude_angle_cal(spa->latitude,spa->solar_declination_angle,spa->solar_hour_angle); 
     /*计算当时方位角*/        
   spa->solar_azimuth_angle =solar_azimuth_angle_cal(spa->latitude,spa->solar_declination_angle,spa->solar_hour_angle,spa->solar_altitude_angle);
     /*计算当天日升角*/     
     spa->solar_hour_angle_sunrise =  solar_hour_angle_sunrise_cal(spa->latitude, spa->solar_declination_angle);
     
     return 0;
}
复制代码

使用距离

复制代码
#include <sunPosition.h>

    spa.year     = 2023;
    spa.month    = 3;
    spa.day    = 3;
    spa.hour     =12;
    spa.minute   =0;
    spa.second   =0;
    spa.latitude =60;
    spa.longitude=110;
    spa.timezone =8;
    
        
    //rt_enter_critical();    
    result = sunpos_cal(&spa);
    /* 退出临界段 */
    //rt_exit_critical();
复制代码

Python程序与测试

复制代码
import math
#计算时角   solar_hour_angle
#太阳真时   apparent_solar_time
#计算高度角  solar_altitude_angle degree
#太阳方位角  solar_azimuth_angle  degree
#太阳高度角  solar_zenith_angle   degree
#纬度   latitude  degree
#经度   altitude  degree
#平太阳时 mean solar time
#北京时间 Beijing time

month=(31,28,31,30,31,30,31,30,31,30,31,30,31)

def mean_solar_time_cal(Beijing_time,longitude_beijing,longitude_location,latitude_locaiton):
    if(latitude_locaiton>=0):
        temp = Beijing_time - 4*(longitude_beijing-longitude_location)
        return temp# 单位:分钟
    else:
        temp = Beijing_time + 4*(longitude_beijing-alongitude_location)
        return temp# 单位:分钟

def time_error(nday):# 可用
    gama_big = 2*math.pi*(nday-1)/365
    ET=0.01719+0.428145*math.cos(gama_big)
    ET=ET-7.3520484*math.sin(gama_big)
    ET=ET-3.349758*math.cos(gama_big*2)
    ET=ET-9.371988*math.sin(gama_big*2)
    return ET

def solar_hour_angle_cal(apparent_solar_time):#时度 -180 180
    return((apparent_solar_time-12*60)*15/60)

#计算赤纬角  solar_declination_angle
# nday
def solar_declination_angle_cal(nday):
    return(23.45 *math.sin(2*math.pi*(284+nday)/365))


def sin_solar_altitude_angle(latitude,solar_declination_angle,solar_hour_angle):
    item1=math.sin(math.pi*latitude/180)*math.sin(math.pi*solar_declination_angle/180)
    item2=math.cos(math.pi*latitude/180)*math.cos(math.pi*solar_declination_angle/180)*math.cos(math.pi*solar_hour_angle/180)
    mysum=item1+item2
    if(mysum<-1):
        mysum=-1
    if(mysum>1):
        mysum=1
    return(mysum)

def solar_altitude_angle_cal(latitude,solar_declination_angle,solar_hour_angle):
    item1=math.sin(math.pi*latitude/180)*math.sin(math.pi*solar_declination_angle/180)
    item2=math.cos(math.pi*latitude/180)*math.cos(math.pi*solar_declination_angle/180)*math.cos(math.pi*solar_hour_angle/180)
    mysum=item1+item2
    if(mysum<-1):#截断误差处理
        mysum=-1
    if(mysum>1):#截断误差处理
        mysum=1
    return(math.asin(mysum)*180/math.pi)

def solar_hour_angle_sunrise_cal(latitude,solar_declination_angle):
    temp = math.sin(math.pi*latitude/180)*math.sin(math.pi*solar_declination_angle/180)
    temp = temp/math.cos(math.pi*latitude/180)/math.cos(math.pi*solar_declination_angle/180)
    if(temp>1): #截断误差处理
        temp=1
    if(temp<-1):#截断误差处理
        temp=-1
    return(-math.acos(-temp)*180/math.pi)


def solar_azimuth_angle_cal(latitude,solar_declination_angle,solar_hour_angle,solar_altitude_angle):
    item1=math.cos(solar_declination_angle*math.pi/180)
    item1=item1*math.sin(solar_hour_angle*math.pi/180)
    item1= item1/math.cos(solar_altitude_angle*math.pi/180)
    item2= math.sin(solar_altitude_angle*math.pi/180)
    item2= item2*math.sin(latitude*math.pi/180)
    item2=item2-math.sin(solar_declination_angle*math.pi/180)
    item2= item2/math.cos(solar_altitude_angle*math.pi/180)
    item2= item2/math.cos(latitude*math.pi/180)
    if(item2>1):
        item2=1
    if(item2<-1):
        item2=-1
    if(item1<=0):
        return(-math.acos(item2)*180/math.pi)
    if(item1>0):
        return(math.acos(item2)*180/math.pi)
    
if __name__ == '__main__':    
#测试数据
    ndays = 59  # 3月1日 
    Beijing_time = 12*60
    longitude_beijing =120
    longitude_location=116.33
    latitude_locaiton =39.933
    mean_solar_time = mean_solar_time_cal(Beijing_time,longitude_beijing,longitude_location,latitude_locaiton)
    print("mean_solar_time {:3.2f}分钟".format(mean_solar_time))
    apparent_solar_time = mean_solar_time +time_error(ndays)
    print("time_error {:3.2f}分钟".format(time_error(ndays)))
    print("apparent_solar_time {:3.2f}分钟".format(apparent_solar_time))

    solar_hour_angle=solar_hour_angle_cal(apparent_solar_time)
    print("solar_hour_angle {:3.2f}°".format(solar_hour_angle))

    solar_declination_angle = solar_declination_angle_cal(ndays)
    print("solar_declination_angle {:3.2f}°".format(solar_declination_angle))

    solar_altitude_angle = solar_altitude_angle_cal(latitude_locaiton,solar_declination_angle,solar_hour_angle)
    print("solar_altitude_angle {:3.2f}°".format(solar_altitude_angle))

    solar_azimuth_angle= solar_azimuth_angle_cal(latitude_locaiton,solar_declination_angle,solar_hour_angle,solar_altitude_angle)
    print("solar_azimuth_angle {:3.2f}°".format(solar_azimuth_angle))

    solar_hour_angle_sunrise=solar_hour_angle_sunrise_cal(latitude_locaiton,solar_declination_angle)
    print("solar_hour_angle_sunrise {:3.2f}°".format(solar_hour_angle_sunrise))
复制代码

 

算例结果:

mean_solar_time 705.32分钟
time_error -13.08分钟
apparent_solar_time 692.24分钟
solar_hour_angle -6.94°
solar_declination_angle -8.67°
solar_altitude_angle 40.97°
solar_azimuth_angle -9.10°
solar_hour_angle_sunrise -82.67°

ESP32 使用arduino代码与测试

头文件:sunrise.h

复制代码
/**
 ****************************************************************************************************
 * @file          sunrise.h
 * @author      辛和
 * @version     V1.0
 * @date        2024-6-24
 * @brief       公用变量
 * @license    辛和
 ****************************************************************************************************
 * @attention
 * MCU:ESP32S3
 *
 * 修改说明
 * V1.0 20240707
 *
 ****************************************************************************************************
 * 
 */

#ifndef __SUNRISE_H
#define __SUNRISE_H

#include "Arduino.h"

/* 变量类型定义*/
typedef struct
{
  //----------------------INPUT VALUES------------------------
  int16_t year;      //输入: 4-digit year,          valid range: 0 to 6000, error code: 1
  int16_t month;     //输入: 2-digit month,         valid range: 1 to 12, error code: 2
  int16_t day;       //输入: 2-digit day,           valid range: 1 to 31, error code: 3
  int16_t hour;      //输入: Observer local hour,   valid range: 0 to 24, error code: 4
  int16_t minute;    //输入: Observer local minute, valid range: 0 to 59, error code: 5
  int16_t second;    //输入: Observer local second, valid range: 0 to 59, error code: 6
  int16_t timezone;  //输入: 当地Observer time zone (negative west of Greenwich)
                     //       valid range: -18   to   18 hours,   error code: 7

  float longitude;  // 输入:当地经度Observer longitude (negative west of Greenwich)
                    //           valid range: -180  to  180 degrees, error code: 8

  float latitude;  // 输入:当地纬度Observer latitude (negative south of equator)
                   //           valid range: -90   to   90 degrees, error code: 9
                   //---------------------OUTPUT VALUES------------------------
  int16_t ndays;   //UTC时间,从元旦(含)计,全年第几天

  float mean_solar_time;           //输出:当地平太阳时, 由hhmmss计算,单位 分钟
  float error_time;                //时差:当地真太阳时与平太阳时之差,单位 分钟
  float apparent_solar_time;       //输出:当地太阳真时    0-24,单位 分钟
  float solar_hour_angle;          //输出:当地太阳时角    -180~180度
  float solar_declination_angle;   //输出:当地赤纬角      -23.45 ~ 23.45度
  float solar_altitude_angle;      //输出:当地太阳高度角  -90~90度
  float solar_azimuth_angle;       //输出:当地太阳方位角  -180~ 180度
  float solar_hour_angle_sunrise;  //输出:当地升起时刻的太阳时角-180~0度
  float sunrise_std_time;          //输出:当地当日日出标准时间    0-24,单位 分钟
} SUN_POS_t;

/* 变量定义 */
extern SUN_POS_t sun_pos_x;

/* 函数声明 */
uint8_t IsLeapYear(int16_t year);
uint32_t Date2Days(uint8_t year, uint8_t month, uint8_t day);
float time_error(double nday);
float solar_declination_angle_cal(float nday);
float solar_hour_angle_sunrise_cal(float latitude, float solar_declination_angle);

float local_mean_solar_time_cal(float std_minutes, int16_t timezone,float longitude_location) ;
float local_std_time_cal(float local_solar_mean_minutes, int16_t timezone,float longitude_location); /*计算设备所在地的平太阳时*/


double solar_altitude_angle_cal(float latitude,float solar_declination_angle,float solar_hour_angle);
double solar_azimuth_angle_cal(float latitude, float solar_declination_angle, float solar_hour_angle, float solar_altitude_angle);

uint8_t sunpos_cal(SUN_POS_t *spa);
#endif
复制代码

源文件 sunrise.c

复制代码
/**
 ****************************************************************************************************
 * @file        sunrise.c
 * @author      辛和
 * @version     V1.0
 * @date        2024-6-24
 * @brief       公用变量
 * @license    辛和
 ****************************************************************************************************
 * @attention
 * MCU:ESP32S3
 *
 * 修改说明
 * V1.0 20240624
 *
 ****************************************************************************************************
 * 
 */
#include "sunrise.h"
#include "Arduino.h"
#include <String.h>

/***********************SUNRISE******************************/
#ifndef PI
#define PI 3.1415926535897932384626433832795028841971
#endif
// 定义一个变量来存储太阳位置特性
SUN_POS_t sun_pos_x;

/*内部函数声明*/

/**
* @brief       判断是否为闰年
* @param       year 公历年,例如2024,输入24
* @retval      0x01 闰年,0x00 平年
*/
uint8_t IsLeapYear(int16_t year) {
  uint8_t Is_Leap_Year;
  Is_Leap_Year = 0x00;
  if ((year % 4 == 0) && (year % 100 != 0))
    Is_Leap_Year = 0x01;  //闰年
  else {
    if (year % 400 == 0) Is_Leap_Year = 0x01;  //闰年
  }
  return Is_Leap_Year;
}
/**
  * @brief  根据年月日,计算是全年第几天
  * @param  
  *         year  年,两位整数
  *         month 月,两位整数
  *         day   日,两位整数
  * @retval 
  *         当年第几天,从元旦计,1~365
  */
uint32_t Date2Days(uint8_t year, uint8_t month, uint8_t day) {
  const uint8_t Leap_Year_month_days[] = { 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
  const uint8_t Ordinary_Year_month_days[] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
  uint8_t Is_Leap_Year;
  uint32_t totalDays;
  uint8_t i;
  if ((year <= 0) || (year > 100)) return 0;
  if ((month <= 0) || (month > 12)) return 0;

  Is_Leap_Year = IsLeapYear(year);

  if (Is_Leap_Year == 0x01) {
    if ((day <= 0) || (day > Leap_Year_month_days[month - 1])) return 0;
  } else {
    if ((day <= 0) || (day > Ordinary_Year_month_days[month - 1])) return 0;
  }

  totalDays = 0;
  for (i = 0; i < month - 1; i++) {
    if (Is_Leap_Year == 0x01) {
      totalDays += Leap_Year_month_days[i];
    } else {
      totalDays += Ordinary_Year_month_days[i];
    }
  }
  totalDays += day;
  return totalDays;
}
/**
  * @brief  计算ET   平太阳时+ET=真太阳时
  * @param 
  *         nday       当年第几天,从元旦计,1~365
  * @retval 
  *         时差,单位:分钟
  */
float time_error(double nday) { /*计算当天的ET*/
  float gama_big;
  float ET;
  gama_big = 2 * PI * (nday - 1) / 365;
  ET = 0.01719f + 0.428145f * cos(gama_big);
  ET = ET - 7.3520484f * sin(gama_big);
  ET = ET - 3.349758f * cos(gama_big * 2);
  ET = ET - 9.371988f * sin(gama_big * 2);
  return ET;
}
/**
  * @brief  计算当天赤纬角度
  * @param 
  *         nday       当年第几天,从元旦计,1~365
  * @retval 
  *         赤纬角,单位:度
  */
float solar_declination_angle_cal(float nday) { /*计算当天赤纬角度*/
  double solar_declination_angle;
  solar_declination_angle = 23.45f * sin(2 * PI * (284 + nday) / 365);
  return solar_declination_angle;
}
/**
  * @brief  计算当天日升高度角度(直接按照4分钟1度换算为真太阳时)
  * @param 
  *         latitude   纬度,单位:度
  *         solar_declination_angle 赤纬角,单位:度
  * @retval 
  *         日出时角,单位:度
  */
float solar_hour_angle_sunrise_cal(float latitude, float solar_declination_angle) { /*计算计算当天日升高度角度*/
  float solar_hour_angle_sunrise;
  /*计算余弦值*/
  solar_hour_angle_sunrise = sin(PI * latitude / 180) * sin(PI * solar_declination_angle / 180);
  solar_hour_angle_sunrise = solar_hour_angle_sunrise / (cos(PI * latitude / 180) * cos(PI * solar_declination_angle / 180));
  if (solar_hour_angle_sunrise > 1)  //#截断误差处理
    solar_hour_angle_sunrise = 1;
  if (solar_hour_angle_sunrise < -1)  //#截断误差处理
    solar_hour_angle_sunrise = -1;
  /*计算角度值*/
  solar_hour_angle_sunrise = (-acos(-solar_hour_angle_sunrise) * 180 / PI);
  return solar_hour_angle_sunrise;
}

/**
  * @brief  直接由标准时间hhmmss计算当地平太阳时,单位分钟
  * @param 
  *         std_minutes        按照时区的标准时间,由当天的hhmmss计算得到 单位:分钟
  *         timezone           时区
  *         longitude_location 当地时间经度,东半球正,单位:度
  * @retval 
  *         当地平太阳时,单位:分钟
  */
float local_mean_solar_time_cal(float std_minutes, int16_t timezone,float longitude_location) { /*计算设备所在地的平太阳时*/
  float mean_solar_time;
  if (longitude_location >= 0) {
    mean_solar_time = std_minutes - 4 * (timezone*15.0f - longitude_location);
  } else {
    mean_solar_time = std_minutes + 4 * (timezone*15.0f  - longitude_location);
  }
  return mean_solar_time;
}
/**
  * @brief  直接由当地平太阳时计算标准时间hhmmss,单位分钟
  * @param 
  *         local_solar_mean_minutes  当地平太阳时,单位:分钟
  *         timezone                  时区
  *         longitude_location        当地时间经度,东半球正,单位:度
  * @retval 
  *         当地标准时间,单位:分钟
  */
float local_std_time_cal(float local_solar_mean_minutes, int16_t timezone,float longitude_location) { /*计算设备所在地的平太阳时*/
  float local_std_time;
  if (longitude_location >= 0) {
    local_std_time = local_solar_mean_minutes + 4 * (timezone*15.0f - longitude_location);
  } else {
    local_std_time = local_solar_mean_minutes - 4 * (timezone*15.0f  - longitude_location);
  }
  return local_std_time;
}
/**
  * @brief  计算太阳高度角度
  * @param 
  *         latitude                纬度 ,北半球正,单位:度
  *         solar_declination_angle 赤纬角,单位:度
  *         solar_hour_angle        当地真太阳时角, 单位:度
  * @retval 
  *         solar_altitude_angle,-90~90,单位:度 
  */
double solar_altitude_angle_cal(float latitude, float solar_declination_angle, float solar_hour_angle) { /*计算太阳高度角度*/
  float item1, item2;
  float total;
  float solar_altitude_angle;
  item1 = sin(PI * latitude / 180) * sin(PI * solar_declination_angle / 180);
  item2 = cos(PI * latitude / 180) * cos(PI * solar_declination_angle / 180) * cos(PI * solar_hour_angle / 180);
  total = item1 + item2;
  if (total < -1)  //#截断误差处理
    total = -1;
  if (total > 1)  //#截断误差处理
    total = 1;
  solar_altitude_angle = asin(total) * 180 / PI;
  return solar_altitude_angle;
}
/**
  * @brief  计算太阳方位角度(先计算高度角)
  * @param  
  *         latitude                纬度 ,北半球正,单位:度
  *         solar_declination_angle 赤纬角,单位:度
  *         solar_hour_angle        当地真太阳时角, 单位:度
  *         solar_altitude_angle    当地太阳高度角, 单位:度
  * @retval 
  *         solar_azimuth_angle,   太阳方位角-180~180,单位:度
  */
double solar_azimuth_angle_cal(float latitude, float solar_declination_angle, float solar_hour_angle, float solar_altitude_angle) { /*计算太阳方位角度*/
  float sin_azimuth;
  float cos_azimuth;
  float azimuth;
  sin_azimuth = cos(solar_declination_angle * PI / 180);
  sin_azimuth = sin_azimuth * sin(solar_hour_angle * PI / 180);
  sin_azimuth = sin_azimuth / cos(solar_altitude_angle * PI / 180);

  cos_azimuth = sin(solar_altitude_angle * PI / 180);
  cos_azimuth = cos_azimuth * sin(latitude * PI / 180);
  cos_azimuth = cos_azimuth - sin(solar_declination_angle * PI / 180);
  cos_azimuth = cos_azimuth / cos(solar_altitude_angle * PI / 180);
  cos_azimuth = cos_azimuth / cos(latitude * PI / 180);

  if (cos_azimuth > 1)
    cos_azimuth = 1;
  if (cos_azimuth < -1)
    cos_azimuth = -1;
  if (sin_azimuth <= 0)
    azimuth = -acos(cos_azimuth) * 180 / PI;
  if (sin_azimuth > 0)
    azimuth = acos(cos_azimuth) * 180 / PI;
  return azimuth;
}

/**
  * @brief  SUN_POS_t合法性输入参数检测
  * @param 
  *         SUN_POS_t*   
  * @retval 
  *         0,合法,其他值非法。
  */
uint8_t validate_inputs(SUN_POS_t *spa) {
  if ((spa->year < -2000) || (spa->year > 6000)) return 1;
  if ((spa->month < 1) || (spa->month > 12)) return 2;
  if ((spa->day < 1) || (spa->day > 31)) return 3;
  if ((spa->hour < 0) || (spa->hour > 24)) return 4;
  if ((spa->minute < 0) || (spa->minute > 59)) return 5;
  if ((spa->second < 0) || (spa->second > 59)) return 6;
  if ((spa->hour == 24) && (spa->minute > 0)) return 5;
  if ((spa->hour == 24) && (spa->second > 0)) return 6;

  if (fabs(spa->timezone) > 18) return 8;
  if (fabs(spa->longitude) > 180) return 9;
  if (fabs(spa->latitude) > 90) return 10;

  return 0;
}

/**
  * @brief  根据输入参数,计算太阳位置信息
  * @param 
  *         SUN_POS_t*   
  * @retval 
  *         0,合法,其他值非法。
  */
// 计算Sun Position
uint8_t sunpos_cal(SUN_POS_t *spa) {
  const uint8_t month_days[] = { 31, 28, 31, 30, 31, 30, 31, 30, 31, 30, 31, 30, 31 };

  uint8_t ret;
  uint16_t ndays;
  float std_minutes;
  float local_mean_solar_time;
  uint16_t index;
  ret = 0;
  ret = validate_inputs(spa);
  if (ret != 0) return ret;
  /*计算从元旦计,天数*/
  spa->ndays = Date2Days(spa->year%100, spa->month, spa->day);
  /*计算当天时差(真太阳时-平太阳时)*/
  spa->error_time = time_error(spa->ndays);
  /*计算当天赤纬角*/
  spa->solar_declination_angle = solar_declination_angle_cal(spa->ndays);
  /*计算当天日升角*/
  spa->solar_hour_angle_sunrise = solar_hour_angle_sunrise_cal(spa->latitude, spa->solar_declination_angle);
  /*计算标准时间,分钟*/
  std_minutes = spa->hour*60 + spa->minute +spa->second/60.0f;
  /*计算当地当天平太阳时*/
  spa->mean_solar_time = local_mean_solar_time_cal(std_minutes, spa->timezone,spa->longitude);
  /*计算当天真太阳时*/
  spa->apparent_solar_time = spa->mean_solar_time+spa->error_time;
  /*计算当时高度角*/
  spa->solar_hour_angle = spa->apparent_solar_time/4.0f -180; // 注意范围: -180 ~ 180
  spa->solar_altitude_angle = solar_altitude_angle_cal(spa->latitude, spa->solar_declination_angle, spa->solar_hour_angle);
  /*计算当时方位角*/
  spa->solar_azimuth_angle = solar_azimuth_angle_cal(spa->latitude, spa->solar_declination_angle, spa->solar_hour_angle, spa->solar_altitude_angle);
  /*当地当日日出标准时间*/
  local_mean_solar_time = (spa->solar_hour_angle_sunrise + 180.0f)*4 - spa->error_time;//平太阳时    
  spa->sunrise_std_time = local_std_time_cal(local_mean_solar_time, spa->timezone,spa->longitude);// 当地标准时间
  
  return 0;
}
sunrise.c
复制代码

主程序 sunrise.ino

复制代码
#include "Arduino.h"
#include "sunrise.h"

void setup() {
  Serial.begin(115200);
  // 2023.2.28
  Serial.println(Date2Days(23, 2, 28));  //59

  Serial.println("solar_declination_angle:");
  Serial.println(solar_declination_angle_cal(59));  //-8.67
  Serial.println("solar_hour_angle_sunrise:");
  Serial.println(solar_hour_angle_sunrise_cal(39.933, -8.67));  //-82.67
  Serial.println("time_error:");
  Serial.println(time_error(59));  //-13.08

  Serial.println("local_mean_solar_time_cal:");
  float mean_time =  local_mean_solar_time_cal(123.11, 8, 130);
  Serial.println(mean_time);
  Serial.println("local_std_time_cal:");
  Serial.println( local_std_time_cal(mean_time , 8, 130)); /*计算设备所在地的平太阳时*/

  Serial.println("-------------------------");
  SUN_POS_t spa;
  spa.year = 23;
  spa.month = 2;
  spa.day = 28;
  spa.hour = 12;
  spa.minute = 0;
  spa.second = 0;
  spa.latitude = 39.933;
  spa.longitude = 116.33;
  spa.timezone = 8;

  sunpos_cal(&spa);
  Serial.println("ndays:");
  Serial.println(spa.ndays);  //59
  Serial.println("solar_declination_angle:");
  Serial.println(spa.solar_declination_angle);
  Serial.println("solar_hour_angle_sunrise:");
  Serial.println(spa.solar_hour_angle_sunrise);
  Serial.println("time_error:");
  Serial.println(spa.error_time);
  Serial.println("solar_altitude_angle:");
  Serial.println(spa.solar_altitude_angle);
  Serial.println("solar_azimuth_angle:");
  Serial.println(spa.solar_azimuth_angle);
  Serial.println("sunrise_std_time:");
  Serial.println(spa.sunrise_std_time);
}

String gpsString;

void loop() {
}
复制代码

 

posted @   辛河  阅读(825)  评论(1编辑  收藏  举报
编辑推荐:
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
点击右上角即可分享
微信分享提示