////////////////////////////////////////////////////////////////////////////
// Model LORENZ.CPP               SIMLIB/C++
//
// Lorenzova rovnice:
//
//      dx1/dt = sigma * (x2 - x1)
//      dx2/dt = (1 + lambda - x3) * x1 - x2
//      dx3/dt = x1 * x2 - b * x3
//
// kde:  sigma=10, lambda=24 a b=2  jsou parametry
// počáteční podmínky:  xi(0) = 1.0
//
// Zdroj: SimPack
//

#include "simlib.h"
#include <stdlib.h>

const double StepPrn = 0.01;    // krok výpisu

struct Lorenz {
  Integrator x1, x2, x3;
  Lorenz(double sigma, double lambda, double b) :
    x1(sigma*(x2 - x1), 1),           // dx1/dt = sigma * (x2 - x1)
    x2((1 + lambda - x3)*x1 - x2, 1), // dx2/dt = (1 + lambda - x3) * x1 - x2
    x3(x1*x2 - b*x3, 1) {}            // dx3/dt = x1 * x2 - b * x3
};

Lorenz L(10, 24, 2);            // model systému

// výstupy:
void Sample() { 
  Print("%6.2f %g %g %g\n", T.Value(), L.x1.Value(), L.x2.Value(), L.x3.Value()); 
}
Sampler S(Sample, StepPrn);

int main(int argc, char *argv[]) {  // experiment
  double maxtime = 30.0;
  if(argc>1) 
    maxtime = atof(argv[1]);        // parametr
  if(maxtime<1.0) {
     _Print("\nPoužití:  %s  [maxtime>1,default=30] \n\n", argv[0]);
     return 1;
  }
  SetOutput("lorenz.dat");
  _Print("# LORENZ - Lorenzova rovnice (maxtime=%g) \n", maxtime);
  Print("# Time x1 x2 \n");
  Init(0, maxtime);             // inicializace
  SetAccuracy(1e-3);            // požadovaná přesnost
  Run();                        // simulace
  return 0;
}

// konec