clear all
close all
load data
lum = ['051';'007';'137';'033';'132'];
tempE = ['131';'147';'157';'010'];

% % %--->>> idea: uso solo il sensore a nord di temperatura (147) meno soggetto
% % %             a sbarellamenti di temperatura
% % tempE = ['147'];

% idea 2 aggiungo un sensore di temperatura interno come ingresso
% indipendente
tempE = [tempE; '141'];

pos_uT = name2pos(tempE,data)
pos_uL = name2pos(lum,data)
pos_y = setdiff([1:length(data.name)],[pos_uT pos_uL])

% tolgo il 102, qcosa nn va secondo me
pos_y = setdiff(pos_y,name2pos(['102'],data))

camp = 20; % 5 minuti
d = iddata(double(data.T(pos_y,1:camp:end)'),...
    double([ data.T(pos_uT,1:camp:end)' data.tR(pos_uL,1:camp:end)']),30*camp)

d1 = d(1:round(size(d)/2));
d2 = d(round(size(d)/2+1):end);

%% n4sid su metà dati e num di stati pari a uscite
tic
m = n4sid(d1,length(pos_y),'focus','simulation')
%m = n4sid(d1,10,'focus','simulation')
toc

%% trova forma equivalente con C = I
Ti = m.C;
T = inv(Ti);
C_ = m.C*T;
B_ = Ti*m.B;
figure(100)
imagesc(B_)
A_ = Ti*m.A*T;
figure(101)
imagesc(A_)

%% simula in catena aperta a partire dalle cond iniziali calcolate
y_sim = sim(m,d,'Init',inv(m.C)*d.OutputData(1,:)')
y_sim = sim(m,d)

%%
for canale = 1:length(pos_y)
    figure(1)
    plot([ y_sim.OutputData(:,canale) d.OutputData(:,canale)])
    axis tight
    legend('sim','vero')
    grid on
    figure(2)
    err = y_sim.OutputData(:,canale)-d.OutputData(:,canale);
    plot(err)
    axis tight
    ax = axis();
    axis([ax(1) ax(2) -4 4]);
    grid on
    figure(3)
    hold off
    for k = 1:length(data.name)
        if k == pos_y(canale)
            plot3(data.p(k,1),data.p(k,2),data.p(k,3),'rx','MarkerSize',10 )
        else
            plot3(data.p(k,1),data.p(k,2),data.p(k,3),'+','MarkerSize',3 )
        end
        text(data.p(k,1)+0.1,data.p(k,2)+0.1,data.p(k,3),data.name(k,:))
        hold on 
    end
    xlabel('x')
    ylabel('y')
    zlabel('z')
    title(sprintf('%s, fit %2.3f',data.name(pos_y(canale),:),100*(1 - norm(err)/norm(d.OutputData(:,canale)-mean(d.OutputData(:,canale))))))
    pause
end

