#include <stdio.h>
#include "ode.h"

/*

Likningene er:
dx             dy                  dz
-- = a(y-x)    -- = cx - y - xz    -- = xy - bz
dt             dt                  dt

La a=10, b=8/3, c=30.


 */

void lorentzEq(double t, const double *y, double *dydt) {
  static double a = 10., b=8./3., c=30.;
  dydt[0] = a*(y[1]-y[0]);
  dydt[1] = c*y[0] - y[1] - y[0]*y[2];
  dydt[2] = y[0]*y[1] - b*y[2];
}

void povVect(double a, double b, double c) {
  fprintf(stdout, "<%e, %e, %e>,\n", a, b, c);
}

int main(void) {
  odeSolver *lorentz = new odeSolver(3, lorentzEq);
  double maxt = 200.;
  int plotEvery = 10;
  int iterations = 0;
  int dots = 0;
  
  lorentz->setTime(0.0);
  lorentz->setStepSize(1e-4);
  lorentz->setAlgorithm(ODE_RK4);
  lorentz->setValue(0, 1);
  lorentz->setValue(1, 1);
  lorentz->setValue(2, 1);

  while (lorentz->getTime()<=maxt) {
    lorentz->calculate();
    if (lorentz->getTime()>10.) { //make sure orbit is attracted.
      if ((iterations%plotEvery)==0){
	//lorentz->display(); //for gnuplot.
	povVect(lorentz->getValue(0), lorentz->getValue(1), 
		lorentz->getValue(2));
	dots++; //count dots.
      }

    }
    iterations++;
  }

  fprintf(stderr, "%d points were calculated.\n", iterations);
  fprintf(stdout, "//%d points in array.\n", dots);

}