I have a problem with plotting two vectors, I want to vary 'z' by a certain amount and 'x' by another amount and plot versus each other, how can i make the matrix dimensions match?

x = 2:1:25;
y = x^2;

z = 1:1:28;
plot(z,x)

[the cyclist]

I'm guessing that you mean

y = x.^2;

See this documentation page for details on the difference.

Answer 1

[madhan ravi]

Read about meshgrid() and ndgrid().

[Ikenna Okoye]

Thanks.

[madhan ravi]

Anytime :)

Answer 2

[the cyclist]

Well, it's really up to you to decide. You have a 24-element vector, and a 28-element vector.

Which of those elements align with each other? The first 24 elements of both? Then you could just drop the last 4 elements of z.

Or maybe the two endpoints align, and you need to linearly interpret both vectors onto the same number of points?

It's really not clear from your question. I think you need to think more carefully about what makes sense.

[Ikenna Okoye]

Im trying to see if there is some variable x that varies, variable y that varies and variable z(x,y), that can take both vectors with different lengths and make them the same size in a code. I want to apply this to a more complex code.

[the cyclist]

Yes, there are an infinite number of ways to make these two vectors the same length. Here are three:

Method 1:

x = [x zeros(1,4)]; % Append 4 zeros to the end of x.

Method 2:

z = z(1:24); % Use only the first 24 elements of z.

Method 3:

x = [x nan(1,4004)]; % Append 4004 NaNs to the end of x.
z = [z nan(1,4000)]; % Append 4000 NaN to the end of z.

(The last one is almost certainly not useful, but who knows?)

My point is that the best method will depend on the relationship between the current elements of x and z. How do they "line up" with each other? That relationship is impossible to understand without more detailed information.

Answer 3

[Ikenna Okoye]

Heres my code: It plots thrust versus altitude and varies 'gamma'. I want to be able to create a carpet plot varying 'gamma', 'ex' and one other parameter i cannot remember right now. It runs fine, but as soon as i vary anything except altitude 'z', the error ''matrix dimentions do not agree'' pops up. The code has three parts: the main thrust code (Propulsion_rev3) and a funciton that it calls for altitude perameters , the code of atmosphere model (atmos0_100) and the function that calls the atmosphere model (Fun_atmos0_100).

%Propulsion_rev3
z = 0:1:100;
[T,P,rho,del,sig,the]= atmos0_100(z);
 %% First stage: constants       
    r_o = 6.371e7;                                                          % earth radius, m
    g_o = 9.81;                                                             % gravity, m/s^2
    D = 100;                                                                % Drag, N, from aero
    m_i = 12000;                                                            % initial wet mass, kg
    P_c = 2.2339e7;                                                         % ambient pressure, Pa
    g = 9.81;                                                               % gravitational acceleration, m/s^2
    f = 6;                                                                  % fuel to air ratio
    I_sp = 363;                                                             % specific impulse, s
    A_e = 4.104;                                                            % exit area, m^2
    % From App.B
    T_c = 3400;                                                             % temp, K
    gamma = 1.21;                                                           % specific heat ratio
    Mol = 13.5;                                                             % molecular mass, kg/kmol
    v_o = 10;                                                               % inicial velocity, m/s
    time = 80;                                                              % time to seperation, s, from performance
    R = (8414/(1+f))*((1/28.97)+(f/Mol));                                   % gas constatnt, J/kg*K
    cp = (R*gamma)/(gamma-1);                                               % specific heat, J/kg*K
    
%% Inputs
    F_req = 1334466;                                                        % required thrust, N (from Performance)

%% Engine characteristics/SSME design point    
% assuming combustion efficiency of 1.0
    ex = 77;
    P_c = 2.2339e+7;                                                        % chamber pressure, Pa
    A_t = (ex.^-1).*A_e;                                                       % throat area, m^2
    
     epsilon = sqrt(gamma./(((gamma+1)./2).^((gamma+1)./(gamma-1))));    
    
