%
% reading data and make figure
%

clear
fclose all ;

len=0.25; % length of debris objects 
wid=0.10; % width of debris objects

% for reading all data 
load('../raw_data/data.mat') ;

%  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  = 
% Data Contents in data.mat

% Three variables with structure arrays are included
%     : bathyData is the bathymetry and topography data etc used in my
%       simulation
%     : FDcase is the data of Floating Debris case
%     : LDcase is the data of Land Debris case

% 1. bathyData -- nx,ny: mesh number in each direction
%              -- X(nx), Y(ny): Coordinate information
%              -- Manning(nx,ny): Manning coefficient used in my simulation 
%              -- landFlag(nx,ny): land=1, water=0
%              -- depth(nx,ny): bathymetry and topography data used in my simulation 

% 2. FDcase
%              -- caseNo: case number (1~10) (nCase)
%              -- TimeAtWaveGauge: Time[s] of wave gauge data (ntWave)
%              -- Time: Time[s] of debris data (ntDebris)
%              -- waveDataAtGauge1(nCase,ntWave): time series of measured wave height at WG1
%              -- waveDataAtGauge2(nCase,ntWave): time series of measured wave height at WG2
%              -- CenterX(nCase,ntDebris), CenterY(nCase,ntDebris): time series of center position of debris object
%              -- CornerX(nCase,ntDebris), CornerY(nCase,ntDebris): time series of corner position of debris object
%              -- RotationAngle(nCase,ntDebris): time series of rotataion anlge of debris 

% 3. LDcase 
%              -- same as FDcase
%  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  = 


%% figure
% exp
X=bathyData.X ;
Y=bathyData.Y ;
topo=bathyData.depth ;
[nx,ny]=size(topo);

