function coriolis % 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. % % by Jim Price. January, 2001. % % This may be considered public domain for educational % purposes. % 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) = {'Coriolis:' }; str2(2) = {' '}; str2(3) = {'Examine the kinematics of a rotating coordinate system by '}; str2(4) = {'showing the path of a projectile as seen in an inertial '}; str2(5) = {'reference frame (blue trajectory, Fig 1) and as seen from a'}; str2(6) = {'frame 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 '}; str2(8) = {'a centripetal acceleration and the motion is purely horizontal.)'}; str2(9) = {'The projectile has an initial horizontal speed Vo in the Y '}; str2(10) = {'direction and is subject to a centripetal acceleration at a'}; str2(11) = {'rate -\Omega^2 r = \Omega x \Omega x r as if it were sliding on a '}; str2(12) = {'parabolic surface Z = 0.5 \Omega^2 r^2/g. The green path '}; str2(13) = {'of Fig 2 (computed in the rotating frame) includes the'}; str2(14) = {'Coriolis acceleration only. This rotating frame solution '}; str2(15) = {'should be identical to the observations of the path' }; str2(16) = {'made from the rotating frame (red, Fig 2) if we'}; str2(17) = {'have fully accounted for the kinematic effects of the'}; str2(18) = {'rotating frame by the device of the Coriolis acceleration. '}; 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) = {' '}; str2(22) = {'Hit any key to continue and to step ahead the integration.'}; str2(23) = {''}; str2(24) = {'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 = 0.; dt = 0.001; kase = menu('choose which case', 'standard Vo, Omega and Xo', 'Vo x 2', ... 'Omega x 2', 'tangential Vo, balanced', 'tangential Vo, unbalanced') % the base case: if kase == 1 rotrate = 0.3; Omega = rotrate*2*pi*[0 0 1]; % set the rotation rate initpos = 0.3; x0 = [0 initpos 0.0000]; initang = 0.; speed0 = 1.5; u0 = speed0*[0 cos(initang*pi/180), sin(initang*pi/180)]; end % larger Vo: if kase == 2 rotrate = 0.3; Omega = rotrate*2*pi*[0 0 1]; % set the rotation rate initpos = 0.3; x0 = [0 initpos 0.0000]; initang = 0.; speed0 = 3.0; u0 = speed0*[0 cos(initang*pi/180), sin(initang*pi/180)]; end % larger rotation: if kase == 3 rotrate = 0.6; Omega = rotrate*2*pi*[0 0 1]; % set the rotation rate initpos = 0.3; x0 = [0 initpos 0.0000]; initang = 0.; speed0 = 1.5; u0 = speed0*[0 cos(initang*pi/180), sin(initang*pi/180)]; end % a case with tangential initial velocity if kase == 4 rotrate = 0.3; Omega = rotrate*2*pi*[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 % a case with approx tangential initial velocity, unbalanced if kase == 5 rotrate = 0.3; Omega = rotrate*2*pi*[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 unbalanced Uo end 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 450 450]) 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 projectile is moving freely') hmt2 = text(-1.5, -1, 2.0, 'and horizontally on a surface,') hmt3 = text(-1.5, -1, 1.8, 'Z = 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 = arrow3([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'); pause delete(hm); delete (hmt1); delete(hmt2); delete(hmt3); delete(hmt4); delete(hmt5); corstuff = 'dV/dt = -g \nabla Z = \Omega x \Omega x r'; 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, 50); 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, 50); 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 observed in the rotating frame (red)', ... 'and as computed in the rotating frame (green)'); htt = title(tstr, 'VerticalAlignment', 'top'); corstuff = 'dV_r/dt = - 2 \Omega x V_r '; text(-0.5, 0., 1.4, corstuff, 'Color', 'g') corstuff = 'V_r(0) = V_0 - \Omega x r' text(-0.5, 0., 1.2, corstuff, 'Color', 'g') harrowv0c = arrow3(uc, x0, 0.03, 0.08, 1, 50); set(harrowv0c, 'Color', 'r') u0pos = uc + x0; harrowv0ct = text(u0pos(1), u0pos(2), u0pos(3),... 'V_{r0}', '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.; while max(t) <= 4. j = j + 1; t(j) = (j-1)*dt; % step ahead the inertial frame variables x = x + dt*u; Omegacr = cross(Omega, x); centforce = cross(Omega, Omegacr); % centripetal force a = [0 0 g]; u = u + dt*(a + centforce); % 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 xc = xc + dt*uc; corforce = -2*cross(Omega, uc); % Coriolis force uc = uc + dt*corforce; % 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 title('an inertial frame (blue) and a rotating frame (red)') 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, '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 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, 'Xr', 'Color', 'r'); text(0, 1, 0, 'Yr', '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) pause % hit a key to continue end % end of if on whether to plot or not end % the loop on time, while z > 0 % the end of coriolis 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