Dev - C++ - Problemas Con Almacenar

#include <math.h>

#include <cmath>

#include <cstdlib>

#include <fstream>

#include <iostream>

using namespace std;

#include "vector.hpp"

using namespace cpl;

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

double b = 8.0 / 3.0;

double r = 28;

Vector f(Vector txyz) {

	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 PrimercruceenB() {

	int B;

	double w = sqrt(b * (r - 1));

	while (txyz[1] < w){ // Cruce con plano  B en x

    	RK4Step(txyz, dt);

    }

    	// Aqui ya tengo el primer valor de x cruzando al plano B

    	if (txyz[2] >= w) { // Cruce con plano B en y

    		B = 1;

    		cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << B << endl;

    	}

    	else {

    		while ( txyz[2] < w) {

    			RK4Step(txyz, dt);

    		}

    		if (txyz[1] >= w) {

    			B = 1;

    			cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << B << endl;

    		}

    	}

}

void despuesdecruceplanoB() {

	Vector txyzP = txyz; //vector presente para comprarar con el siguiente

    		while (txyz[1] >= txyzP[1]) {

    			txyzP = txyz;

    			RK4Step(txyz, dt);

    		}

    		cout <<"x se devuelve del plano B  " << "x = " << txyz[1] << "y = " << txyz[2] << "z = " << txyz[3] << endl; // x comienza devolverse

}

void despuesdecruceplanoA () {

	Vector txyzPra = txyz; //vector presente para comprarar con el siguiente

    				while (txyz[1] <= txyzPra[1]) {

    					txyzPra = txyz;

    					RK4Step(txyz, dt);

					}

					cout <<"x se devuelve del plano A  " << "x = " << txyz[1] << "y = " << txyz[2] << "z = " << txyz[3] << endl; // x comienza devolverse

}

void xcambiadesigno() {

	Vector txyzV = txyz;

    	while (true) { // Decision atractor se queda o no en el mismo plano

    		RK4Step(txyz, dt);

    		if (txyz[1] * txyzV[1] < 0) { // Si x cambia de signo

    			break;

    			}

    		else {

    			txyzV = txyz;

    		}

    	}

}

void PrimercruceenA() {

	int B;

	double w = sqrt(b * (r - 1));

	while (txyz[1] > -w) { // Busco cruce con plano A en x

    				RK4Step(txyz, dt);

    		}

    	if (txyz[2] <= -w) {

    		B = 0;

    		cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << B << endl;

    	}

    	else {

    		while (txyz[2] > -w) {

    			RK4Step(txyz, dt);

    		}

    		if (txyz[1] <= -w) {

    			B = 0;

    			cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << B << endl;

    		}

    	}

}

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;*/

    int skip = 50;

    double dt = 0.001;

    double w = sqrt(b * (r - 1));

    int B;

    int points = 1000;

    double t = 250;

    while (txyz[0] < 50) {

    	RK4Step(txyz, dt);

    }

    while (txyz[0] < t){

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

    	 	Inicio:

    	 		if (txyz[0] >= t) {

    	 			break;

    	 		}

    		PrimercruceenB();

    		// Buscar cuando x comienza a devolverse del plano B

    		despuesdecruceplanoB(); // x comienza a devolverse

    		// Decidir si se queda o no en el mismo plano

    		Vector txyzPr = txyz; //vector presente para comprarar con el siguiente

    		while (txyz[1] <= txyzPr[1]) { // Se sale del while cuando se regresa al plano B

    			txyzPr = txyz;

    			RK4Step(txyz, dt);

    			if (txyz[1] * txyzPr[1] < 0) {

    				cout <<"x cambio de signo" << endl; // x se dirige al siguiente plano 

    				PrimercruceenA();

    				goto Parte2;

				} // Cruce con el plano A

			}

			if (txyz[1] > 0) {

				goto Inicio;

			}

			Parte2:

				if (txyz[0] >= t) {

					break;

				}

			if (txyz[1] < 0) {

				despuesdecruceplanoA ();

				// Decidir si se queda o no en el mismo plano

				Vector txyzPraa = txyz;

				while (txyz[1] >= txyzPraa[1]) {

					txyzPraa = txyz;

					RK4Step(txyz, dt);

					if (txyz[1] * txyzPraa[1] < 0) {

    				cout <<"x cambio de signo" << endl; // x se dirige al siguiente plano 

    				goto Inicio;

    				break;

				} // Cruce con el plano A

				}

				if (txyz[1] < 0) {

					PrimercruceenA();

				}

			}

			goto Parte2;

		 }

}