You forgot to put a failsafe in your code so most likely t(counter) never exceeds P. Try it this way
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
%define input variables
qs=1400; %solar input
qe=240; %earth radiation
ar=0.3; %earth albedo
sigma=5.67*10^(-8); %boltzmann constant
R_e=6371000; %radius of earth(m)
Mu=3.986*10^14; %gravitational parameter for earth
h=1336000; %altitude
alpha=0.46; %absorbance of spacecraft
epsilon=0.77; %emittance of spacecraft
m=1000; %mass of spacecraft
P_sun=1000; %power dissipated in sunlight
P_eclipse=500; %power dissipated in eclipse
A_surf=3.35; %surface area of single face
A_total=20.1; %total surface area
cp=900; %specific heat capacity
Ti=263; %initial temperature
R_R=h/R_e;
orbit_period=2*pi*sqrt((R_e+h)^3/(Mu)); %time period of orbit
P=5*orbit_period
%calculate time in eclipse
te=((2*pi*(h+R_e)^(3/2))/sqrt(Mu))*(asin(R_e/(R_e+h))/pi);
%calculate thermal timeconstant
tt=(m*cp)/(sigma*epsilon*alpha*A_surf);
%calculate view factor
F12=1/(1+R_R)^2;
%calculate input
Qin_sun=(qs*alpha*A_surf)+(ar*qs*alpha*F12*A_surf)+(qe*epsilon*F12*A_surf)+P_sun;
Qin_eclipse=qe*A_surf*epsilon;
%calculate equillibrium temperatures
Teq_sun=(Qin_sun/(epsilon*sigma*A_total))^(1/4);
Teq_ecl=(Qin_eclipse/(epsilon*sigma*A_total))^(1/4);
%given values
dt=1; %steps
t=0:dt:P; %time array
T=zeros(size(t)); %temperature array
counter = 1;
% Demonstration of how to avoid an infinite loop by setting up a failsafe.
% Set up a failsafe
maxIterations = 34000; % Way more than you think it would ever need.
loopCounter = 0;
% Now loop until we obtain the required condition: t(counter) > P.
% If that never happens, the failsafe will kick us out of the loop so we do not get an infinite loop.
while (t(counter)<=P) && counter < maxIterations && counter < numel(t)
if rem(t(counter),orbit_period)<=te
T(counter+1) = T(counter)+dt*(Teq_ecl^4-T(counter)^4)/tt;
else
T(counter+1) = T(counter)+dt*(Teq_sun^4-T(counter)^4)/tt;
end
counter = counter + 1;
fprintf('Iteration #%d. t(%d) = %f.\n', counter, counter, t(counter))
end
% Alert user if we exited normally, or if the failsafe kicked us out to avoid an infinite loop.
if counter < maxIterations
% Then the loop found the condition and exited early, which means normally.
fprintf('Loop exited normally after %d iterations.\n', counter);
else
% Then the loop never found the condition and exited when the number of iterations
% hit the maximum number of iterations allowed, which means abnormally.
fprintf('Loop exited abnormally after iterating the maximimum number of iterations (%d) without obtaining the exit criteria.\n', maxIterations);
end
fprintf('All done after %d iterations.\n', counter)
%Plot the results
figure;
plot(t / orbit_period, T, 'LineWidth', 2);
hold on;
plot([te / orbit_period, te / orbit_period], [min(T), max(T)], '--r', 'LineWidth', 2);
xlabel('Orbits');
ylabel('Temperature (K)');
title('Satellite Temperature Variation Over Time');
legend('Temperature', 'Eclipse Start');
grid on;
hold off;
Actually since t is just a vector
t = 0 : dt : P; %time array
you should be using a for loop instead of a while loop
for counter = 1 : numel(t)
