close all

ModelResults.meanElecPow = nanmean(dl.Current.*dl.Voltage);
ModelResults.maxElecPow = max(dl.Current.*dl.Voltage);

ModelResults.PTO_Eff = nanmean((dl.Current.*dl.Voltage)./(dl.Force.*dl.ElonRate));

t = get(Buoy,'t');
%Get Buoy motion results
x = get(Buoy,'x');
xdot = get(Buoy,'xdot');
xddot = get(Buoy,'xddot');
 
%Get Plate motion results
y = get(Plate,'x');
ydot = get(Plate,'xdot');
yddot = get(Plate,'xddot');

%plotting results

Bforce = get(Buoy,'Forces');
Fforce = get(Plate,'Forces');

xoB = get(Buoy,'x0');
xoF = get(Plate,'x0');

Elongation = zeros(length(t),1);
ElongationRate = zeros(length(t),1);
SpringTension = zeros(length(t),1);
StictionTension = zeros(length(t),1);

s = zeros(length(t),3);
for i = 1:length(t)
    [Elongation(i),ElongationRate(i),s(i,:)] = TetherElongations(xoB,x(i,:)',xdot(i,:)',p_aB,xoF,y(i,:)',ydot(i,:)',p_aF);
    if(Elongation(i) > 0)  %If ram is fully retracted, tether is slack and no force is exerted.
        SpringTension(i) = ppval(Spring_Force_Coeffs,Elongation(i));
        %StictionTension(i) = ppval(Stiction_Force_Coeffs,ElongationRate(i));
    else
        Elongation(i) = 0;
        ElongationRate(i) = 0;
    end
end

ModelResults.MeanPTOElongation = mean(Elongation);
ModelResults.MaxPTOElongation = max(Elongation);
ModelResults.MinPTOElongation = min(Elongation);
ModelResults.StdPTOElongation = std(Elongation);

ModelResults.MeanPTOElongationRate = mean(ElongationRate);
ModelResults.MaxPTOElongationRate = max(ElongationRate);
ModelResults.MinPTOElongationRate = min(ElongationRate);
ModelResults.StdPTOElongationRate = std(ElongationRate);

t_eta  = [t t(end)+t1];
etadot = [0 diff(eta0)./diff(t_eta)];


figure('Position',[520 114 707 984])
ah1 = subplot(3,1,1)
plot(t_eta,etadot)
ah2 = subplot(3,1,2)
plot(t,dl.Current);
xlabel('t')
ylabel('Power (W)');
legend(['MeanPow = ' num2str(ModelResults.meanElecPow) 'W']);


ah3 = subplot(3,1,3)
plot(t,Elongation-mean(Elongation),'k-',t,ElongationRate,t,(BForces{4}(:,3)+BForces{1}(:,3)+BForces{2}(:,3)+PForces{2}(:,3)+PForces{1}(:,3))/10000,'r--')
legend('e','edot','F_{exct}')
xlabel('time (s)');

linkaxes([ah1 ah2 ah3],'x');


%Tabulate Power Balance