clear
close all

load 'procDataNS'
load 'procDataWS'

%To plot PSDs in dB set plotdB to 1, linear scaling plotdB = 0;
plotdB = 0;
%Set the maximum frequency for plotting
fMax = 0.7;
%Set the maximum amplitude for plotting
%if autoScaleY = 1 these will be ignored
autoScaleY = 1;
tensAMax = 2000.0;
buoyAccAMax = 2.0;
buoyZAmax = 5.0;
ESPAMax = 2.0;


%calculate tension statistics
tensNSMean = mean(tensionNS);
tensNSVar = var(tensionNS);
tensNSStd = std(tensionNS);
fprintf('Tension (lbf) without skirt: mean = %6.1f  std = %6.1f  var = %6.1f\n',tensNSMean, tensNSStd, tensNSVar);

tensWSMean = mean(tensionWS);
tensWSVar = var(tensionWS);
tensWSStd = std(tensionWS);
fprintf('Tension (lbf) with skirt:    mean = %6.1f  std = %6.1f  var = %6.1f\n\n',tensWSMean, tensWSStd, tensWSVar)

%plot tension
maxTensIndex = 17000;
time = 0.1 * 1:maxTensIndex;
plot(time,tensionNS(1:maxTensIndex),time,tensionWS(1:maxTensIndex))
title('Tension (lbf) vs Time (sec)');
legend('No Skirt', 'With Skirt')
grid

%plot tension PSD
figure
if plotdB == 1
    semilogy(tensionNSPSD.Frequencies,tensionNSPSD.Data, tensionWSPSD.Frequencies,tensionWSPSD.Data)
    ylabel('Power/frequency (dB lbf^2/Hz)')
else
    plot(tensionNSPSD.Frequencies,tensionNSPSD.Data, tensionWSPSD.Frequencies,tensionWSPSD.Data)
    ylabel('Power/frequency (lbf^2/Hz)')
end
title('Tension PSD')
legend('No Skirt', 'With Skirt')
xlabel('Frequency (Hz)')
limits = axis;
if autoScaleY == 1
    axis([limits(1) fMax limits(3) limits(4)]);
else
    axis([limits(1) fMax limits(3) tensAMax]);
end
grid

%calculate Buoy Acceleration Statistics
buoyAccZNSMean = mean(buoyAccZNS);
buoyAccZNSVar = var(buoyAccZNS);
buoyAccZNSStd = std(buoyAccZNS);
fprintf('Buoy Acceleration Z (g) without skirt: mean = %6.3f  std = %6.3f  var = %6.3f\n',buoyAccZNSMean, buoyAccZNSStd, buoyAccZNSVar);

buoyAccZWSMean = mean(buoyAccZWS);
buoyAccZWSVar = var(buoyAccZWS);
buoyAccZWSStd = std(buoyAccZWS);
fprintf('Buoy Acceleration Z (g) with skirt:    mean = %6.3f  std = %6.3f  var = %6.3f\n\n',buoyAccZWSMean, buoyAccZWSStd, buoyAccZWSVar);

%plot Buoy acceleration 
figure
maxBuoyAccZIndex = 17400;
time = 0.1 * 1:maxBuoyAccZIndex;
plot(time,buoyAccZNS(1:maxBuoyAccZIndex),time,buoyAccZWS(1:maxBuoyAccZIndex))
title('Buoy Acceleration Z (g) vs Time');
legend('No Skirt', 'With Skirt')
grid

%plot Buoy acceleration PSD
figure
if plotdB == 1
    semilogy(buoyAccZNSPSD.Frequencies,buoyAccZNSPSD.Data, buoyAccZWSPSD.Frequencies,buoyAccZWSPSD.Data)
    ylabel('Power/frequency (dB g^2/Hz)')
else
    plot(buoyAccZNSPSD.Frequencies,buoyAccZNSPSD.Data, buoyAccZWSPSD.Frequencies,buoyAccZWSPSD.Data)
    ylabel('Power/frequency (g^2/Hz)')
end
title('Buoy Acceleration Z PSD')
legend('No Skirt', 'With Skirt')
xlabel('Frequency (Hz)')
limits = axis;
if autoScaleY == 1
    axis([limits(1) fMax limits(3) limits(4)]);
