
results_filenm = create_filenm([dir 'FlappingFoil_' name],Hs,Ls,Lc,dt,tf,b)
load([results_filenm '.mat']);

t = get(Buoy,'t');

%Incident Wave profile
figure('Position', [100 100 1000 600])
set(gcf,'color','white');
plot(t,eta0(1:N));
legend('Incident wave profile (m)');

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

%plotting results

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

%Cylinder kinematics

figure('Position', [100 100 1000 600])
set(gcf,'color','white');
subplot(3,1,1),plot(t,x(:,1),t,x(:,3),t,x(:,5));
title('Buoy Motions','FontName','times','FontSize', 12);
legend('surge(m)','heave(m)','pitch(rad)');

subplot(3,1,2),plot(t,xdot(:,1),t,xdot(:,3),t,xdot(:,5));
legend('surge vel (m/s)','heave vel (m/s)','pitch vel (m/s)');

subplot(3,1,3),plot(t,xddot(:,1),t,xddot(:,3),t,xddot(:,5));
legend('surge accel (m/s^2)','heave accel (m/s^2)','pitch accel (m/s^2)');
xlabel('time (s)');


%foil kinematics
figure('Position', [100 100 1000 600])
set(gcf,'color','white');

subplot(3,1,1), plot(t,y(:,1),t,y(:,3),t,y(:,5));
title('Foil Motions','FontName','times','FontSize', 12);
legend('surge','heave(m)','pitch(rad)');

subplot(3,1,2), plot(t,ydot(:,1),t,ydot(:,3),t,ydot(:,5));
legend('surge(m/s)','heave vel (m/s)','pitch vel (m/s)');

subplot(3,1,3), plot(t,yddot(:,1),t,yddot(:,3),t,yddot(:,5));
legend('surge(m/s^2)','heave accel (m/s^2)','pitch accel (m/s^2)');
xlabel('time (s)');


V = zeros(1,N);
AOA = zeros(1,N);
for n = 1:N
[V(n) AOA(n)] = relative_vel(Foil,-n,Buoy,wave);
end

figure
plot(AOA*180/pi,V,'x');
xlabel('Angle of Attack (degrees)');
ylabel('Speed (m/s)');

nu = 1.1E-6;  %kinematic viscosity (m^2/s)plo
figure
plot(AOA*180/pi,V*Lc/nu,'x');
xlabel('Angle of Attack (degrees)');
ylabel('Re_c - Re Number based on chord length');

figure
hist(AOA*180/pi,50);
xlabel('Angle of Attach (degrees)');
ylabel('timestep count');


beta = atan2(-ydot(:,3),ydot(:,1));
theta_optimum = 10*pi/180;
theta_max = 20*pi/180;
th_target = zeros(1,length(beta));
for i = 1:length(beta)
   th_target(i) = theta_target(beta(i),theta_optimum,theta_max);
end
figure
plot(t,y(:,5)*180/pi,t,th_target*180/pi)
xlabel('time (s)');
ylabel('Foil Pitch (degrees)');
legend('Actual','Target');


p_aB = [0.0  0.0  -1.25];  %Tether attachment point on Buoy (local CoordSys).
p_aF = [0.0  0.0  0.0];  %Tether attachment point on Foil (local CoordSys).

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



Elongation = zeros(length(t),1);
ElongationRate = zeros(length(t),1);
SpringTension = zeros(length(t),1);
DampingTension = 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));
        DampingTension(i) = ppval(Damping_Force_Coeffs,ElongationRate(i));
    end
end



%Tether kinematics and Dynamics
figure('Position', [100 100 1000 600])
set(gcf,'color','white');

subplot(2,1,1), plot(t,Elongation)
ylabel('Tether Elongation (m)');

subplot(2,1,2), plot(t,SpringTension/4.45)
ylabel('SpringTension - lbs');
xlabel('time (s)');


figure('Position', [100 100 1000 600])
set(gcf,'color','white');

subplot(2,1,1), plot(t,ElongationRate)
ylabel('Tether Elongation Rate (m/s)');

subplot(2,1,2), plot(t,DampingTension/4.45)
ylabel('DampingTension - lbs');
xlabel('time (s)');





%Compute and display auxillary results
pto_power = ElongationRate.*DampingTension;
propulsion_force = Bforce{end}(:,1);
propulsion_power = propulsion_force.*xdot(:,1);
disp_power = -ydot(:,3).*Fforce{3}(:,3);
disp(['Mean Forward Velocity = ' num2str(mean(xdot(:,1))) ' m/s']);
disp(['Mean Propulsio Force = ' num2str(mean(propulsion_force)/4.45) ' lbs']);
disp(['Mean Propulsion Power = ' num2str(mean(propulsion_force)*mean(xdot(:,1))) ' Watts']);
disp(['Mean PTO Power = ' num2str(mean(pto_power)) ' Watts']);
disp(['Mean Power Dissipated = ' num2str(mean(disp_power)) ' Watts']);

%M(:,1) = t;
%M(:,2) = y(:,1);
%M(:,3) = ydot(:,1);
%M(:,4) = yddot(:,1);
%M(:,5) = y(:,3);
%M(:,6) = ydot(:,3);
%M(:,7) = yddot(:,3);
%M(:,8) = y(:,5);
%M(:,9) = ydot(:,5);
%M(:,10) = yddot(:,5);
%xlswrite([results_filenm '.xls'], M)
