/*
    Code to calculate Lorentz attractor.
    https://en.wikipedia.org/wiki/Lorenz_system

 Algorithm:
 
    A 3d point (x,y,z) (representing a simplfied model of atmospheric convextion) evolves according to
    
      dx/dt = sigma (y-x)
      dy/dt = x(rho-z)-y
      dz/dt = x*y-beta z
      
    where Lorenz choose sigma = 10, rho=28   beta=8/3
 
  Euler integration:
  
     dx/dt = f(x,t)  where the rhs is any function.
     
    We can approximate dx/dt = [ x(t+dt)-x(t) ] / dt   for "small" dt
    
    Hence x(t+dt) = dt*(dx/dt) + x(t) 
    
    But we also have dx/dt = f(x,t) 
    
    hence  x(t+dt) = dt*f(x,t) + x(t)
    
    which is an recursion equation that provides approximate solution to
    the differential equation. 
     
      
  Matlab to visualize:
  data=importdata("lorenz.dat")
  plot3(data(:,1),data(:,2),data(:,3),'-');
  
*/


#include <fstream>
#include <iostream>
using namespace std;

int main()
{
   const double beta=8.0/3.0;
   const double rho=28.0;
   const double sigma=10.0;

   const double tmax=100;
   const double dt=0.01;

   double t=0,x=1,y=1,z=1;

   ofstream fo("lorenz.dat");

   while (t<tmax)
      {
         // Calculate change for x,y,z at this time
         double dx=sigma*(y-x);
         double dy=x*(rho-z)-y;
         double dz=x*y-beta*z;
      
         // Update coordinates using Euler integration
         x+=dx*dt;
         y+=dy*dt;
         z+=dz*dt;
      
         fo << x << " " << y << " " << z << endl;
      
         t+=dt;
      }

   fo.close();
}