planet.c

#include "planet.h"
#include <stdlib.h>
#include <time.h>

#define G 6.67e-11
#define M_SUN 1.989e30
#define MAX_RANGE 4e8
#define EARTH_RADIUS 6371/800

#define NUMBER_OF_PLANETS 10

/// @brief Crée et retourne une planète
/// @param semi_major_axis  km
/// @param mass             kg
/// @param eccentricity     **
/// @param velocity         vec2
/// @param radius           pixels
/// @param color            color
/// @param pos              vec2
/// @return 
planet_t create_planet(double semi_major_axis, double mass, double eccentricity, vec2 velocity, double radius, uint32_t color, vec2 pos, bool clockwise)
{
    return (planet_t){semi_major_axis, mass, eccentricity, velocity, radius, color, pos, clockwise};
}

/// @brief Crée et retourne le systeme solaire
/// @param delta_t
/// @return 
system_t create_system(double delta_t, bool active_fictive_planet)
{
    // Allocation de mémoire pour le tableau de planètes
    planet_t *planets = malloc(NUMBER_OF_PLANETS*sizeof(planet_t));

    //*1. Créer une étoile au centre de notre domaine.
    // Soleil
    planet_t sun = create_planet(0, M_SUN, 0, vec2_create_zero(), 30, COLOR_YELLOW, vec2_create_zero(), true);
    
    system_t result = (system_t){sun, NUMBER_OF_PLANETS, planets};
    
    //*2. On ajoute autant de planètes que l’on désire autour de l’étoile.
    // Mercure
    planet_t mercury = create_planet(57910000, 3.301e23, 0.205630, vec2_create_zero(), EARTH_RADIUS*0.383, MAKE_COLOR(219, 211, 204), vec2_create_zero(), true);
    planets[0] = mercury;
    // Vénus 
    planet_t venus = create_planet(108210000, 4.867e24, 0.006772, vec2_create_zero(), EARTH_RADIUS*0.950, MAKE_COLOR(204, 180, 161), vec2_create_zero(), true);
    planets[1] = venus;
    // Earth
    planet_t earth = create_planet(149598023, 5.973e24, 0.0167086, vec2_create_zero(), EARTH_RADIUS, COLOR_BLUE, vec2_create_zero(), true);
    planets[2] = earth;
    // Mars
    planet_t mars = create_planet(227939366, 6.418e23, 0.0934, vec2_create_zero(), EARTH_RADIUS*0.532, MAKE_COLOR(212, 102, 19), vec2_create_zero(), true);
    planets[3] = mars;

    // Jupiter
    planet_t jupyter = create_planet(778340000, 1.8986e27, 0.04839, vec2_create_zero(), EARTH_RADIUS*11.209, MAKE_COLOR(232, 153, 76), vec2_create_zero(), true);
    planets[4] = jupyter;

    // Saturne
    planet_t saturn = create_planet(1426700000, 5.6846e26, 0.0539, vec2_create_zero(), EARTH_RADIUS*9.4492, MAKE_COLOR(203, 175, 112), vec2_create_zero(), true);
    planets[5] = saturn;

    // Uranus
    planet_t uranus = create_planet(2870700000 , 0.04726e25, 0.04726, vec2_create_zero(), EARTH_RADIUS*4.00723, MAKE_COLOR(191, 210, 215), vec2_create_zero(), true);
    planets[6] = uranus;

    // Neptune
    planet_t neptune = create_planet(4498400000 , 1.0243e26, 0.00859, vec2_create_zero(), EARTH_RADIUS*3.883, MAKE_COLOR(117, 191, 214), vec2_create_zero(), true);
    planets[7] = neptune;

    //* Planètes fictives
    if (active_fictive_planet)
    {
        //* PLanètes fictives
        planet_t tatooine = create_planet(249598023, 9.267e24, 0.532, vec2_create_zero(), EARTH_RADIUS*1.640796488, MAKE_COLOR(145, 118, 122), vec2_create(0, 312073200), false);
        planets[8] = tatooine;
        planet_t naboo = create_planet(349588423, 1.52e25, 0.028, vec2_create_zero(), EARTH_RADIUS*1.90028222, MAKE_COLOR(113, 182, 169), vec2_create(0, -212073200), true);
        planets[9] = naboo;
    }
    
    
    //*3. On définit nos conditions initiales (en particulier la vitesse de chaque planète).
    for (size_t i = 0; i < result.nb_planets; i++)
    {
        init_planet(&result.planets[i], delta_t);
    }
    return result;
}

/// @brief Initialise la vélocité et la position d'une planète
/// @param planet 
/// @param delta_t 
void init_planet(planet_t *planet, double delta_t)
{
    //* Calcul de la position initiale de la planète
    // Calul de la périhélie
    double periapsis = planet->semi_major_axis * (1 - planet->eccentricity);
    planet->pos = vec2_create(periapsis, planet->pos.y);

    //* Calcul de la vélocité initiale de la planète
    // Calcul du vecteur perpendiculaire à la direction de la planète
    // Si la planète tourne dans le sens horaire
    vec2 perp;
    if (planet->clockwise)
    {
        perp = vec2_create(-planet->pos.y, planet->pos.x);
    }
    // Si la planète tourne dans le sens anti-horaire
    else
    {
        perp = vec2_create(planet->pos.y, -planet->pos.x);
    }
    
    // Calcul de la velocité
    vec2 velocity = vec2_mul(sqrt(G * M_SUN * (2 / vec2_norm(planet->pos) - 1 / planet->semi_major_axis)), vec2_normalize(perp));
    // Affectation de la velocité
    planet->velocity = velocity;
    // Affectation de la nouvelle position
    vec2 pos = vec2_add(planet->pos, vec2_mul(delta_t, velocity));
    pos = vec2_add(pos, vec2_mul(delta_t*delta_t/(2*planet->semi_major_axis), vec2_mul(-1, planet->pos)));
    // Affectation de la nouvelle position
    planet->pos = pos;
}

