#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define DELTA   0.1
#define NSEC    1000

typedef struct point {
    double x;
    double y;
} point_t;

typedef struct velocity {
    double x;
    double y;
} velocity_t;

typedef struct planet {
    point_t p;
    velocity_t v;
    double mass;
} planet_t;

const double G = 6.67428e-11;

planet_t sun = { {0.0, 0.0}, {0.0, 0.0}, 330000*5.9742e+24 };

planet_t p1 = { {0.0, 150000000e+3}, {0.0, 0.0}, 5.9742e+24 };
planet_t p2 = { {0.0, 4700+150000000e+3}, {0.0, 0.0}, 7.347673e+22 }

/*planet_t p1 = { {0.0, 0.0}, {0.0, 0.0}, 5.9742e+24 };
planet_t p2 = { {0.0, 4700.0e+3}, {0.0, 0.0}, 7.347673e+22 };

planet_t proiettile = { {0.0, 3200.0e+3}, {0, 0}, 1 };*/

double gfield[2];

double dist(planet_t* p1, planet_t* p2) {
    double tmp1;
    tmp1 = pow(p1->p.x - p2->p.x, 2) + pow(p1->p.y - p2->p.y, 2);
    printf("%f\n", tmp1);
    return sqrt(tmp1);
}

void calc_gfield(double x, double y) {
    double tmp1, tmp2;
    double distquad;

    gfield[0] = 0;
    gfield[1] = 0;

    /* Contributo di sun */

    if (x != sun.p.x && y != sun.p.x) {
        distquad = pow(sqrt(pow(x-sun.p.x, 2)+pow(y-sun.p.y, 2)), 3);
        tmp1 = G*sun.mass;
        tmp2 = tmp1/distquad;

        gfield[0] -= tmp2*(x-sun.p.x);
        gfield[1] -= tmp2*(y-sun.p.y);
    }

    /* Contributo di p1 */

    if (x != p1.p.x && y != p1.p.y) {
        distquad = pow(sqrt(pow(x-p1.p.x, 2)+pow(y-p1.p.y, 2)), 3);
        tmp1 = G*p1.mass;
        tmp2 = tmp1/distquad;

        gfield[0] -= tmp2*(x-p1.p.x);
        gfield[1] -= tmp2*(y-p1.p.y);
    }

    /* Contributo di p2 */

    if (x != p2.p.x && y != p2.p.y) {
        distquad = pow(sqrt(pow(x-p2.p.x, 2)+pow(y-p2.p.y, 2)), 3);
        tmp1 = G*p2.mass;
        tmp2 = tmp1/distquad;

        gfield[0] -= tmp2*(x-p2.p.x);
        gfield[1] -= tmp2*(y-p2.p.y);
    }

}

int main(int argc, char* argv[]) {
    int t;
    FILE* out;

    out = fopen("asd", "w");

    fprintf(out, "#X\tY\n");

    int count = 0;

    /* Il proiettile viaggia auna velocità che gli permette
     * di orbitare intorno alla terra, più forte ed esce dall'orbita
     * più piano e cade sulla terra */

/*    proiettile.v.x = sqrt(G*p1.mass/dist(&proiettile, &p1));
    proiettile.v.y = 0;*/

    p1.v.x = sqrt(G*sun.mass/dist(&p1, &sun));
    p1.v.y = 0;

    p2.v.x = sqrt(G*p1.mass/dist(&p1, &p2))+sqrt(G*sun.mass/dist(&p2, &sun));
    p2.v.y = 0;

    while (count < 100000) {
        double e_field[2];
        double e_vel[2];

        double m_field[2];
        double m_vel[2];

        double tmpv[2];

/*        fprintf(out, "%f\t%f\n", proiettile.p.x, proiettile.p.y);*/
        fprintf(out, "%f\t%f\n", p1.p.x, p1.p.y);
        fprintf(out, "%f\t%f\n", p2.p.x, p2.p.y);

        /* Salvo accelerazione e velocità correnti per la terra */

        calc_gfield(p1.p.x, p1.p.y);
        e_field[0] = gfield[0];
        e_field[1] = gfield[1];

        e_vel[0] = p1.v.x;
        e_vel[1] = p1.v.y;

/*        printf("g: (%f, %f), v: (%f, %f)\n", e_field[0], e_field[1], p1.v.x, p2.v.x);*/

        /* Salvo accelerazione e velocità correnti per la luna */

        calc_gfield(p2.p.x, p2.p.y);
        m_field[0] = gfield[0];
        m_field[1] = gfield[1];

        m_vel[0] = p2.v.x;
        m_vel[1] = p2.v.y;

        /* Calcolo l'incremento della terra */

        p1.v.x += DELTA*e_field[0];
        p1.v.y += DELTA*e_field[1];

        p1.p.x += DELTA*e_vel[0];
        p1.p.y += DELTA*e_vel[1];

        /* Calcolo l'incremento della luna */

        p2.v.x += DELTA*m_field[0];
        p2.v.y += DELTA*m_field[1];

        p2.p.x += DELTA*m_vel[0];
        p2.p.y += DELTA*m_vel[1];

        /* Calcolo l'incremento del proiettile */

/*        calc_gfield(proiettile.p.x, proiettile.p.y);
        tmpv[0] = proiettile.v.x;
        tmpv[1] = proiettile.v.y;

        proiettile.v.x += DELTA*gfield[0];
        proiettile.v.y += DELTA*gfield[1];

        proiettile.p.x += DELTA*tmpv[0];
        proiettile.p.y += DELTA*tmpv[1];*/

        count++;
    }

    return 0;
}