function [pl] = coriolisanim ( n, nend, alpha, beta, X0 ) %*****************************************************************************80 % %% CORIOLISANIM % % Discussion: % % This function is discussed in chapter 2 of the reference. % % Modified: % % 01 December 2009 % % Author: % % Benoit Cushman-Roisin, Jean-Marie Beckers. % % Reference: % % Benoit Cushman-Roisin, Jean-Marie Beckers, % Introduction to Geophysical Fluid Dynamics, % Academic Press, 2010. % % Parameters: % % Input, integer N, the number of time steps for a full period. % % Input, integer NEND, the last time step. % NEND = N / 4 + 1 is recommended. % % Input, real ALPHA, the scaled azimuthal initial velocity V_0/(X_0 Omega). % % Input, real BETA, the scaled radial initial velocity U_0/(X_0 Omega). % % Input, real X0, the initial distance (less then 10 for plotting reasons). % matoct(1) = 1; if matoct(1)==1 fig = figure; end pl = 1; % % Circle to plot % for i=1:101 xc(i)=10*cos(i*pi*2/100); yc(i)=10*sin(i*pi*2/100); end % % Time steps to reach nend/n revolutions with rate Omega % for j=1:nend % omega t ot=(j-1)*pi*2/n; % omega delta t dot=pi*2/n; % for line segments to see rotation of plane x1(1)=7*cos((j-1)*pi*2/n); z1(1)=7*sin((j-1)*pi*2/n); x1(2)=10*cos((j-1)*pi*2/n); z1(2)=10*sin((j-1)*pi*2/n); x2(1)=10*cos((j-1)*pi*2/n+pi/2); z2(1)=10*sin((j-1)*pi*2/n+pi/2); x2(2)=7*cos((j-1)*pi*2/n+pi/2); z2(2)=7*sin((j-1)*pi*2/n+pi/2); x3(1)=10.5; z3(1)=0; x3(2)=12; z3(2)=0; x4(1)=0; z4(1)=10.5; x4(2)=0; z4(2)=12; % solution in absolute axes in function of initial conditions x0=X0; XA(1)=x0*cos(ot)+ beta*x0*sin(ot); YA(1)=alpha*x0*sin(ot); % storage of trajectory in absolute axes XAA(j)=XA(1); YAA(j)=YA(1); % rotation of trajectory inprint on rotating axis for i=1:j-1; xn=xt(i)*cos(dot)-yt(i)*sin(dot); yn=xt(i)*sin(dot)+yt(i)*cos(dot); xt(i)=xn; yt(i)=yn; end % adding of new point to the trajectory xt(j)=XA(1); yt(j)=YA(1); % PLOT if matoct(1)==1 h = plot ( x1, z1, 'b-', ... xc, yc, ... XA, YA, 'o', ... xt, yt, '.', ... XAA, YAA, '-', ... x1(2), z1(2), 'ko', ... x2, z2, 'b-', ... x3, z3, 'b-', ... x4, z4, 'b-' ); axis('equal'); axis('off'); %set(h,'LineWidth',2); % set both lines to 2 pixels wide set(h(3),'LineWidth',3) set(h(5),'LineWidth',2) set(h(2),'LineWidth',2) set(h(1),'LineWidth',2) set(h(6),'LineWidth',6) set(h(7),'LineWidth',2) set(h(8),'LineWidth',2) set(h(9),'LineWidth',2) M(j) = getframe; else plot(xc,yc,x1,z1,'b-',XA,YA,'o',xt,yt,'.',XAA,YAA,'-',x1(2),z1(2),'ko',x2,z2,'b-',x3,z3,'b-',x4,z4,'b-'); axis([-15 15 -15 15]) axis('equal'); axis('off'); hold off pause(0.04) end end if matoct(1)==1 movie(M) end return end