// originally by Andy Sloane
// https://www.a1k0n.net/2011/07/20/donut-math.html
// https://www.a1k0n.net/2021/01/13/optimizing-donut.html

#include <FastLED.h>

#define WIDTH 16
#define HEIGHT 16
#define NUM_LEDS ((WIDTH) * (HEIGHT))

CRGB leds[NUM_LEDS + 1];

uint16_t XY(int8_t x, int8_t y) {
  if (x >= WIDTH || x < 0)  return NUM_LEDS;
  if (y >= HEIGHT || y < 0) return NUM_LEDS;
  if (y & 1)
    return (y + 1) * WIDTH - 1 - x;
  else
    return y * WIDTH + x;
}

void setup() {
  Serial.begin(115200);
  FastLED.addLeds<WS2812B, 3, GRB>(leds, NUM_LEDS);
  FastLED.setCorrection(UncorrectedColor);
  FastLED.setTemperature(UncorrectedTemperature);
  FastLED.setDither(DISABLE_DITHER);
}

#define R(mul,shift,x,y) \
  _=x; \
  x -= mul*y>>shift; \
  y += mul*_>>shift; \
  _ = 3145728-x*x-y*y>>11; \
  x = x*_>>10; \
  y = y*_>>10;

void loop() {
  static int frame; frame++;
  uint32_t ms = millis();
  FastLED.clear();

  int8_t z[NUM_LEDS + 1] = {};
  static int32_t sA = 1024, cA = 0, sB = 1024, cB = 0, _;

  int sj=0, cj=1024;
  for (int32_t j = 0; j < 90; j++) {
    int32_t si = 0, ci = 1024;  // sine and cosine of angle i
    for (int32_t i = 0; i < 324; i++) {
      int32_t R1 = 1, R2 = 2048, K2 = 5120*1024;

      int32_t x0 = R1*cj + R2,
          x1 = ci*x0 >> 10,
          x2 = cA*sj >> 10,
          x3 = si*x0 >> 10,
          x4 = R1*x2 - (sA*x3 >> 10),
          x5 = sA*sj >> 10,
          x6 = K2 + R1*1024*x5 + cA*x3,
          x7 = cj*si >> 10,
          x = WIDTH / 2 + WIDTH / 3 * (cB*x1 - sB*x4)/x6,
          y = HEIGHT / 2 +  HEIGHT / 3 *(cB*x4 + sB*x1)/x6,
          N = (-cA*x7 - cB*((-sA*x7>>10) + x2) - ci*(cj*sB >> 10) >> 10) - x5 >> 7;

      uint16_t xy = XY(x, y);
      int32_t o = x + 80 * y;
      int8_t zz = (x6-K2)>>15;
      if (xy < NUM_LEDS && zz < z[xy]) {
        z[xy] = zz;
        if (N > 0)
          leds[xy] = CRGB(N, N, N);
      }
      R(5, 8, ci, si)  // rotate i
    }
    R(9, 7, cj, sj)  // rotate j
  }
  R(5, 7, cA, sA);
  R(5, 8, cB, sB);
  FastLED.show();
}