clear all;
close all;

g = 9.81

month = {'July2006.txt' 'Aug2006.txt' 'Sept2006.txt' 'Oct2006.txt' 'Nov2006.txt' ...
         'Dec2006.txt' 'Jan2007.txt' 'Feb2007.txt' 'Mar2007-1.txt' 'Mar2007-2.txt' ...
         'Apr2007.txt' 'May2007.txt' 'June2007.txt' 'July2007.txt' 'Aug2007.txt'  ...
         'Sept2007.txt' 'Oct2007.txt' 'Nov2007.txt' 'Dec2007.txt'};

t = [];
Atop = []; 
Ltop = [];
Htop = [];
Dtop = [];

for j = 1:length(month)
  filenm = month{j}
  [A,L,H,D] = textread(filenm,'%*s A%f %*s L%f H%f D%f %*s');
  num_records(j) = length(A);
  
  A = (A-0.004)*10000/0.016;  %Convert to lbs (roughly)
  L = (L-0.004)*10000/0.016;  %Convert to lbs (roughly)
  H = (H-0.004)*10000/0.016;  %Convert to lbs (roughly)
  D = D*10000/0.016;  %Convert to lbs (roughly)

  Atop = [Atop; A-1000];  %Adjust for snubber weight
  Ltop = [Ltop; L-1000];
  Htop = [Htop; H-1000];
  Dtop = [Dtop; D];  %No need to adjust std deviation.

  t = [t length(t)+1:length(t)+length(A)];
end


[n,xout] = hist(Atop,40);
delta = diff(xout);
for i = 1:length(xout)
  idx{i} = find((Atop<=(xout(i)+delta(1)/2)) & (Atop>(xout(i)-delta(1)/2)));
  Ltop_bar(i) = mean(Ltop(idx{i}));
  Htop_bar(i) = mean(Htop(idx{i})); 
  Atop_bar(i) = mean(Atop(idx{i}));
  Dtop_bar(i) = mean(Dtop(idx{i}));

  Atop_max(i) = max(Atop(idx{i})); 
  Htop_max(i) = max(Htop(idx{i})); 
  Dtop_max(i) = max(Dtop(idx{i})); 
end


figure
plot(t,Atop,t,Ltop,t,Htop);
legend('Avg.','Max','Min');
title('MTM-4 Mooring Tensions.  10 minute sample period');
xlabel('Sample Period');
ylabel('Tension (lbs)');



figure
bar(xout,n*10/60,'r');
xlabel('Average Tensions (lbs)');
ylabel('Hours in this condition per year');
title('Histogram of Average Tensions Below Snubber');

figure
bar(xout,Htop_bar);
xlabel('Average Tensions (lbs)');
ylabel('Average Peak Tensions during this condition (lbs)');
title('Average Peak loads versus Average loading condition');

figure
bar(xout,Htop_max);
xlabel('Average Tensions (lbs)');
ylabel('Peak Tensions during this condition (lbs)');
title('Maximum Peak loads versus Average loading condition');

figure
bar(xout,Dtop_max);
xlabel('Average Tensions (lbs)');
ylabel('Max Std Deviation of tensions during this condition (lbs)');
title('Maximum Tension Std Deviation versus Average loading condition');

figure
bar(xout,Atop_bar+4*Dtop_max);
xlabel('Average Tensions (lbs)');
ylabel('(lbs)');
title('Peak Tensions based upon average and std deviation');




%Estimate Spectral Density of tensions that matches the shape of Pierson-Moskowitz
%spectrum.


freq_range = 0.01:.01:.3;  %hz
w = freq_range*2*pi; %rad/sec

Tension_range = min(Ltop):(max(Htop)-min(Ltop))/20:max(Htop);  %Create tension bins

a = 0.0081;
b = 0.74;
g = 9.81;

t = 0:.1:600;
Tension = [];
total = zeros(1,length(Tension_range));

for i = 1:length(xout)
U = Dtop_max(i)*20.0/1430;  %Conversion from std_dev load to wind speed
S = a*g^2*exp(-b*((g/U)./w).^4)./(w.^5);
T = zeros(1,length(t));
for j = 1:length(w)
T = T + (Htop_max(end)-Atop_bar(end))*sqrt(S(j))*cos(w(j).*t+2*pi*rand)/8.5;  %Conversion from sig wave height to tension std dev

end
T = T+Atop_bar(i);

mean_T(i) = mean(T);


[m,Tout] = histc(T,Tension_range);

total = total + n(i)*m;  %Multiply stats for this 10 minute simulted record by the number of records observed with this average tension.

end


figure
bar(Tension_range,total/10/3600);


paperwidth = 8.5;
paperheight = 11.0;

h_summary = figure('PaperPositionMode','manual','PaperSize',[paperwidth paperheight],'PaperPosition',[0 0 paperwidth paperheight],'Units','inches','Position',[2 -1 paperwidth paperheight],'Visible','on'); 

left = 1.5; bottom = 7.5; width = 6.0; height = 2.0;
h_TimeDomain  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);

left = 1.5; bottom = 4.5; width = 6.0; height = 2.0;
h_Histogram = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);

subplot(h_TimeDomain)
plot(t,Atop,t,Ltop,t,Htop);
legend('Avg.','Max','Min');
title('MTM-4 Mooring Tensions. 1 Year Record, Results per 10 minute sample period');
xlabel('Sample Period');
ylabel('Tension (lbs)');

subplot(h_Histogram)