int N_BESSEL = 5;
int N_ECCENTRIC_ANOMALY = 5;
float EARTH_ECCENTRICITY = 0.01671;
float EARTH_SEMI_MAJOR_AXIS = 1.0;
float MARS_ECCENTRICITY = 0.0934;
float MARS_SEMI_MAJOR_AXIS = 1.52400;

float earth_angle = 0;
float mars_angle = 0;
float earth_velocity = 0.005;
PVector sun;
float mars_circle_extent = 150;
float planet_extent = 10;

void setup() {
  size(1000, 1000);
  pixelDensity(1);
  sun = new PVector(width * 0.5, height * 0.5);
}

void draw() {
  if (keyPressed == true) {
    return;
  }

  background(255, 255, 255);
  noStroke();

  // Draw the Sun
  noStroke();
  fill(0xFF, 0xA5, 0x00);
  circle(sun.x, sun.y, planet_extent);

  float radius_scale = width * 0.3;


  // Draw Earth
  float t_earth = true_anomaly(eccentric_anomaly(earth_angle, EARTH_ECCENTRICITY, N_ECCENTRIC_ANOMALY), EARTH_ECCENTRICITY);
  float r_earth = calc_orbit(t_earth, EARTH_ECCENTRICITY, EARTH_SEMI_MAJOR_AXIS);
  PVector earth = new PVector(
    sun.x + radius_scale * r_earth * sin(t_earth),
    sun.y + radius_scale * r_earth * cos(t_earth)
    );
  noStroke();
  fill(0, 0, 255);
  circle(earth.x, earth.y, planet_extent);
  earth_angle += PI * earth_velocity;

  // Draw Mars
  float t_mars = true_anomaly(eccentric_anomaly(mars_angle, MARS_ECCENTRICITY, N_ECCENTRIC_ANOMALY), MARS_ECCENTRICITY);
  float r_mars = calc_orbit(t_mars, MARS_ECCENTRICITY, MARS_SEMI_MAJOR_AXIS);
  PVector mars = new PVector(
    sun.x + radius_scale * r_mars * sin(t_mars),
    sun.y + radius_scale * r_mars * cos(t_mars)
    );
  noStroke();
  fill(255, 0, 0);
  circle(mars.x, mars.y, planet_extent);
  mars_angle += PI * earth_velocity * 365.2422 / 686.980;

  // Draw the direction of view from Earth
  PVector sight = new PVector(earth.x - mars.x, earth.y - mars.y);
  sight.normalize();
  stroke(0, 255, 0);
  line(
    mars.x - sight.x * 1000, mars.y - sight.y * 1000,
    mars.x + sight.x * 1000, mars.y + sight.y * 1000
    );

  // Draw the orbit of Mars based on Earth
  PVector mars_circle_pos = new PVector(width * 0.9, height * 0.1);
  noFill();
  stroke(0, 0, 0);
  circle(mars_circle_pos.x, mars_circle_pos.y, mars_circle_extent);
  noStroke();
  fill(0, 0, 255);
  circle(mars_circle_pos.x, mars_circle_pos.y, planet_extent);
  noStroke();
  fill(255, 0, 0);
  circle(
    mars_circle_pos.x + mars_circle_extent * 0.5 * sight.x,
    mars_circle_pos.y + mars_circle_extent * 0.5 * sight.y,
    planet_extent
    );
}

float bessel(float x, float alpha, int N) {
  if (N < 0) {
    throw new RuntimeException("n must be non-negative");
  }
  if (alpha < 0) {
    throw new RuntimeException("alpha must be non-negative");
  }

  float sum = 0;
  for (int m = 0; m <= N; m++) {
    float term =
      (pow(-1, m) * pow(x / 2, alpha + 2 * m)) /
      (factorial(m) * factorial(alpha + m));
    sum += term;
  }
  return sum;
}

float factorial(float n) {
  if (n < 0) {
    throw new RuntimeException("n must be non-negative");
  }
  if (n == 0 || n == 1) {
    return 1;
  }
  float result = 1;
  for (int i = 2; i <= n; i++) {
    result *= i;
  }
  return result;
}

// M: mean anomaly
// e: eccentricity
// N: number of terms in the series
float eccentric_anomaly(float M, float e, int N) {
  if (e < 0 || e >= 1) {
    throw new RuntimeException("e must be in the range [0, 1)");
  }
  if (N < 0) {
    throw new RuntimeException("n must be non-negative");
  }

  float s = 0;
  for (int n = 1; n <= N; n++) {
    s = (bessel(n * e, n, N_BESSEL) * sin(n * M)) / n;
  }
  return M + 2 * s;
}

// E: eccentric anomaly
// e: eccentricity
float true_anomaly(float E, float e) {
  if (e < 0 || e >= 1) {
    throw new RuntimeException("e must be in the range [0, 1)");
  }

  float beta = e / (1 + sqrt(1 - e * e));
  return E + 2 * atan((beta * sin(E)) / (1 - beta * cos(E)));
}

// a: semi-major axis
// e: eccentricity
// theta: true anomaly
float calc_orbit(float theta, float e, float a) {
  if (e < 0 || e >= 1) {
    throw new RuntimeException("e must be in the range [0, 1)");
  }
  if (theta < 0) {
    throw new RuntimeException("theta must be non-negative");
  }
  if (a <= 0) {
    throw new RuntimeException("a must be positive");
  }

  float r = (a * (1 - e * e)) / (1 + e * cos(theta));
  return r;
}