function HW1_Matlab_2022

close all  % Closes any open figure windows


%% Problem 7
% Create a vector of 501 points starting at 0 and ending at 2 pi
disp('Solution to Problem 7: 501 equally spaced points from 0 to 2 pi')
t = linspace(0,2*pi,501)

% For those who prefer loops
% for i=1:501
%     t(i) = 2*pi*(i-1)/500;
% end

% [1 POINT]

%% Problem 8
% Calculate and plot the 'bunny' function
% Create x and y
t = 2*t;   %This creates a vector from 0 to 4 pi
for i=1:length(t)
        if (t(i)<2*pi)
            x(i) = 4*cos(t(i))*sin(t(i))^3;
            y(i) = 4*cos(t(i))^2*sin(t(i))^2;
        else
            x(i) = 0.75*sin(t(i));
            y(i) = 0.75*(cos(t(i))-1);
        end
end
% Plot the graph
plot1 = plot(x,y,'LineWidth',2);
axes1 = gca;
ylabel('y');
xlabel('x');
title('Plot for problem 8');
box(axes1,'on');
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
% [2 POINTS + 1 point credit for labeling the plot] 

%% Problem 9

% Write a function to calculate the 'hailstone sequence' 
y1 = collatz(7);   % Hailstone sequence starting at 7

disp('Hailstone sequence starting at 7');
disp(y1);
disp('Compare the the hand calculated sequence ')
disp([7,22,11,34,17,52,26,13,40,20,10,5,16,8,4,2,1])

% Now calculate the sequences starting at 54 and 97, and plot them.
y2 = collatz(54);
y3 = collatz(97);
figure
plot(y2,'LineWidth',2,'Displayname','I_1=54');
hold on
plot(y3,'LineWidth',2,'Displayname','I_1=97');

axes1 = gca;
ylabel('I_n');
xlabel('n');
title('Hailstone sequence');
box(axes1,'on');
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.2665983606557377 0.7 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);
% [5 POINTS + 1 point for labeling plot]


%% Problem 10

% Find and plot the solution for the first set of constants
nm = 2;   % No. female mosquitoes per human
a = 0.2;  % No. of infected bites by a single mosquito that produce infection per day
c = 0.5;  % Transmission efficiency from human to mosquito
r = 0.05; % Average recovery rate of humans
mu = 0.4; % Mortality rate of mosquitoes

tstart = 0;
tstop = 200;  % Run for 200 days.

initial_w = [0;0.01]; % Initial values of [h;m] at time t=0
[times,sols] = ode45(@(t,w) malaria(t,w,nm,a,c,r,mu),[tstart,tstop],initial_w); % Find the solution

% Plot
figure
  plot1 = plot(times,sols(:,1),'LineWidth',2,'Color',[1,0,0],'DisplayName','Humans');
  hold on
  plot1 = plot(times,sols(:,2),'LineWidth',2,'Color',[0,0,1],'DisplayName','Mosquitoes');
  axes1 = gca;
  ylabel('% infected');
  xlabel('Year');
  title('% Infected (2 female mosquitoes per human)');
  box(axes1,'on');
  set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
    legend1 = legend(axes1,'show');
    set(legend1,...
        'Position',[0.2665983606557377 0.7 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);


nm = 1; % Change the number of mosquitoes per human and get a new solution

[times,sols] = ode45(@(t,w) malaria(t,w,nm,a,c,r,mu),[tstart,tstop],initial_w);

% Now plot the results.
figure
  plot1 = plot(times,sols(:,1),'LineWidth',2,'Color',[1,0,0],'DisplayName','Humans');
  hold on
  plot1 = plot(times,sols(:,2),'LineWidth',2,'Color',[0,0,1],'DisplayName','Mosquitoes');
  axes1 = gca;
  ylabel('% infected');
  xlabel('Year');
  title('% Infected (1 female mosquitoe per human)');
  box(axes1,'on');
  set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
    legend1 = legend(axes1,'show');
    set(legend1,...
        'Position',[0.2665983606557377 0.7 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);
  
%    [5 POINTS + 1 point for labeling plot]


end


%% Differential equation function for problem 10
% The function returns a column vector containing the time
% derivatives of the unknown vector w at time t
% nm, a, c, r, mu are constants.
function dwdt = malaria(t,w,nm,a,c,r,mu)
   h = w(1); m = w(2); % h and m are the fractions of infected humans and mosquitoes, respectively
   
   dhdt = nm*a*m*(1-h) - r*h; % Time derivative of h
   dmdt = a*c*h*(1-m) - mu*m; % Time derivative of m
   
   dwdt = [dhdt;dmdt]; % Store the time derivatives in a column vector

end


%% Function for problem 9
% Calculates the hailstone sequence starting at 'start_integer'
function hailstone_sequence = collatz(start_integer)
    hailstone_sequence(1) = start_integer;  % Store the start
    next = start_integer;   % next will contain the next number in the sequence
    count = 1;              % Count the number of numbers that have been calclated
    while (next>1)          % Continue looping until next=1
        if (mod(next,2)>0)  % Check if next/2 has a remainder using the mod function.   
            next = 3*next+1; % Has remainder (i.e. odd number) multiply by 3 and add 1
        else
            next = next/2;  % No remainder; halve the number
        end
        count = count + 1;  % We calculated one more number
        hailstone_sequence(count) = next;  % Store it
    end
end


