north_coriolis.m

function [t,r] = coriolis(L,m)
  % 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,F,m,R)
    Fx=F(1); Fy=F(2); Fz=F(3); % set each Force component
    g=9.81; % acceleration due to gravity m/s^2
    we=2*pi/23.934/3600; % rotation of Earth (each day is 23.934 hours long)
    ddr=zeros(6,1); % initialize ddr

    ddr(1) = r(2); % x North(+) South (-)
    ddr(2) = 2*we*r(2).*sin(r(1)/R)-r(6).*cos(r(1)/R)+Fx/m; % dx/dt 
    ddr(3) = r(4); % y West (+) East (-)
    ddr(4) = -2*we*(r(2).*sin(r(1)/R)-r(6).*cos(r(1)/R))+Fy/m; % dy/dt
    ddr(5) = r(6); % z altitude
    ddr(6) = -2*we*r(4).*cos(r(1)/R)+Fz/m-g; % dz/dt
  end

  R=6378.1e3; % radius of Earth in m
  % Applied force is 100 N in Northern direction and z-direction equals force of gravity
  % (neutrally buoyant travel)
  F= [ 100 0 9.81*m]; %Applied force in N
  [t,r]=ode45(@(t,r) myode(t,r,[100 0 9.81*m],m, R),[0 200], [L*pi/180*R 0 0 0 10 0]);
  figure()
  % Plot Coriolis effect for deviation Westward
  plot(r(:,3),1e-3*(r(:,1)-r(1,1)))
  xlabel('West (+m)','Fontsize',18)
  ylabel('North (+km)','Fontsize',18)
  text(0,0, 'Storrs, CT')
  title('Coriolis acceleration with north force')

  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