code a solar still in matlab

% Solar still simulation code
% Adapted from L.W. Finger, "The Design of Solar Still", Barlow
% Press, San Francisco, 1972

clear; clc;

%% Input parameters
y = 10; % Depth of water [cm]
A = 1; % Area of still [m^2]
H = 0.3; % Height of sloping cover [m]
n = (y+H)/(2*H); % Ratio of water depth to height
Ta = 35; % Ambient temperature [deg C]
Ts0 = 80; % Initial temperature of water [deg C]

%% Calculation
Gsc = 1367; % Solar constant [W/m^2]
solar_time = 0:0.5:24; % Solar time [hours]
Ts = Ts0*ones(size(solar_time)); % Initialize output temperature vector
for i = 1:length(solar_time)
    % Solar radiation
    B = 2*pi*(solar_time(i)-12)/24; % Hour angle [rad]
    delta = 23.45*sin(2*pi*(284+solar_time(i))/365); % Declination angle [deg]
    phi = 35.5; % Latitude of still [deg] 
    Hs = acosd(-tand(phi)*tand(delta)); % Hour angle of sunrise/sunset [deg]
    if B <= -Hs || B >= Hs % Sun not in sky
        Qs = 0;
    else
        cos_theta = sind(phi)*sind(delta) + cosd(phi)*cosd(delta)*cosd(B);
        I = Gsc*cosd(phi)*cosd(delta)*sind(Hs)*(sind(Hs)-sind(phi)*sind(delta*cosd(B)));
        Qs = I*A*cos_theta; % Solar radiation on still [W]
    end
    
    % Convective and conductive heat loss
    h = 2.35 + 5.5*n; % Average film coefficient [W/m^2*K]
    Ql = h*A*(Ts(i)-Ta); % Convective heat loss [W]
    Qk = 19.4*A*(Ts(i)-Ta)/y; % Conductive heat loss [W]
    
    % Energy balance
    Qe = Qs - Ql - Qk; % Energy input to still [W]
    m = Qe/(880*(Ts(i)-Ta)); % Mass of water evaporated [kg/s]
    Ts(i+1) = Ts(i) - m*2256/(1000*A); % Final temperature of water [deg C]
end

%% Plot
plot(solar_time,Ts(1:end-1),'LineWidth',2)
xlabel('Solar time [hours]')
ylabel('Temperature of water [deg C]')
title('Temperature of water in a solar still')
grid on
1777 chars
50 lines