function rotation
%
%  Rotation;  Demonstrate the kinematic effects of a rotating 
%  reference frame by computing the motion of a particle in 3-d
%  as seen from an interial and rotating reference frames.
%  The projectile is launched upwards and with a horizontal
%  veclocity; the rotating reference frame (red) shows the 
%  position as seen from a rotating reference frame that
%  rotates at a steady rate. In Fig 2 the trajectory as seen from
%  the rotating reference frame is computed by including 
%  the cetrifugal and Coriolis accelerations.  The solution
%  computed in this rotating frame should be exactly the same  
%  as the position seen from the rotating frame if the centrifugal
%  and Coriolis accelerations are as required to account for
%  the kinematic effects of the rotating frame. To see the 
%  Coriolis accleration only (more like the motion of fluids 
%  on the rotating earth) try the script 'Coriolis.m'. 
%  The initial velocity, rotation rate etc. are here hardwired
%  to reasonable values; to make them variables, remove 
%  the percent sign in front of the lines with 'input'.
%
%  by Jim Price.  January, 2001.
%
%  This may be considered public domain for educational
%  purposes. 
%  This code was written on the assumption that the monitor would
%  have a resolution of 1024x768.  If the monitor has less, then
%  all best are off. 

clear
clear memory
close all


path(path, 'c:/jpsource/matextras')

% set default graphics things
set(0,'DefaultTextFontSize',12)
set(0,'DefaultAxesFontSize',12)
set(0,'DefaultAxesLineWidth',1.0)
set(0,'DefaultLineLineWidth',1.0)


str2(1) =  {'Rotation:' };
str2(2) =  {'  '};
str2(3) =  {'Examine the kinematics of a rotating coordinate system'};
str2(4) =  {'by showing the path of a free projectile as seen in an '};
str2(5) =  {'inertial reference frame (blue, Fig 1) and as seen from a'};
str2(6) =  {'frame rotating at a rate \Omega (red, Fig 1 and Fig 2).'}; 
str2(7) =  {'The projectile is launched with speed Vo in the '}; 
str2(8) =  {'Y direction and with an elevation of 60 deg. It is'}; 
str2(9) =  {'then accelerated downward at a rate g = -1. '}; 
str2(10) = {'The green trajectory of Fig 2 is the trajectory '}; 
str2(11) = {'computed in the rotating frame by including the Coriolis'};
str2(12) = {'and centrifugal accelerations.  This solution is'};
str2(13) = {'identical to the observations of the inertial trajectory' };
str2(14) = {'made from a rotating frame, thus showing that we have'};
str2(15) = {'been able to account for the kinematic effects of the'};
str2(16) = {'rotating frame. The + symbols are the trajectory'};
str2(17) = {'in 3-d; the . symbols are the projection onto z = 0 '};
str2(18) = {''};
str2(19) = {' Hit any key to continue and to step ahead.'};
str2(20) = {''};
str2(21) = {'Jim Price, January, 2001. '};
hf3 = figure(10);
clf
set(hf3,'Position',[50 50 620 620])
set(gca,'Visible','off')
text(-0.1, 0.50, str2,'FontSize',12,'HorizontalAlignment','left')
pause


% t = zeros(2000,1); xs = zeros(2000,3); xcs = xs; xrs = xs;
hpn = 0; k = 0; 

g = -1;
dt = 0.001;

% the initial position and velocity in the inertial frame

rotrate = 0.2;
initpos = 0.3;
initang = 60.0;
speed0 = 1.0;

% to have the particle at rest in the inertial frame try:
%  speed0 = 0.001; g = 0.; initpos = 1.0;           

% initpos = input(...
%   'Enter the initial position along y axis, 0 - 1  ')
% initang = input('Enter the elevation angle, 0 - 90 (deg)  ')
% speed0 = input('Enter the initial speed  ')
% rotrate = input('Enter rotation rate, cycles/unit time =  ')

%  compute the vectors of the IC

u0 = speed0*[0 cos(initang*pi/180), sin(initang*pi/180)];
x0 = [0 initpos 0.0001];
Omega = rotrate*2*pi*[0 0 1.0];   %  set the rotation rate 


x = x0;
u = u0;

%  initial position in the rotating frame
xr = x0;

%  initial position and velocity in the Coriolis frame
xc = x0;
uc = u0 - cross(Omega, x0);

j = 1;
t(1) = 0.;
ang = zeros(1,100);

figure(1)
clf reset

