Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Hi. First of all sorry for asking the same question again and again but I am really close to the solution (probably).

Please do not mind long lines because most of them are constants and loops.

This code tries to solve 6 ODEs with 6 state variables [horizontal position (x1 and x2), altitude (x3), the true airspeed (x4), the heading angle (x5) and the mass of the aircraft (x6)] and 3 control inputs [engine thrust (u1), the bank angle (u2) and the flight path angle (u3)] by using RK4.

Different flight maneuvers are performed for the specified time intervals.

Velocities.m, Cruise_Vel.m, Des_Vel.m, Thr_cl.m, Thr_cr.m, Thr_des.m, fuel_cl.m, fuel_cr.m, fuel_des.m,den.m,den_cr.m,den_des.m,drag.m,drag_cr.m,drag_des.m,lift.m,lift_cr.m,lift_des.m are functions in seperate tabs.

File for equations:

function dydt = f_1(t,x,p,Cl,C_D,f,u1,u2,u3)
dydt = zeros(6,length(t));
% Aircraft Properties
W = .44225E+06;                             % .44225E+03 tons = .44225E+06 kg 
S = .51097E+03;                             % Surface Area [m^2]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]
% Initial conditions
dydt(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
dydt(1) = x(4) * cos(x(5)) * cos(u3);       
dydt(2) = x(4) * sin(x(5)) * cos(u3);       
dydt(3) = x(4) * sin(u3);               
dydt(4) = -C_D*S*p*(x(4)^2)/(2*x(6))-g0*sin(u3)+u1/x(6); 
dydt(5) = -Cl*S*p*x(4)./(2*x(6))*sin(u2);
dydt(6) = -f;
end

Code in the main.m is:

function main
W = .44225E+06;                             % .44225E+03 tons = .44225E+06 kg 
% Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
T = 0:60;
C =[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];

x(4) = Velocities(3608.9,5249.3);                 % Changing speed [m/s] line 10 %%%%%%%%%%%%%%%%%%%%
x(5) = 0;                                         % Changing head angle [deg]
f = fuel_cl(3608.9,5249.3);                       % Changing fuel flow [kg/min] 

u1 = Thr_cl(3608.9,5249.3);                       % Changing thrust [N]
u2 = 0;                                           % Changing bank angle [deg]
u3 = 5;                                           % Changing flight path angle [deg]
V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
C_D = drag(3608.9,5249.3,x(4));                   % Changing drag coefficient
Cl = lift(3608.9,5249.3,x(4));                    % Changing lift coefficient
p = den(3608.9,5249.3);                           % Changing density [kg/m3]

[t1,x1] = ode45(@(t,x)f_1(t,x,p,Cl,C_D,u1,u2,u3),T,C);

% Perform cruise flight for t=60 minutes.
T = 60:3660;
C = x1(end,:);
x(4) = Cruise_Vel(5249.3);                        % Changing speed [m/s]
x(5) = 0;                                         % Changing head angle [deg]
f = fuel_cr(5249.3);                              % Changing fuel flow [kg/min] 

u1 = Thr_cr(5249.3);                              % Changing thrust [N]
u2 = 0;                                           % Changing bank angle [deg]
u3 = 0;                                           % Changing flight path angle [deg]
V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,x(4));                       % Changing drag coefficient
Cl = lift_cr(5249.3,x(4));                        % Changing lift coefficient
p = den_cr(5249.3);                               % Changing density [kg/m3]

[t2,x2] = ode45(@(t,x)f_1(t,x,p,Cl,C_D,u1,u2,u3),T,C);

% Turn with β=30 bank angle until heading is changed by η=270◦.
T =3660:3720;
C = x2(end,:);
x(4) = Cruise_Vel(5249.3);                        % Changing speed [m/s]
x(5) = 0:30:270;                                  % Changing head angle [deg]
f = fuel_cr(5249.3);                              % Changing fuel flow [kg/min] 

u1 = Thr_cr(5249.3);                              % Changing thrust [N]
u2 = 30;                                          % Changing bank angle [deg]
u3 = 0;                                           % Changing flight path angle [deg]
V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,x(4));                       % Changing drag coefficient
Cl = lift_cr(5249.3,x(4));                        % Changing lift coefficient
p = den_cr(5249.3);                               % Changing density [kg/m3]

[t3,x3] = ode45(@(t,x)f_1(t,x,p,Cl,C_D,u1,u2,u3),T,C);

% Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
T = 3720:3780;
C = x3(end,:);
x(4) = Des_Vel(5249.3,3608.9);                    % Changing speed [m/s]
x(5) = 270;                                       % Changing head angle [deg]
f = fuel_des(5249.3,3608.9);                      % Changing fuel flow [kg/min] 

