MATLAB: Missile Lauch friction gamma factor – Output deviating from desired output, please help!

earth curvaturefactorfrictiongammagravitylaunchmissilepepeno

Goal:
I'm trying to implement friction into this script of a launch of a missile.
The friction force should be related to the velocity like: Ff = -gamma|v|v.
Subsequently we also have the gravity force acting on the missile: Fg = mg.
The values of ad should be as follows, but I am unable to get this desired output:
ad = [300332, 302763, 303903, 303771, 302384, 299763, 295927, 290896, 284693];
for degrees of launch 37, 39, 41, 43, 45, 47, 49, 51, 53 respectively.
Script:
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
Output:
ad = [200358, 203662, 205972, 207257, 207539, 206800, 205056, 202324, 198615];
for degrees of launch 37, 39, 41, 43, 45, 47, 49, 51, 53 respectively.
Any help / suggestions are immensly appreciated!

Best Answer

Step 1: Draw a free body diagram that shows the body with all of the forces that act on it and the coordinate frame(s) that are being used.
Step 2: For a point-mass in 2 dimensions, you should be able to use the free body diagram (by inspection) to ensure that you are calculating the force components correctly; sum of the forces in X = m*Ax, sum of forces in Y = m*Ay.
Are forces defined in body or inertial axes? (The aerodynamic force is usually in body axes, but the gravity force is defined in Inertial axes)
Show us this drawing and we can work from there.
Related Question