%  Coriolis;  Demonstrate the Coriolis force
%  by computing the motion of a particle 
%  as seen from a rotating coordinate system.
%  This is a companion to rotation.m.  The
%  plots are shown in 3-d, as in rotation.m,
%  even though the trajectories here are 2-d.
%  This uses a very simple time-stepping
%  scheme (2nd order RK) that is not highly 
%  accurate, but good enough for the graphics
%  that are the purpose here. (Would be nice to
%  reduce the flickering, but ...)
%  
%  by Jim Price. January, 2001. 
%  revised December 2002.
%
%  This is public domain for all personal, educational
%  purposes. 
%
%  This uses a funtion arrow3h to plot arrows that is a slighltly 
%  modified version of 
%

clear
clear memory
close all
format compact

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

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


str2(1) =  {'Coriolis:' };
str2(2) =  {'  '};
str2(3) =  {'   Examine the kinematics of a rotating coordinate system by showing'};
str2(4) =  {'the path of a particle as seen in an inertial reference frame'};
str2(5) =  {'(blue trajectory, Fig 1) and as seen from a reference frame'};
str2(6) =  {'that is rotating at a rate \Omega (red trajectory, Figs 1 and 2).'};
str2(7) =  {'(This differs from rotation.m in that there is assumed to be a'};
str2(8) =  {'centripetal force = -\Omega^2r and the motion is horizontal.)'};
str2(9) =  {'The particle is given an initial horizontal speed Vo and is then'}; 
str2(10) = {'subject to a centripetal force as if it were sliding on a ' }; 
str2(11) = {'frictionless parabolic surface \Phi = 0.5 \Omega^2 r^2/g.'};
str2(12) = {'   Fig. 2 shows the path as seen from the rotating reference frame'}; 
str2(13) = {'(the red path) and as computed in the rotating frame using the '};
str2(14) = {'Coriolis force only.  This rotating frame solution '};
str2(15) = {'should be identical to the observations of the path made' };
str2(16) = {'from the rotating frame (compare red and green paths) if we'};
str2(17) = {'have fully accounted for the kinematic effects of the'};
str2(18) = {'rotating frame by the Coriolis force. '};
str2(19) = {'   The circular motion at a rate 2\Omega observed and computed '};
str2(20) = {'in the rotating frame is often termed an `inertial motion`. '};
str2(21) = {'   The Foucault pendulum cases have an additional radial force'};
str2(22) = {'equal to freq^2, where freq is in units of Omega.  Vo is radial.'};
str2(23) = {'    '};
str2(24) = {'    Hit any key to continue.      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

hpn = 0; k = 0; 
g = 0.;
dt = 0.002;  %  the time step in units 1/Omega;  keep it small

kase = menu('choose which case', 'tangential Vo, balanced', ...
    'tangential Vo, Vo x 1.5', ...
    'tangential Vo, radial impulse',  'tangential Vo, Omega x 2',  ...
   'Foucault pendulum, freq = O.', 'Foucault pendulum, freq = 6*Omega')   
  
%  a balanced Vo
if kase == 1
Omega = [0 0 1];   %  set the rotation rate 
initang = 0;
initpos = 1.0;
x0 = [0 initpos 0.0000];
centa = Omega(3)^2*initpos;
speed0 = sqrt(initpos*centa);
u0 = speed0*[-1 0 0];
end

%  initial velocity unbalanced because too large:
if kase == 2
Omega = [0 0 1];   %  set the rotation rate       
initpos = 1.0;
x0 = [0 initpos 0.0000];
initang = 0.;
centa = Omega(3)^2*initpos;
speed0 = sqrt(initpos*centa);
u0 = 1.5*speed0*[-1 0 0];   %   note the extra speed
end

%  initial velocity unbalanced because of radial impulse
if kase == 3
Omega = [0 0 1];   %  set the rotation rate 
initang = 0;
initpos = 1.0;
x0 = [0 initpos 0.0000];
centa = Omega(3)^2*initpos;
speed0 = sqrt(initpos*centa);
u0 = speed0*[-1 0.5 0];    %  note the extra (radial) Uo
end

%  larger (than balanced) rotation:
if kase == 4
Omega = [0 0 1];   %  set the rotation rate       
initpos = 1.0;
x0 = [0 initpos 0.0000];
initang = 0.;
centa = Omega(3)^2*initpos;
speed0 = sqrt(initpos*centa);
u0 = speed0*[-1 0 0];
Omega = 2.*Omega;   %  increase the rotation rate       
end

%  set some default parameters for Foucault pendulum cases
freqfac = 0.; 
plotfrq = 100.;

%  Foucault pendulum, freq = 0. (an inertial motion)
if kase == 5
Omega = [0 0 1];   %  set the rotation rate 
initang = 0;
initpos = 0.;
x0 = [0 initpos 0.0000];
centa = Omega(3)^2*initpos;
speed0 = sqrt(initpos*centa);
u0 = [0 2 0];    %  note the extra unbalanced Uo
end

%  Foucault pendulum, higher freq (like a merry-go-round)
if kase == 6
freqfac = 6.;
plotfrq = round(plotfrq/freqfac);
Omega = [0 0 1];   %  set the rotation rate vector
initang = 0;
initpos = 0.;
x0 = [0 initpos 0.0000];
centa = Omega(3)^2*initpos;
speed0 = sqrt(initpos*centa);
u0 = [0 freqfac  0];    %  note the  unbalanced Uo
end


%  set the ICs:
%  the inertial frame
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.;

figure(1)
clf reset

set(gcf,'position', [25 100 470 470])
set(gca,'xColor', [0 0 1], 'ycolor', [0 0 1], 'zcolor', [0 0 1])
xlabel('X'); ylabel('Y'), zlabel('Z')
axis([-1.5 1.5 -1.5 1.5 0 1])
view(30., 60.)
grid

% define and draw an equipotential surface
alpha = 0.2;
xessize = 1.3
for mx = 1:21
   for my = 1:21
      xes(mx) = xessize*(mx-11)/10.;
      yes(my) = xessize*(my-11)/10.;
      res = sqrt(xes(mx)^2 + yes(my)^2);
      zes(my, mx) = alpha*res^2/2;
   end
end
hold on
hm = mesh(xes, yes, zes);
hmt1 = text(-1.5, -1, 2.2, 'The particle is moving horizontally');
hmt2 = text(-1.5, -1, 2.0, 'and without friction on a surface,');
hmt3 = text(-1.5, -1, 1.8, '\Phi = 0.5 \Omega^2 r^2/g, that could be');
hmt4 = text(-1.5, -1, 1.6, 'the free surface of a fluid in ');
hmt5 = text(-1.5, -1, 1.4, 'solid body rotation at the rate \Omega ');
harrowo = arrow3h([0 0 1], [0 0 0], 0.03, 0.08, 1, 50);
set(harrowo, 'color', 'r');
harrowot = text(0., 0., 1.1, '\Omega', 'color', 'r');
if j == 1 httt = text(-0.5, -1.7, 0., 'hit any key to continue'); end

pause
delete(hm); delete (hmt1); delete(hmt2); delete(hmt3); delete(hmt4);
delete(hmt5); delete(httt);


hict = text(-2, 0., 1.3, 'initial condition on position and velocity', 'Color', 'k')

hold on

u0p = u0;
s0 = norm(u0)    %  just for plotting's sake
if s0 >= 1.5
    u0p = u0*(1.5/s0);
end

harrowv0 = arrow3h(u0p, x0, 0.03, 0.08, 1, 50);
u0pos = u0p + x0;
harrowv0t = text(u0pos(1), u0pos(2), u0pos(3)+0.1, 'V_0', 'color', 'b');
hball = plot3( x(1), x(2), x(3), '.k', 'markersize', 16);

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

%  add a circle at the initial radius
dth = 2*pi/500;
for j=1:500
    th = dth*j;
    xcr(j) = initpos*cos(th);
    ycr(j) = initpos*sin(th);
end
plot(xcr, ycr, 'g:', 'linewidth', 0.6)

httt = text(1., -1., 0., 'hit any key to continue');
drawnow


delete(httt); delete(harrowv0); delete(harrowv0t); 
delete(harrowot); delete(hict); 

%  add some text that will stay of Fig 1
corstuff = 'dV/dt = -g \nabla \Phi '; 
if kase >=6
corstuff = 'dV/dt = -g \nabla \Phi - \nu^2 e_r'; 
end
text(-2, 0., 1.3, corstuff, 'Color', 'b')
corstuff = 'V(0) = V_0';
text(-2, 0., 1.0, corstuff, 'Color', 'b')

figure(2)
clf reset
hold on

set(gcf, 'position', [500 100 470 470])
xlabel('X`'); ylabel('Y`'); zlabel('Z')
clear tstr
tstr = str2mat('as observed in the rotating frame (red)', ...
   'and as computed in the rotating frame (green)');
htt = title(tstr, 'VerticalAlignment', 'top');

corstuff = 'dV`/dt = - 2 \Omega x V` '; 
if kase >=6
corstuff = 'dV`/dt = - 2 \Omega x V`  - \nu^2 e_r'; 
end
text(-2, 0., 1.0, corstuff, 'Color', 'g')
corstuff = 'V`(0) = V_0 - \Omega x X`';
text(-2, 0., 0.7, corstuff, 'Color', 'g')

harrowv0c = arrow3h(uc, x0, 0.03, 0.08, 1, 50);
set(harrowv0c, 'Color', 'r')
u0pos = uc + x0;
harrowv0ct = text(u0pos(1), u0pos(2), u0pos(3),...
   'V`_{0}', 'color', 'r');

axis([-1.5 1.5 -1.5 1.5 0 1])
grid
view(30., 60.)
drawnow

jp = 0;
j = 0;
t(1) = 0.;

%  the time loop starts here:

maxt = 2*pi;
while max(t) <= maxt
   
j = j + 1;
t(j) = (j-1)*dt;
xis(j) = x(1);
yis(j) = x(2);
xcs(j) = xc(1);
ycs(j) = xc(2);
      
%  step ahead the inertial frame variables   
xh = x + 0.5*dt*u;
Omegacr = cross(Omega, x);
centforceh = cross(Omega, Omegacr);  %  the centripetal force 
uh = u + 0.5*dt*(1. + freqfac^2)*centforceh;

x = x + 0.5*dt*(u + uh);
Omegacr = cross(Omega, xh);
centforce = cross(Omega, Omegacr);  %  the centripetal force 
u = u + 0.5*dt*(1. + freqfac^2)*(centforce + centforceh);

%  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);