u1 = Thr_des(5249.3,3608.9);                      % Changing thrust [N]
u2 = 0;                                           % Changing bank angle [deg]
u3 = 4;                                           % Changing flight path angle [deg]
V_ver = x(4)*sin(u3);                             % Changing vertical speed [m/s]
C_D = drag_des(5249.3,3608.9,x(4));               % Changing drag coefficient
Cl = lift_des(5249.3,3608.9,x(4));                % Changing lift coefficient
p = den_des(5249.3,3608.9);                       % Changing density [kg/m3]

[t4,x4] = ode45(@(t,x)f_1(t,x,p,Cl,C_D,u1,u2,u3),T,C);

% Complete a 360◦ turn (loiter) at level flight.
T = 3780:3840;
C = x4(end,:);
x(4) = Cruise_Vel(3608.9);                          % Changing speed [m/s]
lon = [270 300 360 60 120 180 240 270];
x(5) = wrapTo360(lon);                            % Changing head angle [deg]
f = fuel_cr(3608.9);                              % Changing fuel flow [kg/min] 

u1 = Thr_cr(3608.9);                              % Changing thrust [N] 
u2 = 0;                                           % Changing bank angle [deg]
u3 = 0;                                           % Changing flight path angle [deg]
V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
C_D = drag_cr(3608.9,x(4));                       % Changing drag coefficient
Cl = lift_cr(3608.9,x(4));                        % Changing lift coefficient
p = den_cr(3608.9);                               % Changing density [kg/m3]
[t5,x5] = ode45(@(t,x)f_1(t,x,p,Cl,C_D,u1,u2,u3),T,C);

% Descent to h3=800 [m] with κ=4.5◦ flight path angle.
T = 3840:3900;
C = x5(end,:);
x(4) = Des_Vel(3608.9,2624.67);                   % Changing speed [m/s] 
x(5) = 270;                                       % Changing head angle [deg]
f = fuel_des(3608.9,2624.67);                     % Changing fuel flow [kg/min] 

u1 = Thr_des(3608.9,2624.67);                     % Changing thrust [N] 
u2 = 0;                                           % Changing bank angle [deg]
u3 = 4.5;                                         % Changing flight path angle [deg]
V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
C_D = drag_des(3608.9,2624.67,x(4));              % Changing drag coefficient
Cl = lift_des(3608.9,2624.67,x(4));               % Changing lift coefficient
p = den_des(3608.9,2624.67);                      % Changing density [kg/m3]

[t6,x6] = ode45(@(t,x)f_1(t,x,p,Cl,C_D,u1,u2,u3),T,C);

t = [t1;t2;t3;t4;t5;t6];
x = [x1;x2;x3;x4;x5;x6];
tot=cell2mat(f);                % Total fuel consumption during mission [kg/min]
Tot_fuel=sum(tot);
figure(1)
plot3(x(1),x(2),x(3));          % 3D position graph
figure(2)
plot(t,x(4));                   % Vtas − Time graph
figure(3)
plot(t,V_ver(:));               % V_vertical − Time graph
figure(4)
plot(t,x(5));                   % Heading − Time graph
figure(5)
plot(t,x(6));                   % Mass − Time graph
figure(6)
plot(t,u1(:));                  % Thrust − Time graph
figure(7)
plot(t,u2(:));                  % Bank Angle − Time graph
figure(8)
plot(t,u3(:));                  % Flight Path Angle − Time graph
fprintf('Total fuel consumption during mission is %.2f [kg]',Tot_fuel*tend/60);
end

And when I run:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.
Error in main (line 10)
x(4) = Velocities(3608.9,5249.3);                 % Changing speed [m/s]

The problem is, functions I use output multiple values and I try to store them in x(4), x(5), f, u1, C_D,Cl,p by using arrays. It is mandatory to update x(4), x(5), f, u1, C_D,Cl,p with changing altitude ( for example 3608.9 feet to 5249.3 feet )

The question is, how can I feed x(4), x(5), f, u1, C_D,Cl,p with these updating values without using arrays, so I won't face with size problems?

Thanks.

For example:

x(4) = Velocities(3608.9,5249.3);                 % Changing speed [m/s] line 10 %%%%%%%%%%%%%%%%%%%%
   C_D = drag(3608.9,5249.3,x(4));                     % Changing drag coefficient

and Velocities.m is:

 function [Vtas_cl] = Velocities(H1,H2)
