clear all
close all
load data
lum = ['051';'007';'137';'033';'132'];
tempE = ['131';'147';'157';'010'];
pos_uT = name2pos(tempE,data)
pos_uL = name2pos(lum,data)
pos_y = setdiff([1:length(data.name)],[pos_uT pos_uL])

camp = 10; % 5 minuti
d = iddata(double(data.T(pos_y,1:camp:end)'),...
    double([ data.T(pos_uT,1:camp:end)' data.aR(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')
toc

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

%%
for canale = 1:56
    plot([ y_sim.OutputData(:,canale) d.OutputData(:,canale)])
    pause
end

