WoWs Penetration calc

pen.m

clear all; clc
% 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 = input('Shell Weight in Kg: '); % SHELL WEIGHT
D = input('Shell Diameter in Meters: '); % SHELL DIAMETER
c_D = input('Shell Drag Coefficient: ');  % SHELL DRAG
V_0 = input('Shell Starting Velocity in m/s: '); % SHELL MUZZLE VELOCITY
K = input('Shell KRUPP Value: '); % SHELL KRUPP
angle = input('Gun Elevation Angle in Degrees: ');
AT = input('Armor Plate Thickness in mm: ');

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:angle/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);
      
        y = y + dt;

    end
    
    v_t(i) = (v_y^2 + v_x^2)^0.5; % Total Velocity Index
    p_athit = C * v_t(i)^1.1 * W^0.55 / (D*1000)^0.65; % PENETRATION FORMULA
    IA = atan(abs(v_y)/abs(v_x)); % IMPACT ANGLE ON BELT ARMOR
    IA2 = 0;
    armor(i) = cos(IA)*p_athit;
    armor2(i) = cos(IA2)*p_athit;
    distance(i) = x;
    
end

for i = 1:length(alpha)
    v_2(i) = v_t(i) * (1-exp((AT-p_athit)/AT)); %Shell Velocity After Penetration of Inital Armor
    v_x2 = cos(alpha(i))*v_2;
    v_y2 = sin(alpha(i))*v_2;
    y2 = 0;
    x2 = 0;
    
    while y2 >= 0 

        x2 = x2 + dt*v_x2;
        y2 = y2 + dt*v_y2;
        
        T = T_0 - L*y2;
        p = p_0*(1-L*y2/T_0).^(a*M/(R*L));
        rho = p*M/(R*T);
        
        v_x2 = v_x2 - dt*k*rho*(cw_1*v_x2.^2+cw_2*v_x2);
        v_y2 = v_y2 - dt.*a - dt.*k.*rho.*(cw_1.*v_y2.^2+cw_2.*abs(v_y2)).*sign(v_y2);
      
        y2 = y2 + dt;

    end
    
    p_athit2 = C * v_2(i)^1.1 * W^0.55 / (D*1000)^0.65;
    IA1 = atan(abs(v_y2)/abs(v_x2));
    armor3(i) = cos(IA1)*p_athit2;
    if armor3(i) < 0
        armor3(i) = 0;
    end
end

subplot(2,2,1);
plot(distance,armor,'-r')
hold on
plot(distance,armor2,'-b')
plot(distance,armor3)
str = sprintf('Penetration Sweep for Shell: %.0f mm, %.0f kg, %.0f m/s MV,Drag Constant %.3f, At %.0f Degress Turret Elevation\nWith Shell Velocity After %.0f mm Armor Plate Penetration',D*1000,W,V_0,c_D,angle,AT);
title(str)
ylabel('Penetration [mm]')
xlabel('Distance [m]')
grid 
set(gca,'FontSize',10,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',10,'fontWeight','bold')
legend('With Impact Angle and 0 Degrees Armor Angle','At 90 Degrees Impact Angle', 'After Armor Plate Penetration')
subplot(2,2,2);
plot(distance,v_t)
hold on
plot(distance,v_2)
str = sprintf('Velocity Sweep for Shell: %.0f mm, %.0f kg, %.0f m/s MV,Drag Constant %.3f, At %.0f Degrees Turret Elevation\nWith Shell Velocity After %.0f mm Armor Plate Penetration',D*1000,W,V_0,c_D,angle,AT);
title(str)
ylabel('Velocity [m/s]')
xlabel('Distance [m]')
grid 
set(gca,'FontSize',10,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',10,'fontWeight','bold')
legend('Shell Velocity with Respect to Distance','Shell Velocity After Armor Plate')