% Newton-Raphson method for finding exit mach number    
    b = 2.*(gamma-1)./(gamma+1);
    c = ((gamma+1)./(gamma-1)).*((ex).^(b));
    d = 2./(gamma-1);
    
    Me_0 = 7.5;
    original = (c .* (Me_0 .^ b)) - (Me_0 .^ 2) - d;
    derivative = (b .* c .* (Me_0 .^ (b - 1))) - (2 .* Me_0);
    
     Me_1 = Me_0 - (original ./ derivative);
    
    while abs(original) > 10^-2
        Me_0 = Me_1;
        
        original = (c .* (Me_0 .^ b)) - (Me_0 .^ 2) - d;
        derivative = (b .* c .* (Me_0 .^ (b - 1))) - (2 .* Me_0);
        Me_1 = Me_0 - (original ./ derivative);
        
        fprintf('Iteration: M_e = %.20f, error = %.20f \n', Me_1, original);
    end
    
    mf_f = 1-exp(-(sqrt(g_o*r_o))./(I_sp*v_o*((1-D)/F_req)));               % fuel mass fraction
    mf = mf_f *m_i;                                                         % mass of fuel, kg,

T_e = T_c./((1 + gamma.*R.*Me_1.^2)./(2.*cp));
P_e = P_c.*(1+((gamma-1)./2).*(Me_1.^2)).^(gamma./(1-gamma));               % exit static pressure, Pa (pg.144, Mattingly)
P_te = P_e./(1+((gamma-1)./(Me_1.^2))).^(-gamma./(gamma-1));                % exit total pressure, Pa
 
% using isentropic relations
% T_e = ((P_e./P_te).*T_c.^(gamma./(gamma-1))).^(1./(gamma./(gamma-1)));
    v_e = Me_1.*sqrt(gamma.*R.*T_e);
    a = (2*gamma)/(gamma-1);
    coeff = gamma/(gamma-1);
    Pr = ((2/(gamma+1))*(1+((gamma-1)/2))^(-gamma/(gamma-1)));
    
        m_dot = ((P_c*A_t)/sqrt(T_c))*sqrt((2/R)*coeff*((Pr^(2/gamma))-(Pr^((gamma+1)/gamma))));
        F_ava =((m_dot.*v_e) + (P_e-P).*(A_t.*ex));
        S = (m_dot.*f)./F_ava;                                              % specific fuel consumption, kg/s*N
        f_t = (S.*F_ava).*time;
        C_f = F_ava/(P_c*A_t);
        R_thrusttoweight = (F_ava)/(f_t*9.81);
    %C_f = epsilon*sqrt(a*(1-(P_e/P_c)^coeff))+(P_e/P_c - P/P_c)*ex;
    
%% engine sizing   
    D_e = 0.00357.*F_ava+14.48;                                              % engine diameter, m
    m_e = F_ava./(g*0.0006098.*F_ava+13.44);                                 % engine mass, kg
    L_e = 0.000030*F_ava+327.7;                                              % engine length, m
    Lstar = 1;                                                               % characteristic length, m
    A_c = A_t*(8*((2*sqrt(A_t/pi))^(-0.6)) + 1.25);                          % average combustion chamber cx area, m^2
    L_c = Lstar*(A_t/A_c);                                                   % combusiton chamber length, m            
    V_c = Lstar*A_t;                                                         % combustion chamber volume, m^3
    d_e  = sqrt((4.*ex.*A_t)/pi);                                            % exit diameter, m

%% Tank sizing
V_f = (f_t/70.8)*(1/7);                                                      % fuel volume, m^3
V_o = (f_t/1141)*(6/7);                                                      % oxidizer volume, m^3
V_t = V_f + V_o;                                                             % total propellant volume, m^3

%% Plots/outputs

    fprintf('m_dot = %f [kg/s] \n',m_dot)
    fprintf('v_e = %f [m/s] \n',v_e)
    fprintf('V_c = %f [m^3] \n,',V_c)
    fprintf('A_c = %f [m^2] \n', A_c)
    fprintf('m_e = %f [kg] \n',m_e)
    fprintf('L_e = %f [m] \n', L_e)
    fprintf('D_e = %f [m]',D_e)
    fprintf('L_c = %f [m] \n', L_c)
    fprintf('f_t = %f [kg]\n', f_t)
