Index exceeds the number of array elements (5).

Hi, I need help how to solve this issue at line67 if(Af(i)<0). When I change for i =[1,2,3,4,5,6] to for i=1:m-1, the problems still persist. Please enlighten me.

clc; 
close all; 
clear all;

% Constants
rGEO  = 42164e3;            % [m] Radius of the GEO orbit    
rE    = 6378e3;             % [m] Radius of Eearth
mu    = 3.986e14;           % {m^3/kg] Gravitational parameter of Earth

sigma = 5.66961e-8;         % [W/m^2/K^4] Stefan-Boltzmann constant
S     = 1372;               % [W/m^2] Solar constant

% Material parameters
c_AL = 896;             % [J/kg/K] Heat capacity of Aluminium
r_AL = 2.7e3;           % [kg/m^3] Density of Aluminium
a_AL = 0.15;            % [-] Absorption of Aluminium
e_AL = 0.05;            % [-] Emissivity of Aluminium    
k_AL = 229;             % [W/m/K] Heat transfer of Aluminium    

% Satellite configuration
e_ik = 0.2;         % [-] View factor
l    = 2;           % [m] Length of plate
d    = 0.02;        % [m] Thickness of plate
mk   = l^2*d*r_AL;  % [kg] Mass of plate

% GEO Orbit
theta = deg2rad(23.14);             % [rad] Sun height in summer
omega_s = sqrt(mu/rGEO)/rGEO;      % [rad/s] Orbital angular rotation speed of the satellite

% Simulation parameters
dt = 60;                % [s]  Step size
n  = 4;                 % [-]  Number of orbits
tE = n*2*pi/omega_s;           % [s]  Duration of simulation for 90days

% -------------------------------------------------------------------------
%% Initialization
phi = 0.1;                  % [rad] Anomaly of the satellite
t = 0;                      % [s] Start time
m=1;                        % Counter [-]
Tsp(m,:)=273*ones(1,6);     % [K] Initial temperature for simulation in spring 
Tsu(m,:)=273*ones(1,6);     % [K] Initial temperature for simulation in summer 
% -------------------------------------------------------------------------
%%Eclipse period
alpha=2*asin(rE/rGEO);
% -------------------------------------------------------------------------
%% Main loop
while t<tE
    %% Calculation of temperature with the Euler met
    m=m+1;
    %%The projected areas for the nodes with respect to the incident Sunlight are
Af(1)=l^2*cos(phi);                           %Spring node 1
Af(2)=0;                                      %Spring node 2
Af(3)=l^2*sin(phi);                           %Spring node 3
Af(4)=-l^2*sin(phi);                          %Spring node 4
Af(5)=0;                                      %Spring node 5
As(6)=-l^2*cos(phi);                          %Spring node 6

As(1)=l^2*cos(phi)*cos(theta);                %Summer node 1
As(2)=l^2*sin(theta);                         %Summer node 2
As(3)=l^2*sin(phi)*cos(theta);                %Summer node 3
As(4)=-l^2*sin(phi)*cos(theta);               %Summer node 4
As(5)=0;                                      %Summer node 5
As(6)=-l^2*cos(phi)*cos(theta);               %Summer node 6
    
%Applying Loop over all nodes
for i =[1,2,3,4,5,6]
%For negative areas both spring and summer, we need to make them 0 instead
%of absolute values
if(Af(i)<0)
    Af(i)=0;
end
if(As(i)<0)
    As(i)=0;
end
%Emitted energy from Sun in summer
Qs_sun = S * a_AL * As(i);
%Emitted Energy from sun in spring
if (phi>(pi-alpha/2) && phi <(pi+alpha/2))
    Qf_sun=0;
else
    Qf_sun = S * a_AL*Af(i);
end
%Heat emission is conducted on both sides
%Heat emission from all plates are proportional to T to the power of 4
Qf_em = -e_AL * sigma * 2 * l^2 * Tsp(m-1,i)^4;                                     %[W] Heat emission during spring
Qs_em = -e_AL * sigma * 2 * l^2 * Tsu(m-1,i)^4;                                     %[W] Heat emission during summer
Qf_ab =0.0;
Qs_ab =0.0;
Qf_HT =0.0;
Qs_HT =0.0;
for k=[1,2,3,4,5,6]
    %Summation of the absorbtion from all other plates
    if (i ~= k)
        Qf_ab = + Qf_ab + a_AL * e_ik * (e_AL * sigma * Tsp(m-1,k)^4) * l^2;
        Qs_ab = + Qs_ab + a_AL * e_ik * (e_AL * sigma * Tsu(m-1,k)^4) * l^2;
    end 
    %Summation of the heat conduction
    if (k + i ~= 7)
        Qf_HT = Qf_HT - k_AL * d * (Tsp(m-1,i) - Tsp(m-1,k));
        Qs_HT = Qs_HT - k_AL * d * (Tsu(m-1,i) - Tsu(m-1,k));
    end
end
%Total temperature of spring
        Tsp(m,i) = Tsp(m-1,i) + (Qf_em+Qf_ab+Qf_HT+Qf_sun)*dt/(mk*c_AL);
%Total temperature of summer
        Tsu(m,i) = Tsu(m-1,i) + (Qs_em+Qs_ab+Qs_HT+Qs_sun)*dt/(mk*c_AL);
end
%There are two loops of variables i and t. It is needed to end them
%respectively. First loop is i loop

phi = phi + omega_s * dt;

if phi > 2*pi
    phi = phi - 2*pi;
end
t = t + dt;
end %while
% -------------------------------------------------------------------------
%% Plot of temperature curves
figure('Position',[100 100 1000 800],'Name','Node-model');

subplot(2,1,1)
hold on;
plot(Tsp(:,1),'k'); 
plot(Tsp(:,2),'b');
plot(Tsp(:,3),'y');
grey=plot(Tsp(:,4)); set(grey,'Color',[0.5 0.5 0.5]);
plot(Tsp(:,5),'g');
plot(Tsp(:,6),'r');
title('Temperature Profile of 6 nodes on Spring');xlabel('Time [min]');ylabel('Temperature [K]')
grid on; legend('Node1','Node2','node3','Node4','Node5','Node6');

subplot(2,1,2)
hold on;
plot(Tsu(:,1),'k');
plot(Tsu(:,2),'b');
plot(Tsu(:,3),'y');
grey=plot(Tsu(:,4)); set(grey,'Color',[0.5 0.5 0.5]);
plot(Tsu(:,5),'g');
plot(Tsu(:,6),'r');
title('Temperature Profile of 6 nodes on Summer');xlabel('Time [min]');ylabel('Temperature [K]')
grid on;legend('Node1','Node2','Node3','Node4','Node5','Node6');

Answer 1

[Tommy]

Likely a typo in line 56:

 As(6)=-l^2*cos(phi);                          %Spring node 6
% ^ should be Af?

As it is, Af only has 5 values.

[Hafizuddin Bin Mohd Lowhim]

Yes thank you! now solved!