function EN40_HW1_2019

 close all
 
%% Problem 7
% Solution with a loop

for i=1:501
    t(i) = -2*pi + (i-1)*4*pi/500;
end

t = -2*pi:pi/500:2*pi;  % Direct dot notation


% [1 POINT]

%% Problem 8
% Create vectors of the Fresnel integrals (the Cornu spiral)
x = fresnelc(t);
y = fresnels(t);

% Plot the graph
plot(x,y,'LineWidth',2,'Color',[1 0 0]);
axes1 = gca;
ylabel('y');
xlabel('x');
title('A Cornu sprial');
box(axes1,'on');
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');

% [2 POINTS + 1 point credit for labeling the plot] 


%% Problem 9 drawing and modifying a black and white image with a matrix
% 9.1 Create an image
M = ones(101);
M(1:10,1:25) = 0; % This creates the rectangle
M(:,51) = 0;      % Vertical black line (the column 51 can be 50 or 52) 
M(51,:) = 0;      % Horizontal black line (the row 51 can be 50 or 52)
figure
imshow(M,'Initialmagnification',200);
title('Original Image')

% [2 POINTS]

%9.2 Mirror image about diagonal
Mnew = transpose(M);
figure
imshow(Mnew,'Initialmagnification',200);
title('Mirrored image');

% [1 POINT]

%9.3 Rotate by 90 degrees counterclockwise
[nr,nc] = size(M);
for i=1:nr
    for j = 1:nc
        Mnew(nr+1-i,j) = M(j,i);
    end
end
figure
imshow(Mnew,'Initialmagnification',200);
title('Rotated Image');
% Other approaches are fine too - there are many solutions
% [2 POINTS]

%% Problem 10 Potts model
% This is all copied from the HW assignment
M = 2*randi(2,101)-3;
M(:,1) = 0;
M(1,:) = 0;
M(101,:) = 0;
M(:,101) = 0;
figure
for i=1:100
   M = update_matrix(M);
   imshow((M+1)./2,'Initialmagnification',200)
   pause(0.1)  
end
% Some people seem to get weirdly symmetric
% figures for the evolution - I dont understand
% how that can happen but the codes look correct
% [3 POINTS - for correctly implementing function - see below]

    
 %% Problem 11 Simulation oscillations in the glycolysis reaction   

 % Values of the constants
 v0 = 1;
 k1 = 1;
 k2 = 1.01;
 
 % Initial conditions
 S0 = 1.;
 P0 = 0.8;
 initial_w = [S0,P0];
 
 % Start time and end time
 tstart = 0;
 tstop = 100.;
 
 [times,sols] = ode45(@(t,w) glycolysis(t,w,v0,k1,k2),[tstart,tstop],initial_w);

  figure
  plot1 = plot(times,sols);
  set(plot1(1),'DisplayName','Fructose 6 phospate');
  set(plot1(2),'DisplayName','ADP');
  axes1 = gca;
  ylabel('Concentrations (mols/cc)');
  xlabel('Time (min)');
  title('Higgins Selkov (k_2=1.01 min^{-1})');
  box(axes1,'on');
  set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
  legend(axes1,'show');
  
   k2 = 1.075;
   [times,sols] = ode45(@(t,w) glycolysis(t,w,v0,k1,k2),[tstart,tstop],initial_w);

  figure
  plot1 = plot(times,sols);
  set(plot1(1),'DisplayName','Fructose 6 phospate');
  set(plot1(2),'DisplayName','ADP');
  axes1 = gca;
  ylabel('Concentrations (mols/cc)');
  xlabel('Time (min)');
  title('Higgins Selkov (k_2=1.075 min^{-1})');
  box(axes1,'on');
  set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
  legend(axes1,'show');
  
   % [3 POINTS + 1 POINT for labeling plots]
  
  
end

%%   Functions

function B = update_matrix(A)
% Function for the Potts model simulation
    [nrows,ncols] = size(A);
    B = A;
    for ii=1:nrows*ncols
            i = randi(nrows-2)+1;
            j = randi(nrows-2)+1;
            e = 1-A(i,j)*sum(sum(A(i-1:i+1,j-1:j+1)));
            if (e>=0) B(i,j) = -A(i,j); end
    end


end

function dwdt = glycolysis(t,w,v0,k1,k2)
%  Function for the glycolysis differential equation
   S = w(1);   P = w(2);
   dSdt = v0 - k1*S*P^2;
   dPdt = k1*S*P^2 - k2*P;
   dwdt = [dSdt;dPdt];
end