set(gcf,'position', [25 100 450 450])
set(gca,'xColor', [0 0 1], 'ycolor', [0 0 1], 'zcolor', [0 0 1])
xlabel('X'); ylabel('Y'), zlabel('Z')
title('an inertial frame (blue) and a rotating frame (red)')
axis([-1 1 -0.5 1.5 0 1])
view(30., 60.)
grid

corstuff = 'dV/dt = -g'; 
text(-0.5, 0., 1.8, corstuff, 'Color', 'b')
corstuff = 'V(0) = V_0'
text(-0.5, 0., 1.5, corstuff, 'Color', 'b')

hold on

harrowv0 = arrow3(u0, x0, 0.03, 0.08, 1, 20)
u0pos = u0 + x0;
harrowv0t = text(u0pos(1), u0pos(2), u0pos(3), 'V_0', 'color', 'b')

harrowo = arrow3([0 0 1], [0 0 0], 0.03, 0.08, 1, 20)
set(harrowo, 'color', 'r')
harrowot = text(0., 0., 1.1, '\Omega', 'color', 'r')

drawnow

figure(2)
clf reset
hold on

set(gcf, 'position', [500 100 450 450])
xlabel('Xr'); ylabel('Yr'); zlabel('Z')

tstr = str2mat('as seen from the rotating frame (red)', ...
   'and as computed in the rotating frame (green)');
htt = title(tstr, 'VerticalAlignment', 'top'); 

corstuff = 'dV_r/dt = -(g + 2 \Omega x V_r + \Omega x \Omega x r)'; 
text(-0.5, 0., 1.8, corstuff, 'Color', 'g')
corstuff = 'V_r(0) = V_0 - \Omega x r'
text(-0.5, 0., 1.5, corstuff, 'Color', 'g')

harrowv0c = arrow3(uc, x0, 0.03, 0.08, 1, 20);
set(harrowv0c, 'Color', 'r')
u0pos = uc + x0;
harrowv0ct = text(u0pos(1), u0pos(2), u0pos(3),...
   'V_{r0}', 'color', 'r')
axis([-1 1 -0.5 1.5 0 1])
grid
view(30., 60.)
drawnow

jp = 0;
j = 0;

while x(3) >= 0
   
j = j + 1;
t(j) = (j-1)*dt;
      
%  step ahead the inertial frame variables   
x = x + dt*u;
a = [0 0 g];
u = u + dt*a;
   
%  compute the position as seen from the rotating frame
angle = t(j)*Omega(3);
xr(1) = x(1)*cos(angle) + x(2)*sin(angle);
xr(2) = x(2)*cos(angle) - x(1)*sin(angle);
xr(3) = x(3);

%  time step the position and speed in the Coriolis frame
xc = xc + dt*uc;
Omegacr = cross(Omega, xc);
centforce = -cross(Omega, Omegacr);  %  centrifugal accel
corforce = -2*cross(Omega, uc);      %  Coriolis accel
uc = uc + dt*(a + corforce + centforce);

%  plot some things, occasionally  ************plotting

plotfrq = 100;
if mod(j,plotfrq) == 1
   
jp = jp + 1

if jp == 2
   figure(1)
   delete(harrowv0);
   delete(harrowv0t);
   delete(harrowo);
   delete(harrowot);
   figure(2)
   delete(harrowv0c);
   delete(harrowv0ct);
end

figure(1)
hold on


if exist('htim') == 1
  delete(htim);
end
htim = text(-0.8, -1.2, ['time = ', num2str(t(j),3)]);

if exist('hball') == 1
   delete(hball)
end
hball = plot3( x(1), x(2), x(3), '.k', 'markersize', 16);

plot3(x(1), x(2), x(3), '+b')
 
%  plot the rotating frame
xa = 1*cos(angle);
ya = 1*sin(angle);
xb = -ya;
yb = xa;

%  save data and plot the surface position in the inertial frame
xcn(jp,:) = x;
if jp >=1
for n=1:jp
    ang(n) = -(jp-n)*plotfrq*dt*Omega(3);
end
end

if exist('hppr') == 1
delete(hppr); delete(hppr1);
end

for n=1:jp
xr1(n,1) = xcn(n,1)*cos(ang(n)) + xcn(n,2)*sin(ang(n));
xr1(n,2) = xcn(n,2)*cos(ang(n)) - xcn(n,1)*sin(ang(n));
xr1(n,3) = xcn(n,3);
end
hppr = plot3(xr1(:,1), xr1(:,2), 0*xr1(:,2), 'r.');
hppr1 = plot3(xr1(:,1), xr1(:,2), xr1(:,3), 'r+');


if exist('hpp') == 1
   delete(hpp); delete(ht1); delete(ht2);
