coriolis.m

function [t,r] = coriolis(L)
  % In class we ran this function using L=41.8084 (the latitude of Storrs, CT). This
  % function takes the latitude (L) in degrees and mass (m) and calculates the trajectory
  % of a particle with a 100 N load directed North. The initial conditions are set as L
  % (in radians) * radius of Earth, 0 m/s initial x-velocity, 0 m E-W position (add to
  % -72.261319 degrees for longitude of Storrs), 0 m/s initial E-W velocity, 10 m initial
  % altitude, 0 m/s initial z-velocity [L*pi/180*R 0 0 0 10 0], the force is given as 100
  % N North, 0 West and 9.81*m z (neutrally buoyant) [100 0 9.81*m]
  %
  % the output of myode is ddr=[dx/dt d2x/dt2 dy/dt d2y/dt2 dz/dt d2z/dt2]' and the input
  % to myode is r=[x dx/dt y dy/dt z dz/dt]'
  % using ode23 solver solves for r as a function of time, here we solve from 0 to 200 s
  % r(:,1) = x (the north-south position from 0 to 200 s)
  % r(:,3) = x (the West-East position from 0 to 200 s)
  % r(:,5) = x (the altitude from 0 to 200 s)

  % define ordinary differential equation in terms of 6 first order ODE's
  function ddr = myode(t,r,R,L)
    g=9.81; % acceleration due to gravity m/s^2
    l=10; % 10 m long cable
    we=2*pi/23.934/3600; % rotation of Earth (each day is 23.934 hours long)
    ddr=zeros(4,1); % initialize ddr

    ddr(1) = r(2); % x North(+) South (-)
    ddr(2) = 2*we*r(4).*sin(L)-g/l*r(1); % dx/dt 
    ddr(3) = r(4); % y West (+) East (-)
    ddr(4) = -2*we*(r(2).*sin(L))-g/l*r(3); % dy/dt
  end

  R=6378.1e3; % radius of Earth in m
  L=L*pi/180;
  [t,r]=ode45(@(t,r) myode(t,r,R,L),[0 30000], [1 0 0 0 ]);
  figure()
  z=-sqrt(10^2-r(:,1).^2-r(:,3).^2);
  figure(1)
  we=2*pi/23.934/3600; % rotation of Earth (each day is 23.934 hours long)
  
  plot(t,tan(we*sin(L)*t),t,-r(:,3)./r(:,1),'.')
  xlabel('time (s)','Fontsize',18)
  ylabel('-y/x','Fontsize',18)
  % Plot Coriolis effect path
  figure(2)
  title('Path at 0 hr, 4.1 hrs, 8.3 hrs','Fontsize',24) 
  N=length(t);
  i1=[1:100]; i2=[floor(N/2):floor(N/2)+100]; i3=[N-100:N];
  plot3(r(i1,1),r(i1,3),z(i1))
  hold on
  plot3(r(i2,1),r(i2,3),z(i2),'k-')
  plot3(r(i3,1),r(i3,3),z(i3),'g-')
  xlabel('X (m)','Fontsize',18)
  ylabel('Y (m)','Fontsize',18)
  zlabel('Z (m)','Fontsize',18)
  title('Coriolis acceleration Foucalt Pendulum')

%  figure()
%  % Plot Eotvos effect for deviation upwards
%  plot(1e-3*(r(:,1)-r(1,1)),r(:,5))
%  xlabel('North (+km)','Fontsize',18)
%  ylabel('Altitude (+m)','Fontsize',18)
%  title('Eotvos effect with north force')
%
end