Dev - C++ - #include vector.hpp

#include <math.h>

#include <cmath>

#include <cstdlib>

#include <fstream>

#include <iostream>

#include <iomanip>

#include <iterator>

#include <sstream>

#include <string>

#include <sstream>

#include <string>

#include <vector>

#include <algorithm>

#include <math.h>

#include <cmath>

#include <cstdlib>

#include <fstream>

#include <iostream>

#include <iomanip>

#include<iterator>

#include <stdio.h>

#include <cctype>

using namespace std;

#include "vector.hpp"

using namespace cpl;

//#include "Vector.hpp"

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

double b = 8.0 / 3.0;

double r = 28;

int t = 400;

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

	char B;

	char Bi;

	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

                Bi = '3';

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

        }

        else {

               cout << " y < w" << endl;

                while ( txyz[2] < w) {

                        RK4Step(txyz, dt);

                }

                if (txyz[1] >= w) {

                        Bi = '3';

                       cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << 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 PrimercruceenA() {

	char Bi;

	char A;

	Vector u;

	//double g;

	//int h;

	//const int tamano_entrada = 16;

	//double delta;

	//double in[tamano_entrada] = {20.5722, 16.9518, 11.4379, 17.7182, 13.1959, 16.9518, 16.6336, 13.1959, 20.5722, 20.5722, 16.9518, 16.6336, 13.1959, 20.5722, 16.9518, 16.6336};

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

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

		u = txyz;

		RK4Step(txyz, dt);

	}

	/*cout << u[3] << " " << txyz[3] << endl;

				delta = in[h] - u[3];

           		 h = h + 1;

				u[3] = u[3] + delta;*/

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

		Bi = '2';

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

	}

	else {

		cout << " y es > w" << endl;

		while (txyz[2] > -w) {

			RK4Step(txyz, dt);

		}

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

			Bi = '2';

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

		}

	}

}

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

}

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

	int n = 6;

	int i, c;

	int h = 0;

	int j = 0;

	const int tamano_entrada = 14;

	char sim[tamano_entrada];

	double delta;

	double equis[tamano_entrada + 1];

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

	double out[tamano_entrada-1];

	//Vector u = txyz;

	char Bi;

	char A;

	Vector u;

	//double in[tamano_entrada];

	if (n==6) {

		//double in[] = {20.5722, 20.5722, 20.5722, 16.9518, 16.6336, 13.1959, 20.5722, 16.9518, 16.6336, 13.1959, 16.9518, 11.4379, 17.7182, 13.1959, 16.9518, 11.4379, 20.7618, 17.7182, 13.1959};

		 double in[] = {17.1772, 12.7452, 19.7595, 18.2817, 15.8382, 14.2816, 17.6427, 14.2049, 17.6369, 14.549, 16.9111, 10.6945, 22.135, 22.135};//, 14.2816, 17.6427, 14.5687, 17.1772, 13.2075, 18.8481, 16.9111, 10.6945};

		 initialize();

		while(txyz[0] < 50) {

			RK4Step(txyz, dt);

		}

		while (h < tamano_entrada + n - 2) {

			PrimercruceenB();

			cout << txyz[3] << endl; // 14.4911 que me genera 322 

			despuesdecruceplanoB();

			// Aqui ya se hizo la perturbacion del sistema

			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

                                 goto Parte2;

						}

				}

				if (txyz[1] > 0) {

					goto Inicio;

					break;

				}

			Parte2:

				{

			//cout << "ENTRO A PARTE 2" << endl;

				if (h > tamano_entrada + n - 2) {

					break;

				}

			//PrimercruceenA(); // Aqui

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

				//u = txyz;

				RK4Step(txyz, dt);

			}

			cout << txyz[3] << endl;

			delta = in[h] - txyz[3];

			cout << "delta es " << delta << endl;

			h = h + 1;

			txyz[3] = txyz[3] + delta;

			cout << "zmean es " << txyz[3] << endl;

			//RK4Step(u, dt);

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

				Bi = '2';

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

			}

            equis[j] = txyz[1];

            cout << "x es " << equis[j] << endl;

            j = j + 1;

			despuesdecruceplanoA ();

			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;

				} // Cruce con el plano A

			}

			if (txyz[1] < 0) {

                goto Parte2;

			}

		}

		Inicio:

			{

				if (h > tamano_entrada + n - 2) {

				break;

				}

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

					RK4Step(txyz, dt);

    			}

				cout << txyz[3] << endl;

				delta = in[h] - txyz[3];

				h = h + 1;

				cout << "delta es " << delta << endl;

				txyz[3] = txyz[3] + delta;

				cout << "zmean es " << txyz[3] << endl;

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

                	Bi = '3';

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

        		}

				equis[j] = txyz[1];

				cout << "x es " << equis[j] << endl;

				j = j + 1;

				despuesdecruceplanoB();

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

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

                        txyzPr2 = txyz;

                        RK4Step(txyz, dt);

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

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

                                goto Parte2;

						}

				}

				if (txyz[1] > 0) {

				goto Inicio;

			}

		}

		}

		for (i=n-1; i< tamano_entrada + n - 1; i++) {

		if (equis[i] < 0) {

			sim[i] = '2';

		}

		if (equis[i] > 0) {

			sim[i] = '3';

		}

		cout << sim[i];

	}

	}

	cout << endl;

	char desim = sim[n-1];

	cout << "desim es " << desim << endl;

	for (i=n; i<tamano_entrada; i++) {

		int l = 0;

		if (sim[i] == desim) {

			desim = sim[i];

			out[l] = 1;

		}

		if (sim[i] != desim) {

			desim = sim[i];

			out[l] = 0;

		}

		cout << out[l];

		l = l + 1;

	}

}

