#include "PlanetarySystem.h"

#include <GL/glut.h>
#include <math.h>
#include <stdlib.h>

#include "CelestialObject.h"
#include "Vector2.h"

const uint32_t SCREEN_WIDTH = 1000;
const uint32_t SCREEN_HEIGHT = 1000;

const double G = 6.67430 * 1E-11;

// Eccentricity : https://fr.wikipedia.org/wiki/Excentricit%C3%A9_orbitale
const double MERCURY_ECCENTRICITY = 0.20563069;
const double VENUS_ECCENTRICITY = 0.00677323;
const double EARTH_ECCENTRICITY = 0.01671022;
const double MARS_ECCENTRICITY = 0.09341233;

// Mass : https://promenade.imcce.fr/fr/pages5/557.html
const double SUN_MASS = 1.988900 * 1E30;
const double MERCURY_MASS = 0.33018 * 1E24;
const double VENUS_MASS = 4.8685 * 1E24;
const double EARTH_MASS = 5.9736 * 1E24;
const double MARS_MASS = 0.64185 * 1E24;

// Semi-major axis : https://www.le-systeme-solaire.net/demi-grand-axe.html
const double MERCURY_SEMI_MAJOR_AXIS = 579.09227 * 1E8;
const double VENUS_SEMI_MAJOR_AXIS = 108.208475 * 1E9;
const double EARTH_SEMI_MAJOR_AXIS = 149.598262 * 1E9;
const double MARS_SEMI_MAJOR_AXIS = 227.943824 * 1E9;

PlanetarySystem *planetary_system_create() {
    PlanetarySystem *planetary_system = (PlanetarySystem *)malloc(sizeof(PlanetarySystem));

    planetary_system->objects_length = 6;
    planetary_system->objects = (CelestialObject **)malloc(sizeof(PlanetarySystem *) * planetary_system->objects_length);
    planetary_system->zoom_factor = 1;
    planetary_system->reference_frame_object_index = 0;

    planetary_system->objects[0] = celestial_object_create(SUN_MASS, 0, 0, 50, 0xFFFFFF);
    planetary_system->objects[1] = celestial_object_create(MERCURY_MASS, MERCURY_SEMI_MAJOR_AXIS, MERCURY_ECCENTRICITY, 10, 0xDBCECA);
    planetary_system->objects[2] = celestial_object_create(VENUS_MASS, VENUS_SEMI_MAJOR_AXIS, VENUS_ECCENTRICITY, 20, 0x8B7D82);
    planetary_system->objects[3] = celestial_object_create(EARTH_MASS, EARTH_SEMI_MAJOR_AXIS, EARTH_ECCENTRICITY, 20, 0x6b93d6);
    planetary_system->objects[4] = celestial_object_create(MARS_MASS, MARS_SEMI_MAJOR_AXIS, MARS_ECCENTRICITY, 12, 0xBC2732);
    planetary_system->objects[5] = celestial_object_create(7.34767309 * 1E22, 384.399 * 1E6, 0.0549, 5, 0x8A2BE2);
    planetary_system->objects[5]->current_position.x += planetary_system->objects[3]->current_position.x;

    return planetary_system;
}

CelestialObject *planetary_system_get_star(PlanetarySystem *planetary_system) {
    return planetary_system->objects[0];
}

Vector2 calculate_gravitational_acceleration(PlanetarySystem *planetary_system, uint32_t object_index) {
    Vector2 a = vector2_create_zero();

    for (uint32_t i = 0; i < planetary_system->objects_length; i += 1) {
        if (i == object_index)
            continue;
        Vector2 r = vector2_substract(planetary_system->objects[i]->current_position, planetary_system->objects[object_index]->current_position);
        double a_scalar = G * planetary_system->objects[i]->mass * pow(pow(vector2_norm(r), 2), -1);
        a = vector2_add(a, vector2_multiply(vector2_normalize(r), a_scalar));
    }

    return a;
}

