%% SL 5
clear;
N = 6;
F = 1;
for n = 1:N
    F = F*n;
end

prod(1:N)

%% SL 9
clear;
N = 100; % calculate up to 100
k = 1;   % start with 1
s = 0;   % allocation for final sum
while k <= N
   s = s + k; % add current value to the sum
   k = k + 1; % move to the next value
end

sum(1:N)

%% SL 13
%% script fills a matrix
close all; clear;
M = 100;
N = 20;
A = nan(M, N);
for m = 1:M
    for n = 1:N
        A(m,n) = (m*n)/4 + m/(2*n);
    end
end
pcolor(A); % plot

%% SL 14
close all, clear;
M = 100;
N = 20;
m = (1:M).';
n = 1:N;
A = (m.*n)/4 + m./(2*n);
pcolor(A)

%% SL 15
%% script visualizes guitar string
close all; clear;
N   = 101;
om0 = 2*pi;
x   = linspace(-pi/2, pi/2, N);
t   = linspace(0, 4*pi, N);
I0  = cos(x);
for n = 1:N
    I = I0*exp(-1j*om0*t(n));
    figure(1);
    plot(x, real(I));
    axis([x(1) x(end) -1 1]);
    pause(0.1);
end

%% SL XX
r = 2*rand(8, 1) - 1;
Pos = []; % allocation
Neq = []; % allocation
for n = 1:length(r)
    if r(n) >= 0
        Pos(end+1) = r(n)
    else
        Neq(end+1) = r(n)
    end
end
% OR
Pos = r(r >= 0);
Neq = r(r < 0);

%% SL XX
clear, clc
s1 = 4;
s2 = 3;
type = 'hyp';
switch type
    case 'hyp'
        c = sqrt(s1^2 + s2^2);
    case 'leg'
        a = sqrt(max(s1, s2)^2 - min(s1, s2)^2);
    otherwise
        c = 0; a = 0;
end

%% SL
r = 2*rand(100, 1) - 1;
pos = 0;
neg = 0;
for n = 1:length(r)
    if r(n) >= 0
        pos = pos + 1;
    else
        neg = neg + 1;
    end
end

pos = sum(r >= 0);
neg = sum(r < 0);

clear, clc
z = rand(1) + 1j*rand(1) - 1/2 * (1 + 1j);
if real(z) >= 0 % I. or IV. quadrant
    if imag(z) >= 0
        disp('I. quadrant');
    else
        disp('IV. quadrant');
    end
else % II. or III. quadrant
    if imag(z) >= 0
        disp('II. quadrant');
    else
        disp('III. quadrant');
    end
end

clear, clc, close all
limit = 1000; % set the limit
F = [1 1]; % initial values
while F(end) + F(end-1) < limit
   F(end+1) = F(end) + F(end-1);
end
figure(1), plot(F, '-o')
% spiral length
sl = pi/2 * sum(F);
% golden ratios
phi = F(2:end) ./ F(1:end-1);
figure(2), plot(phi, '-or')

%% density  of prime  numbers
clear, clc, close all
N = 7;
nPrimes = nan(1, N); % alocation
orders = 10.^(1:N);
for n = 1:N
   nPrimes(n) = ...
     length(primes(orders(n)));
end
figure, plot(nPrimes)
figure, loglog(orders, nPrimes)

%% density  of prime  numbers
% improved performance
% primes is called only once
clear, clc, close all
N = 7;
nPrimes = nan(1, N); % alocation
orders = 10.^(1:N);
P = primes(orders(end));
for n = 1:N
   nPrimes(n) = sum(P < orders(m));
end
figure, plot(nPrimes)
figure, loglog(orders, nPrimes)

clear, clc
nStudents = 40;
points = randi(100, nStudents, 1);
A = 0; B = 0; C = 0;
for iS = 1:nStudents
   if points(iS) >= 80
      A = A+1;
   elseif points(iS) >= 35
      B = B+1;
   else
      C = C+1;
   end
end

A2 = sum(points >= 80);
B2 = sum(points >= 35 & points < 80);
C2 = sum(points < 35);

grade = 4;
switch grade
   case 1
      reward = 1000;
   case 2
      reward = 500;
   case 3
      reward = 200;
   case 4
      reward = 50;
   otherwise
      reward = 0;
end
disp(reward)

v = nan(4, 1);
while any(isnan(v))
   z = rand(1) + 1i * rand(1) ...
       - 1/2 * (1 + 1i);
   if real(z) >= 0
      if imag(z) >= 0
         quadrant = 1;
      else
         quadrant = 4;
      end
   else
      if imag(z) >= 0
         quadrant = 2;
      else
         quadrant = 3;
      end
   end
   v(quadrant) = z;
end

mask = [1 1; -1 1; -1 -1; 1 -1];
numbers = mask .* rand(4, 2);
v = 0.5 * complex(numbers(:, 1), numbers(:, 2));

%%
n1 = 0;
err = 1e-6;
estPi1 = 0;
while abs(estPi1 - pi/4) > err
   estPi1 = estPi1 + ...
      (-1)^n1 * 1/(2*n1+1);
   n1 = n1+1;
end % n1 = 250001

%%
n2 = 0;
estPi2 = 0;
while abs(estPi2 - pi/8) > err
   estPi2 = estPi2 + ...
      1/((4*n2+1) * (4*n2+3));
   n2 = n2+1;
end % n2 = 62500

%%
n3 = 0;
estPi3 = 0;
argumentValues = [1/8, 1/57, 1/239];
coefs = [6, 2, 1];
while abs(estPi3 - pi/4) > err
   estPi3 = 0;
   for iValue = 1:3
      x = argumentValues(iValue);
      n = 0:n3;
      y = sum((-1).^n .* x.^(2*n+1) ./ ...
         (2*n+1));
      estPi3 = estPi3 + coefs(iValue) * y;
   end
   n3 = n3 + 1;
end % n3 = 3

%%
grade = 4;
reward = [1000 500 200 50 0];
reward(grade == 1:5)

%%
phi = angle(z);
quadrant = ceil(2*phi/pi + 4*(phi < 0))