#ifndef CPL_VECTOR_HPP

#define CPL_VECTOR_HPP

#include <complex>

#include <iostream>

namespace cpl {

    class Vector {

        public:

            Vector ( int dim = 1 );

            Vector ( const Vector &dv );

            ~Vector() { delete [] v; }

            int dimension() const { return dim; }

            void resize ( const int );

            const double operator[] ( const int i ) const { return v[i]; }

            double &operator[] ( const int i ) { return v[i]; }

            Vector &operator = ( const Vector &dv );

            Vector &operator += ( const Vector &dv );

            Vector &operator -= ( const Vector &dv );

            Vector &operator *= ( double d );

            Vector &operator /= ( double d );

            double abs();

            double norm();

            double dot ( const Vector &dv );

            friend std::ostream &operator<< ( std::ostream &os, const Vector &dv );

        private:

            int dim;

            double *v;

    };

    inline Vector operator + ( const Vector &dv ) {

        return dv;

    }

    extern Vector operator - ( const Vector &dv );

    extern Vector operator * ( const Vector &dv, double d );

    extern Vector operator * ( double d, const Vector &dv );

    extern Vector operator / ( const Vector &dv, double d );

    extern Vector operator + ( const Vector &v1, const Vector &v2 );

    extern Vector operator - ( const Vector &v1, const Vector &v2 );

    class ComplexVector {

        public:

            ComplexVector ( int dim = 1 );

            ComplexVector ( const ComplexVector &cv );

            ~ComplexVector() { delete [] v; }

            int dimension() const { return dim; }

            const std::complex<double> operator[] ( const int i ) const { return v[i]; }

            std::complex<double> &operator[] ( const int i ) { return v[i]; }

            ComplexVector &operator = ( const ComplexVector &cv );

        private:

            int dim;

            std::complex<double> *v;

    };

    class FFT {

        public:

            FFT() { N = 0; f = 0; inverse = false; }

            void transform ( ComplexVector &data );

            void inverseTransform ( ComplexVector &data );

            Vector power ( ComplexVector &data );

        private:

            int N;

            ComplexVector *f;

            bool inverse;

            void bitReverse();

            void DanielsonLanczos ( int n );

    };

    extern void solveTridiagonal (

        ComplexVector &a, ComplexVector &b, ComplexVector &c,

        ComplexVector &r, ComplexVector &u );

    extern void solveTridiagonalCyclic (

        ComplexVector &a, ComplexVector &b, ComplexVector &c,

        std::complex<double> alpha, std::complex<double> beta,

        ComplexVector &r, ComplexVector &x );

}  /* end namespace cpl */

