基于微控制器的日出/设置算法实现
Posted
技术标签:
【中文标题】基于微控制器的日出/设置算法实现【英文标题】：Microcontroller based sunrise/set algorithm implementation 【发布时间】：2017-06-12 19:39:35 【问题描述】：两个问题：
我正在尝试实现一种算法来确定日出和日落时间，并将纬度、经度和日期作为微控制器上的输入。我知道有算法可以做到这一点，但我无法找到有关如何将其放在微控制器上的任何细节。这可能吗？
其次，如果我要实现这个，我大约需要多少内存以及我必须使用什么样的控制器？
请告诉我。
【问题讨论】：
是的，这是可能的。但是在缩小第二个问题的答案之前，您需要考虑更多的输入。例如，计算需要多精确？需要多快/多长时间执行一次计算？地理定位将如何提供给微控制器？日出/日落时间将如何提供给用户？这是安全关键系统的一部分还是一个爱好项目？ +/-15 分钟的公差就可以了。计算需要每天执行一次 - 想法是在日出和日落时触发输出。地理位置将作为用户输入提供（通过触摸屏或智能手机）。这不是一种爱好，也不是一个安全关键系统——只是想点子:) 我怀疑接口需求比日出/日落需求更重要。我的意思是，如果微控制器可以管理触摸屏或蓝牙无线电与智能手机通信，那么它可能可以处理算法。要考虑的另一件事是您的功率预算。设备是否由电池供电？ 微控制器与此有什么关系？算法就是算法，无论您的目标是什么，您都可以以速度或大小或两者的准确性为代价来调整算法。 具体目标（架构和其他资源）使答案更加模糊。初学者的目标架构是什么（添加一个或几个标签）。编译器和语言在优化和/或实现中起着重要作用。内存不是这里唯一的因素，只是其中之一。 【参考方案1】：嗯，在 MCU 上没有太多内存，而且 FPU 并不总是存在，所以你需要自己做很多事情。有几种计算方法。
日出日历
这是最简单的，但需要存储日历数据（可以在程序存储器中，因此不需要 RAM）。这种方法仅限于单一地理位置，因此如果位置随时间变化，则它不可用。欲了解更多信息，请参阅：
Calculate whether it is close to dawn/dusk based on sunrise/sunset?开普勒方程
如果编码正确，这将工作很长时间（例如 1000 年）。但是这种方法需要 FPU 计算。如果不存在，您需要自己编写代码（例如：double,sin,cos,+,-,*,/）
开普勒方程将为您提供地球和太阳在日心坐标系中的实际相对位置。如果您添加地球的每日自转和本地 NEH（东北高度/alt）变换，您将获得任何地理位置和时间的方位坐标中的太阳位置。然后只需根据日出/日落几何限制检查海拔高度。欲了解更多信息，请参阅：
How to draw sky chart? Is it possible to make realistic n-body solar system simulation in matter of size and mass? calculate the time when the sun is X degrees below/above the Horizon Understanding 4x4 homogenous transform matrices这种方法需要 RAM 来存储实际坐标（每个位置 3 个双精度 + 一个温度）、转换 3 个双精度 4x4 变换矩阵（地球、NEH、温度）和少量迭代变量。我估计 512 字节就足够了。
直接方程
您需要为此谷歌。有一些近似方程可以直接为您提供任何地理位置的太阳方位角坐标。它比 #2 简单一些，但需要更高级的测角函数，例如 arctan,acos，并且不确定是否还需要一些双曲函数。如果您没有 FPU 支持，那么这种方法可能难以以所需的精度实现。
[备注]
为了提高准确性，您需要将几何太阳的位置转换为视觉位置。有表多少仰角偏移...
【讨论】：
对于（1）；通过日历数据，您是指特定位置的日出和日落时间吗？一年来，这个数据集在我的电脑上大约有 4 Kb。 @genesis 您可以将表格存储在程序存储器中，而不是现在的 RAM 中三次，所以它只是几个 2D 点......如果你没有 FPU，它也可以在整数算术上完成，请参阅How can i produce multi point linear interpolation?【参考方案2】：日出/日落时间可以使用互联网上可用的算法计算，这通常需要三角函数，这反过来意味着支持浮点。
如果没有 FPU，那么您将不得不使用编译器库。或查找表，an approximation or even CORDIC：
如果你：
没有 FPU 或不想使用浮点数。 不想将编译器库用于三角函数。 需要近似值无法满足的精度。 无法使用固定表，因为时区可能会改变。然后我会使用 CORDIC，它似乎很自然地适用于定点数。
此示例使用来自Ed Williams's website 的算法，其中还列出了example of this 算法。我选择使用 Q10.21 格式的定点数，因为它似乎可以提供尽可能高的精度而不会出现算术溢出。因此目标应该是 32 位的。可以改变宽度，但是会影响精度。
我正在使用John Burkardt 的自定义版本的 CORDIC 函数，这些函数经过修改以使用定点算术并减少乘法运算的数量。他在他的 cmets 中提到，他参考了 Wikipedia 上的 MATLAB 示例，这有助于理解算法的工作原理。
我已经针对NOAA sunrise/sunset calculator 对几个时区（+08:00、+04:00、-03:00）进行了几轮测试，这个示例的准确度似乎为 ±1 分钟。 虽然时区被指定为参数，但示例中使用的时间（和 DST 设置）取自您的操作系统。
关于定点运算的一些注意事项：
可以使用常规 +/- 运算符完成加法和减法。 乘法和除法必须分别使用 IFMUL 和 IFDIV 宏来完成。 模除法可以用 % 完成，但我添加了一个 IFMOD 宏以与除法运算保持一致。main.c:
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "fixedpoint.h"
#include "cordic_helper.h"
#define ZENITH FROMFLOAT(-.83f)
/*
Computes the sunrise/sunset.
Based on the algorithm described here:
- https://***.com/questions/7064531/sunrise-sunset-times-in-c
- https://edwilliams.org/sunrise_sunset_algorithm.htm
- https://www.edwilliams.org/sunrise_sunset_example.htm
@param dayOfYear The day of year (1 - 365).
@param lat The latitude, in decimal degrees.
@param lng The longitude, in decimal degrees.
@param localOffset The timezone offset in minutes.
@param daylightSavings Whether DST is in effect. 1 = DST in effect, 0 = no DST.
@param sunrise Whether to compute the sunrise. Otherwise, the sunset will be computed.
@return The sunrise/sunset, in minutes since 00:00 of local time. */
static int calculateSunriseSunset(int dayOfYear, fixedfloat_t lat, fixedfloat_t lng, int localOffset, int daylightSavings, int sunrise);
static fixedfloat_t sunLocalHourAngle(fixedfloat_t lat, fixedfloat_t sinDec, fixedfloat_t cosDec);
/**
* Computes the day of year (where 1 = January 1st).
* Explained here: https://astronomy.stackexchange.com/questions/2407/calculate-day-of-the-year-for-a-given-date
* @param year The year
* @param month The month within the year (1-12)
* @param day The day of the month (1-31)
*/
static int computeDayOfYear(int year, int month, int day)
//1. first calculate the day of the year
int N1 = 275 * month / 9;
int N2 = (month + 9) / 12;
int N3 = (1 + ((year - 4 * (year / 4) + 2) / 3));
int N = N1 - (N2 * N3) + day - 30;
return N;
static int calculateSunriseSunset(int dayOfYear, fixedfloat_t lat, fixedfloat_t lng, int localOffset, int daylightSavings, int sunrise)
fixedfloat_t sinDec;
fixedfloat_t cosDec;
fixedfloat_t cosH;
fixedfloat_t H;
fixedfloat_t time;
int utc;
fixedfloat_t t;
fixedfloat_t M;
fixedfloat_t L;
fixedfloat_t RA;
int Lquadrant;
int RAquadrant;
int r;
fixedfloat_t temp;
//2. convert the longitude to hour value and calculate an approximate time
const fixedfloat_t lngHour = IFDIV(lng, FROMINT(15));
t = FROMINT(dayOfYear) + IFDIV((FROMINT(sunrise ? 6 : 18) - lngHour), FROMINT(24));
//3. calculate the Sun's mean anomaly
// from, https://physics.stackexchange.com/questions/80034: M = nt + M0, where n = 0.9856, M0 = 3.289
//M = (0.9856f * t) - 3.289f;
M = IFMUL(FROMFLOAT(0.9856f), t) - FROMFLOAT(3.289f);
//4. calculate the Sun's true longitude
//L = fmodf(M + (1.916f * sin((PI / 180)*M)) + (0.020f * sin(2 * (PI / 180) * M)) + 282.634f, 360.0f);
// Break up this expression to avoid arithmetic overflows.
temp = M + (IFMUL(FROMFLOAT(1.916f), sin(RADIANS(M)))) + IFMUL(FROMFLOAT(0.020f), sin(RADIANS(M) << 1));
if (TOINT(temp) + 283 >= 360)
temp -= FROMINT(360); // Avoid overflow.
L = IFMOD(temp + FROMFLOAT(282.634f), FROMINT(360));
//5a. calculate the Sun's right ascension
//RA = fmodf(180 / PI * atan(0.91764f * tan((PI / 180)*L)), 360.0f);
RA = IFMOD(DEGREES(atan(IFMUL(FROMFLOAT(0.91764f), tan(RADIANS(L))))), FROMINT(360));
//5b. right ascension value needs to be in the same quadrant as Lval
Lquadrant = TOINT(IFDIV(L, FROMINT(90))) * 90;
RAquadrant = TOINT(IFDIV(RA, FROMINT(90))) * 90;
RA = RA + FROMINT(Lquadrant - RAquadrant);
//5c. right ascension value needs to be converted into hours
RA = IFDIV(RA, FROMINT(15));
//6. calculate the Sun's declination
//sinDec = 0.39782f * sin((PI / 180) * L);
sinDec = IFMUL(FROMFLOAT(0.39782f), sin(RADIANS(L)));
cosDec = cos(asin(sinDec));
//7a. calculate the Sun's local hour angle
//cosH = (sin((PI/180)*ZENITH) - (sinDec * sin((PI/180)*lat))) / (cosDec * cos((PI/180)*lat));
cosH = sunLocalHourAngle(lat, sinDec, cosDec);
/*
if (cosH > 1)
the sun never rises on this location (on the specified date)
if (cosH < -1)
the sun never sets on this location (on the specified date)
*/
//7b. finish calculating H and convert into hours
H = sunrise ? FROMINT(360) - DEGREES(acos(cosH)) : DEGREES(acos(cosH));
H = IFDIV(H, FROMINT(15));
//8. calculate local mean time of rising/setting
//time = H + RA - (0.06571 * t) - 6.622;
time = H + RA - IFMUL(FROMFLOAT(0.06571), t) - FROMFLOAT(6.622);
//9. adjust back to UTC
utc = (TOINT(time - lngHour) % 24 * 60) + TOINT(IFMUL(FRAC(time - lngHour), FROMINT(60)));
//10. convert UTC value to local time zone of latitude/longitude
r = utc + localOffset + daylightSavings * 60;
return r < 0 ? r + (24 * 60) : (r % (24 * 60));
static fixedfloat_t sunLocalHourAngle(fixedfloat_t lat, fixedfloat_t sinDec, fixedfloat_t cosDec)
//cosH = (sin((PI/180)*ZENITH) - (sinDec * sin((PI/180)*lat))) / (cosDec * cos((PI/180)*lat));
fixedfloat_t sinValue, cosValue;
sincos(RADIANS(lat), &sinValue, &cosValue);
return IFDIV(sin(RADIANS(ZENITH)) - IFMUL(sinDec, sinValue), IFMUL(cosDec, cosValue));
int main(int argc, char *argv[])
time_t now;
float lat;
float lng;
int localOffset;
struct tm *tm;
int localSunriseT, localSunsetT;
if (argc != 4)
printf("%s <latitude> <longitude> <timezone offset in minutes>\nExample: %s 1.3364926804464332 103.74412865673372 480\n", argv[0], argv[0]);
return EINVAL;
lat = (float) atof(argv[1]);
lng = (float) atof(argv[2]);
localOffset = atoi(argv[3]);
printf("Coordinates: %f,%f\nTimezone Offset: %d\n", lat, lng, localOffset);
now = time(NULL);
tm = localtime(&now);
// Get the time of day.
//const int dayOfYear = computeDayOfYear(tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday);
const int dayOfYear = tm->tm_yday + 1; // If you have a source for the day in the year (1 to 365, inclusive).
const int dst = tm->tm_isdst;
// OR, set it manually:
/* const int dayOfYear = computeDayOfYear(2021, 1, 1);
const int dst = 0; */
localSunriseT = calculateSunriseSunset(dayOfYear, FROMFLOAT(lat), FROMFLOAT(lng), localOffset, dst, 1);
localSunsetT = calculateSunriseSunset(dayOfYear, FROMFLOAT(lat), FROMFLOAT(lng), localOffset, dst, 0);
printf("Sunrise: %02d:%02d\nSunset: %02d:%02d\n", localSunriseT / 60, localSunriseT % 60, localSunsetT / 60, localSunsetT % 60);
return 0;
定点.h：
typedef int fixedfloat_t;
// Must be double the length of fixedfloat_t.
typedef long long int fixedfloat_double_t;
/* Fixed-point format: Q10.21
Whole numbers: 0-1024
Bit 31 indicates sign */
#define FRAC_PLACES 21
#define FORMAT (1 << FRAC_PLACES)
#define FRAC_MASK (FORMAT - 1)
#define FRAC(a) ((a) & FRAC_MASK)
#define FROMFLOAT(a) ((fixedfloat_t)(FORMAT * (a)))
#define TOFLOAT(a) (((float) (a)) / FORMAT)
#define FROMINT(a) ((fixedfloat_t)((a) << FRAC_PLACES))
#define TOINT(a) ((int)((a) >> FRAC_PLACES))
//#define IFADD(a, b) ((a) + (b)) // Just for illustration
//#define IFSUB(a, b) ((a) - (b)) // Just for illustration
#define IFMUL(a, b) ((fixedfloat_t)((((fixedfloat_double_t)(a)) * (b)) >> FRAC_PLACES))
#define IFDIV(a, b) ((fixedfloat_t)((((fixedfloat_double_t)(a)) << FRAC_PLACES) / (b)))
#define IFMOD(a, b) ((a) % (b))
#define IFFLOORF(a) ((a) & ~FRAC_MASK)
#define IFABS(a) ((a) < 0 ? -(a) : (a))
#define PI FROMFLOAT(3.141592653589793)
#define TWOPI FROMFLOAT(2 * 3.141592653589793)
#define MIN(a, b) ((a) < (b) ? (a) : (b))
// Order was changed and some shifting is done, to avoid overflows.
#define RADIANS(a) IFDIV(IFMUL((a), PI >> 1), FROMINT(180 >> 1))
#define DEGREES(a) IFMUL(IFDIV((a), PI >> 1), FROMINT(180 >> 1))
cordic_helper.h：
/* CORDIC Helper functions. Presently developed around the CORDIC implementation from https://people.sc.fsu.edu/~jburkardt/cpp_src/cordic/cordic.html */
#define ITERATIONS 20
fixedfloat_t ifsin(fixedfloat_t a);
fixedfloat_t ifcos(fixedfloat_t a);
void ifsincos(fixedfloat_t a, fixedfloat_t *sinout, fixedfloat_t *cosout);
fixedfloat_t iftan(fixedfloat_t a);
fixedfloat_t ifacos(fixedfloat_t a);
fixedfloat_t ifasin(fixedfloat_t a);
#define ifatan(a) ifatan2(FROMINT(1), a)
fixedfloat_t ifatan2(fixedfloat_t x, fixedfloat_t y);
#define sin(x) ifsin(x)
#define cos(x) ifcos(x)
#define sincos(x, y, z) ifsincos(x, y, z)
#define tan(x) iftan(x)
#define acos(x) ifacos(x)
#define asin(x) ifasin(x)
#define atan(x) ifatan(x)
#define atan2(x, y) ifatan2(x, y)
cordic_helper.c：
#include "fixedpoint.h"
#include "cordic.h"
#include "cordic_helper.h"
fixedfloat_t ifsin(fixedfloat_t a)
fixedfloat_t cosVal;
fixedfloat_t sinVal;
sincos(a, &sinVal, &cosVal);
return sinVal;
fixedfloat_t ifcos(fixedfloat_t a)
fixedfloat_t cosVal;
fixedfloat_t sinVal;
sincos(a, &sinVal, &cosVal);
return cosVal;
void ifsincos(fixedfloat_t a, fixedfloat_t *sinout, fixedfloat_t *cosout)
cossin_cordic(a, ITERATIONS, cosout, sinout);
fixedfloat_t iftan(fixedfloat_t a)
fixedfloat_t cosVal;
fixedfloat_t sinVal;
sincos(a, &sinVal, &cosVal);
return IFDIV(sinVal, cosVal);
fixedfloat_t ifacos(fixedfloat_t a)
return arccos_cordic(a, ITERATIONS);
fixedfloat_t ifasin(fixedfloat_t a)
return arcsin_cordic(a, ITERATIONS);
fixedfloat_t ifatan2(fixedfloat_t x, fixedfloat_t y)
return arctan_cordic(x, y, ITERATIONS);
cordic.h：
/* Define if a fixed value of N will be used. If so, then define KPROD_VALUE, which determines the K value for correction.
The value should be selected from the kprod array, with the index being the value of N.
If N is greater than KPROD_LENGTH, choose the last value. */
#define FIXED_N 1
#define KPROD_VALUE FROMFLOAT(0.60725293500924945172)
fixedfloat_t arccos_cordic(fixedfloat_t t, fixedfloat_t n);
fixedfloat_t arcsin_cordic (fixedfloat_t t, int n );
fixedfloat_t arctan_cordic (fixedfloat_t x, fixedfloat_t y, int n );
void cossin_cordic(fixedfloat_t beta, int n, fixedfloat_t *c, fixedfloat_t *s);
cordic.c:
# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>
# include "fixedpoint.h"
# include "cordic.h"
static fixedfloat_t angle_shift(fixedfloat_t alpha, fixedfloat_t beta);
# define ANGLES_LENGTH 60
static const fixedfloat_t angles[ANGLES_LENGTH] =
FROMFLOAT(7.8539816339744830962E-01),
FROMFLOAT(4.6364760900080611621E-01),
FROMFLOAT(2.4497866312686415417E-01),
FROMFLOAT(1.2435499454676143503E-01),
FROMFLOAT(6.2418809995957348474E-02),
FROMFLOAT(3.1239833430268276254E-02),
FROMFLOAT(1.5623728620476830803E-02),
FROMFLOAT(7.8123410601011112965E-03),
FROMFLOAT(3.9062301319669718276E-03),
FROMFLOAT(1.9531225164788186851E-03),
FROMFLOAT(9.7656218955931943040E-04),
FROMFLOAT(4.8828121119489827547E-04),
FROMFLOAT(2.4414062014936176402E-04),
FROMFLOAT(1.2207031189367020424E-04),
FROMFLOAT(6.1035156174208775022E-05),
FROMFLOAT(3.0517578115526096862E-05),
FROMFLOAT(1.5258789061315762107E-05),
FROMFLOAT(7.6293945311019702634E-06),
FROMFLOAT(3.8146972656064962829E-06),
FROMFLOAT(1.9073486328101870354E-06),
FROMFLOAT(9.5367431640596087942E-07),
FROMFLOAT(4.7683715820308885993E-07),
FROMFLOAT(2.3841857910155798249E-07),
FROMFLOAT(1.1920928955078068531E-07),
FROMFLOAT(5.9604644775390554414E-08),
FROMFLOAT(2.9802322387695303677E-08),
FROMFLOAT(1.4901161193847655147E-08),
FROMFLOAT(7.4505805969238279871E-09),
FROMFLOAT(3.7252902984619140453E-09),
FROMFLOAT(1.8626451492309570291E-09),
FROMFLOAT(9.3132257461547851536E-10),
FROMFLOAT(4.6566128730773925778E-10),
FROMFLOAT(2.3283064365386962890E-10),
FROMFLOAT(1.1641532182693481445E-10),
FROMFLOAT(5.8207660913467407226E-11),
FROMFLOAT(2.9103830456733703613E-11),
FROMFLOAT(1.4551915228366851807E-11),
FROMFLOAT(7.2759576141834259033E-12),
FROMFLOAT(3.6379788070917129517E-12),
FROMFLOAT(1.8189894035458564758E-12),
FROMFLOAT(9.0949470177292823792E-13),
FROMFLOAT(4.5474735088646411896E-13),
FROMFLOAT(2.2737367544323205948E-13),
FROMFLOAT(1.1368683772161602974E-13),
FROMFLOAT(5.6843418860808014870E-14),
FROMFLOAT(2.8421709430404007435E-14),
FROMFLOAT(1.4210854715202003717E-14),
FROMFLOAT(7.1054273576010018587E-15),
FROMFLOAT(3.5527136788005009294E-15),
FROMFLOAT(1.7763568394002504647E-15),
FROMFLOAT(8.8817841970012523234E-16),
FROMFLOAT(4.4408920985006261617E-16),
FROMFLOAT(2.2204460492503130808E-16),
FROMFLOAT(1.1102230246251565404E-16),
FROMFLOAT(5.5511151231257827021E-17),
FROMFLOAT(2.7755575615628913511E-17),
FROMFLOAT(1.3877787807814456755E-17),
FROMFLOAT(6.9388939039072283776E-18),
FROMFLOAT(3.4694469519536141888E-18),
FROMFLOAT(1.7347234759768070944E-18) ;
#ifdef FIXED_N
#ifndef KPROD_VALUE
#error KPROD_VALUE must be defined (refer to the KPROD table, for N iterations)
#endif
#else
# define KPROD_LENGTH 33
static const fixedfloat_t kprod[KPROD_LENGTH] =
FROMFLOAT(0.70710678118654752440), //1
FROMFLOAT(0.63245553203367586640),
FROMFLOAT(0.61357199107789634961),
FROMFLOAT(0.60883391251775242102),
FROMFLOAT(0.60764825625616820093),
FROMFLOAT(0.60735177014129595905),
FROMFLOAT(0.60727764409352599905),
FROMFLOAT(0.60725911229889273006),
FROMFLOAT(0.60725447933256232972),
FROMFLOAT(0.60725332108987516334), //10
FROMFLOAT(0.60725303152913433540),
FROMFLOAT(0.60725295913894481363),
FROMFLOAT(0.60725294104139716351),
FROMFLOAT(0.60725293651701023413),
FROMFLOAT(0.60725293538591350073),
FROMFLOAT(0.60725293510313931731),
FROMFLOAT(0.60725293503244577146),
FROMFLOAT(0.60725293501477238499),
FROMFLOAT(0.60725293501035403837),
FROMFLOAT(0.60725293500924945172), //20
FROMFLOAT(0.60725293500897330506),
FROMFLOAT(0.60725293500890426839),
FROMFLOAT(0.60725293500888700922),
FROMFLOAT(0.60725293500888269443),
FROMFLOAT(0.60725293500888161574),
FROMFLOAT(0.60725293500888134606),
FROMFLOAT(0.60725293500888127864),
FROMFLOAT(0.60725293500888126179),
FROMFLOAT(0.60725293500888125757),
FROMFLOAT(0.60725293500888125652), //30
FROMFLOAT(0.60725293500888125626),
FROMFLOAT(0.60725293500888125619),
FROMFLOAT(0.60725293500888125617) ;
#endif
/******************************************************************************/
static fixedfloat_t angle_shift (fixedfloat_t alpha, fixedfloat_t beta )
/******************************************************************************/
/*
Purpose:
ANGLE_SHIFT shifts angle ALPHA to lie between BETA and BETA+2PI.
Discussion:
The input angle ALPHA is shifted by multiples of 2 * PI to lie
between BETA and BETA+2*PI.
The resulting angle GAMMA has all the same trigonometric function
values as ALPHA.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
19 January 2012
Author:
John Burkardt
Parameters:
Input, double ALPHA, the angle to be shifted.
Input, double BETA, defines the lower endpoint of
the angle range.
Output, double ANGLE_SHIFT, the shifted angle.
*/
return (alpha < beta) ? beta - IFMOD( beta - alpha, TWOPI) + TWOPI : beta + IFMOD( alpha - beta, TWOPI);
/******************************************************************************/
fixedfloat_t arccos_cordic (fixedfloat_t t, fixedfloat_t n )
/******************************************************************************/
/*
Purpose:
ARCCOS_CORDIC returns the arccosine of an angle using the CORDIC method.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
19 January 2012
Author:
John Burkardt
Reference:
Jean-Michel Muller,
Elementary Functions: Algorithms and Implementation,
Second Edition,
Birkhaeuser, 2006,
ISBN13: 978-0-8176-4372-0,
LC: QA331.M866.
Parameters:
Input, double T, the cosine of an angle. -1 <= T <= 1.
Input, int N, the number of iterations to take.
A value of 10 is low. Good accuracy is achieved with 20 or more
iterations.
Output, double ARCCOS_CORDIC, an angle whose cosine is T.
Local Parameters:
Local, double ANGLES(60) = arctan ( (1/2)^(0:59) );
*/
fixedfloat_t angle;
int i;
int j;
int shiftamount;
int sigma;
int sign_z2;
fixedfloat_t theta;
fixedfloat_t x1;
fixedfloat_t x2;
fixedfloat_t y1;
fixedfloat_t y2;
theta = FROMINT(0);
x1 = FROMINT(1);
y1 = FROMINT(0);
for ( j = 1; j <= n; j++ )
sign_z2 = ( y1 < 0 ) ? -1 : 1;
sigma = ( t <= x1 ) ? +sign_z2 : -sign_z2;
angle = (j <= ANGLES_LENGTH) ? angles[j - 1] : angle >> 1;
shiftamount = j - 1; // Index of the root to compute, where 3 = square root.
for ( i = 1; i <= 2; i++ )
// sigma * poweroftwo, where poweroftwo = 0.5, 0.25...
x2 = x1 - ((sigma < 0 ? -y1 : y1) >> shiftamount);
y2 = ((sigma < 0 ? -x1 : x1) >> shiftamount) + y1;
x1 = x2;
y1 = y2;
theta = theta + ((sigma < 0 ? -angle : angle) << 1); // 2 * sigma * angle
t = t + ((t >> shiftamount) >> shiftamount); // t + t * poweroftwo * poweroftwo
return theta;
/******************************************************************************/
fixedfloat_t arcsin_cordic (fixedfloat_t t, int n )
/******************************************************************************/
/*
Purpose:
ARCSIN_CORDIC returns the arcsine of an angle using the CORDIC method.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
19 January 2012
Author:
John Burkardt
Reference:
Jean-Michel Muller,
Elementary Functions: Algorithms and Implementation,
Second Edition,
Birkhaeuser, 2006,
ISBN13: 978-0-8176-4372-0,
LC: QA331.M866.
Parameters:
Input, double T, the sine of an angle. -1 <= T <= 1.
Input, int N, the number of iterations to take.
A value of 10 is low. Good accuracy is achieved with 20 or more
iterations.
Output, double ARCSIN_CORDIC, an angle whose sine is T.
Local Parameters:
Local, double ANGLES(60) = arctan ( (1/2)^(0:59) );
*/
fixedfloat_t angle;
int i;
int j;
int shiftamount;
int sigma;
int sign_z1;
fixedfloat_t theta;
fixedfloat_t x1;
fixedfloat_t x2;
fixedfloat_t y1;
fixedfloat_t y2;
theta = FROMINT(0);
x1 = FROMINT(1);
y1 = FROMINT(0);
for ( j = 1; j <= n; j++ )
sign_z1 = ( x1 < 0 ) ? -1 : 1;
sigma = ( y1 <= t ) ? +sign_z1 : -sign_z1;
angle = ( j <= ANGLES_LENGTH ) ? angles[j-1] : angle >> 1;
shiftamount = j - 1; // Index of the root to compute, where 3 = square root.
for ( i = 1; i <= 2; i++ )
// sigma * poweroftwo, where poweroftwo = 0.5, 0.25...
x2 = x1 - ((sigma < 0 ? -y1 : y1) >> shiftamount);
y2 = ((sigma < 0 ? -x1 : x1) >> shiftamount) + y1;
x1 = x2;
y1 = y2;
theta = theta + ((sigma < 0 ? -angle : angle) << 1); // 2 * sigma * angle
t = t + ((t >> shiftamount) >> shiftamount);
return theta;
/******************************************************************************/
fixedfloat_t arctan_cordic (fixedfloat_t x, fixedfloat_t y, int n )
/******************************************************************************/
/*
Purpose:
ARCTAN_CORDIC returns the arctangent of an angle using the CORDIC method.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
15 June 2007
Author:
John Burkardt
Reference:
Jean-Michel Muller,
Elementary Functions: Algorithms and Implementation,
Second Edition,
Birkhaeuser, 2006,
ISBN13: 978-0-8176-4372-0,
LC: QA331.M866.
Parameters:
Input, double X, Y, define the tangent of an angle as Y/X.
Input, int N, the number of iterations to take.
A value of 10 is low. Good accuracy is achieved with 20 or more
iterations.
Output, double ARCTAN_CORDIC, the angle whose tangent is Y/X.
Local Parameters:
Local, double ANGLES(60) = arctan ( (1/2)^(0:59) );
*/
fixedfloat_t angle;
int j;
int shiftamount;
int sigma;
int sign_factor;
fixedfloat_t theta;
fixedfloat_t x1;
fixedfloat_t x2;
fixedfloat_t y1;
fixedfloat_t y2;
x1 = x;
y1 = y;
/*
Account for signs.
*/
if ( x1 < 0 && y1 < 0 )
x1 = -x1;
y1 = -y1;
if ( x1 < 0 )
x1 = -x1;
sign_factor = -1;
else if ( y1 < 0 )
y1 = -y1;
sign_factor = -1;
else
sign_factor = 1;
theta = FROMINT(0);
for ( j = 1; j <= n; j++ )
sigma = (y1 <= 0) ? 1 : -1;
angle = ( j <= ANGLES_LENGTH ) ? angles[j-1] : angle >> 1;
shiftamount = j - 1; // Index of the root to compute, where 3 = square root.
// sigma * poweroftwo, where poweroftwo = 0.5, 0.25...
x2 = x1 - ((sigma < 0 ? -y1 : y1) >> shiftamount);
y2 = ((sigma < 0 ? -x1 : x1) >> shiftamount) + y1;
theta = theta - (sigma < 0 ? -angle : angle);
x1 = x2;
y1 = y2;
theta = sign_factor < 0 ? -theta : theta;
return theta;
/******************************************************************************/
void cossin_cordic (fixedfloat_t beta, int n, fixedfloat_t *c, fixedfloat_t *s )
/******************************************************************************/
/*
Purpose:
COSSIN_CORDIC returns the sine and cosine of an angle by the CORDIC method.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
19 January 2012
Author:
Based on MATLAB code in a Wikipedia article.
C++ version by John Burkardt
Parameters:
Input, double BETA, the angle (in radians).
Input, int N, the number of iterations to take.
A value of 10 is low. Good accuracy is achieved with 20 or more
iterations.
Output, double *C, *S, the cosine and sine of the angle.
Local Parameters:
Local, double ANGLES(60) = arctan ( (1/2)^(0:59) );
Local, double KPROD(33). KPROD(j) = product ( 0 <= i <= j ) K(i)
where K(i) = 1 / sqrt ( 1 + (1/2)^(2i) ).
*/
fixedfloat_t angle;
fixedfloat_t c2;
int j;
int shiftamount;
fixedfloat_t s2;
int sigma;
int sign_factor;
fixedfloat_t theta;
/*
Shift angle to interval [-pi,pi].
*/
theta = angle_shift ( beta, -PI);
/*
Shift angle to interval [-pi/2,pi/2] and account for signs.
*/
if ( theta < - (PI >> 1) ) // theta < -0.5 * pi
theta = theta + PI;
sign_factor = -1;
else if ( (PI >> 1) < theta ) // 0.5 * pi < theta
theta = theta - PI;
sign_factor = -1;
else
sign_factor = 1;
/*
Initialize loop variables:
*/
*c = FROMINT(1);
*s = FROMINT(0);
angle = angles[0];
/*
Iterations
*/
for ( j = 1; j <= n; j++ )
sigma = (theta < 0) ? -1 : 1;
shiftamount = j - 1; // Index of the root to compute, where 3 = square root.
// sigma * poweroftwo, where poweroftwo = 0.5, 0.25...
c2 = *c - ((sigma < 0 ? -*s : *s) >> shiftamount);
s2 = ((sigma < 0 ? -*c : *c) >> shiftamount) + *s;
*c = c2;
*s = s2;
/*
Update the remaining angle.
*/
theta = theta - (sigma < 0 ? -angle : angle);
/*
Update the angle from table, or eventually by just dividing by two.
*/
angle = (ANGLES_LENGTH < j + 1) ? angle >> 1 : angles[j];
/*
Adjust length of output vector to be [cos(beta), sin(beta)]
KPROD is essentially constant after a certain point, so if N is
large, just take the last available value.
*/
#ifndef FIXED_N
if ( 0 < n )
*c = IFMUL(*c, kprod [ MIN ( n, KPROD_LENGTH ) - 1 ]);
*s = IFMUL(*s, kprod [ MIN( n, KPROD_LENGTH ) - 1 ]);
#else
*c = IFMUL(*c, KPROD_VALUE);
*s = IFMUL(*s, KPROD_VALUE);
#endif
/*
Adjust for possible sign change because angle was originally
not in quadrant 1 or 4.
*/
if (sign_factor < 0)
*c = -*c;
*s = -*s;
return;
/******************************************************************************/
【讨论】：
以上是关于基于微控制器的日出/设置算法实现的主要内容，如果未能解决你的问题，请参考以下文章