Code covered by the BSD License  

Highlights from
a Coriolis tutorial

image thumbnail

a Coriolis tutorial

by

 

16 Feb 2003 (Updated )

Scripts that demonstrate aspects of rotating reference frames, Coriolis force and geostrophy.

rotation_1.m
%  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 particle 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 (relevent to the motion of fluids 
%  on the rotating earth) see 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 bets are off. 

clear
clear memory
close all


path(path, 'c:/jpsource/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) =  {'Rotation:' };
str2(2) =  {'  '};
str2(3) =  {'Examine the kinematics of a rotating coordinate system'};
str2(4) =  {'by showing the path of a free particle 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 particle 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 (nearly)'};
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.'};
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

initpos = 0.3;
initang = 70.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 = [0 0 1];   %  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 470 470])
set(gca,'xColor', [0 0 1], 'ycolor', [0 0 1], 'zcolor', [0 0 1])
xlabel('X'); ylabel('Y'), zlabel('Z')
htit = title('inertial frame IC for position and velocity')
axis([-1 1 -0.5 1.5 0 1])
view(30., 60.)
grid

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

hold on

harrowv0 = arrow3h(u0, x0, 0.03, 0.08, 1, 20);
u0pos = u0 + x0;
harrowv0t = text(u0pos(1), u0pos(2), u0pos(3), '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, 20);
set(harrowo, 'color', 'r')
harrowot = text(0., 0., 1.1, '\Omega', 'color', 'r');



figure(2)
clf reset
hold on

set(gcf, 'position', [500 100 470 470])
xlabel('X`'); ylabel('Y`'); zlabel('Z')
htit2 = title('rotating frame IC for position and velocity')


corstuff = 'dV`/dt = -(g + 2 \Omega x V` + \Omega x \Omega x X`)'; 
text(-0.8, -0.4, 1.8, corstuff, 'Color', 'g')
corstuff = 'V`_0 = V_0 - \Omega x X`'
text(-0.8, -0.4, 1.5, corstuff, 'Color', 'g')

harrowv0c = arrow3h(uc, x0, 0.03, 0.08, 1, 20);
set(harrowv0c, 'Color', 'r')
u0pos = uc + x0;
harrowv0ct = text(u0pos(1), u0pos(2), u0pos(3),...
   'V`_{0}', 'color', 'r')
hball = plot3( x(1), x(2), x(3), '.k', 'markersize', 16);
axis([-1 1 -0.5 1.5 0 1])
grid
view(30., 60.)
if j == 1
httt = text(0.3, 0., 0., 'hit any key to continue');
drawnow
pause
end

figure(1)
delete(htit)
tstr = str2mat('  ', ...
   ' an inertial frame (blue) and a rotating frame (red)');
htt = title(tstr, 'VerticalAlignment', 'top'); 

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

drawnow

jp = 0;
j = 0;

while x(3) >= 0
   
j = j + 1;
t(j) = (j-1)*dt;
      
%  step ahead the inertial frame variables 
a = [0 0 g];
uh = u + 0.5*dt*a;
x = x + 0.5*dt*(u + uh);
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
xch = xc + 0.5*dt*uc;
corforceh = -2*cross(Omega, uc);      %  Coriolis force
Omegacr = cross(Omega, xc);
centforceh = -cross(Omega, Omegacr);  %  centrifugal force
uch = uc + 0.5*dt*(corforceh + centforceh + a) ;

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

%  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.4, -1.2, ['time/(2\pi/\Omega) = ', ...
           num2str(t(j)/(2*pi/norm(Omega)),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, '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 
htimc = text(-0.4, -1.2, ['time/(2\pi/\Omega) = ', ...
           num2str(t(j)/(2*pi/norm(Omega)),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.)
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)

end  %  end of if on whether to plot or not
end    %  the loop on time, while z > 0

Contact us