void planetary_system_update(PlanetarySystem *planetary_system, double interval) {
    for (uint32_t i = 1; i < planetary_system->objects_length; i += 1) {
        CelestialObject *object = planetary_system->objects[i];
        Vector2 current_position = object->current_position;
        Vector2 new_position;

        if (vector2_is_zero(object->previous_position)) {
            CelestialObject *star = planetary_system_get_star(planetary_system);

            if (i == 5) {
                star = planetary_system->objects[3];
            }

            double periapsis_velocity_scalar = sqrt((G * star->mass * (1 + object->eccentricity)) / (object->semi_major_axis * (1 - object->eccentricity)));

            if (i == 5) {
                periapsis_velocity_scalar += 30290.322245;
            }

            Vector2 r = vector2_normalize(vector2_create(current_position.y, -current_position.x));
            Vector2 periapsis_velocity = vector2_multiply(r, periapsis_velocity_scalar);

            new_position = vector2_add(current_position, vector2_multiply(periapsis_velocity, interval));
            Vector2 a = calculate_gravitational_acceleration(planetary_system, i);
            new_position = vector2_add(new_position, vector2_multiply(a, 0.5 * pow(interval, 2)));
        } else {
            new_position = vector2_substract(vector2_multiply(current_position, 2), object->previous_position);
            Vector2 a = calculate_gravitational_acceleration(planetary_system, i);
            new_position = vector2_add(new_position, vector2_multiply(a, pow(interval, 2)));
        }

        object->previous_position = object->current_position;
        object->current_position = new_position;

        int32_t length = object->points_length;

        if (length > 0 && vector2_norm(vector2_substract(object->points[0], object->current_position)) < 1E9 * 1.5)
            continue;
        for (int32_t j = (length == 200) ? length - 1 : length; j >= 1; j -= 1) {
            object->points[j] = object->points[j - 1];
        }

        object->points[0] = object->current_position;

        if (length < 200)
            object->points_length += 1;
    }
}

Vector2 scale_position(PlanetarySystem *planetary_system, Vector2 object_position, Vector2 reference_frame_position) {
    Vector2 unscaled_position = vector2_substract(object_position, reference_frame_position);
    unscaled_position = vector2_multiply(unscaled_position, planetary_system->zoom_factor);

    Vector2 scaled_position = vector2_multiply(unscaled_position, 1.0 / (300 * 1E9));
    scaled_position = vector2_fit_canvas(scaled_position, SCREEN_WIDTH, SCREEN_HEIGHT);
    return scaled_position;
}

void draw_disc(Vector2 position, uint32_t radius) {
    glBegin(GL_POLYGON);
    for (uint32_t i = 0; i < 360; i += 1) {
        double theta = i * 3.1415 / 180;
        double x = position.x + radius * cos(theta);
        double y = position.y + radius * sin(theta);
        glVertex2f(x, y);
    }
    glEnd();
}

void planetary_system_draw(PlanetarySystem *planetary_system) {
    CelestialObject *reference_frame = planetary_system->objects[planetary_system->reference_frame_object_index];

    for (uint32_t i = 0; i < planetary_system->objects_length; i += 1) {
        CelestialObject *object = planetary_system->objects[i];
        Vector2 scaled_position = scale_position(planetary_system, object->current_position, reference_frame->current_position);

        uint32_t color = object->drawing_color;
        glColor3ub((color & 0xFF0000) >> 16, (color & 0x00FF00) >> 8, (color & 0x0000FF) >> 0);

        if (planetary_system->zoom_factor < 1) {
            draw_disc(scaled_position, object->drawing_disc_radius * planetary_system->zoom_factor);
        } else {
            draw_disc(scaled_position, object->drawing_disc_radius);
        }

        if (i != 5) {
            glLineWidth(4.0);

            for (int32_t j = 0; j < object->points_length - 1; j += 1) {
                glBegin(GL_LINES);
                Vector2 aaa = scale_position(planetary_system, object->points[j], reference_frame->current_position);
                glVertex2f(aaa.x, aaa.y);
                Vector2 bbb = scale_position(planetary_system, object->points[j + 1], reference_frame->current_position);
                glVertex2f(bbb.x, bbb.y);
                glEnd();
            }
        }
    }
}