%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%
Vcl_1 = 335;             % Standard calibrated airspeed [kt]
Vcl_2 = 172.3;           % Standard calibrated airspeed [kt] -> [m/s] (To find the  Mach transition altitude)
Vcl_2_in_knots = 335;    % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cl  = 0.86;            % Standard calibrated airspeed [kt]
K = 1.4;                 % Adiabatic index of air
R = 287.05287;           % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;           % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15;             % Standard atmospheric temperature at MSL [K]
g0 = 9.80665;            % Gravitational acceleration [m/s2]
a0= 340.294;             % Speed of Sound [m/s]
Vd_CL1 = 5;              % Climb speed increment below 1500 ft (jet)
Vd_CL2 = 10;             % Climb speed increment below 3000 ft (jet)
Vd_CL3 = 30;             % Climb speed increment below 4000 ft (jet)
Vd_CL4 = 60;             % Climb speed increment below 5000 ft (jet)
Vd_CL5 = 80;             % Climb speed increment below 6000 ft (jet)
CV_min = 1.3;            % Minimum speed coefficient 
Vstall_TO = .14200E+03;  % Stall speed at take-off [KCAS]
CAS_climb  = Vcl_2;   
Mach_climb = M_cl;
delta_trans = (((1+((K-1)/2)*(CAS_climb/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_climb^2)^(K/(K-1)))-1);    % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0);      % Temperature ratio at the transition altitude
H_p_trans_climb = (1000/0.348/6.5)*(T0*(1-teta_trans));   %  Transition altitude for climb [ft]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% constants
H_climb = H1:H2; %%%%%%%%%%% Input %%%%%%%%%%%%%%%%%%%%%%%%%
Vnom_climb_jet = zeros(1, length(H_climb));
for k1 = 1:length(H_climb)
    
    if (0<=H_climb(k1)&&H_climb(k1)<1500)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;
elseif (1500<=H_climb(k1)&&H_climb(k1)<3000)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
        
elseif (3000<=H_climb(k1)&&H_climb(k1)<4000)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL3; 
        
elseif (4000<=H_climb(k1)&&H_climb(k1)<5000)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL4;
        
elseif (5000<=H_climb(k1)&&H_climb(k1)<6000)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL5;
        
elseif (6000<=H_climb(k1)&&H_climb(k1)<10000)
        Vnom_climb_jet (k1)= min(Vcl_1,250);
        
elseif (10000<=H_climb(k1)&&H_climb(k1)<=H_p_trans_climb)
        Vnom_climb_jet(k1) = Vcl_2_in_knots;
        
elseif (H_p_trans_climb<H_climb(k1))
        Vnom_climb_jet(k1) = M_cl;
    end
Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514;   % [kn] -> [m/s]
H_climb (k1)= H_climb(k1) * 0.3048;         % [feet] -> [m]
K = 1.4;                                    % Adiabatic index of air
R = 287.05287;                              % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;                              % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0;                                 % Value of the real temperature T in ISA conditions [K]
T0 = 288.15;                                % Standard atmospheric temperature at MSL [K]
P0 = 101325;                                % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]
p0 = 1.225;                                 % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
T(k1) = T0 + deltaT + Bt * H_climb(k1);         % Temperature [K]
P (k1)= P0*((T(k1)-deltaT)/T0).^((-g0)/(Bt*R));     % Pressure [Pa]
p (k1)= P(k1) ./ (R*T(k1));                             % Density [kg/m^3]
% Using above values:
Vtas_cl(k1) = (2*P(k1)/visc/p(k1)*((1 + P0/P(k1)*((1 + visc*p0*Vcas_cl(k1)*Vcas_cl(k1)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
% Output
end
end

Answer 1

[KSSV]

This error occurs when you try to save ore number of elements in LHS than it is initialized for.

A = zeros(2) ; 
A(1,:) = rand(1,2) ;  % no error, as said two elements are being saved 
%
B = zeros(3) ; 
B(1,:) = rand(1,4) ;  % error, as B should have three elements but you are sabing 4 elements

Unable to perform assignment because the size of the left side is 1-by-3 and the size of the right side is 1-by-4.

In the same way check your initialization. Read about debugging the code and check the dimensions.

If you are not sure of the dimensions, initialize them into cell and store.

B = cell(2,1) ;
B{1} = rand(1,3) ;
B{2} = rand(1,4) ; 

[Turgut Ataseven]

Thanks for the response.

But the question is more about how should I modify the Velocities.m ( last codes in the question) so that it is possible to transfer outputs to main.m without using arrays? How can I update variables without arrays, in other words?

[Stephen23]

"How can I update variables without arrays, in other words?"

Basically everything in MATLAB is an array: https://www.mathworks.com/help/matlab/matlab_prog/fundamental-matlab-classes.html

As KSSV wrote, the issue seems to be your use of indexing or the choice of array type, not the use of arrays per se. Using arrays is what MATLAB is design for, so trying to use MATLAB "without arrays" is .... really trying to shoot yourself in the foot before running a race.