% SHELL CONSTANTS %

C = 0.5561613; % PENETRATION
a = 9.81; % GRAVITY
T_0 = 288; % TEMPERATURE AT SEA LEVEL
L = 0.0065; % TEMPERATURE LAPSE RATE
p_0 = 101325; % PRESSURE AT SEA LEVEL
R = 8.31447; % UNIV GAS CONSTANT
M = 0.0289644; % MOLAR MASS OF AIR

% SHELL CONSTANTS %

W = 55; % SHELL WEIGHT
D = 0.152; % SHELL DIAMETER
c_D = 0.321;  % SHELL DRAG
V_0 = 950; % SHELL MUZZLE VELOCITY
K = 2216; % SHELL KRUPP

cw_1 = 1; % QUADRATIC DRAG COEFFICIENT
cw_2 = 100+1000/3*D; % LINEAR DRAG COEFFICIENT

C = C * K/2400; % KRUPP INCLUSION
k = 0.5 * c_D * (D/2)^2 * pi / W; % CONSTANTS TERMS OF DRAG

alpha = [0:0.001:15/360*2*pi]; % ELEV. ANGLES 0...15
dt = 0.1; % TIME STEP

for i = 1:length(alpha) % for each alpha angle do:

    v_x = cos(alpha(i))*V_0;
    v_y = sin(alpha(i))*V_0;
    y = 0;
    x = 0;

    while y >= 0; % follow flight path until shell hits ground again

        x = x + dt*v_x;
        y = y + dt*v_y;

        T = T_0 - L*y;
        p = p_0*(1-L*y/T_0)^(a*M/(R*L));
        rho = p*M/(R*T);

        v_x = v_x - dt*k*rho*(cw_1*v_x^2+cw_2*v_x);
        v_y = v_y - dt*a - dt*k*rho*(cw_1*v_y^2+cw_2*abs(v_y))*sign(v_y);

        t = t + dt;

    end

    v_total = (v_y^2 + v_x^2)^0.5;
    p_athit = C * v_total^1.1 * W^0.55 / (D*1000)^0.65; % PENETRATION FORMULA
    IA = atan(abs(v_y)/abs(v_x)); % IMPACT ANGLE ON BELT ARMOR

    armor(i) = cos(IA)*p_athit;
    distance(i) = x;

end

% INSERT CODE FOR PLOTS HERE %
% EX.: plot(distance,armor,'-m')