else
    axis([limits(1) fMax limits(3) buoyAccAMax]);
end
grid

%Calculate buoy displacement statistics
buoyZNSMean = mean(buoyZNS);
buoyZNSVar = var(buoyZNS);
buoyZNSStd = std(buoyZNS);
fprintf('Buoy Z (m) without skirt: mean = %6.3f  std = %6.3f  var = %6.3f\n',buoyZNSMean, buoyZNSStd, buoyZNSVar);

buoyZWSMean = mean(buoyZWS);
buoyZWSVar = var(buoyZWS);
buoyZWSStd = std(buoyZWS);
fprintf('Buoy Z (m) with skirt:    mean = %6.3f  std = %6.3f  var = %6.3f\n\n',buoyZWSMean, buoyZWSStd, buoyZWSVar);

%plot buoy displacements
figure
maxBuoyZIndex = 17400;
time = 0.1 * 1:maxBuoyZIndex;
plot(time,buoyZNS(1:maxBuoyZIndex),time,buoyZWS(1:maxBuoyZIndex))
title('Buoy Z (m) vs Time');
legend('No Skirt', 'With Skirt')
grid

%plot buoy displacement PSD
figure
if plotdB == 1
    semilogy(buoyZNSPSD.Frequencies,buoyZNSPSD.Data, buoyZWSPSD.Frequencies,buoyZWSPSD.Data)
    ylabel('Power/frequency (dB m^2/Hz)')
else
    plot(buoyZNSPSD.Frequencies,buoyZNSPSD.Data, buoyZWSPSD.Frequencies,buoyZWSPSD.Data)
    ylabel('Power/frequency (m^2/Hz)')
end
title('Buoy Z PSD')
legend('No Skirt', 'With Skirt')
xlabel('Frequency (Hz)')
limits = axis;
if autoScaleY == 1
    axis([limits(1) fMax limits(3) limits(4)]);
else
    axis([limits(1) fMax limits(3) buoyAccAMax]);
end
axis([limits(1) fMax limits(3) buoyZAMax]);
grid

%calculate ESP Acceleration Statistics
ESPAccZNSMean = mean(ESPAccZNS);
ESPAccZNSVar = var(ESPAccZNS);
ESPAccZNSStd = std(ESPAccZNS);
fprintf('ESP Acceleration Z (g) without skirt: mean = %6.3f  std = %6.3f  var = %6.3f\n',ESPAccZNSMean, ESPAccZNSStd, ESPAccZNSVar);

ESPAccZWSMean = mean(ESPAccZWS);
ESPAccZWSVar = var(ESPAccZWS);
ESPAccZWSStd = std(ESPAccZWS);
fprintf('ESP Acceleration Z (g) with skirt:    mean = %6.3f  std = %6.3f  var = %6.3f\n\n',ESPAccZWSMean, ESPAccZWSStd, ESPAccZWSVar);

%plot ESP acceleration 
figure
maxESPAccZIndex = 17400;
time = 0.1 * 1:maxESPAccZIndex;
plot(time,ESPAccZNS(1:maxESPAccZIndex),time,ESPAccZWS(1:maxESPAccZIndex))
title('ESP Acceleration Z (g) vs Time');
legend('No Skirt', 'With Skirt')
grid

%plot ESP acceleration PSD
figure
if plotdB == 1
    semilogy(ESPAccZNSPSD.Frequencies,ESPAccZNSPSD.Data, ESPAccZWSPSD.Frequencies,ESPAccZWSPSD.Data)
    ylabel('Power/frequency (dB g^2/Hz)')
else
    plot(ESPAccZNSPSD.Frequencies,ESPAccZNSPSD.Data, ESPAccZWSPSD.Frequencies,ESPAccZWSPSD.Data)
    ylabel('Power/frequency (g^2/Hz)')
end    
title('ESP Acceleration Z PSD')
legend('No Skirt', 'With Skirt')
xlabel('Frequency (Hz)')
limits = axis;
axis([limits(1) fMax limits(3) ESPAMax]);
grid