Demo entry 6344529

Matlab

   

Submitted by anonymous on Jan 13, 2017 at 18:39
Language: Matlab. Code size: 3.3 kB.

clear;
status0=[
    4000*cosd(-5);%V
    4000*sind(-5);
    0;
    0;%P
    70000;
    0;
];
optimal = odeset('AbsTol',1e-12,'RelTol',1e-6);
[t,hypersonic] = ode45(@evoluate,[0 90],status0,optimal);
figure(1);
plot(hypersonic(:,4),hypersonic(:,5));
grid on;
title('trace of vehicle');
figure(2)
plot(t,sqrt(hypersonic(:,1).^2+hypersonic(:,2).^2+hypersonic(:,3).^2));
grid on;
title('velocity-t');
%乌鲁木齐 87:36E 43:46N
%北纬30.67度,东经104.06度
function d_status = evoluate(t, status)

    mass = 1500;%质量
    u = 398603e9;%地心引力常数
    we = 0;%2*pi/24/3600;%地球自传速度
    A0 = 2.356;%发射角
    B0 = 0.7637;%发射点地理纬度
    Sm = 14;%特征面积

    V = status(1:3);
    P = status(4:6);

    q = 0.5*Rho(P(2))*(V'*V);

    We = [we*cos(B0)*cos(A0); we*sin(B0); -we*cos(B0)*sin(A0)];

    R0 = 6376500;

    Ro = [0; R0; 0];

    r = norm(P+Ro);

    Psi = 0;
    Phi = atan(V(2)/V(1));
    Beta = 0;
    Alpha = control(Mah(P(2)));%to be decided

    G_v = [cos(Phi)*cos(Psi), -sin(Phi), cos(Phi)*sin(Psi);
           sin(Phi)*cos(Psi), cos(Psi),  sin(Phi)*sin(Psi);
           -sin(Psi),         0,         cos(Psi)];
    
    Cx = 0.04083 - 0.004167*norm(V)/Mah(P(2)) + 0.001*Alpha;%linear curve fitting
        %阻力系数
    Cy_a = 0.2 - 0.0475*norm(V)/Mah(P(2)) + 0.06193*Alpha;%linear curve fitting
        %升力系数
    R_f = [-Cx*q*Sm*Alpha; Cy_a*q*Sm*Alpha; -Cy_a*q*Sm*Beta];

    g_r = -u/r^2;

    a = [We(1)^2-we^2,      We(1)*We(2),     We(1)*We(3);
         We(1)*We(2),       We(2)^2-we^2,    We(2)*We(3);
         We(1)*We(3),       We(2)*We(3),     We(3)^2-we^2];

    b = [0,                 -2*We(3),        2*We(2);
         2*We(3),           0,               -2*We(1);
         -2*We(2),          2*We(1),         0];
    
    dV = G_v*R_f/mass + g_r/r*(P+Ro) - a*(P+Ro) - b*V;
    dP = V;

    d_status = [dV;dP];

end

function rho = Rho(Height)
    H = Height / (1 + Height / 6356766) / 1000;
    if Height <= 11019.1
        W = 1 - H / 44.3308;
        rho = W^4.2559;
    elseif Height <= 20063.1
        W = exp((14.9647 - H) / 6.3416);
        rho = 1.5898 * 0.1 * W;
    elseif Height <= 32161.9
        W = 1 + (H - 24.9021) / 221.552;
        rho = 3.2722 * 0.01 * W^-35.1629;
    elseif Height <= 47350.1
        W = 1 + (H - 39.7499) / 89.4107;
        rho = 3.2618 * 0.001 * W^-13.2001;
    else%if Height <= 71802
        W = 1 - (H - 59.439) / 88.2218;
        rho = 2.528 * 0.0001 * W^11.2001;
    end
    rho = rho * 1.225;
end
function aT = Mah(Height)
    H = Height / (1 + Height / 6356766) / 1000;
    if Height <= 11019.1
        W = 1 - H / 44.3308;
        T = 288.15 * W;
    elseif Height <= 20063.1
        T = 216.65;
    elseif Height <= 32161.9
        W = 1 + (H - 24.9021) / 221.552;
        T = 221.552 * W;
    elseif Height <= 47350.1
        W = 1 + (H - 39.7499) / 89.4107;
        T = 250.35 * W;
    else%if Height <= 71802
        W = 1 - (H - 59.439) / 88.2218;
        T = 247.021 * W;
    end
    aT = 20.0468 * sqrt(T);
end
function a = control(Ma)
    if Ma >= 7
        a = 10*exp(-0.964*(Ma-7))/57.3;
    elseif Ma<7 && Ma>3
        a = (20 + 5  * (5-Ma))/57.3;
    else
        a = (45 - 5 * Ma)/57.3;
    end
end

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).