%  step ahead the position and speed in the Coriolis frame
xch = xc + 0.5*dt*uc;
corforceh = -2*cross(Omega, uc);      %  Coriolis force
Omegacr = cross(Omega, xc);
centforceh = cross(Omega, Omegacr);  %  centripetal + centrifugal force
uch = uc + 0.5*dt*(corforceh + (freqfac^2)*centforceh) ;

xc = xc + 0.5*dt*(uc + uch);
corforce = -2*cross(Omega, uch);      %  Coriolis force
Omegacr = cross(Omega, xch);
centforce = cross(Omega, Omegacr);  %  centripetal + centrifugal force
uc = uc + 0.5*dt*((corforceh + corforce) + ...
    (freqfac^2)*(centforceh + centforce)) ;


%  plot some things, occasionally  
if mod(j,plotfrq) == 1;
   
jp = jp + 1;
if jp == 2
   figure(2)
   delete(harrowv0c);
   delete(harrowv0ct);
end

figure(1)
hold on
clear tstr
tstr = str2mat(' an inertial frame (blue) and a rotating frame (red)');
htt = title(tstr, 'VerticalAlignment', 'top');

if exist('htim') == 1
   delete(htim);
end
tnd = t(j)/(2*pi/Omega(3));
htim = text(-0.8, -2.8, ['time/(2\pi/\Omega) = ', num2str(tnd,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 projectile 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), 0*xr1(:,2), '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, 'X', 'Color', 'r'); 
ht2 = text(xb, yb, 0, 'Y', '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 
tnd = t(j)/(2*pi/Omega(3));
htimc = text(-0.8, -2.8, ['time/(2\pi/\Omega) = ', num2str(tnd,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, 'X', 'Color', 'r'); 
text(0, 1, 0, 'Y', 'Color', 'r')
view(30., 60.)
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)


end  %  end of if on whether to plot or not

end    %  the loop on time

%  plot some saved data

figure(3)
clf reset

ri = abs(xis.^2 + yis.^2);
rc = abs(xcs.^2 + ycs.^2);
ts = t/(2*pi/Omega(3));

subplot(2,1,1)
plot(ts, ri, ts, rc, 'r')
xlabel('time, rotation periods')
ylabel('radial distance')
legend('inertial frame', 'rotating frame')
legend boxoff

subplot(2,1,2)
plot(ts, yis, ts, ycs, 'r')
xlabel('time, rotation periods')
ylabel('Y distance')
legend('inertial frame', 'rotating frame')
legend boxoff
%
%  the end of coriolis