SOLACOOKA

function solar_cooker()
% Models the internal energy of the surface of a pot and the water inside
% of it over time, heated from the sun's solar energy.
    insolation = 1000; % W/m^2
    k = 55; % W/m/K average of cast iron conductivity
    radius_pot = 0.3; %m inner radius where water is at
    height_water = 0.2; %m
    height_pot = 0.3; %m
    density_pot = 1300; % kg / m^3
    d = 0.01; %m
    mass_pot = pi * (radius_pot + d)^2 * height_pot * density_pot;
    c_pot = 500; % J / (kg * K)
    % assuming the pot is a cylinder, short and stout with the lid on
    SA_pot_in = 2 * pi * radius_pot * height_pot + 2 * pi * radius_pot^2;
    SA_pot_out = 2 * pi * (radius_pot + d) * height_pot...
        + 2 * pi * (radius_pot + d)^2;
    mass_water = pi * radius_pot^2 * height_water; %density is 1 kg/m^3
    c_water = 4186; % J /(kg * K)
    Tenv = 293; % K
    Twater = 273 + 10; % CHANGE LATA
    e = 0.97; % reflection efficiency
    emmisivity = 0.65; % of cast iron
    sigma = 5.67e-8; %w/m^2/K Stephen Boltzmann's constnat
    h = k / d; % heat transef coefficient of pot

    function res = change_U_pot(t, params)
        % Params: internal energy of pot, internal energy of water
        temp_pot = params(1) / mass_pot / c_pot;
        temp_water = params(2) / mass_water / c_water;
        % Assume that the surface area of the cooker is 1 m^2
        % 1 - emmisivity percent of the radiated energy is reflected from
        % the pot. Therefore, only a certain percentage, corresponding to
        % emmisivity is actually absorbed into the pot and available for us
        % of conduction and convection and then radiation outward
        radiation = insolation * e * emmisivity;
        % convection depends on the outside surface area in contact with the
        % environment
        convection = h * SA_pot_out * (Tenv - temp_pot);
        % conduction depends on the surface area of pot touching the inside
        conduction = k * SA_pot_in / d * (temp_pot - temp_water)...
             * (temp_water < 373);
%          fprintf('Water temp: %d\nConduction amount: %d\n\n', temp_water,...
%              conduction)
        % whatever energy is absorbed into the pot is also emmitted through
        % radiation outwards
        emmission = emmisivity * sigma * SA_pot_out * temp_pot^4;
        res = [radiation - convection - conduction - emmission; ...
            conduction];
    end

    [t, u] = ode45(@change_U_pot, [0, 1000], [Tenv * mass_pot * c_pot, ...
        Twater * mass_water * c_water]);
    %disp(u(:, 1) / mass_pot / c_pot)
    %disp(u(:, 2) / mass_pot / c_pot)
    figure(1)
    plot(t, u(:, 2) / mass_water / c_water - 273)
    figure(2)
    plot(t, u(:, 1) / mass_pot / c_pot - 273)
end