% E dipole formula for harmonic source in a whole-space (unbound media).
% Formula in Ward and Hohman (eqn 2.40), which is attributed to Banos.
%
% S.H. Ward and G.W. Hohmann. Electromagnetic theory for geophysical
% applications. In M.N. Nabighian, editor, Electromagnetic methods 
% in applied geophysics, pages 131-312. Soc. Expl. Geophys.,
% 01 edition, 1988.
%
% Matlab movie script by: 
% Kerry Key
% Scripps Institution of Oceanography
% kkey@ucsd.edu
%
% Saves 60 frames to an uncompressed avi movie file 
% that is really large (120 MB).  You can recompress the movie using 
% any number of free compression programs to a few
% megabytes.  Another option would be to turn on the avi compression in
% Matlab, but this isn't available on the Mac OS X version of Matlab. 
% 
%

    clear all
% Some parameters:
    filename = 'edipolemovie';
    sig = 3.3;
    w   = 2*pi*.25;
    II  = sqrt(-1);
    MU  = 4*pi*10^-7;

% Grid of points 
    x   = linspace(-5000,5000,200);
    y   = linspace(-5000,5000,200);
    z   = 0;
    [X Y Z ] =  meshgrid(x,y,z);
     
% compute E fields:
    R  = sqrt(X.^2 + Y.^2 + Z.^2);
    k  = sqrt( -II*w*MU*sig);
    
    f1 =  1./(4*pi*sig*R.*R.*R).*exp(-II.*k.*R);
    f2 = -k.*k.*R.*R + 3*II.*k.*R + 3;
    f3 =  k.*k.*R.*R - II.*k.*R - 1;

    ex = f1 .* ( f2.*X.*X./(R.*R) + f3) ;
    ey = f1 .*   f2.*X.*Y./(R.*R) ;
    ez = f1 .*   f2.*X.*Z./(R.*R) ;
    ax = abs(ex);
    ay = abs(ey);
    az = abs(ez);
    px = atan2(imag(ex),real(ex));
    py = atan2(imag(ey),real(ey));
    pz = atan2(imag(ez),real(ez));
    
% Polarization ellipse pmax and pmin:
    dp = py-px; 
    alpha = 1/2*atan2(2*ax.*ay.*cos(dp),(ax.^2-ay.^2));
    pmax = abs(ay.*exp(II*dp).*sin(alpha) + ax.*cos(alpha));
    pmin = abs(ay.*exp(II*dp).*cos(alpha) - ax.*sin(alpha));    
    
 % Create an avi movie file:
    mov = avifile(filename,'fps',24);
    
 % Open a figure:   
    fig=figure('position',[66 253 887 839],'color','w'); 
    set(fig,'DoubleBuffer','on');
    set(gcf,'renderer','zbuffer')
 
  % Colorplot of pmax:
    surf(Y/1000,X/1000,log10(pmax));
    shading flat; view(2)
    cb = colorbar;
    set(get(cb,'ylabel'),'string','log10(Pmax), V/ Am')
    ca = caxis;
    xlabel('y (km)');
    ylabel('x (km)');
    axis equal;
    axis tight;
    hold on;
    
  % add skin depth scale bar:
    del = pi*500*sqrt(2*pi/(w*sig)); % half wavelength 
    ndel = [ 2];  % plot up half and full wavelength
    hl = plot3([0 del*ndel]/1000, [ 0 0],[100 100]/1000,'w','linewidth',2);
    for i=0:ndel
        hdel(i+1) = plot3([del*i del*i]/1000, [ -100 100]/1000,[100 100]/1000,'w','linewidth',2);
        %hdeltext = text(2*pi*del*i,-200,num2str(i),'color','w');
    end

  % define a cycle in time with 60 steps      
    t = linspace(0,2*pi/w,60);
    t = t(1:end-1);

   % add quiver arrows that trace out the polarization ellipses in time
    dx = 5;  % skip values since the grid plotted above is too dense for the arrow plots
    offset = 20; % offset so that scaling log10 works ok
    xdist = (x(dx+1)-x(1))/1000;  
    for it = 1:length(t)
        sx = cos(w*t(it)+px(1:dx:end,1:dx:end)); 
        sy = cos(w*t(it)+py(1:dx:end,1:dx:end));
       % plot the arrows: 
        q1 = quiver(Y(1:dx:end,1:dx:end)/1000,X(1:dx:end,1:dx:end)/1000, ....
            xdist*sy.*(ay(1:dx:end,1:dx:end))./(pmax(1:dx:end,1:dx:end)),...
            xdist*sx.*(ax(1:dx:end,1:dx:end))./(pmax(1:dx:end,1:dx:end)), ...
            0,'k');
        set(q1,'tag','arrow')
      % pie chart for time increment:
        hp = pie([(it)/length(t) (length(t)-it+.001)/length(t)],{'' ''});
        % now move the pie chart to a better location
        xd = get(hp(1),'xd');
        yd = get(hp(1),'yd');
        xd = xd+3.5;
        yd = -yd-3.5;
        zd = 100*ones(size(xd));
        set(hp(1),'xd',xd,'yd',yd,'zd',zd,'facecolor','w');
        xd = get(hp(3),'xd');
        yd = get(hp(3),'yd');
        xd = xd+3.5;
        yd = -yd-3.5;
         zd = 100*ones(size(xd));
        set(hp(3),'xd',xd,'yd',yd,'zd',zd,'facecolor','k');
        caxis(ca);
        box on;
        drawnow;
        F = getframe(gcf);
        mov = addframe(mov,F);
        delete(findobj('tag','arrow'))
        delete(hp)
    end % loop over time
    mov = close(mov);
% 

