Dev - C++ - Ayuda con un codigo de Lorentz

#include <cmath>

#include <cstdlib>

#include <fstream>

#include <iostream>

using namespace std;

#include "vector.h"

using namespace cpl;

double sigma = 10;                 //  Lorenz model constants in textbook

double b = 8.0 / 3.0;

double r = 25;

Vector f(Vector txyz) {            //  The Lorenz equations

    double t = txyz[0];

    double x = txyz[1];

    double y = txyz[2];

    double z = txyz[3];

    Vector f(4);

    f[0] = 1;

    f[1] = - sigma * x + sigma * y;

    f[2] = - x * z + r * x - y;

    f[3] = x * y - b * z;

    return f;

}

void RK4Step(                      //  4th order Runge-Kutta

    Vector& y,                     //  extended solution vector

    double h)                      //  step size

{

    Vector k1 = h * f(y);

    Vector k2 = h * f(y + 0.5 * k1);

    Vector k3 = h * f(y + 0.5 * k2);

    Vector k4 = h * f(y + k3);

    y += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;

}

Vector txyz(4);                    //  global variable to hold t,x,y,z

void initialize() {                //  initial conditions in textbook

    txyz[0] = 0.0;

    txyz[1] = 0.0;

    txyz[2] = 1.0;

    txyz[3] = 0.0;

}

double dt = 0.001;                 //  time step for integration

void findNextCrossing() {          //  find next Poincare section point

    Vector txyzOld = txyz;         //  save the old point

    while (true) {

        RK4Step(txyz, dt);         //  take a step

        // check whether y changes sign, i.e., crosses y = 0

        if (txyz[2] * txyzOld[2] < 0)

            break;                 // stop stepping

        else

            txyzOld = txyz;

    }

    //  use linear interpolation like in projectile.cpp to find intersection

    double r = txyzOld[2] / txyz[2];

    for (int i = 0; i < 4; i++)

        txyz[i] = (txyzOld[i] + r * txyz[i]) / (r + 1);

}

int main(int argc, char *argv[]) {

    cout << " Trajectory and Poincare Section for the Lorenz Attractor\n"

         << " using 4th order Runge-Kutta with time step dt = " << dt << "\n"

         << " sigma = " << sigma << ", b = " << b << ", r = " << r << endl;

    initialize();

    cout << " initial conditions: x = " << txyz[1] << "\t"

         << ", y = "  << txyz[2] << "\t" << ", z = "  << txyz[3] << endl;

    //  transient trajectory

    double t = 50;

    string fileName = "transient.data";

    cout << " Integrating to time t = " << t << "\n"

         << " trajectory in file " << fileName << endl;

    ofstream dataFile(fileName.c_str());

    dataFile << txyz[0] << "\t" << txyz[1] << "\t"

             << txyz[2] << "\t" << txyz[3] << "\n";

    int step = 0;

    int skip = 5;

    while (txyz[0] < t) {

        RK4Step(txyz, dt);

        if (++step % skip != 0)    // record every skip steps

           continue;

        dataFile << txyz[0] << "\t" << txyz[1] << "\t"

                 << txyz[2] << "\t" << txyz[3] << "\n";

    }

    dataFile.close();

    //  Poincare section

    fileName = "section.data";

    int points = 1000;

    cout << " Finding " << points << " Poincare section points at y = 0\n"

         << " section data in file " << fileName << endl;

    dataFile.open(fileName.c_str());

    for (int point = 0; point < points; point++) {

        findNextCrossing();

        dataFile << txyz[0] << "\t" << txyz[1] << "\t"

                 << txyz[2] << "\t" << txyz[3] << "\n";

    }

}

#include <iostream>

#include <math.h>

#ifdef __APPLE__

#include <OpenGL/OpenGL.h>

#include <GLUT/glut.h>

#else

#include <glut.h>

#endif

void myinit() {

	// set the background color

	glClearColor(0.0f, 0.0f, 0.0f, 1.0f);

	// set the foreground (pen) color

	glColor4f(1.0f, 1.0f, 1.0f, 1.0f);

}

void mydisplay() {

	// clear the screen

	glClear(GL_COLOR_BUFFER_BIT);

	// swap the buffers

	glutSwapBuffers();

}

void mykey(unsigned char mychar, int x, int y) {

	// exit the program when the Esc key is pressed

	if (mychar == 27) {

		exit(0);

	}

}

int main (int argc, char **argv) {

	// initialize GLUT

	glutInit(&amp argc, argv);

	// set up our display mode for color with alpha and double buffering

	glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);

	// create a 400px x 400px window

	glutInitWindowSize(400, 400);

	glutCreateWindow("Strange Attractors in C++ and OpenGL Tutorial");

	// register our callback functions

	glutDisplayFunc(mydisplay);

	glutKeyboardFunc(mykey);

	// call our initialization function

	myinit();

	// start the program

	glutMainLoop();

    return 0;

}