/* lorenz.c: skeletal ODE displayer */

#include <stdio.h>
#include <math.h>
#include "swgraph.h"

#define real double

// Put in your own 3D ODEs here, if you like
// These are the Lorenz equations
#define xdot(X,Y,Z) a*((Y)-(X))
#define ydot(X,Y,Z) ((c-(Z))*(X)-(Y))
#define zdot(X,Y,Z) ((X)*(Y)-b*(Z))

real x,y,z;       // State
real a,b,c;       // Parameters of simulated system */
real t,timestep;  // Time and the integrator time step

  // The following structs are convenient ways to store coordinates
  //
typedef enum {Square, Regular} FrameType;
typedef struct
  {
  real hlo,hhi,vlo,vhi;
  } FrameRect,*FrameRectPtr;
  // This describes the range of Lorenz variables
FrameRect TrajRange = {-45.,45.,0.,50.};
  // This is the piece of screen real estate that we plot to
FrameRect TrajFrame = {.2,.8,.4,1.};

main(argc,argv)
int argc;
char **argv;
  {
  int i = 0,num;

  x = -20.;   // Initial conditions
  y = -3.;
  z = 10.;
  num = 100000;  // Number of steps to integrate
  a   = 10.;     // Lorenz system parameters
  b   = 2.667;
  c   = 28.;
  timestep = .004;  // Integrator time step

     // Initialize the graphics driver, get a window to draw in
  initgraph(setdisplaydevice(argc,argv));
     // Clear the window
  clear();
     // Puts axes and draws labels
  template(Square,&TrajRange,&TrajFrame,"x + y",
	(char*)0,(char*)0,"z",(char*)0,(char*)0,"Lorenz Equations");
     // Set the coordinate range for the variables
  SetWindow(&TrajRange,&TrajFrame);

  for(i=1000;i;i--) rungekutta(); // Exorcise those transients
  move(x+y,z);
  pendown();  // Start drawing
  for(i=0;i<num;i++)
	{
  	rungekutta(); // Integrate the ODEs one timestep
	move(x+y,z);
	}
  getchar(); // Hit <RETURN> to exit the simulation
  penup();   // Stop drawing
  exitgraph(); // Clean up and exit
  }

  // This is a 4th order Runge-Kutta integrator
rungekutta()
  {
  real x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
  real x0,y0,z0;
  real x01,y01,z01,timetemp;

  timetemp = timestep/6.;

  x1 = xdot(x,y,z);
  y1 = ydot(x,y,z);
  z1 = zdot(x,y,z);
  x01 = x+timestep*x1*.5;
  y01 = y+timestep*y1*.5;
  z01 = z+timestep*z1*.5;
  x2 = xdot(x01,y01,z01);
  y2 = ydot(x01,y01,z01);
  z2 = zdot(x01,y01,z01);
  x01 = x+timestep*x2*.5;
  y01 = y+timestep*y2*.5;
  z01 = z+timestep*z2*.5;
  x3 = xdot(x01,y01,z01);
  y3 = ydot(x01,y01,z01);
  z3 = zdot(x01,y01,z01);
  x01 = x+timestep*x3;
  y01 = y+timestep*y3;
  z01 = z+timestep*z3;
  x4 = xdot(x01,y01,z01);
  y4 = ydot(x01,y01,z01);
  z4 = zdot(x01,y01,z01);
  x += (x1+2.*x2+2.*x3+x4)*timetemp;
  y += (y1+2.*y2+2.*y3+y4)*timetemp;
  z += (z1+2.*z2+2.*z3+z4)*timetemp;
  t += timestep;
  }

   // Establish drawing region and coordinate scaling
SetWindow(range,frame)
FrameRectPtr range,frame;
  {
   // Say which fractional part of the window you will draw in
  sqrwindow(frame->hlo,frame->vlo,frame->hhi-frame->hlo);
   // Associate floating point ranges with that area
  scale(range->hlo,range->hhi,range->vlo,range->vhi);
  }

   // Non-essential code: Draw axes and label with text
   //
template(type,range,frame,hlabel,hsub,hsup,vlabel,vsub,vsup,Title)
FrameRectPtr range,frame;
char *hlabel,*hsub,*hsup,*vlabel,*vsub,*vsup,*Title;
FrameType type;
  {
  char string[50];
  real xmin = range->hlo;
  real xmax = range->hhi;
  real ymin = range->vlo;
  real ymax = range->vhi;

  switch(type)
	{
  case Square:
	sqrwindow(frame->hlo,frame->vlo,frame->hhi-frame->hlo);
	break;
  case Regular:
	window(frame->hlo,frame->vlo,frame->hhi,frame->vhi);
	break;
	}
  scale(0.,1.,0.,1.);

  sprintf(string,"%g",ymin);
  move(-.06,0.); orientlabel(2,1,string,(char*)0,(char*)0);
  move(-.06,.5); orientlabel(2,1,vlabel,vsub,vsup);
  sprintf(string,"%g",ymax);
  move(-.06,1.); orientlabel(2,1,string,(char*)0,(char*)0);

  sprintf(string,"%g",xmin);
  move(0.,-.05); orientlabel(1,2,string,(char*)0,(char*)0);
  move(.5,-.05);orientlabel(1,2,hlabel,hsub,hsup);
  sprintf(string,"%g",xmax);
  move(1.,-.05); orientlabel(1,2,string,(char*)0,(char*)0);

  move(.5,1.1); orientlabel(1,0,Title,(char*)0,(char*)0);

  scale(xmin,xmax,ymin,ymax);
  border();
  axes(0.,0.,.1*(xmax-xmin),.1*(ymax-ymin),.01,.01);
  }