function [t,x,d,te,z,e] = NetworkSimRun(tF,Dat)
%Deterministic simulation of a single smart feeder
%[t,x,d,te,z,e] = NetworkSimRun(tF,Dat)
%execute run with default initial conditions
%inputs:
%tf: final time
%Dat: perturbed network paramaters; Dat = [X2,Ps,Tm]
%outputs:
%t: time trajectory
%x: continuous state trajectory
%te: event time trajectory
%z: dynamic parameter trajectory
%e: event trajectory






%get network description from external file
%   p0 = nominal static parameters
%   z0 = initial discrete state
%   St = state modular structure
%   Co = module interconnection
%   fnames = names of ode, event, Jacobian functions
%   xalg = list of continuous states
%   options = solver options
[p0,z0,St,Co,fnames,xalg,oevents,FSA,FSAsup] = getNetworkPars;

%get network initial state and adjusted parameters
options = odeset('Jacobian',@(t,x,d)jfun(t,x,d,z0,p0,St,Co,fnames),'events',@(t,x,d)efun(t,x,d,z0,p0,St,fnames,length(oevents)));
[x0,d0,p0,z0] = getNetworkX0(-4,p0,z0,St,Co,fnames,options);


%perturb the parameter vector according to Dat
p = perP(p0,Dat);


%initialize trajectory arrays and current values
t = 0; %continuous time
x = x0; %state
d = d0; %state derivative
te = 0; %event time
e = {''}; %event
z = [z0{:}]'; %dynamic parameters
qSup = ones(length(FSAsup),1); %supervisor state
qSys = ones(length(FSA),1); %system state

%initialize current values
tc = t(end);
xc = x(:,end); %state
dc = d(:,end); %state derivative
zc = z0; %dynamic parameters
qSupc = qSup(:,end); %supervisor state
qSysc = qSys(:,end); %system state

%compute the instantanous events and state reset resulting from static
%parameter change
options = odeset('Jacobian',@(t,x,d)jfun(t,x,d,zc,p,St,Co,fnames),'events',@(t,x,d)efun(t,x,d,zc,p,St,fnames,length(oevents)));
[xc,dc,enew] = instantChange(xc,dc,zc,p0,zc,p,St,Co,fnames,xalg,oevents,options);
t = [t tc];
x = [x xc];
d = [d dc];
flagi = ~isempty(enew);
while tc<tF
    %execute events
    while flagi==1||flagi==2
        %update FSAs
        if flagi==1
            if ~isempty(enew)
                qSupnew = updateFSA(FSAsup,enew,qSupc);
                qSysnew = updateFSA(FSA,enew,qSysc);
                qSupc = qSupnew(:,end);
                qSysc = qSysnew(:,end);
                
                tenew = tc*ones(1,length(enew));
                znew = z(:,end)*ones(1,length(enew));
                
                te = [te tenew];
                e = [e enew];
                z = [z znew];
                qSup = [qSup qSupnew];
                qSys = [qSys qSysnew];
            end
            %compute allowed controllable event (assumed that controllable
            %events are not observable)
            Gam = [];
            for i = 1:length(FSAsup)
                Gam = [Gam;FSAsup{i}(qSupc(i)).events(FSAsup{i}(qSupc(i)).con == 1)];
            end
        else
            Gam = enew;
        end
        %update the dynamic parameters

        [zcp,znew,qSysnew,enew] = IO(xc,p,zc,qSysc,Gam,FSA);

        if ~isempty(enew)
            qSupnew = updateFSA(FSAsup,enew,qSupc);
            
            qSupc = qSupnew(:,end);
            qSysc = qSysnew(:,end);
            
            tenew = tc*ones(1,length(enew));
            
            te = [te tenew];
            e = [e enew];
            z = [z znew];
            qSup = [qSup qSupnew];
            qSys = [qSys qSysnew];
            
            %check for more instant events
            options = odeset('Jacobian',@(t,x,d)jfun(t,x,d,zcp,p,St,Co,fnames),'events',@(t,x,d)efun(t,x,d,zcp,p,St,fnames,length(oevents)));
            [xc,dc,enew] = instantChange(xc,dc,zc,p,zcp,p,St,Co,fnames,xalg,oevents,options);
            zc = zcp;
            t = [t tc];
            x = [x xc];
            d = [d dc];
        end
        if flagi == 1
            flagi = ~isempty(enew);
        elseif flagi == 2
            flagi = 1;
        end
    end
    Gamsys = [];
    for i = 1:length(FSA)
        Gamsys = [Gamsys;FSA{i}(qSysc(i)).events];
    end
    if ~isempty(intersect(Gamsys,{'tf' 'cf'}))
        enew = intersect(Gamsys,{'tf' 'cf'});     
        flagi = 2;
    else
        sol = ode15i(@(t,x,d)odefun(t,x,d,zc,p,St,Co,fnames),[tc tF],xc,dc,options);
        [xnew,dnew] = deval(sol,sol.x);
        
        t = [t sol.x];
        x = [x xnew];
        d = [d dnew];
        
        tc = t(end);
        xc = x(:,end);
        dc = d(:,end);
        
        %check if system event occured
        enew = oevents(sol.ie);
        if ~isempty(intersect(Gamsys,{'tcstop' 'ttstop'}))
            flagi = 2;
        else
            flagi = ~isempty(enew);
        end
    end
