function prob5()

clf

% Compute coefficients in the Legendre series
A = compute_coeffs();

% --------------------
% Numerical parameters
% --------------------
M = 80;
N = 160;
nmax = 12;   % Terms in the series - A0 through A12


% Phi coordinates
pvec = linspace(0,pi,M);

% r coordinates
rvec = linspace(0,1,N);

% --------------------
% Compute the solution
% --------------------

% Solution for vertical slices
soln = zeros(N,M);
for j = 1:N   % r loop
    r = rvec(j);
    y = zeros(size(pvec));
    for n = 0:nmax
        k = n + 1;
        error('TODO : Compute the series solution using coefficients computed above');
        % Use 'my_legendre', defined below.        
    end
    soln(j,:) = y(:);
end

plot_solution(M,N,soln,1);
hold on;

% Solution for horizontal plate at the equator
soln = zeros(N,2*M);
for j = 1:N   % r loop
    r = rvec(j);
    y = 0;
    for n = 0:nmax
        k = n + 1;
        error('TODO : Compute the series solution for phi=pi/2');
    end
    soln(j,:) = y;
end

plot_solution(2*M,N,soln,2);

% Color axis between 0 and 1
caxis([0,1]);
h = colorbar;
set(h,'color','k','fontweight','bold');

% Make the background black
% set(gcf,'color','k');

% Make the sphere a sphere and not an ellipsoid
daspect([1,1,1]);
axis off;
view(3)

% Don't chop off sphere when rotating
set(gca,'clipping','off');

shg   % Bring current figure window to the top

end


function P = my_legendre(n,x)

% Use this function to compute the Legendre functions
% from the book.
y = legendre(n,x);
P = reshape(y(1,:),size(x));

end


function A = compute_coeffs()

error('TODO : Provide coefficients A0 through A12');

A = zeros(13,1);
% You may use your analytic expressions from WolframAlpha, or 
% compute them numerically using 'integral'.  


end

% --------------------------------------------------
% Plot the solution - NOTHING TO DO BELOW THIS LINE!
% --------------------------------------------------

function plot_solution(M,N,s,flag)

if flag == 1
    % Plot 5 vertical slices at locations in thetav
    thetav = linspace(0,2*pi,5);
    plot_slices(M,N,s,thetav);
elseif flag == 2
    % Plot horizontal plate at the equator
    plot_horizontal(M,N,s);
end

end

    
function [pout,hout] = plot_slices(M,N,s,thetav)

% phi coordinates
phiv = linspace(0,pi,M);

% r coordinates
rv = linspace(0,1,N);

[phim,rm] = meshgrid(phiv,rv);

% Plot vertical slices at constant theta values in thetev
p = zeros(size(thetav));
hv = [];
for i = 1:length(thetav)
    th = thetav(i);
    [xem,yem,zem] = map2cart(rm,th,phim);
    p(i) = patch(surf2patch(xem,yem,zem));   
    patch_graphics(p(i),s);
    hold on;
    
    % Add contour lines 
    flag = 1;  % phi varies, th is constant
    h = create_contourlines(rv,phiv,s,th,flag);
    hv = vertcat(h(:),hv);
    hold on;
end

if nargout > 0
    pout = p;
    hout = hv;
end

end

function [pout,hout] = plot_horizontal(M,N,s)

% theta coordinates at nodes
thetav = linspace(0,2*pi,M);

% r coordinates at nodes
rv = linspace(0,1,N);

[thetam,rm] = meshgrid(thetav,rv);

% Plot at constant phi = pi/2;
[xem,yem,zem] = map2cart(rm,thetam,pi/2);

p = patch(surf2patch(xem,yem,zem));
patch_graphics(p,s);

% Add contour lines
hold on;
flag = 2;    % theta varies, phi = pi/2 is constant
h = create_contourlines(rv,thetav,s,pi/2,flag);

if nargout > 0
    pout = p;
    hout = h;
end

end

function patch_graphics(p,s)

set(p,'cdata',s);
set(p,'edgecolor','none');
set(p,'facecolor','interp');
set(p,'facealpha',1);

end

% Mapping to spherical coordinates
function [xp,yp,zp] = map2cart(r,th,phi)

xp = r.*sin(phi).*cos(th);
yp = r.*sin(phi).*sin(th);
zp = r.*cos(phi);

end

function h = create_contourlines(rv,av,s,b_constant,flag)

% Input : 
%        rv  : is a vector of radial values
%        av  : is a vector of angular values that vary (theta or phi)
%        s   : solution; size = length(rv) x length(av)
% b_constant : Plays the role of angle that does NOT vary
%      flag  : flag=1 phi varies, theta is b_constant
%              flag-2 theta varies, phi is b_constant
%
% Output : 
%          h :  Vector of handles to line objects (segments making up the
%               contour line).
% 
% NOTE : Since we are plotting on both vertical and horizonal slices,
% we need to specify which angles varies (av) and which is constant
% (b_constant)
% 

smin = min(s(:));
smax = max(s(:));
if abs(smin-smax) < 1e-8
    h = [];
    % Don't add any contour lines if solution is constant
    return
end

cv = linspace(0,1,21);  % Add 21 contour lines
c = contourc(av,rv,s,cv);

st_idx = 1;
line_cnt = 0;
while (1)
    next_length = c(2,st_idx);
    line_cnt = line_cnt + 1;
    
    pdata = c(1,st_idx+1:st_idx+next_length);
    rdata = c(2,st_idx+1:st_idx+next_length);   
    if flag == 1
        % Plotting vertical slices - pdata acts like phi
        theta = b_constant;
        phi = pdata;
    elseif flag == 2
        % Plotting horizontal slices - pdata acts like theta
        theta = pdata;
        phi = b_constant;
    end    
    [xdata, ydata, zdata] = map2cart(rdata,theta,phi);
    
    h(line_cnt) = line('XData',xdata,'YData',ydata,'ZData',zdata);
    st_idx = st_idx + next_length + 1;
    if (st_idx > length(c))
        break;
    end
end

set(h,'color','k');

end