#include "library.h"

const double max_steps = 1000,
             win_width = 600, 
             win_height = 600,
             range_min_real = -2,
             range_max_real = 1,
             range_min_imag = -1.5, 
             range_max_imag = 1.5;

struct complex
{ double r, i; };

complex make_c(double r0, double i0)
{ complex x;
  x.r = r0;
  x.i = i0;
  return x; }

complex add(complex a, complex b)
{ complex x;
  x.r = a.r + b.r;
  x.i = a.i + b.i;
  return x; }

complex sub(complex a, complex b)
{ complex x;
  x.r = a.r - b.r;
  x.i = a.i - b.i;
  return x; }

complex mul(complex a, complex b)
{ complex x;
  x.r = a.r * b.r - a.i * b.i;
  x.i = a.r * b.i + a.i * b.r;
  return x; }

complex div(complex a, complex b)
{ double bottom = b.r * b.r + b.i * b.i;
  if (bottom == 0)
    cout << "Look out! Division by zero\n";
  complex x;
  x.r = (a.r * b.r + a.i * b.i) / bottom;
  x.i = (a.i * b.r - a.r * b.i) / bottom;
  return x; }

double magnitude(complex a)
{ return sqrt(a.r * a.r + a.i * a.i); }

ostream & operator<<(ostream & out, complex a)
{ out << a.r;
  if (a.i > 0)
    out << "+" << a.i << "i";
  else if (a.i < 0)
    out << a.i << "i";
  return out; }

double r_to_x(double r)
{ double p = (r - range_min_real) / (range_max_real - range_min_real);
  return p * win_width; }

double x_to_r(double x)
{ double p = x / win_width;
  return range_min_real + p * (range_max_real - range_min_real); }

double i_to_y(double i)
{ double p = (i - range_min_imag) / (range_max_imag - range_min_imag);
  return win_height - p * win_height; }

double y_to_i(double y)
{ double p = (win_height - y) / win_height;
  return range_min_imag + p * (range_max_imag - range_min_imag); }

struct point
{ double x, y; };

point make_p(double x, double y)
{ point p;
  p.x = x;
  p.y = y;
  return p; }

point complex_to_point(complex a)
{ point p;
  p.x = r_to_x(a.r);
  p.y = i_to_y(a.i);
  return p; }

complex point_to_complex(point p)
{ complex c;
  c.r = x_to_r(p.x);
  c.i = y_to_i(p.y);
  return c; }

void move_to(complex c)
{ point p = complex_to_point(c);
  move_to(p.x, p.y); }

void draw_to(complex c)
{ point p = complex_to_point(c);
  draw_to(p.x, p.y); }

void draw_point(complex c)
{ point q = complex_to_point(c);
  draw_point(q.x, q.y); }

void draw_axes()
{ set_pen_width(1);
  set_pen_color(color::grey);
  complex c;
  double axis_step = (range_max_real - range_min_real) / 10;
  for (c.r = range_min_real; c.r <= range_max_real; c.r += axis_step)
  { c.i = range_min_imag;
    move_to(c);
    c.i = range_max_imag;
    draw_to(c); }
  axis_step = (range_max_imag - range_min_imag) / 10;
  for (c.i = range_min_imag; c.i <= range_max_imag; c.i += axis_step)
  { c.r = range_min_real;
    move_to(c);
    c.r = range_max_real;
    draw_to(c); }
  set_pen_color(color::black);
  move_to(make_c(0.0, range_min_imag));
  draw_to(make_c(0.0, range_max_imag));
  move_to(make_c(range_min_real, 0.0));
  draw_to(make_c(range_max_real, 0.0)); }

point click_to_point()
{ wait_for_mouse_click();
  point p;
  p.x = get_click_x();
  p.y = get_click_y();
  return p; }

complex click_to_complex()
{ return point_to_complex(click_to_point()); }

int explore(complex c)
{ int steps = 0;
  complex z = c;
  move_to(z);
  while (true)
  { if (steps == max_steps)
      return -1;
    if (magnitude(z) > 5)
      return steps;
    z = add(mul(z, z), c);
    draw_to(z);
    steps += 1; } }

void main()
{ make_window(win_width, win_height);
  draw_axes();
  set_pen_width(1);
  complex c;
  bool first_time = true;
  while (true)
  { complex newc = click_to_complex();
    if (! first_time)
    { set_pen_color(color::light_grey);
      explore(c); }
    c = newc;
    cout << "c = " << c << "\n";
    set_pen_color(color::black);
    int s = explore(c);
    if (s < 0)
      cout << "stable\n";
    else
      cout << s << " steps to explode\n";
    first_time = false;
    draw_axes(); } }