end
te = [te tF];
e = [e {''}];
z = [z [zc{:}]'];
qSup = [qSup qSupc];
qSys = [qSys qSysc];


figure(1)

Vm = p{1}(3);
Vp = p{1}(4);
V2 = sqrt(x(3,:).^2 + x(4,:).^2);
figure(1)
subplot(2,3,1)
plot(t,V2)
line([t(1) t(end)],[Vm Vm],'color','r','linestyle','--')
line([t(1) t(end)],[Vp Vp],'color','r','linestyle','--')
title('Voltage Sensor')
ylabel('|V2|')
xlabel('time (sec)')
set(gca,'xlim',[0 tF])

subplot(2,3,2)
plot(t,x(14,:))
title('Tap Timer')
ylabel('ySS')
xlabel('time (sec)')
set(gca,'xlim',[0 tF],'ylim',[0 p{2}(1)])

subplot(2,3,3)
nT = z0{1}(1);
nTs = p{1}(5);
TU = 2*nTs + nT;
TL = -2*nTs + nT;
plot(te,z(1,:),'-o')
line([t(1) t(end)],[TL TL],'color','r','linestyle','--')
line([t(1) t(end)],[TU TU],'color','r','linestyle','--')
title('Transformer Tap')
ylabel('nT')
xlabel('time (sec)')
set(gca,'xlim',[0 tF])


Qs = x(4,:).*x(9,:) - x(3,:).*x(10,:);
Qc = x(4,:).*x(11,:) - x(3,:).*x(12,:);
Qerr = Qs + Qc;
subplot(2,3,4)

Qm = p{1}(6);
Qp = p{1}(7);
plot(t,Qerr)
line([t(1) t(end)],[Qm Qm],'color','r','linestyle','--')
line([t(1) t(end)],[Qp Qp],'color','r','linestyle','--')
title('Power Sensor')
ylabel('Q_{err}')
xlabel('time (sec)')
set(gca,'xlim',[0 tF])

subplot(2,3,5)
plot(t,x(16,:))
title('Capacitor Timer')
ylabel('ySS')
xlabel('time (sec)')
set(gca,'xlim',[0 tF],'ylim',[0 p{2}(1)])

subplot(2,3,6)
nCs = p{1}(8);
nC = z0{5}(1);
CU = 2*nCs + nC;
CL = -2*nCs + nC;
plot(te,z(7,:),'-o')
line([t(1) t(end)],[CL CL],'color','r','linestyle','--')
line([t(1) t(end)],[CU CU],'color','r','linestyle','--')
title('Capacitor Bank')
ylabel('nC')
xlabel('time (sec)')
set(gca,'xlim',[0 tF])






