% Plot Orbit of Displacements
% To run this m-file, you need "wave_dat" from "grfft.f"

explanation= 'Plot Orbit of Displacements; You need "wave_dat" from "grfft.f" to run this.'

clear
close

fig=figure('Units','centimeters','Position',[0.0 0.0 21.0 21.0]...
,'PaperUnits','centimeters','PaperPosition',[0.0 0.0 21.0 21.0]);

load t071_d2;
wave_dat=t071_d2;
[nt,n]=size(wave_dat)
amax=270

delt=wave_dat(2,1)-wave_dat(1,1)
tf=0.0;
te=delt*nt;

for i=1:nt,
t(i)=wave_dat(i,1);
z_dat(i)=wave_dat(i,2);
x_dat(i)=wave_dat(i,3);
y_dat(i)=wave_dat(i,4);
end

yrange=[-amax amax];
ax=[tf te yrange];
subplot(3,1,1), plot(t,z_dat), axis(ax)
 title('u-d'), ylabel('displacement (cm)')
subplot(3,1,2), plot(t,x_dat), axis(ax)
 title('n-s'), ylabel('displacement (cm)')
subplot(3,1,3), plot(t,y_dat), axis(ax)
 title('e-w'), ylabel('displacement (cm)')
xlabel('time (sec)')

explanation='Click twice with the mouse on the screen to select a time interval'

t1=ginput(1);
t1(2)=t1(1);
subplot(3,1,1), hold on, plot(t1,yrange,'r-')
subplot(3,1,2), hold on, plot(t1,yrange,'r-')
subplot(3,1,3), hold on, plot(t1,yrange,'r-')

t2=ginput(1);
t2(2)=t2(1);
subplot(3,1,1), hold on, plot(t2,yrange,'r-')
subplot(3,1,2), hold on, plot(t2,yrange,'r-')
subplot(3,1,3), hold on, plot(t2,yrange,'r-')

print

explanation='Zoom in on the selected time interval'
trange=sort([t1(1),t2(1)])
k=find((trange(1)<=t) & (t<=trange(2)));
ud=z_dat(k);
ns=x_dat(k);
ew=y_dat(k);
t=t(k);
ax=[trange,yrange];
subplot(3,1,1), plot(t,ud), axis(ax)
 title('u-d'), ylabel('displacement (cm)')
subplot(3,1,2), plot(t,ns), axis(ax)
 title('n-s'), ylabel('displacement (cm)')
subplot(3,1,3), plot(t,ew), axis(ax)
 title('e-w'), ylabel('displacement (cm)')
xlabel('time (sec)')

%print
explanation='Enter Return'
pause

% Orbit of the displacement

amaxu=max(ud);
aminu=min(ud);
if amaxu < -aminu
  amaxu=-aminu;
end
amaxn=max(ns);
aminn=min(ns);
if amaxn < -aminn
  amaxn=-aminn;
end
amaxe=max(ew);
amine=min(ew);
if amaxe < -amine
  amaxe=-amine;
end

amax=amaxu;
if amaxn > amax
   amax=amaxn;
end
if amaxe > amax
   amax=amaxe;
end
amax=1.1*amax;

subplot(2,2,1), plot(ns,ud,'y-');
axis([-amax,amax,-amax,amax]), xlabel('S<-->N (cm)'), ylabel('D<-->U (cm)')
n=size(t)
text(ns(1),ud(1),'* start')

subplot(2,2,2), plot(ew,ud,'m-');
axis([-amax,amax,-amax,amax]), xlabel('W<-->E (cm)'), ylabel('D<-->U (cm)')
text(ew(1),ud(1),'* start')

subplot(2,2,3), plot(ew,ns,'c-');
axis([-amax,amax,-amax,amax]), xlabel('W<-->E (cm)'), ylabel('S<-->N (cm)')
text(ew(1),ns(1),'* start')

% 3-D trajectory of the displacement
subplot(2,2,4), plot3(ns,ew,ud,'r-');
axis([-amax,amax,-amax,amax,-amax,amax])
xlabel('S-N'), ylabel('W-E'), zlabel('D-U')
text(ns(1),ew(1),ud(1),'* start')
set(gca,'box','on')

print
explanation='Enter Return'
pause

% Large 3-D trajectory of the displacement 
subplot(1,1,1)
plot3(ns,ew,ud,'r');
axis([-amax,amax,-amax,amax,-amax,amax])
set(gca,'box','on')
xlabel('S-N'), ylabel('W-E'), zlabel('D-U')
title('3-D trajectory of displacement (unit: cm)')
text(ns(1),ew(1),ud(1),'* start')
hold on
plt=plot3(0,0,0,'.','erasemode','xor','markersize',24);
for k=1:length(ew)
  set(plt,'xdata',ns(k),'ydata',ew(k),'zdata',ud(k))
  drawnow
end
hold off
s='for k=1:length(ew)';
s=[s,'set(plt,''xdata'',ns(k),''ydata'',ew(k),''zdata'',ud(k)),'];
s=[s,'drawnow,end'];
echo off

%print





