/* hopf_ode.c  bifurcation, parameter p */
/* z = (x,y) */
/* dz/dt = g(z,p) = gx(x,y,p) = px-y-x(x^2+y^2) */
/*                = gy(x,y,p) = x+py-x(x^2+y^2) */
/* solve dz/dt for various values of p, e.g. -0.5, +0.5 */

#include <stdio.h>
#include <stdlib.h>
#include <GL/glut.h>

double gx(double x, double y, double p)
{
  return p*x - y - x*(x*x+y*y);
}

double gy(double x, double y, double p)
{
  return x + p*y - y*(x*x+y*y);
}

static double x, y, p, dt;
static int np, n, ninc;
static double xp[1000], yp[1000];

static void display(void)
{
  int i;

  glClear(GL_COLOR_BUFFER_BIT);
  glLoadIdentity();

  glPointSize(2.0);
  glColor3f(1.0, 0.0, 0.0);
  glBegin(GL_POINTS);
  for(i=0; i<np; i++)
  {
      glVertex3f(xp[i], yp[i], 0.0);
  }
  glEnd();
  glFlush(); 
} /* end display */

static void physics(void)
{
  double dx, dy, dr;
  char c;

  dx = gx(x, y, p);
  dy = gy(x, y, p);
  dr = sqrt(dx*dx+dy*dy);
  x = x + dx*dt/dr;
  y = y + dy*dt/dr;
  ninc++;
  if(ninc>0)
  {
    xp[np] = x;
    yp[np] = y;
    ninc = 0;
  }
  if(np>=999)
  {
    scanf("%c",c);
    exit(0);
  }
  np++;
  glutPostRedisplay();
} /* end physics */

void init()
{
  x = 0.4;
  y = 0.4;
  p = -0.1;
  n = 1000;
  dt = 0.01;
  xp[0] = x;
  yp[0] = y;
  np = 1;
  ninc = 0;
  printf("initial x=%g, y=%g, p=%g, dt=%g, steps=%d \n", x, y, p, dt, n);
} /* end init */

static void myReshape(int w, int h)
{
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(-0.5, 0.5, -0.5, 0.5, -1.0, 1.0);
  glMatrixMode(GL_MODELVIEW);
} /* end myReshape */

int main(int argc, char * argv[])
{
  printf("Hopf_ode.c running \n");
  init();
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGB);
  glutInitWindowSize(500, 500);
  glutCreateWindow(argv[0]);
  glutReshapeFunc(myReshape);
  glutDisplayFunc(display);
  glutIdleFunc(physics);
  glClearColor(1.0,1.0,1.0,1.0);
  glutMainLoop();
  return 0;
} /* end hopf_ode */