function roc = rocketPlot(mission, rocket)
arguments
    mission Mission
    rocket Rocket
end
    deltaXLE = rocket.rear.finsDeltaXFreeChord;
    r = rocket.diameter/2;
    height = rocket.rear.finsHeight;
    C1 = rocket.rear.finsRootChord;
    C2 = rocket.rear.finsFreeChord;

    Lcent = rocket.lengthCenter;
    Lnos = rocket.parafoil.noseLength;

    Xle1 = Lcent + Lnos - rocket.rear.finsAxialDistance - C1;
    
    Daft = rocket.rear.boatFinalDiameter;
    Laft = rocket.rear.boatLength;

    xNos = Lnos - Lnos*cos(linspace(0, pi/2, 50));
    if strcmp(rocket.parafoil.noseType, 'KARMAN')
        theta = @(x) acos( 1 - ( (2*x)./Lnos ) );
        Karman = @(x) ( r/sqrt(pi) ) * sqrt( theta(x) - ( sin(2*theta(x))./ 2 ) );
        yNos = Karman(xNos);
    elseif strcmp(rocket.parafoil.noseType, 'HAACK')
        theta = @(x) acos( 1 - ( (2*x)./Lnos ) );
        Haack = @(x) ( r/sqrt(pi) ) * sqrt( theta(x) - ( sin(2*theta(x))./ 2 ) + (1/3)*sin(theta(x)).^3 );
        yNos = Haack(xNos);
    elseif strcmp(rocket.parafoil.noseType, 'OGIVE')
        rho = (r^2 + Lnos^2)/(2*r);
        Ogive = @(x) sqrt(rho^2 - (Lnos - x).^2) + r - rho;
        yNos = Ogive(xNos);
    elseif strcmp(rocket.parafoil.noseType, 'POWER')
        power = rocket.parafoil.nosePower;
        Power = @(x) r * (x/Lnos).^(power);
        yNos = Power(xNos);
    elseif strcmp(rocket.parafoil.noseType, 'MHAACK')
        p = rocket.parafoil.nosePMod;
        c = rocket.parafoil.noseCMod;
        xMod = @(x, p) (x/Lnos).^(p)*Lnos;
        thetaMod = @(x, p) acos( 1 - ( (2*xMod(x, p))./Lnos ) );
        haackSeriesMod = @(x, p, C) ( r/sqrt(pi) ) * sqrt( thetaMod(x, p) ...
            - ( sin(2*thetaMod(x, p))./ 2 ) + C*sin(thetaMod(x, p)).^3 );
        yNos = haackSeriesMod(xNos, p, c);
    end

    if strcmp(rocket.rear.boatType, 'OGIVE') % [-] Boat-tail shape. 0: Cone, 1: Tangent Ogive
        [xBoat, yBoat] = computeTangentBoatPoints(2*r, Daft, Laft);
    else
        xBoat = [0 Laft];
        yBoat = [r Daft/2];
    end

    %%% figure begin
    roc = figure();
    hold on

    %%% NOSECONE
    plot(xNos, yNos, 'k');
    plot(xNos, -yNos, 'k');
    plot([Lnos Lnos], [-r r], 'k');

    %%% CENTERBODY
    plot([Lnos Lnos+Lcent], [r r], 'k');
    plot([Lnos Lnos+Lcent], [-r -r], 'k');
    plot([Lnos+Lcent Lnos+Lcent], [-r r], 'k');

    %%% BOAT-TAIL PLOT
    plot(xBoat+Lnos+Lcent, yBoat, 'k');
    plot(xBoat+Lnos+Lcent, -yBoat, 'k');
    plot([Lnos+Lcent+Laft Lnos+Lcent+Laft], [-Daft/2 Daft/2], 'k');

    %%% FINS PLOT
    % top
    plot([Xle1 Xle1+deltaXLE], [r r+height],'k'); 
    plot([Xle1+C1 Xle1+deltaXLE+C2], [r r+height],'k');
    plot([Xle1+deltaXLE Xle1+deltaXLE+C2], [r+height r+height],'k');
    plot([Xle1 Xle1+C1], [r r], 'k');
    % bottom
    plot([Xle1 Xle1+deltaXLE], [-r -r-height],'k');
    plot([Xle1+C1 Xle1+deltaXLE+C2], [-r -r-height],'k');
    plot([Xle1+deltaXLE Xle1+deltaXLE+C2], [-r-height, -r-height],'k');
    plot([Xle1 Xle1+C1], [-r -r], 'k');

    % plot([Xle1 Xle1+deltaXLE], [r+height r+height], 'k--')
    % plot([Xle1+deltaXLE Xle1+deltaXLE], [r+height r],'k--')

    %%% BAYS
    lengths = rocket.absolutePositions.values + 2*Lnos;
    bays = rocket.absolutePositions.keys;
    xline(0,'r-.', 'Label', 'NOSE', ...
        'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ...
        'LabelHorizontalAlignment', 'right');
    for i = 1:length(lengths)
        xline(lengths(i), '--', 'Label', bays(i), ...
        'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ...
        'LabelHorizontalAlignment', 'right');
    end
    boatPlotLength = rocket.length - rocket.rear.boatLength;
    xline(boatPlotLength,'r-.', 'Label', 'BOAT-TAIL', ...
        'LabelOrientation', 'aligned', 'LabelVerticalAlignment', 'bottom', ...
        'LabelHorizontalAlignment', 'right');
    
    axis equal
    % set(gca, 'xlim', [-100, Lnos+Lcent+Laft+100]);

    %%% TITLE
    title(string(strrep(mission.name, '_', ' ')));
end
