
function atten_movie

% Plots E-field as function of depth and illustrates
% skin depth effect. Animated as movie

clear all
close all

mu = 4*pi*0.0000001
sigma = 0.01
freq = [10,40]
omega = 2*pi*freq
c=3e4
tc=0.2
rc=-(1-tc)

xg =[-2.5,2.5]
zg =[0,0]

% Set up array of depths
iz=0
for zdum =-5000:10:5000;
   iz=iz+1;
   z(iz) = zdum;
end
nz = length(z)

%k(1) = (1-j)*sqrt((mu*sigma*omega(1))/2)
%k(2) = (1-j)*sqrt((mu*sigma*omega(2))/2)

k_air(1)  =  omega(1)/c;
k_air(2)  =  omega(2)/c;
real_k(1) =  sqrt((mu*sigma*omega(1))/2)
imag_k(1) = -sqrt((mu*sigma*omega(1))/2)

real_k(2) =  sqrt((mu*sigma*omega(2))/2)
imag_k(2) = -sqrt((mu*sigma*omega(2))/2)

%Start time step loop

% With 2 frequencies, need to account for separate frequencies

n=60;
M=moviein(n);
for it=1:n
  t=(it)/(freq(1)*n);
  for iz=1:nz
     if z(iz)>= 0
        ex_t(1,iz) = tc*cos(imag_k(1)*z(iz)+omega(1)*t)*exp(-real_k(1)*z(iz));
        ex_t(2,iz) = tc*cos(imag_k(2)*z(iz)+omega(2)*t)*exp(-real_k(2)*z(iz));
        ex_i(1,iz) = NaN;
        ex_i(2,iz) = NaN;
        ex_r(1,iz) = NaN;
        ex_r(2,iz) = NaN;
        
        
     else
        ex_i(1,iz) = cos(-k_air(1)*z(iz)+omega(1)*t);
        ex_i(2,iz) = cos(-k_air(2)*z(iz)+omega(2)*t);
        ex_r(1,iz) = rc*cos(k_air(1)*z(iz)+omega(1)*t);
        ex_r(2,iz) = rc*cos(k_air(2)*z(iz)+omega(2)*t);
        
        ex_t(1,iz) = NaN;
        ex_t(2,iz) = NaN;

        %realex_t(1,iz) = cos(imag_k(1)*z(iz)+omega(1)*t);
        %realex_t(2,iz) = cos(imag_k(2)*z(iz)+omega(2)*t);
     end
  end
  %subplot(1,2,1)
  plot(-1+ex_t(1,:),-z*0.001,'b');
  hold on
  plot(1+ex_t(2,:),-z*0.001,'b');
  
  % Plot sum of I+R in air
  %plot(-1+ex_i(1,:)+ex_r(1,:),-z*0.001,'b');
  %plot(1+ex_i(2,:)+ex_r(2,:),-z*0.001,'b');
  
  % Plot I and R separately
  plot(-1+ex_r(1,:),-z*0.001,'r');
  plot( 1+ex_r(2,:),-z*0.001,'r');
  plot(-1+ex_i(1,:),-z*0.001,'b');
  plot( 1+ex_i(2,:),-z*0.001,'b');
  
  
  plot(xg,zg,':')
  hold off
  axis([-2.1,2.1, -5,5])
  title('MT plane wave: f=10 Hz and 40 Hz: rho = 100 ohm.m')
  xlabel ('E_x')
  ylabel ('depth (km)')
  
  %subplot(1,2,2)
  %plot(realex_t(2,:),-z*0.001,'b');
  %hold on
  %plot(xg,zg,':')
  %hold off
  %axis([-1.1,1.1, -5,5])
  %title('MT plane wave: f=10 Hz : rho = 10 ohm.m')
  %xlabel ('E_x')
  %ylabel ('depth (km)')

  
  M(:,it)=getframe;
end
  
movie(M,50)