// Programa para encontrar la orbita de los tres cuerpos usando matrices

//--- Librerias ---
#include <iostream>
#include <cmath>
#include <fstream>
#include <cstdio>
//-----------------

//--- Directiva de using ---
using namespace std;
//--------------------------

// --- CONSTANTE GLOBAL ---
const double pi = acos(-1.0);
const double GMs = 4 * pi * pi; // u.a.^3 / año^2 (para el Sol)[cite: 2]

// Estructura matricial para almacenar estados de los 3 cuerpos en 3D (Fila: Cuerpo, Columna: Coordenada x,y,z)
struct Matriz3x3 {
    double m[3][3];
};

struct CCeleste {
    double r[3]; // r[0]->x, r[1]->y, r[2]->z[cite: 1, 2]
    double v[3]; // v[0]->vx, v[1]->vy, v[2]->vz[cite: 1, 2]
    double m_cuerpo;
    double Periodo;

    void EulerCromer(double Dt, CCeleste Planeta[], int j, int Np)
    {
        double ax = 0.0;
        double ay = 0.0;
        double az = 0.0;

        // Suma de fuerzas gravitacionales de los demás cuerpos (Ecuación de N cuerpos)[cite: 1]
        for (int k = 0; k < Np; k++) {
            if (k == j) continue;

            double rx = r[0] - Planeta[k].r[0];
            double ry = r[1] - Planeta[k].r[1];
            double rz = r[2] - Planeta[k].r[2];
            double r2 = rx * rx + ry * ry + rz * rz;
            double r3 = r2 * sqrt(r2);

            // GM_k en función de la masa del Sol[cite: 2]
            double GMk = GMs * (Planeta[k].m_cuerpo / Planeta[0].m_cuerpo);

            ax -= (GMk * rx) / r3;
            ay -= (GMk * ry) / r3;
            az -= (GMk * rz) / r3;
        }

        v[0] += ax * Dt;
        v[1] += ay * Dt;
        v[2] += az * Dt;

        r[0] += v[0] * Dt;
        r[1] += v[1] * Dt;
        r[2] += v[2] * Dt;
    }
};

