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

int step = 0;
int max_steps = 30000;

double h = 0.001;

double x = 0.0;
double y = 0.0;
double z = 0.0;

double sigma = 10.0;
double r = 28.0;
double b = 8.0 / 3.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;

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() {
  
  x = 0.001;
  y = 0.001;
  z = 0.001;

  FILE *fp = fopen("lorenz_out","w");

  metafl("XWIN");  

  disini();
  winmod("NOHOLD");  

  //  axis3d(2.0,2.0,2.0);
  graf3d(-50,50,-50,10,-50,50,-50,10,0,50,0,10);

  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);

  
    x += (m_0x + 2 * m_1x + 2 * m_2x + m_3x) / 6.0;
    y += (m_0y + 2 * m_1y + 2 * m_2y + m_3y) / 6.0;
    z += (m_0z + 2 * m_1z + 2 * m_2z + m_3z) / 6.0;

    sphe3d(x, y, z, 0.1, 10, 10);
    fprintf(fp,"%f %f %f\n", x, y, z);
  }
  
  fclose(fp);

  //symfil("XWIN","DELETE");
  disfin();
  
}