Skip to content

merge #9

Merged
merged 2 commits into from Apr 25, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
Binary file added lecture_01/octave-workspace
Binary file not shown.
Binary file added lecture_02/octave-workspace
Binary file not shown.
Binary file modified lecture_20/octave-workspace
Binary file not shown.
Binary file modified lecture_21/octave-workspace
Binary file not shown.
Binary file added lecture_22/octave-workspace
Binary file not shown.
6 changes: 6 additions & 0 deletions lecture_23/.ipynb_checkpoints/lecture_23-checkpoint.ipynb
@@ -0,0 +1,6 @@
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
32 changes: 32 additions & 0 deletions lecture_23/boussinesq_lookup.m
@@ -0,0 +1,32 @@
function sigma_z=boussinesq_lookup(q,a,b,z)
% function that determines stress under corner of an a by b rectangular platform
% z-meters below the platform. The calculated solutions are in the fmn data
% m=fmn(:,1)
% in column 2, fmn(:,2), n=1.2
% in column 3, fmn(:,2), n=1.4
% in column 4, fmn(:,2), n=1.6

fmn= [0.1,0.02926,0.03007,0.03058
0.2,0.05733,0.05894,0.05994
0.3,0.08323,0.08561,0.08709
0.4,0.10631,0.10941,0.11135
0.5,0.12626,0.13003,0.13241
0.6,0.14309,0.14749,0.15027
0.7,0.15703,0.16199,0.16515
0.8,0.16843,0.17389,0.17739];

m=a/z;
n=b/z;
if n < 1.3
f=fmn(:,2);
elseif n > 1.5
f=fmn(:,4);
else
f=fmn(:,3);
end
[~,i_fit]=sort(abs(m-fmn(:,1)));
x=fmn(i_fit(1:4),1);
y=f(i_fit(1:4));
f_out = Newtint(x,y,m);
sigma_z=q*f_out;
end
63 changes: 63 additions & 0 deletions lecture_23/coriolis.m
@@ -0,0 +1,63 @@
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

29 changes: 29 additions & 0 deletions lecture_23/eulode.m
@@ -0,0 +1,29 @@
function [t,y] = eulode(dydt,tspan,y0,h,varargin)
% eulode: Euler ODE solver
% [t,y] = eulode(dydt,tspan,y0,h,p1,p2,...):
% uses Euler's method to integrate an ODE
% input:
% dydt = name of the M-file that evaluates the ODE
% tspan = [ti, tf] where ti and tf = initial and
% final values of independent variable
% y0 = initial value of dependent variable
% h = step size
% p1,p2,... = additional parameters used by dydt
% output:
% t = vector of independent variable
% y = vector of solution for dependent variable
if nargin<4,error('at least 4 input arguments required'),end
ti = tspan(1);tf = tspan(2);
if ~(tf>ti),error('upper limit must be greater than lower'),end
t = (ti:h:tf)'; n = length(t);
% if necessary, add an additional value of t
% so that range goes from t = ti to tf
if t(n)<tf
t(n+1) = tf;
n = n+1;
end
y=zeros(n,length(y0));
y(1,:)=y0;
for i = 1:n-1 %implement Euler's method
y(i+1,:) = y(i,:) + dydt(t(i),y(i,:),varargin{:})'*(t(i+1)-t(i));
end