#endif  /* CPL_VECTOR_HPP */

#include "Vector.hpp"

namespace cpl {

    Vector::Vector ( int dim ) {

        v = new double [this->dim = dim];

        for ( int i = 0; i < dim; i++ ) v[i] = 0;

    }

    Vector::Vector ( const Vector &dv ) {

        v = new double [dim = dv.dim];

        for ( int i = 0; i < dim; i++ ) v[i] = dv.v[i];

    }

    void Vector::resize ( const int dimension ) {

        delete [] v;

        v = new double [dim = dimension];

        for ( int i = 0; i < dim; i++ ) v[i] = 0;

    }

    Vector &Vector::operator = ( const Vector &dv ) {

        if ( this != &dv ) {

            if ( dim != dv.dim ) {

                delete [] v;

                v = new double [dim = dv.dim];

            }

            for ( int i = 0; i < dim; i++ ) v[i] = dv[i];

        }

        return *this;

    }

    Vector &Vector::operator += ( const Vector &dv ) {

        for ( int i = 0; i < dim; i++ ) v[i] += dv[i];

        return *this;

    }

    Vector &Vector::operator -= ( const Vector &dv ) {

        for ( int i = 0; i < dim; i++ ) v[i] -= dv[i];

        return *this;

    }

    Vector &Vector::operator *= ( double d ) {

        for ( int i = 0; i < dim; i++ ) v[i] *= d;

        return *this;

    }

    Vector &Vector::operator /= ( double d ) {

        for ( int i = 0; i < dim; i++ ) v[i] /= d;

        return *this;

    }

    Vector operator - ( const Vector &dv ) {

        int dim = dv.dimension();

        Vector temp ( dim );

        for ( int i = 0; i < dim; i++ ) temp[i] = -dv[i];

        return temp;

    }

    Vector operator * ( const Vector &dv, double d ) {

        int dim = dv.dimension();

        Vector temp ( dim );

        for ( int i = 0; i < dim; i++ ) temp[i] = dv[i] * d;

        return temp;

    }

    Vector operator * ( double d, const Vector &dv ) {

        int dim = dv.dimension();

        Vector temp ( dim );

        for ( int i = 0; i < dim; i++ ) temp[i] = dv[i] * d;

        return temp;

    }

    Vector operator / ( const Vector &dv, double d ) {

        int dim = dv.dimension();

        Vector temp ( dim );

        for ( int i = 0; i < dim; i++ ) temp[i] = dv[i] / d;

        return temp;

    }

    Vector operator + ( const Vector &v1, const Vector &v2 ) {

        int dim = v1.dimension();

        Vector temp ( dim );

        for ( int i = 0; i < dim; i++ ) temp[i] = v1[i] + v2[i];

        return temp;

    }

    Vector operator - ( const Vector &v1, const Vector &v2 ) {

        int dim = v1.dimension();

        Vector temp ( dim );

        for ( int i = 0; i < dim; i++ ) temp[i] = v1[i] - v2[i];

        return temp;

    }

    double Vector::abs() {

        return std::sqrt ( norm() );

    }

    double Vector::norm() {

        double sum = 0;

        for ( int i = 0; i < dim; i++ ) sum += v[i] * v[i];

        return sum;

    }

    double Vector::dot ( const Vector &dv ) {

        double sum = 0;

        for ( int i = 0; i < dim; i++ ) sum += v[i] * dv[i];

        return sum;

    }

    std::ostream &operator<< ( std::ostream &os, const Vector &dv ) {

        for ( int i = 0; i < dv.dim; i++ ) {

            os << dv.v[i];

            if ( i < dv.dim - 1 )

                os << '\t';

            else

                os << '\n';

        }

        return os;

    }

    // ComplexVector implementation

    ComplexVector::ComplexVector ( int dim ) {

        v = new std::complex<double> [this->dim = dim];

        for ( int i = 0; i < dim; i++ ) v[i] = 0.0;

    }

    ComplexVector::ComplexVector ( const ComplexVector &cv ) {

        v = new std::complex<double> [dim = cv.dim];

        for ( int i = 0; i < dim; i++ ) v[i] = cv.v[i];

    }