end

hpp = plot3([xa 0 xb], [ya 0 yb], [ 0 0 0], 'r', 'LineWidth', 1.4);
plot3(x(1), x(2), 0, '.b')
ht1 = text(xa, ya, 0, 'Xr', 'Color', 'r'); 
ht2 = text(xb, yb, 0, 'Yr', 'Color', 'r');
view(30., 60.)

hpa2 = plot3([x(1) x(1)], [x(2) x(2)], [0 x(3)],...
   'b-', 'LineWidth', 0.5);

drawnow
pause(0.1)
delete(hpa2);


figure(2)
hold on

if exist('htimc') == 1
  delete(htimc);
end 
htimc = text(-0.8, -1.2, ['time = ', num2str(t(j),3)]);

if exist('hballr') == 1
   delete(hballr)
end
hballr = plot3(xr(1), xr(2), xr(3), '.k', 'markersize', 16)

plot3(  xr(1), xr(2), x(3), '+r', 'markersize', 8)
plot3( xr(1), xr(2), 0., '.r', 'markersize', 12)
plot3( xc(1), xc(2), x(3), '+g', ...
    xc(1), xc(2), 0., '.g')
plot3([0 0 1], [1 0 0], [0 0 0], 'r', 'LineWidth', 1.4)  %  the extra axis
text(1, 0, 0, 'Xr', 'Color', 'r'); 
text(0, 1, 0, 'Yr', 'Color', 'r')
view(30., 60.)
axis([-1 1 -0.5 1.5 0 1])
set(gca, 'xcolor', [1 0 0], 'ycolor', [1 0 0], 'zcolor', [1 0 0])

hpa = plot3([xr(1) xr(1)], [xr(2) xr(2)], [0 xr(3)],...
   'r-', 'LineWidth', 0.5)
drawnow
pause(0.1)
delete(hpa)


pause    %  hit a key to continue

end  %  end of if on whether to plot or not

end    %  the loop on time, while z > 0

%
function hp3 = arrow3(v,x0,radius,l,scale,ntet,c)
%  arrow3( V , X0 , R , L , Scale , N )
%
%          DRAW A 3-D ARROW (as a segment plus a cone)
%
% V   vector to be represented as an arrow
% X0  point where vector start -            (default [0 0 0])
% R   arrow width  (cone radius) (in units of V)
%                                           (default 0.2)
% L   arrow length (cone height) (in units of V)
%     if L>1 the segment is not plotted     (default 0.3)
% Scale is to scale the vector              (default 1)
% N   is the resolution (number of lines to draw a cone)
%                                           (default 12)
%
if nargin<2 x0=[0 0 0]; end
if nargin<3 radius=0.2; end
if nargin<4 l=0.3; end
if nargin<5 scale=1; end
if nargin<6 ntet=12; end
if nargin<7 c=[1 1 1]; end
%create circle normal to vector v
V=norm(v);
salpha=v(3)/V;calpha=sqrt(v(1)*v(1)+v(2)*v(2))/V;
sbeta=0;cbeta=1;
if calpha~=0,sbeta=v(2)/V/calpha;cbeta=v(1)/V/calpha;end
tet=(0:pi/ntet:2*pi)';ct=radius*V*cos(tet);st=radius*V*sin(tet);
x(:,1)=+ct*salpha*cbeta+st*sbeta;
x(:,2)=+ct*salpha*sbeta-st*cbeta;
x(:,3)=-ct*calpha;
ntet2=2*ntet;
%graphic tools
v=v*scale;x=x*scale;
p=x0+v;
%b=axis;d(1:3)=b(2:2:6)-b(1:2:5);d=d/max(d);c=d;
for i=1:3,x(:,i)=x0(i)+(1-l)*v(i)+c(i)*x(:,i);end
if l<1, 
hp3 = plot3(x(:,1),x(:,2),x(:,3),'b-',[p(1)*ones(ntet,1) x(1:2:ntet2,1)]',...
   [p(2)*ones(ntet,1) x(1:2:ntet2,2)]',[p(3)*ones(ntet,1) x(1:2:ntet2,3)]',...
   'b-',[x0(1) p(1)],[x0(2) p(2)],[x0(3) p(3)],'b-');
else
hp3 = plot3(x(:,1),x(:,2),x(:,3),'b-',[p(1)*ones(ntet,1) x(1:2:ntet2,1)]',...
   [p(2)*ones(ntet,1) x(1:2:ntet2,2)]',[p(3)*ones(ntet,1) x(1:2:ntet2,3)]',...
   'b-');
end