plotElectron(varargin)

function f = plotElectron(varargin)
% Plot electrons on conducting body after solution finishes.
% Copyright 2010 The MathWorks, Inc.

if nargin == 0
    % return function handle to be used as a plotting function
    
    f = @(x,itervals,flag) plotElectron(x,itervals,flag);
    
else
    % plots with specified arguements
    
    f = myElectronPlot(varargin{:});
    
end
end

function stop = myElectronPlot(x,itervals,flag, varargin)

if ~exist('flag','var')
    myElectronPlot([],[],'init');
    flag = 'done';
end

switch flag
    case 'iter'
        % Make updates to plot or guis as needed
    case 'interrupt'
        % Probably no action here. Check conditions to see
        % whether optimization should quit.
    case 'init'
        tic
    case 'done'
        % Cleanup of plots, guis, or final plot
        toc
        E = length(x)/3;
        % These two lines set up standard polar coordinates
        R=0:.0025:1;
        TH= 2*pi*(0:.0025:1);
        % This line finds the maximum  radius value for each angle TH
        Rmax = diag((abs(cos(TH)) + abs(sin(TH))) ./ (1 + abs(sin(TH).*cos(TH))));
        % The equation to solve is
        % z = -abs(x)  abs(y) (pyramid)
        % x^2 + y^2 +(z+1)^2 = 1 (sphere centered at [0,0,-1])
        % Solve for x and y, get the equation
        % x^2 + y^2 = x + y  x*y
        % Set x = R*cos(TH), y = R*sin(TH)
        % R^2 = R(cos(TH) + sin(TH))  R^2*sin(TH)*cos(TH)
        % R = (cos(TH) + sin(TH)) / (1 + sin(TH)*cos(TH)) modulo signs
        
        % These two lines set up polar coordinates, then scale the maximum radius
        % at each angle TH by the value of Rmax(TH)
        % This is done in a vectorized, efficient calculation
        X=R'*cos(TH) * Rmax;
        Y=R'*sin(TH) * Rmax;
        
        Z1 = - abs(X) - abs(Y); % pyramid
        Z2 = -1 - real(sqrt(1 - X.^2 - Y.^2)); % sphere
        
        surf(X,Y,Z1,'LineStyle','none'); % plot pyramid
        hold on
        surf(X,Y,Z2,'LineStyle','none'); % plot sphere
        
        set(gcf,'Color','w') % white background
        view(-44,18)
        for i = 1:E
            plot3(x(3*i-2),x(3*i-1),x(3*i),'r.','MarkerSize',25);
        end
        rotate3d
        axis equal
        drawnow
        hold off
    otherwise
        % not used
end
stop = false;
end