%   fprintf('m_f = %f [kg]\n', m_f)
 
    figure(1)
    plot(z,F_ava,'linewidth',1.25)
    title('Avalible thrust vs. Altitude')
    xlabel('Altitude (km)')
    ylabel('F_ava (N)')
    
    figure(2)
    plot(z,S,'linewidth',1.25)
    title('Specific fuel consimption vs. Altitiude')
    xlabel('Altitude (km)')
    ylabel('S (kg/s*N)')
    
    figure(3)
    plot(z,C_f,'linewidth',1.25)
    title('Thrust Coefficient vs. Altitiude')
    xlabel('Altitude (km)')
    ylabel('C_f')
    
%% Trade Studies
% plot(F_ava,gamma)  % vary gamma from 1.2-1.4

%% Second stage: Constants
r_o = 6.371e7;                                                          % earth radius, m
                                                             % gravity, m/s^2
    D = 100;                                                                % Drag, N, from aero
    m_i = 12000;                                                            % initial wet mass, kg
    P_c = 2.2339e7;                                                         % ambient pressure, Pa
    g = 9.81;                                                               % gravitational acceleration, m/s^2
    f2 = 6;                                                                  % fuel to air ratio
    I_sp = 363;                                                             % specific impulse, s
    A_e2 = 4.104;                                                            % exit area, m^2
    % From App.B
    T_c2 = 3400;                                                             % temp, K
    gamma2 = 1.21;                                                           % specific heat ratio
    Mol2 = 13.5;                                                             % molecular mass, kg/kmol
    R2 = (8414/(1+f2))*((1/28.97)+(f2/Mol2));                                   % gas constatnt, J/kg*K
    cp2 = (R*gamma2)/(gamma2-1);      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fun_atmos0_100
prompt = 'For SI select 1 for imperial select 2?';
unit = input(prompt);
if unit==1
prompt = 'What is your height value in kilometers?';
z = input(prompt);
else unit==2
prompt = 'What is your height value in feet';
zi = input(prompt);
z = convlength([zi],'ft','km');
end  
[T,P,rho,del,sig,the]= atmos0_100(z)
fprintf('For altitude of   z = %.3f [km]     ',z);

fprintf('        SI                    Imperial                          Ratios   \n');
fprintf('------------------------------------------------------------------------- \n');
fprintf('P =   %.3f [Pa]                                                del = %.5f  \n',P,del);
fprintf('T =   %.3f [K]                                                 the = %.5f  \n',T,the);
fprintf('rho = %f [kg/m^3]                                              sig = %f  \n',rho,sig);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%atmos0_100
function[T,P,rho,del,sig,the]= atmos0_100(z)
global r R_o
%% Constants
% Sea level conditions
P_o   = 101325; % (Pa)
rho_o = 1.225; % (Kg/m^3)
T_o   = 288.15; % (K)
mu_o  = 1.789*10^-5; % (Kg/m/s)
nu_o  = 1.46*10^-5; % (m^2/s)
tc_o  = 0.02596; % (W/m/K)
R_o   = 287.053; % (J/Kg/K)
g_o   = 9.80665; % (m/s^2)
Cp    = 1005; % (J/Kg/K)
r_o   = 6356.766; % (km) 
r     = 1.4;
%% Process
for i = 1:length(z)    
  
h(i) = (r_o*z(i))/(r_o+z(i));

if h(i) == 0
    P(i) = P_o;
    T(i) = T_o;
  rho(i) = rho_o;
  
elseif 0<h(i) && h(i)<=11
    P(i) = P_o*(T_o/(T_o-6.5*h(i)))^(34.1632/-6.5);
    T(i) = T_o-6.5*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = z(i)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 11<h(i) && h(i)<=20
    T(i) = 216.65;
    P(i) = 22632.06*exp(-34.1632*(h(i)-11)/T(i));
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-11)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 20<h(i) && h(i)<=32
    P(i) = 5474.889*(216.65/(216.65+(h(i)-20)))^34.1632;
    T(i) = 196.65+h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-20)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
 
