//************************************************************************ 
// Lorenz model with the classical RK4 method
//  x'=sigma(y-x)
//  y'=rx-y-xz,
//  z'=xy-bz,
//  x=x(t), y=y(t), z=z(t),  T>t>0, x(0)=x_0, y(0)=y_0, z(0)=z_0
//************************************************************************
/* 
compile with	
gcc -o lorenz lorenz.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/

#include "/usr/local/dislin/dislin.h"
#include <stdio.h>
#include <math.h>


int step = 0.;
int max_steps = 20000.;   // 2000.0 2000.0 5.000 10000.0

double h = 0.005;
// initial conditions
double x = 0.0;
double y = 1.0;
double z = 0.0;
// parameters
double sigma = 10.0;
double b = 8.0 / 3.0;
//------- control parameter -------------
double r = 26.0; // 0.5 3.0 16.0 26.0
//---------------------------------------

double m_0x = 0.0;
double m_1x = 0.0;
double m_2x = 0.0;
double m_3x = 0.0;
double m_0y = 0.0;
double m_1y = 0.0;
double m_2y = 0.0;
double m_3y = 0.0;
double m_0z = 0.0;
double m_1z = 0.0;
double m_2z = 0.0;
double m_3z = 0.0;

// define r.h.s.
double funct_x(double x, double y, double z) {
  return -sigma*x + sigma * y;
}

double funct_y(double x, double y, double z) {
  return r*x-y-x*z;
}

double funct_z(double x, double y, double z) {
  return -b*z+x*y;
}

main() {
      
  double zp = funct_z(x,y,z);
  double zp_alt=zp;
  double z_min_alt = 0.0;
 
  winsiz (800, 800);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini ();
  pagera ();
  hwfont ();
  texmod ("on");


  // 3d plot (x,y,z)
  axspos(140,2000);
  axslen(1500,1500);
  titlin ("3D Phase space", 4);
  //graf3d(0,1,0,1,0,1,0,1,0,0.2,0,0.1); //for stable origin
  // graf3d(0,3,0,1,0,3,0,1,0,3,0,1); //for stable c+, c-
  //graf3d(-15,15,-15,5,-20,20,-20,10,0,30,0,15); //for unstable limit circle
 graf3d(-20,20,-20,10,-50,50,-50,25,0,50,0,25); //for strange attractor
  title();
  endgrf();

  // Lozenz map
  axspos(1700,1800);
  axslen(800,800);
  titlin ("Lorenz map", 4);
  graf(25.0,45.0,25.0,10.0,25.0,45.0,25.0,10.0);
  title();
  endgrf();

    
  for ( step = 0; step < max_steps; step++ ) {

    m_0x = h * funct_x(x, y, z);
    m_0y = h * funct_y(x, y, z);
    m_0z = h * funct_z(x, y, z);
    
    m_1x = h * funct_x(x + 0.5 * m_0x, y + 0.5 * m_0y, z + 0.5 * m_0z);
    m_1y = h * funct_y(x + 0.5 * m_0x, y + 0.5 * m_0y, z + 0.5 * m_0z);
    m_1z = h * funct_z(x + 0.5 * m_0x, y + 0.5 * m_0y, z + 0.5 * m_0z);

    m_2x = h * funct_x(x + 0.5 * m_1x, y + 0.5 * m_1y, z + 0.5 * m_1z);
    m_2y = h * funct_y(x + 0.5 * m_1x, y + 0.5 * m_1y, z + 0.5 * m_1z);
    m_2z = h * funct_z(x + 0.5 * m_1x, y + 0.5 * m_1y, z + 0.5 * m_1z);

    m_3x = h * funct_x(x + m_2x, y + m_2y, z + m_2z);
    m_3y = h * funct_y(x + m_2x, y + m_2y, z + m_2z);
    m_3z = h * funct_z(x + m_2x, y + m_2y, z + m_2z);

    zp=funct_z(x,y,z);
    zp_alt = zp;
    x += (m_0x + 2.0 * m_1x + 2.0 * m_2x + m_3x) / 6.0;
    y += (m_0y + 2.0 * m_1y + 2.0 * m_2y + m_3y) / 6.0;
    z += (m_0z + 2.0 * m_1z + 2.0 * m_2z + m_3z) / 6.0;
    zp=funct_z(x,y,z);
 //--------------------------------------------------------  
    setrgb(0,0,0);
    axspos(140,2000);
    axslen(1500,1500);
    //graf3d(0,1,0,1,0,1,0,1,0,0.2,0,0.1); //for stable origin
   // graf3d(0,3,0,1,0,3,0,1,0,3,0,1); //for stable c+, c-
    //graf3d(-15,15,-15,5,-20,20,-20,10,0,30,0,15); //for unstable limit circle
    graf3d(-20,20,-20,10,-50,50,-50,25,0,50,0,25); //for strange attractor
    setrgb(0,0,1);
    sphe3d(x, y, z, 0.15, 10.0, 10.0); // 0.005 0.01 0.1 0.15
    endgrf();

// find zz=min(z) and plot the Lorenz map zz_(n+1)(zz_n)
if ( (zp_alt > 0.0) && (zp < 0.0) ) {
	if (z_min_alt) {
          setrgb(0,0,0);
          axspos(1700,1800);
          axslen(800,800);
          graf(25.0,45.0,25.0,10.0,25.0,45.0,25.0,10.0);
	  setrgb(1,0,0);
          hsymbl(12);
	  rlsymb(21,z_min_alt,z);
	}
	z_min_alt = z;
      }
  endgrf();

  }

  disfin();
  
}