figure(1)
clf
% bath/topo data 
contour(X,Y,topo',[0 0]) ;
shading flat ;
axis equal tight; hold on ;
set(gca,'Ydir','reverse') ;
xlabel('X [m]') ; 
ylabel('Y [m]') ; 
title('Initial Debris Locations')
hold on

% plot initial debris position for FD case
num=FDcase.caseNo ;
for i=1:length(num)
    x=FDcase.CenterX ;
    y=FDcase.CenterY ;
    t=FDcase.RotationAnlge ;
    % this function(func_vertex) is calculated the position of debris
    % vertex using the center of gravity, rotation angle, width and length. 
    % I stored this function in the same folder as the data. 
    [xd,yd]=func_vertex(x(i,1),y(i,1),-t(i,1),len,wid) ;
    plot(xd,yd) ;
end

% plot initial debris position for LD case
num=LDcase.caseNo ;
for i=1:length(num)
    x=LDcase.CenterX ;
    y=LDcase.CenterY ;
    theta=LDcase.RotationAnlge ;
    [xd,yd]=func_vertex(x(i,1),y(i,1),-theta(i,1),len,wid) ;
    plot(xd,yd) ;
end

axis([0 4 -Inf Inf])
print -djpeg100 initial_locations.jpg

%%
% make figure for debris trajectory and time series of wave height in each
% experiment case 
% the figures "FD_case01-06.png", "FD_case07-10.png", "LD_case01-06.png", "LD_case07-10.png"
% are based on the following program. 

% for FD case

time=FDcase.Time ;
x=FDcase.CenterX ;
y=FDcase.CenterY ;
theta=FDcase.RotationAnlge ;
%plot as time series of x-postion, y-postion, and rotation angle

dt_interp=0.01;
t_interp=[0:dt_interp:330];
x_interp=zeros(length(t_interp),length(num));
y_interp=x_interp;
theta_interp=x_interp;

figure(2)
clf
subplot(3,1,1)
hold on
for i=1:length(num)
    t_c=time(i,:) ;
    x_c=x(i,:);
    x_interp(:,i)=interp1(t_c,x_c,t_interp);
    plot(t_c,x_c)
    legend_str{i}=['Trial ' num2str(i)];
end
x_mean=mean(x_interp,2);
x_std=std(x_interp,1,2);
plot(t_interp,x_mean,'r','LineWidth',4)
plot(t_interp,x_mean+x_std,'r--','LineWidth',2)
plot(t_interp,x_mean-x_std,'r--','LineWidth',2)
xlabel('Time (sec)')
ylabel('X-Postion (m)')
title('X-Position Time Series of FD Debris')
legend_str{i+1}='Mean Position';
legend_str{i+2}='Mean +/- StDev';
legend(legend_str)
data=[t_interp' x_mean x_std x_interp];
save FD_Xposition.txt data -ascii

subplot(3,1,2)
hold on
for i=1:length(num)
    t_c=time(i,:) ;
    y_c=y(i,:);
    y_interp(:,i)=interp1(t_c,y_c,t_interp);
    plot(t_c,y_c)
end
y_mean=mean(y_interp,2);
y_std=std(y_interp,1,2);
plot(t_interp,y_mean,'r','LineWidth',4)
plot(t_interp,y_mean+y_std,'r--','LineWidth',2)
plot(t_interp,y_mean-y_std,'r--','LineWidth',2)
xlabel('Time (sec)')
ylabel('Y-Postion (m)')
title('Y-Position Time Series of FD Debris')
legend(legend_str)
data=[t_interp' y_mean y_std y_interp];
save FD_Yposition.txt data -ascii

subplot(3,1,3)
hold on
for i=1:length(num)
    t_c=time(i,:) ;
    theta_c=theta(i,:);
    theta_interp(:,i)=interp1(t_c,theta_c,t_interp);
    plot(t_c,theta_c)
end
theta_mean=mean(theta_interp,2);
theta_std=std(theta_interp,1,2);
plot(t_interp,theta_mean,'r','LineWidth',4)
plot(t_interp,theta_mean+theta_std,'r--','LineWidth',2)
plot(t_interp,theta_mean-theta_std,'r--','LineWidth',2)
xlabel('Time (sec)')
ylabel('Rotation Angle (deg)')
title('Debris Rotation Angle Time Series of FD Debris')
legend(legend_str)
data=[t_interp' theta_mean theta_std theta_interp];
save FD_rotation.txt data -ascii

print -djpeg100 FD_postion_time_series.jpg

%plot 2D trajectories
figure(3)
clf
% bath/topo data 
contour(X,Y,topo',[0 0]) ;
shading flat ;
axis equal tight; hold on ;
set(gca,'Ydir','reverse') ;
xlabel('X [m]') ; 
ylabel('Y [m]') ; 
title('FD Case Debris Locations')
hold on

for i=1:length(num)
    plot(x(i,:),y(i,:),'Linewidth',0.5); % plot debris trajectory
end
plot(x_mean,y_mean,'r','Linewidth',4)
dnt_plot=round(5/dt_interp);
for n=1:dnt_plot:find(t_interp>80,1)
    [xd,yd]=func_vertex(x_mean(n),y_mean(n),-theta_mean(n),len,wid) ;
    plot(xd,yd,'r','Linewidth',0.5) ;
end

axis([0 4 -Inf Inf])
print -djpeg100 FD_trajectories.jpg


% for LD case

time=LDcase.Time ;
x=LDcase.CenterX ;
y=LDcase.CenterY ;
theta=LDcase.RotationAnlge ;
%plot as time series of x-postion, y-postion, and rotation angle

dt_interp=0.01;
t_interp=[0:dt_interp:320];
x_interp=zeros(length(t_interp),length(num));
y_interp=x_interp;
theta_interp=x_interp;

figure(4)
clf
subplot(3,1,1)
hold on
for i=1:length(num)
    t_c=time(i,:) ;
    x_c=x(i,:);
    x_interp(:,i)=interp1(t_c,x_c,t_interp);
    plot(t_c,x_c)
end
x_mean=mean(x_interp,2);
x_std=std(x_interp,1,2);
plot(t_interp,x_mean,'r','LineWidth',4)
plot(t_interp,x_mean+x_std,'r--','LineWidth',2)
plot(t_interp,x_mean-x_std,'r--','LineWidth',2)
xlabel('Time (sec)')
ylabel('X-Postion (m)')
title('X-Position Time Series of LD Debris')
legend(legend_str)
data=[t_interp' x_mean x_std x_interp];
save LD_Xposition.txt data -ascii

subplot(3,1,2)
hold on
for i=1:length(num)
    t_c=time(i,:) ;
    y_c=y(i,:);
    y_interp(:,i)=interp1(t_c,y_c,t_interp);
    plot(t_c,y_c)
end
y_mean=mean(y_interp,2);
y_std=std(y_interp,1,2);
plot(t_interp,y_mean,'r','LineWidth',4)
plot(t_interp,y_mean+y_std,'r--','LineWidth',2)
plot(t_interp,y_mean-y_std,'r--','LineWidth',2)
xlabel('Time (sec)')
ylabel('Y-Postion (m)')
title('Y-Position Time Series of LD Debris')
legend(legend_str)
data=[t_interp' y_mean y_std y_interp];
save LD_Yposition.txt data -ascii

subplot(3,1,3)
hold on
for i=1:length(num)
    t_c=time(i,:) ;
    theta_c=theta(i,:);
    theta_interp(:,i)=interp1(t_c,theta_c,t_interp);
    plot(t_c,theta_c)
end
theta_mean=mean(theta_interp,2);
theta_std=std(theta_interp,1,2);
plot(t_interp,theta_mean,'r','LineWidth',4)
plot(t_interp,theta_mean+theta_std,'r--','LineWidth',2)
plot(t_interp,theta_mean-theta_std,'r--','LineWidth',2)
xlabel('Time (sec)')
ylabel('Rotation Angle (deg)')
title('Debris Rotation Angle Time Series of LD Debris')
legend(legend_str)
data=[t_interp' theta_mean theta_std theta_interp];
save LD_rotation.txt data -ascii

print -djpeg100 LD_postion_time_series.jpg

%plot 2D trajectories
figure(5)
clf
% bath/topo data 
contour(X,Y,topo',[0 0]) ;
shading flat ;
axis equal tight; hold on ;
set(gca,'Ydir','reverse') ;
xlabel('X [m]') ; 
ylabel('Y [m]') ; 
title('LD Case Debris Locations')
hold on

for i=1:length(num)
    plot(x(i,:),y(i,:),'Linewidth',0.5); % plot debris trajectory
end
plot(x_mean,y_mean,'r','Linewidth',4)
dnt_plot=round(10/dt_interp);
for n=1:dnt_plot:find(t_interp>280,1)
    [xd,yd]=func_vertex(x_mean(n),y_mean(n),-theta_mean(n),len,wid) ;
    plot(xd,yd,'r','Linewidth',0.5) ;
end

axis([0 4 -Inf Inf])
print -djpeg100 LD_trajectories.jpg