    ComplexVector &ComplexVector::operator = ( const ComplexVector &cv ) {

        if ( this != &cv ) {

            if ( dim != cv.dim ) {

                delete [] v;

                v = new std::complex<double> [dim = cv.dim];

            }

            for ( int i = 0; i < dim; i++ ) v[i] = cv[i];

        }

        return *this;

    }

    // FFT implementation

    void FFT::transform ( ComplexVector &data ) {

        N = data.dimension();

        f = &data;

        bitReverse();

        for ( int n = 1; n < N; n *= 2 )

            DanielsonLanczos ( n );

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

            ( *f ) [i] /= std::sqrt ( double ( N ) );

    }

    void FFT::inverseTransform ( ComplexVector &data ) {

        inverse = true;

        transform ( data );

        inverse = false;

    }

    void FFT::bitReverse() {

        int j = 1;

        for ( int i = 1; i < N; ++i ) {

            if ( i < j ) {

                std::complex<double> temp = ( *f ) [i - 1];

                ( *f ) [i - 1] = ( *f ) [j - 1];

                ( *f ) [j - 1] = temp;

            }

            int k = N / 2;

            while ( k < j ) {

                j -= k;

                k /= 2;

            }

            j += k;

        }

    }

    void FFT::DanielsonLanczos ( int n ) {

        const double pi = 4 * atan ( 1.0 );

        std::complex<double> W ( 0, pi / n );

        W = inverse ? std::exp ( -W ) : std::exp ( W );

        std::complex<double> W_j ( 1, 0 );

        for ( int j = 0; j < n; ++j ) {

            for ( int i = j; i < N; i += 2 * n ) {

                std::complex<double> temp = W_j * ( *f ) [n + i];

                ( *f ) [n + i] = ( *f ) [i] - temp;

                ( *f ) [i] += temp;

            }

            W_j *= W;

        }

    }

    Vector FFT::power ( ComplexVector &data ) {

        Vector P ( 1 + N / 2 );

        P[0] = std::norm ( data[0] ) / double ( N );

        for ( int i = 1; i < N / 2; i++ )

            P[i] = ( std::norm ( data[i] ) + std::norm ( data[N - i] ) ) / double ( N );

        P[N / 2] = std::norm ( data[N / 2] ) / double ( N );

        return P;

    }

    // Solving tridiagonal complex matrices

    void solveTridiagonal (

        ComplexVector &a, ComplexVector &b, ComplexVector &c,

        ComplexVector &r, ComplexVector &u ) {

        int n = a.dimension();

        ComplexVector gamma ( n );

        std::complex<double> beta = b[0];

        u[0] = r[0] / beta;

        for ( int j = 1; j < n; j++ ) {

            gamma[j] = c[j - 1] / beta;

            beta = b[j] - a[j] * gamma[j];

            u[j] = ( r[j] - a[j] * u[j - 1] ) / beta;

        }

        for ( int j = n - 2; j >= 0; j-- )

            u[j] -= gamma[j + 1] * u[j + 1];

    }

    void solveTridiagonalCyclic (

        ComplexVector &a, ComplexVector &b, ComplexVector &c,

        std::complex<double> alpha, std::complex<double> beta,

        ComplexVector &r, ComplexVector &x ) {

        int n = a.dimension();

        ComplexVector bb ( n ), u ( n ), z ( n );

        std::complex<double> gamma = -b[0];

        bb[0] = b[0] - gamma;

        bb[n - 1] = b[n - 1] - alpha * beta / gamma;

        for ( int i = 1; i < n - 1; i++ )

            bb[i] = b[i];

        solveTridiagonal ( a, bb, c, r, x );

        u[0] = gamma;

        u[n - 1] = alpha;

        for ( int i = 1; i < n - 1; i++ )

            u[i] = 0.0;

        solveTridiagonal ( a, bb, c, u, z );

        std::complex<double> fact = x[0] + beta * x[n - 1] / gamma;

        fact /= 1.0 + z[0] + beta * z[n - 1] / gamma;

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

            x[i] -= fact * z[i];

    }

}  /* end namespace cpl */