/// @brief Met à jour l'affichage
/// @param ctxt 
/// @param system 
void show_system(struct gfx_context_t *ctxt, system_t *system, double scale)
{
    // Afficher l'étoile
    draw_full_circle(
            ctxt,
            500,
            500,
            system->star.radius / scale,
            system->star.color
            );
    // Afficher chaque planète
    for (uint32_t i = 0; i < system->nb_planets; i++)
    {
        // Afficher la planète
        vec2 pos = vec2_create(system->planets[i].pos.x/MAX_RANGE, system->planets[i].pos.y/MAX_RANGE);
        pos = vec2_mul(1/scale, pos);
        coordinates coord = vec2_to_coordinates(pos, ctxt->width, ctxt->height);
        draw_full_circle(
            ctxt,
            coord.column,
            coord.row,
            system->planets[i].radius / scale,
            system->planets[i].color
            );
    }
}

/// @brief Applique la force de planet_b sur planet_a
/// @param planet_a planète recevante
/// @param planet_b
/// @return 
vec2 compute_force(planet_t *planet_a, planet_t *planet_b)
{
    // Calcul de la distance entre les deux planètes
    double distance = vec2_distance(planet_a->pos, planet_b->pos);
    // Calcul de la force gravitationnelle
    double force = G * planet_a->mass * planet_b->mass / (distance * distance);
    // Calcul de la force résultante
    vec2 direction = vec2_normalize(vec2_sub(planet_b->pos, planet_a->pos));
    vec2 force_vector = vec2_mul(force, direction );
    return force_vector;
}

/// @brief Applique la force de tout les objets du système sur une planète
/// @param planet 
/// @param system 
/// @return 
vec2 compute_all_forces(planet_t *planet, system_t *system)
{
    vec2 result = vec2_create_zero();
    // *Calcule de la force résultante
    // Avec chaque planètes
    for (size_t i = 0; i < system->nb_planets; i++)
    {
        planet_t other_planet = system->planets[i];
        if (other_planet.mass != planet->mass)
        {
            vec2 force = compute_force(planet, &other_planet);
            result = vec2_add(result, force);
        }
    }
    // Avec le soleil
    vec2 force = compute_force(planet, &system->star);
    result = vec2_add(result, force);
    return result;
}

/// @brief Met à jour l'entiereté du système
/// @param system 
/// @param delta_t 
void update_system(system_t *system, double delta_t)
{
    // *Calculer la nouvelle position de chaque planète
    for (size_t i = 0; i < system->nb_planets; i++)
    {
        planet_t *planet = &system->planets[i];
        // Calcul de la force résultante
        //*     2. On calcule la force résultante sur chacune des planètes.
        vec2 force = compute_all_forces(planet, system);
        // Calcul de l'accélération
        vec2 acceleration = vec2_mul(1/planet->mass, force);
        // Calcul de la vitesse
        vec2 velocity = vec2_add(planet->velocity, vec2_mul(delta_t, acceleration));
        //*     3. On calcule la prochaine position de chacune des planètes
        // Calcul de la position
        vec2 pos = vec2_add(planet->pos, vec2_mul(delta_t, velocity));
        //*     4. On met à jour la position des planètes de notre système.
        // Affectation de la position
        planet->pos = pos;
        // Affectation de la vitesse
        planet->velocity = velocity;
    }
}

/// @brief Libère la mémoire du système
/// @param system 
void free_system(system_t *system)
{
    free(system->planets);
}

/// @brief Ajoute aléatoirement des étoiles suivant une graine donnée
/// @param ctxt 
/// @param quantity 
/// @param seed 
void put_stars(struct gfx_context_t *ctxt, int quantity, int seed)
{
    srand(seed);
    for (int i = 0; i < quantity; i++)
    {
        int x = rand() % ctxt->width;
        int y = rand() % ctxt->height;
        gfx_putpixel(ctxt, x, y, COLOR_WHITE);
    }
}

/// @brief Ajoute une planète aléatoire dans le système
/// @param system 
/// @param delta_t 
void add_random_planet_in_system(system_t *system, double delta_t)
{
    system->planets = realloc(system->planets,(system->nb_planets+1)*sizeof(planet_t));

    double min, max = 0;
    //* On génère des valeurs aléatoires pour notre planète
    // Le demi grand axe
    min = 1e8; max = 4e10;
    double semi_major_axis = fmod(rand(), (max - min + 1)) + min;
    // La masse
    min = 3e23; max = M_SUN;
    double mass = fmod(rand(), (max - min + 1)) + min;
    // L'excentricité
    min = 0.0001; max = 0.9999; 
    double eccentricity = fmod(rand(), (max - min + 1)) + min;
    // Le rayon
    min = 0.3; max = 2; 
    double radius = fmod(rand(), (max - min + 1)) + min;
    // La couleur
    int red = rand() % 255;
    int green = rand() % 255;
    int blue = rand() % 255;
    // La position
    vec2 pos = vec2_create(0, 0);
    // Le sens de rotation
    bool clockwise = rand() % 2;
    // On crée notre planète
    planet_t p = create_planet(semi_major_axis, mass, eccentricity, vec2_create_zero(), EARTH_RADIUS*radius, MAKE_COLOR(red, green, blue), pos, clockwise);
    init_planet(&p, delta_t);
    system->planets[system->nb_planets] = p;
    system->nb_planets+=1;
}