function [q,t_temp] = BECK(ihcp,r)
%BECK Funkce implementujici Beckovu sekvencni metodu pro odhad tepelneho
%toku
q = [];
options = optimset('Display','iter');
U0 = ihcp.eq.phi(ihcp.eq.domain.p(1,:),ihcp.eq.domain.p(2,:));
for M = 1:size(ihcp.t,2)-r
    fprintf('BSA: Krok %i z %i.\n',M,size(ihcp.t,2)-r);
    if M > 1
        t_temp = ihcp.t(1:M);
        ihcp.eq.beta{ihcp.ind_edge} = @(x,y,t,T) vect2fun(q,t_temp,t);
        if ihcp.eq.linear
            [U,~] = THETA(ihcp.eq,t_temp,[],0.5,false,5);
        else
            [U,~] = ODE(ihcp.eq,t_temp,[],false);
        end
    else
        U = U0;
    end
    q0 = 0;
    if M > 1
        q0 = q(M-1);
    end
    q = [q,fminsearch(@(qM)wrapper(qM,M,ihcp,r,U(:,end)),q0,options)];
end
end

function S = wrapper(qM,M,ihcp,r,U0)
    t_temp = ihcp.t(M:M+r-1);
    q_temp = qM*ones(1,r);
    ihcp.eq.beta{ihcp.ind_edge} = @(x,y,t,T)vect2fun(q_temp,t_temp,t);
    if ihcp.eq.linear
        [T,~] = THETA(ihcp.eq,t_temp,U0,0.5,false,0);
    else
        [T,~] = ODE(ihcp.eq,t_temp,U0,false);
    end
    [T,~] = INTERPOLATE(ihcp.eq,T,t_temp,ihcp.location,false,'');
    S = sum(sum((ihcp.T(:,M:M+r-1)-T).^2));
end

function value = vect2fun(vector,segments,t)
    segments = segments(:);
    segments = segments';
    smaller = find(segments<=t);
    index = smaller(end);
    equal = find(segments==t);
    if equal == size(segments,2)
        index = size(segments,2)-1;
    end
    value = vector(index);
end

