punto7periodica

clear all;
clc;
format long;
windmill = xlsread('windmill.xlsx');
for i = 1:size(windmill,1)
   m(i) = i/12;
end


cont = size(windmill,1);
b = windmill;

for i = 0:59
   for j = 1:4
      if j == 1
         a(i+1,j) = 1;
      end
      if j == 2
         a(i+1,j) = cos(2*pi*i/12);
      end
      if j == 3
         a(i+1,j) = sin(2*pi*i/12);
      end
      if j == 4
         a(i+1,j) = cos(4*pi*i/12);
      end
   end
end

at = transpose(a);

newa = mtimes(at, a);
newb = mtimes(at, b);

n = size(newa, 1);

for j = 1 : n - 1
    if abs(newa(j,j)) < eps; error('Pivote cero encontrado'); end;
    for i = j + 1 : n
        mult = newa(i,j)/newa(j,j);
        for k = j + 1 : n
            newa(i,k) = newa(i,k) - mult * newa(j,k);
        end
        newb(i) = newb(i) - mult * newb(j);
    end
end

for i = n : -1 : 1
    for j = i + 1 : n
        newb(i) = newb(i) - newa(i,j)*x(j);
    end
    x(i) = newb(i) / newa (i, i);
end

x = transpose(x);
ax = mtimes(a, x);
gamma = minus(b, ax);

ec = 0;
for i = 1 : size(gamma)
   ec = ec + (gamma(i))^2;
end

rmec = sqrt(ec/(cont));


for i = 1 : cont
   d(i) = x(1) + x(2)*cos(2*pi*i/12) + x(3)*sin(2*pi*i/12) + x(4)*cos(4*pi*i/12);
end

plot(m, windmill, '*',m,d);