Microcontroller based sunrise/set algorithm implementation

Question

Two questions:

I am trying to implement an algorithm to determine sunrise and sunset times, with latitude, longitude and date as inputs on a microcontroller. I know there are algorithms to do this, but I wasn't able to find any detail on how to put this on a microcontroller. Is this possible?

Second, if I am to implement this, how much memory would I need approximately and what kind of controller would I have to use?

Please let me know.

Answer 1

Sunrise/sunset timings can be computed with algorithms available on the Internet, which generally require trigonometry functions, which in turn implies support for floating-point.

If there is no FPU, then you will have to use a compiler library. Or lookup tables, an approximation or even CORDIC:

If you:

• Do not have a FPU or do not want to use floating-point.
• Do not want to use compiler libraries for trigonometry.
• Need accuracy that cannot be satisfied with an approximation.
• Cannot use fixed tables because the timezone may change.

Then I would use CORDIC, which seems to lend itself naturally to fixed-point numbers.

This example uses the algorithm from Ed Williams's website, which also has an example of this algorithm in action listed. I chose to use fixed-point numbers in Q10.21 format, as it seems to give the highest-possible accuracy without arithmetic overflows. Thus the target should be 32-bit. It is possible to change the width, but the accuracy will be affected.

I am using custom versions of the CORDIC functions by John Burkardt, which were modified to use fixed-point arithmetic and to reduce the number of multiplication operations. He mentioned in his comments that he referred to the MATLAB example on Wikipedia, making it useful for understanding how the algorithm works.

I have done a few rounds of testing with a few timezones (+08:00, +04:00, -03:00) against the NOAA sunrise/sunset calculator, and the accuracy of this example seems to be ±1 minute. While the timezone is specfied as an argument, the time (and DST setting) used in the example is taken from your OS.

Some notes regarding the fixed-point operations:

• Addition and substraction can be done using regular +/- operators.
• Multiplication and division must be done with the IFMUL and IFDIV macros respectively.
• Modulo division can be done with %, but I added a IFMOD macro for consistency with the division operation.

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://stackoverflow.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;
}

fixedpoint.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;
}
/******************************************************************************/

Answer 2

Well on MCU there is not much memory and also FPU is not always present so you need to do much on your own. There are few ways how to compute this.

1. Sunrise calendar

   This is the simplest but require to store the calendar data (can be in program memory so no need for RAM). This approach is limited to single geo location so if the location change with time then it is not usable. For more info see:

   • Calculate whether it is close to dawn/dusk based on sunrise/sunset?

2. Kepler's equation

   This will work for a large amount of time (for example 1000 years) when coded right. But this approach needs FPU computations. If not present you need to code it yourself (things like: double,sin,cos,+,-,*,/)

   Kepler's equation will give you actual relative position of Earth and Sun in Heliocentric coordinate system. If you add Earth's daily rotation and local NEH (north east height/alt) transformations you will obtain Suns position in azimutal coordinates for any geolocation and time. Then it is simply matter of checking elevation against sunrise/sunset geometric limits. For more info see:

   • 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

   This approach need RAM to store actual coordinates (3x double per position + one temp), transform 3x double 4x4 transform matrices (earth,NEH,temp) and few iteration variables. I estimate 512 Byte would be enough for this.

3. Direct equation

   You need to google for this. There are out there some approximation equations that will give you directly the azimutal coordinates of Sun for any geolocation. It is a bit simpler then #2 but needs more advanced goniometric functions like arctan,acos and not sure if not also some hyperbolic functions. In case you do not have FPU support then this approach may be hard to implement in desired precision.

[Notes]

To improve accuracy you need to convert geometric sun's position to visual. There are tables how much to shift elevation angle ...