clear;
close all;
clc;
tic
m = 500;
aT = [];
ad = [];
gamma = 1.0e-5
for i = [37:2:53]
    r = 6371 * 10^3;
    G = 6.674 * 10^-11;
    M = 5.972 * 10^24;
    g = (G * M)/(r^2);
    theta0 = i;
    ax = 0;
    ay = r;
    v0 = 2000;
    vx0 = v0*cosd(theta0);
    vy0 = v0*sind(theta0);
    x = 0;
    y = r;
    vx = vx0;
    vy = vy0;
    T = 0;
    dt = 0.01;
    at = 0;
    landed = 0;
    z = 1;
    while landed == 0
        z = z + 1;
        T = T + dt;
        xo = x;
        yo = y;
        x = x + vx * dt;
        y = y + vy * dt;
        d = sqrt(x^2 + y^2);
        alpha = atand(x/y);
        g = (G*M)/(d^2);
        gy = cosd(alpha) * g;
        gx = sind(alpha) * g;
        FgY = m * gy;
        FgX = m * gx;
        FricY = gamma*abs(vy)*vy;
        FricX = gamma*abs(vx)*vx;
        FnetY = FgY - FricY;
        FnetX = FgX - FricX;
        vy = vy - (gy * dt) - (FnetY / m)*dt; % This should substract gravity and friction from the velocity, not sure if this is correct.
        vx = vx - (gx * dt) - (FnetX / m)*dt; % Same for this line of code
        v = vx/sin(alpha);
        ax = [ax, x];
        ay = [ay, y];
        if d < r
            landed = 1;
        end
    end

    aT = [aT, T];
    distance = (alpha/360) * 2 * pi * r;
    ad = [ad, distance];

    fprintf('Checked for degree: %.0f\n', i)

end
toc