elseif 32<h(i) && h(i)<=47
    P(i) = 868.0187*(228.65/(228.65+2*r*(h(i)-32)))^(34.1632/(2*r));
    T(i) = 139.05+2*r*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-32)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 47<h(i) && h(i)<=51
    T(i) = 270.65;
    P(i) = 110.9063*exp(-34.1632*(h(i)-47)/T(i));
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-47)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 51<h(i) && h(i)<=71
    P(i) = 66.93884*(270.65/(270.65-2*r*(h(i)-51)))^(34.1632/(-2*r));
    T(i) = 413.45-2*r*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-51)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 71<h(i) && h(i)<=85
    P(i) = 3.956420*(214.65/(214.65-2*(h(i)-71)))^(34.1632/-2);
    T(i) = 356.65-2*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-71)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 85<z(i) && z(i)<=91
    T(i) = 186.8673;
    P(i) = exp(2.159582*10^(-6)*z(i)^3 + -4.836957*10^(-4)*z(i)^2 + -.14251928*z(i)  +  13.47530);
  rho(i) = exp(-3.322622*10^(-6)*z(i)^3 + 9.111460*10^(-4)*z(i)^2  +  -.2609971*z(i)  +  5.944694);
  eps(i) = (z(i)-85)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
else 
    T(i) = 263.1905-76.3232*sqrt(1-((z(i)-91)/-19.9429)^2);
    P(i) = exp(3.304895*10^-5*z(i)^3  +  -.009062730*z(i)^2  +  .6516698*z(i)  +  -11.03037);
  rho(i) = exp(2.873405*10^-5*z(i)^3  +  -.008492037*z(i)^2  +  .6541179*z(i)  +  -23.6201);
  eps(i) = (z(i)-85)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
end
end

[the cyclist]

This is now almost too much detail. :-)

Can you give an example of something you are trying to change here, and then give the full text of the error message you are getting as a result?

It would be easier for users here to help if you removed all parts of the code that are not directly related to the problem you are having. Then you could post a much smaller piece of code that shows the error.

Also, it would have been better if you had posted this as a comment rather than an answer, since it is not an answer.

[Ikenna Okoye]

Gotcha and gotcha, didnt know it was in the answers section.

% im trying to vary 'gamma' dimentions wont agree with 'F_ava', avalible thrust
gamma = 1.2:0.05:1.4;  
ex = 60:4:76;
    
T_e = T_c./((1 + gamma.*R.*Me_1.^2)./(2.*cp));
P_e = P_c.*(1+((gamma-1)./2).*(Me_1.^2)).^(gamma./(1-gamma));               % exit static pressure, Pa (pg.144, Mattingly)
P_te = P_e./(1+((gamma-1)./(Me_1.^2))).^(-gamma./(gamma-1));                % exit total pressure, Pa
 
% using isentropic relations
% T_e = ((P_e./P_te).*T_c.^(gamma./(gamma-1))).^(1./(gamma./(gamma-1)));
    v_e = Me_1.*sqrt(gamma.*R.*T_e);
    a = (2*gamma)/(gamma-1);
    coeff = gamma/(gamma-1);
    Pr = ((2./(gamma+1)).*(1+((gamma-1)./2)).^(-gamma./(gamma-1)));
    
        m_dot = ((P_c*A_t)./sqrt(T_c)).*sqrt((2/R).*coeff.*((Pr.^(2./gamma))-(Pr.^((gamma+1)./gamma))));
        F_ava =((m_dot.*v_e)+ (P_e-P).*(A_t.*ex));

error message:

            Error in Pro_rev3 (line 74)
                    F_ava =((m_dot.*v_e)+ (P_e-P).*(A_t.*ex)); Matrix dimensions must agree.

Most of the equaions depend on 'gamma' and that is the variable i want to vary. I want to do the same with 'ex' expansion ratio also to make a carpet plot similar to the image attatched.

[the cyclist]

So, you are talking about a pretty major change in your code. You want to take something that was written as a function of one variable, and change it to now be a function two (or maybe three?) variables.

madhan ravi's idea that you will need to use meshgrid is likely part of the solution. Given input of two vectors, that function will create a 2-d grid of all possible combinations of the elements of them.

Then you would use the output of that to build the functions of those two variables. You will probably need to rewrite a lot of lines of code.