32F通过经纬度坐标与日期计算当天日出日落时间
本帖最后由 DebugLab 于 2024-9-11 10:56 编辑(1)先计算出从格林威治时间公元2000年1月1日到计算日天数days;
(2)计算从格林威治时间公元2000年1月1日到计算日的世纪数t, 则t＝(days+UTo／360)／36525;（一个世纪的分数日）
(3)计算太阳的平黄径 ： L=280.460+36000.770×t;
(4)计算太阳的平近点角 ：G＝357.528+35999.050×t
(5)计算太阳的黄道经度 ：λ＝L+1.915×sinG+0.020xsin(2G);
(6)计算地球的倾角 ε＝23.4393-0.0130×t;
(7)计算太阳的偏差 δ＝arcsin(sinε×sinλ);
(8)计算格林威治时间的太阳时间角GHA： GHA=UTo-180-1.915×sinG-0.020×sin(2G) +2.466×sin(2λ)-0.053×sin(4λ)
(9)计算修正值e： e=arcos{[ sinh-sin(Glat)sin(δ)]/cos(Glat)cos(δ)}
(10)计算新的日出日落时间 ：UT＝UTo-(GHA+Long±e); 其中“+”表示计算日出时间,“-”表示计算日落时间;
(11)比较UTo和UT之差的绝对值,如果大于0.1°即0.007小时,把UT作为新的日出日落时间值,重新从第(2)步开始进行迭代计算,如果UTo和UT之差的绝对值小于0.007小时,则UT即为所求的格林威治日出日落时间;
(12)上面的计算以度为单位,即180°=12小时,因此需要转化为以小时表示的时间,再加上所在的时区数Zone,即要计算地的日出日落时间为 ：T=UT/15+Zone
上面的计算日出日落时间方法适用于小于北纬60°和南纬60°之间的区域,如果计算位置为西半球时,经度Long为负数。
计算结果和查询网站的结果会有1~2分钟左右的误差
不足之处是STC32_FPMU_LARGE库无法运算asin函数，所以无法使用STC32_FPMU_LARGE库运算
同时由于keil与单片机无法运行double双精度数据。只能运算float，所以输出数值略有误差。但是不影响时钟分钟数据。
适合用于路灯控制，与太阳能灯照明。降低非必要时间照明！
打印输出数据与网上查询的日出日落时间对比
#include <stdio.h>
#include <math.h>
#define PI 3.1415926
//返回全局变量 sunrize sunset
double sunrizetime,sunsettime;
void CalSunRise(double year,double month,double day,double jingdu,double weidu)
{
double rssj;//日升时间
double rlsj;//日落时间
double ztsj;//中天时间
double rszztsj;//日升至中天时间
double ztzrlsj;//中天至日落时间
double jsc;//均时差（小时）
double shiqu;//时区8
double JD1;//儒略日中间值
double JD2;//儒略日
double G4;//参考A
double H4;//参考B
double D8,D9,D10,D11,D12,D13,D14,D15,D16,D17,D18,D19,D20,D21,D22;
double F10,F11,F12,F13,F14,F15,F16,F17,F18,F19,F20,F21,F22;
double D25,D26,F26,H26;
double K8,K9,K10,K11,K12,K13,K14,K15,K16,K17,K18,K19,K20,K21,K22;
double M10,M11,M12,M13,M14,M15,M16,M17,M18,M19,M20,M21,M22;
double rzsjdyb;//日照时间的一半
double weiduhudu;//纬度弧度
double xzlzxz;//修正量正弦值
double R8,R9,R10,R11,R12,R13,R14,R15,R16,R17,R18,R19,R20,R21,R22;
double T10,T11,T12,T13,T14,T15,T16,T17,T18,T19,T20,T21,T22;
//jingdu=DFM(116,23,29);//天安门
//weidu=DFM(39,54,20);
shiqu=8;
day=day 0.5;
if(year>0)
{
year=year;
}
else
{
year=year-1;
}
if(month>2)
{
year=year;
month=month;
}
else
{
year=year-1;
month=month 12;
}
weiduhudu=weidu/180*PI;
G4=floor(year/100);
H4=2-G4 floor(G4/4);
JD1=floor(365.25*(year 4716)) floor(30.6001*(month 1)) day H4-1524.5;
if(JD1>=2299160.5)
{
JD2=JD1;
}
else
{
JD2=JD1-H4;
}
D8=JD2-jingdu/360;
D9=(D8-2451545)/36525;
D10=84381.448-46.815*D9-0.00059*D9*D9 0.001813*D9*D9*D9;
F10=D10/3600*PI/180;
D11=280.46645 36000.76983*D9 0.0003032*D9*D9;
F11=D11*PI/180;
D12=357.5291 35999.0503*D9 - 0.0001559*D9*D9 -0.00000048*D9*D9*D9;
F12=D12*PI/180;
D13=0.016708617 - 0.000042037*D9 - 0.0000001236*D9*D9;
F13=D13;
D14=(1.9146 - 0.004817*D9 -0.000014*D9*D9)*sin(F12) (0.019993 - 0.000101*D9)*sin(2*F12) 0.00029*sin(3*F12);
F14=D14*PI/180;
D15=D11 D14;
F15=D15*PI/180;
D16=D12 D14;
F16=D16*PI/180;
D17= 1.000001018*(1-F13*F13) / (1 F13*cos(F16));
F17=D17;
D18= 125.04 - 1934.136*D9;
F18=D18*PI/180;
D19=D15 - 0.00569 -0.00478*sin(D18);
F19=D19*PI/180;
F20=asin(sin(F10)*sin(F15));
D20=F20*180/PI;
D21=D10/3600 0.00256*cos(F18);
F21=D21/180*PI;
F22=asin(sin(F21)*sin(F19));
D22=F21/PI*180;
D25=(tan(F10/2))*(tan(F10/2));//sin/cos
D26=D25*sin(2*F11)-2*F13*sin(F12) 4*F13*D25*sin(F12)*cos(2*F11)-(D25*D25/2)*sin(4*F11)-(5/4*F13*F13)*sin(2*F12);
F26=D26*180/PI;
H26=F26*4;
//================
//xzlzxz=sin((16/F17 34)*-1/60*PI/180)；//34大气折射修正 单位：分
xzlzxz=sin((16/F17 34)*-1/60*PI/180);
rzsjdyb=acos((xzlzxz-sin(F22)*sin(weiduhudu))/(cos(F22)*cos(weiduhudu)))/PI*12;
K8=D8-rzsjdyb/24;
K9=(K8-2451545)/36525;
K10=84381.448-46.815*K9-0.00059*K9*K9 0.001813*K9*K9*K9;
M10=K10/3600*PI/180;
K11=280.46645 36000.76983*K9 0.0003032*K9*K9;
M11=K11*PI/180;
K12=357.5291 35999.0503*K9 - 0.0001559*K9*K9 -0.00000048*K9*K9*K9;
M12=K12*PI/180;
K13=0.016708617 - 0.000042037*K9 - 0.0000001236*K9*K9;
M13=K13;
K14=(1.9146 - 0.004817*K9 -0.000014*K9*K9)*sin(M12) (0.019993 - 0.000101*K9)*sin(2*M12) 0.00029*sin(3*M12);
M14=K14*PI/180;
K15=K11 K14;
M15=K15*PI/180;
K16=K12 K14;
M16=K16*PI/180;
K17= 1.000001018*(1-M13*M13) / (1 M13*cos(M16));
M17=K17;
K18= 125.04 - 1934.136*K9;
M18=K18*PI/180;
K19=K15 - 0.00569 -0.00478*sin(K18);
M19=K19*PI/180;
M20==asin(sin(M10)*sin(M15));
K20=M20*180/PI;
K21=K10/3600 0.00256*cos(M18);
M21=K21/180*PI;
M22=asin(sin(M21)*sin(M19));
K22=M22/PI*180;
//============================================
R8=D8 rzsjdyb/24;
R9=(R8-2451545)/36525;
R10=84381.448-46.815*R9-0.00059*R9*R9 0.001813*R9*R9*R9;
T10=R10/3600*PI/180;
R11=280.46645 36000.76983*R9 0.0003032*R9*R9;
T11=R11*PI/180;
R12=357.5291 35999.0503*R9 - 0.0001559*R9*R9 -0.00000048*R9*R9*R9;
T12=R12*PI/180;
R13=0.016708617 - 0.000042037*R9 - 0.0000001236*R9*R9;
T13=R13;
R14=(1.9146 - 0.004817*R9 -0.000014*R9*R9)*sin(T12) (0.019993 - 0.000101*R9)*sin(2*T12) 0.00029*sin(3*T12);
T14=R14*PI/180;
R15=R11 R14;
T15=R15*PI/180;
R16=R12 R14;
T16=R16*PI/180;
R17= 1.000001018*(1-T13*T13) / (1 T13*cos(T16));
T17=R17;
R18= 125.04 - 1934.136*R9;
T18=R18*PI/180;
R19=R15 - 0.00569 -0.00478*sin(R18);
T19=R19*PI/180;
T20=asin(sin(T10)*sin(T15));
R20=T20/PI*180;
R21=R10/3600 0.00256*cos(T18);
T21=R21/180*PI;
T22=asin(sin(T21)*sin(T19));
R22=T22/180*PI;
rszztsj=acos((xzlzxz-sin(M22)*sin(weiduhudu))/(cos(M22)*cos(weiduhudu)))/PI*12;
ztzrlsj=acos((xzlzxz-sin(T22)*sin(weiduhudu))/(cos(T22)*cos(weiduhudu)))/PI*12;
//? jsc;//均时差（小时）
jsc=H26/60;
//? ztsj;//中天时间
ztsj=12-jsc-(jingdu-shiqu*15)*4/60;
rssj=ztsj-rszztsj;//日升时间=中天时间-日升至中天时间
rlsj=ztsj ztzrlsj;//日落时间=中天时间 中天至日落时间
sunrizetime=rssj;
sunsettime=rlsj;
}
void printtime(double t)
{
printf("%d:%d:%d\n",(int)(floor(t)),(int)floor((t-floor(t))*60),(int)floor((t*60-floor(t*60))*60));
}
double DFM(double d,double f,double m)
{
return d f/60 m/3600;
}
void main()
{
CalSunRise(2013,11,11,DFM(116,23,29),DFM(39,54,20));
printtime(sunrizetime);
printtime(sunsettime);
}#include <stdio.h>
#include <time.h>
#include <math.h>
#define RADEG ( 180.0 / M_PI )
#define sind(x)sin((x)*DEGRAD)
#define cosd(x)cos((x)*DEGRAD)
#define tand(x)tan((x)*DEGRAD)
#define DEGRAD ( M_PI / 180.0 )
#define atand( x ) (RADEG*atan(x))
#define asind(x) (RADEG*asin(x))
#define acosd(x) (RADEG*acos(x))
#define atan2d(y,x) (RADEG*atan2(y,x))
#define days_since_2000_Jan_0(y,m,d) \
(367*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530)
void sunriset( int year, int month, int day, double lon,
double lat, double *rise, double *set );
//year, month, date = calendar date, 1801-2099 only.
//Eastern longitude positive, Western longitude negative
//Northern latitude positive, Southern latitude negative
//The longitude value is critical in this function!
//*rise = where to store the rise time
//*set= where to store the settime
void sunpos( double d, double *lon, double *r );
//Computes the Sun's ecliptic longitude and distance at an instant given in d
//number of days since 2000 Jan 1
//The Sun's ecliptic latitude is not computed, it's always very near 0.
void sun_RA_dec( double d, double *RA, double *dec, double *r );
//Compute Sun's ecliptic longitude add slant
double revolution( double x );
//Reduce angle to within 0..360 degrees
double rev180( double x );
//Reduce angle to within -180..+180 degrees
double GMST0( double d );
//This function computes the Greenwich Mean Sidereal Time at 0h UT
//GMST is then the sidereal time at Greenwich at any time of the day
//GMST = (GMST0) + UT * (366.2422/365.2422)
int main()
{
//get system time
time_t t;
t = time( NULL );
struct tm *local = localtime( &t );
int year = local->tm_year + 1900;
int month = local->tm_mon+1;
int day = local->tm_mday;
int hour = local->tm_hour;
int minute;
//get coordinate
double lon, lat;
printf("Eastern longitude positive, Western longitude negative\n");
printf("Northern latitude positive, Southern latitude negative\n");
printf("such as:114.31667 30.51667\n");
printf("please input longitude and latitude:\n");
//this coordinate is Wuhan
scanf("%lf %lf", &lon, &lat);
double sunrise, sunset;
sunriset( year, month, day, lon, lat,
&sunrise, &sunset );
//calculate time zone
struct tm *gmt = gmtime( &t );
int zone = hour - gmt->tm_hour + (day - gmt->tm_mday)*24;
//calcute sunrise and sunset time
hour = sunrise+zone;
minute = (sunrise - hour+zone)*60;
if( hour > 24 ) hour -= 24;
else if( hour < 0 ) hour += 24;
printf("sunrise:%2d:%2d\n", hour, minute);
hour = sunset+zone;
minute = (sunset - hour+zone)*60;
if( hour > 24 ) hour -= 24;
else if( hour < 0 ) hour += 24;
printf("sunset: %2d:%2d\n", hour, minute);
return 0;
}
void sunriset( int year, int month, int day, double lon, double lat, double *rise, double *set )
{
double sr; //sun's distance
double sRA; //sun's ecliptic longitude
double sdec; //slant of the sun
double t; //offset
double d = days_since_2000_Jan_0( year, month, day ) + 0.5 - lon / 360.0;
//Compute d of 12h local mean solar time, from 2000 Jan 1 to now
double sidtime = revolution( GMST0(d) + 180.0 + lon );
//Compute local sidereal time of this moment
double altit = -35.0/60.0;
//Compute Sun's RA + Decl at this moment
sun_RA_dec( d, &sRA, &sdec, &sr );
//Compute time when Sun is at south - in hours UT
double tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
//Compute the Sun's apparent radius, degrees
double sradius = 0.2666 / sr;
//Do correction to upper limb, if necessary
altit -= sradius;
//Compute the diurnal arc that the Sun traverses to reach
//the specified altitude altit:
double cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
( cosd(lat) * cosd(sdec) );
if ( cost >= 1.0 )
t = 0.0; //Sun always below altit
else if ( cost <= -1.0 )
t = 12.0; //Sun always above altit
else
t = acosd(cost)/15.0; //The diurnal arc, hours
//Store rise and set times - in hours UT
*rise = tsouth - t;
*set= tsouth + t;
}
void sunpos( double d, double *lon, double *r )
{
//Compute mean elements
double M = revolution( 356.0470 + 0.9856002585 * d ); //anomaly
double w = 282.9404 + 4.70935E-5 * d; //sun's ecliptic longitude
double e = 0.016709 - 1.151E-9 * d; //earth's offset
//Compute true longitude and radius vector
double E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) ); //eccentric anomaly
double x = cosd(E) - e;
double y = sqrt( 1.0 - e*e ) * sind(E); //orbit coordinate
*r = sqrt( x*x + y*y ); //Solar distance
double v = atan2d( y, x ); //True anomaly
*lon = v + w; //True solar longitude
if ( *lon >= 360.0 )
*lon -= 360.0; //Make it 0..360 degrees
}
void sun_RA_dec( double d, double *RA, double *dec, double *r )
{
double lon;
//Compute Sun's ecliptical coordinates
sunpos( d, &lon, r );
//Compute ecliptic rectangular coordinates (z=0)
double x = *r * cosd(lon);
double y = *r * sind(lon);
/* Compute obliquity of ecliptic (inclination of Earth's axis) */
double obl_ecl = 23.4393 - 3.563E-7 * d;
/* Convert to equatorial rectangular coordinates - x is unchanged */
double z = y * sind(obl_ecl);
y = y * cosd(obl_ecl);
/* Convert to spherical coordinates */
*RA = atan2d( y, x );
*dec = atan2d( z, sqrt(x*x + y*y) );
}
double revolution( double x )
{
return( x - 360.0 * floor( x / 360 ) );
}
double rev180( double x )
{
return( x - 360.0 * floor( x / 360 + 0.5 ) );
}
double GMST0( double d )
{
//Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 deg
//L = M + w, as defined in sunpos().
//Any decent C compiler will add the constants at compile time
double sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) +
( 0.9856002585 + 4.70935E-5 ) * d );
return sidtim0;
}
用FPU计算的吗
页:
[1]