//--- Funcion principal ---
int main()
{
    const int Np = 3; // Sol (0), Tierra (1), Luna (2)[cite: 1]

    string Nombre[Np] = {"Sol", "Tierra", "Luna"};
    double Masa[Np];

    // Masas en kg[cite: 2]
    Masa[0] = 1.989E30; // Sol
    Masa[1] = 5.972E24; // Tierra
    Masa[2] = 7.342E22; // Luna

    CCeleste Planeta[Np];

    for (int j = 0; j < Np; j++) {
        Planeta[j].m_cuerpo = Masa[j];
    }

    // --- MATRICES DE CONDICIONES INICIALES SEGÚN LA GUÍA ---[cite: 2, 3]
    Matriz3x3 R_CM; // Posiciones iniciales (Sol y CM Tierra-Luna)
    Matriz3x3 V_CM; // Velocidades iniciales en au/día

    // Sol (0)
    R_CM.m[0][0] =  4.504795855674601E-3;
    R_CM.m[0][1] =  7.732544746890742E-4;
    R_CM.m[0][2] =  2.685039985573271E-4;
    V_CM.m[0][0] = -3.524457445683394E-7;
    V_CM.m[0][1] =  5.177637780672221E-6;
    V_CM.m[0][2] =  2.229113252400401E-6;

    // Tierra en CM Tierra-Luna (1)
    R_CM.m[1][0] =  1.205197129002194E-1;
    R_CM.m[1][1] = -9.258323035784091E-1;
    R_CM.m[1][2] = -4.015377611192672E-1;
    V_CM.m[1][0] =  1.681126759384403E-2;
    V_CM.m[1][1] =  1.748308899617618E-3;
    V_CM.m[1][2] =  7.582041930175932E-4;

    // Luna en CM Tierra-Luna (2)
    R_CM.m[2][0] = -8.081772358351251E-4;
    R_CM.m[2][1] = -1.994630037441996E-3;
    R_CM.m[2][2] = -1.087262721620868E-3;
    V_CM.m[2][0] =  6.010848314829119E-4;
    V_CM.m[2][1] = -1.674454691500606E-4;
    V_CM.m[2][2] = -8.556208109904861E-5;

    // --- TRANSFORMACIÓN DE MATRICES AL SISTEMA ABSOLUTO (INERCIAL) ---[cite: 3]
    double m2 = Masa[1];
    double m3 = Masa[2];
    double factor_r3 = m3 / (m2 + m3);
    double factor_r2 = m2 / (m2 + m3);

    Matriz3x3 R_Abs;
    Matriz3x3 V_Abs;

    // Sol Absoluto
    for (int k = 0; k < 3; k++) {
        R_Abs.m[0][k] = R_CM.m[0][k];
        V_Abs.m[0][k] = V_CM.m[0][k] * 365.25; // Conversión de au/día a au/año[cite: 3]
    }

    // Tierra y Luna Absolutas usando operaciones matriciales por componentes[cite: 3]
    for (int k = 0; k < 3; k++) {
        // r_2 = r_2' - [m3/(m2+m3)] * r_3'
        R_Abs.m[1][k] = R_CM.m[1][k] - factor_r3 * R_CM.m[2][k];
        V_Abs.m[1][k] = (V_CM.m[1][k] - factor_r3 * V_CM.m[2][k]) * 365.25;

        // r_3 = r_2' + [m2/(m2+m3)] * r_3'
        R_Abs.m[2][k] = R_CM.m[1][k] + factor_r2 * R_CM.m[2][k];
        V_Abs.m[2][k] = (V_CM.m[2][k] + factor_r2 * V_CM.m[2][k]) * 365.25;
    }

    // Asignar a los objetos celestes
    for (int j = 0; j < Np; j++) {
        Planeta[j].m_cuerpo = Masa[j];
        for (int k = 0; k < 3; k++) {
            Planeta[j].r[k] = R_Abs.m[j][k];
            Planeta[j].v[k] = V_Abs.m[j][k];
        }
    }

    ofstream archivos[Np];
    ofstream archivosLz[Np];

    for (int i = 0; i < Np; i++) {
        string nombre = Nombre[i] + ".dat";
        archivos[i].open(nombre);

        string nombreLz = Nombre[i] + "_Lz.dat";
        archivosLz[i].open(nombreLz);

        if (!archivos[i] || !archivosLz[i]) {
            cout << "Error al crear los archivos de salida." << endl;
            return 1;
        }
    }

    double dt = 0.001; // Paso temporal recomendado en la guía[cite: 2]
    int N_pasos = 1000;

    //===========================================================
    // SIMULACION DE LAS ORBITAS Y CALCULO DEL MOMENTO ANGULAR
    //===========================================================

    for (int j = 0; j < Np; j++) {
        double t = 0.0;
        archivos[j] << t << " " << Planeta[j].r[0] << " " << Planeta[j].r[1] << " " << Planeta[j].r[2] << "\n";

        double Lz = Planeta[j].m_cuerpo * (Planeta[j].r[0] * Planeta[j].v[1] - Planeta[j].r[1] * Planeta[j].v[0]);
        double dAdt = 0.5 * (Planeta[j].r[0] * Planeta[j].v[1] - Planeta[j].r[1] * Planeta[j].v[0]);

        archivosLz[j] << t << " " << Lz << " " << dAdt << "\n";
    }

    for (int k = 1; k <= N_pasos; k++) {
        double t = k * dt;

        // Actualizar todos los cuerpos mediante el método matricial Euler-Cromer
        for (int j = 0; j < Np; j++) {
            Planeta[j].EulerCromer(dt, Planeta, j, Np);
        }

        for (int j = 0; j < Np; j++) {
            double Lz = Planeta[j].m_cuerpo * (Planeta[j].r[0] * Planeta[j].v[1] - Planeta[j].r[1] * Planeta[j].v[0]);
            double dAdt = 0.5 * (Planeta[j].r[0] * Planeta[j].v[1] - Planeta[j].r[1] * Planeta[j].v[0]);

            archivosLz[j] << t << " " << Lz << " " << dAdt << "\n";
            archivos[j] << t << " " << Planeta[j].r[0] << " " << Planeta[j].r[1] << " " << Planeta[j].r[2] << "\n";
        }
    }

    for (int j = 0; j < Np; j++) {
        Planeta[j].Periodo = N_pasos * dt;
        archivos[j].close();
        archivosLz[j].close();

        cout << "Cuerpo: " << Nombre[j] << endl;
        cout << "Tiempo simulado: " << Planeta[j].Periodo << " años" << endl;
        cout << "-----------------------------------" << endl;
    }

    //===========================================================
    // GRAFICAS CON GNUPLOT
    //===========================================================

    FILE* gnuplotPipe = popen("gnuplot", "w");

    if (gnuplotPipe != nullptr) {
        fprintf(gnuplotPipe, "set terminal pngcairo size 800,800 font 'Arial,12'\n");
        fprintf(gnuplotPipe, "set grid\n");
        fprintf(gnuplotPipe, "set zeroaxis lt 1 lc rgb 'gray' lw 1\n");

        for (int j = 0; j < Np; j++) {
            string archivoPNG = Nombre[j] + "_orbita3D.png";
            string archivoDAT = Nombre[j] + ".dat";

            fprintf(gnuplotPipe, "set output '%s'\n", archivoPNG.c_str());
            fprintf(gnuplotPipe, "set title 'Orbita Tridimensional de %s'\n", Nombre[j].c_str());
            fprintf(gnuplotPipe, "set xlabel 'Posicion X (u.a.)'\n");
            fprintf(gnuplotPipe, "set ylabel 'Posicion Y (u.a.)'\n");
            fprintf(gnuplotPipe, "set zlabel 'Posicion Z (u.a.)'\n");
            // Gráfica 3D solicitada en la guía usando splot[cite: 3]
            fprintf(gnuplotPipe, "splot '%s' using 2:3:4 with lines lw 2 lc rgb '#1f77b4' title 'Trayectoria'\n", archivoDAT.c_str());

            cout << "Grafica '" << archivoPNG << "' generada." << endl;
        }

        fprintf(gnuplotPipe, "unset output\n");
        pclose(gnuplotPipe);

        cout << "\nTodas las graficas se han generado correctamente." << endl;
    } else {
        cerr << "No se pudo conectar con Gnuplot." << endl;
    }

    return 0;
}