% Colin Joye
% 6.641, PS #13
% 5/13/03T
% Plots the time growth of a convective instability and plots the 
% corresponding dispersion relation.

% Problem P13.4c
k  = [-3:0.01:3]; % wave number
v  = 1;     % phase velocity
u  = 1;     % translation velocity
wc2= 1;     % cutoff frequency
wa = (k*u) + sqrt((k*v).^2 - 3*wc2);
wb = (k*u) - sqrt((k*v).^2 - 3*wc2);
wc = (k*u) + sqrt((k*v).^2 - 1*wc2);
wd = (k*u) - sqrt((k*v).^2 - 1*wc2);
clf(1)
hold on; grid on;
plot([0 0],[k(1) k(end)],'g');
plot(k,real(wa),'b',k,imag(wa),':r');
plot(k,real(wb),'b',k,imag(wb),':r');
plot(k,real(wc),'b',k,imag(wc),':r');
plot(k,real(wd),'b',k,imag(wd),':r');
hold off;
title('Dispersion Relation, P13.4c -- Colin Joye');
xlabel('k [normalized]');
ylabel('w_r_e_a_l (blue), w_i_m_a_g (red) [normalized]');

% Problem P13.4g
% sketch spatial ksee1 = -ksee2 for v0=0 
t = 0;     % time (look at space dependence)
x = [0:0.01:20]; % spatial vector
d = 12;     % membrane separation
k = 1;     % wave number
v = 1;     % Voltage applied
u = 1;     % tranlation velocity
w = 7;     % driven wave frequency
wc2=1;     % = e0*V0^2/sm/d^3; cufoff freq.
vp2=2;     % phase velocity squared
ksee_ = 1; % drive magnitude
%
a = u*w / (u^2 - vp2);
b = sqrt(w^2*vp2 - wc2*(u^2 - vp2))/(u^2 - vp2)^2;
ksee1 = ksee_*( cos(w*t-a*x).*cos(b*x) - a/real(b)*sin(w*t-a*x).*sin(b*x) );
ksee2 = -ksee1;
figure(2)
clf(2)
hold on; grid on;
plot(x,real(ksee1)+d/2,'b',x,real(ksee2)-d/2,'b');

% Plot again with v0 large enough for the wave mode to grow
b = 4+i*1e-1;
% need b slightly imaginary for growing waves
ksee1 = ksee_*( cos(w*t-a*x).*cos(b*x) - a/real(b)*sin(w*t-a*x).*sin(b*x) );
ksee2 = -ksee1;
plot(x,real(ksee1)+d/2,'r',x,real(ksee2)-d/2,'r');
hold off;
title('Membrane response, V0 = 0 (blue) and k complex (red) -- Colin Joye');
xlabel('Space [normalized]');
ylabel('Deflection: V0=0, k is real (blue), k complex (red)');
%