function HW1_Matlab_2020

close all  % Closes any open figure windows


%% Problem 7
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

% Create x and y
a = [1,2,5];
for i=1:length(t)
    for j=1:3
        x(i,j) = cos(t(i))*(abs(cos(t(i))))^(a(j)-1);
        y(i,j) = sin(t(i))*(abs(sin(t(i))))^(a(j)-1);
    end
end

% Plot the graph
plot1 = plot(x,y,'LineWidth',2);
for j=1:length(a)
    set(plot1(j),'DisplayName',['a = ',num2str(a(j))]);
end
axes1 = gca;
ylabel('y');
xlabel('x');
title('Plot for problem 8');
box(axes1,'on');
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.2665983606557377 0.708333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);
% [2 POINTS + 1 point credit for labeling the plot] 

%% Problem 9

figure
M = zeros(101,101,3);
% Blue
M(1:101,1:101,3) = 1;
% Green
M(25:75,25:75,2) = 1;
M(25:75,25:75,1) = 0;
M(25:75,25:75,3) = 0;
% Red
M(27:73,27:73,1) = 1;
M(27:73,27:73,2) = 0;
M(27:73,27:73,3) = 0;
imshow(M,'InitialMagnification',400)
title('Image for problem 9');

% [2 POINTS]
blueno = imread('blueno4.jpg');
figure
imshow(blueno);
title('Blueno');
[m,n,k] = size(blueno);
pinkno = blueno;
bluenohsv = rgb2hsv(blueno);
for i=1:m
    for j=1:n
        if (bluenohsv(i,j,1) > 0.44 && bluenohsv(i,j,1)<0.55 )
            pinkno(i,j,1) = 242;
        end
    end
end

figure
imshow(pinkno)
title('Pinkno');

%% Optional problem
bluenpinkno = zeros(m,2*n,3);
bluenpinkno(1:m,1:n,1:3) = double(blueno)/256;
for i=1:m
    for j=1:n
        bluenpinkno(i,2*n+1-j,1:3) = double(blueno(i,j,1:3))/256;
        if (bluenohsv(i,j,1) > 0.44 && bluenohsv(i,j,1)<0.55 )
            bluenpinkno(i,2*n+1-j,1) = 242/256;
        end
    end
end
figure
imshow(bluenpinkno);
title('Blueno and pinkno');



% [3 POINTS]

%% Problem 10

  lambda = 1;
  initial_I = 0.05;
  time_int = [0,30];
  mu = 0.85;
  [time_vals,sol_vals] = icantbelieveitsnotode45(@(t,I) infectionrate(t,I,lambda,mu),time_int,initial_I);
  figure
  plot(time_vals,sol_vals,'DisplayName','\mu=0.85','Color',[0 0.447 0.741],...
     'LineWidth',2);
  hold on
  mu = 0.5;
  [time_vals,sol_vals] = icantbelieveitsnotode45(@(t,I) infectionrate(t,I,lambda,mu),time_int,initial_I);
  plot(time_vals,sol_vals,'DisplayName','\mu=0.5','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
  mu = 0.1;
  [time_vals,sol_vals] = icantbelieveitsnotode45(@(t,I) infectionrate(t,I,lambda,mu),time_int,initial_I);
  plot(time_vals,sol_vals,'DisplayName','\mu=0.1','Color',[0.929 0.694 0.125],...
     'LineWidth',2);

 axes1 = gca;
 xlabel({'Time (days)'},'FontName','Times New Roman');
 title('SI Infection Model');
 ylabel({'% Infected population'},'FontName','Times New Roman');
 set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');
 legend1 = legend(axes1,'show');
 set(legend1,...
    'Position',[0.6665983606557377 0.7 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [5 POINTS]

%% Problem 11

alpha1 = 0.005;
alpha2 = -40/(10^5);
alpha3 = 5.2/1000;
beta1 = 0.031;
beta2 = -1/(10^13);

tstart = 1850;
tstop = 2150;

initial_w = [1.15;0.21];
[times,sols] = ode45(@(t,w) gdp(t,w,alpha1,alpha2,alpha3,beta1,beta2),[tstart,tstop],initial_w);

figure
  yyaxis left
  plot1 = plot(times,sols(:,1),'LineWidth',2);

  axes1 = gca;
  ylabel('Population (Billions)');
  xlabel('Year');
  title('Population and GDP');
  box(axes1,'on');
  set(axes1,'FontSize',16,'XGrid','on','YGrid','on');

  yyaxis right
  plot2 = plot(times,sols(:,2),'LineWidth',2);
  ylabel('GDP (Trillions)');

% Optional problem
  [times,sols] = icantbelieveitsnotode45(@(t,w) gdp(t,w,alpha1,alpha2,alpha3,beta1,beta2),[tstart,tstop],initial_w);

  figure
  yyaxis left
  plot1 = plot(times,sols(:,1),'LineWidth',2);

  axes1 = gca;
  ylabel('Population (Billions)');
  xlabel('Year');
  title('Population and GDP');
  box(axes1,'on');
  set(axes1,'FontSize',16,'XGrid','on','YGrid','on');

  yyaxis right
  plot2 = plot(times,sols(:,2),'LineWidth',2);
  ylabel('GDP (Trillions)');
  
  
  
  % [4 POINTS]
  

end

function dIdt = infectionrate(t,I,lambda,mu)

   dIdt = (lambda-mu)*I - lambda*I^2;


end

function dwdt = gdp(t,w,alpha1,alpha2,alpha3,beta1,beta2)
   p = w(1); g = w(2);
   
   dpdt = alpha1*p  + alpha2*g*p+ alpha3*g;
   dgdt = beta1*g + beta2*p*g;
   
   dwdt = [dpdt;dgdt];

end

function [time_vals,sol_vals] = icantbelieveitsnotode45(f,time_int,initial_w)
% ODE solver function for problem 10
    N = 500;   % 500 points in solution
    nvars = length(initial_w);  % Number of unknowns
    time_vals = zeros(N,1);     % Storage for time values in solution
    sol_vals = zeros(N,nvars);  % Storage for solution
    time_vals(1) = time_int(1); % Start time
    sol_vals(1,1:nvars) = initial_w(1:nvars); % Initial solution

    dtime = (time_int(2)-time_int(1))/(N-1); % Time step size

    for i=2:N
        time_vals(i) = time_vals(i-1) + dtime; % Update time
        dwdt0 = f(time_vals(i-1),sol_vals(i-1,:)'); % Slope at time t
        dwdt = f(time_vals(i-1)+dtime/2,sol_vals(i-1,:)' + (dtime/2)*dwdt0); % Estimate of slope at time t+dt/2
        sol_vals(i,:) = sol_vals(i-1,:) + dtime*dwdt'; % Update the solution
    end

end



