function stokes_soln(N,M)

global plot_soln example

if nargin == 0
    N = 20;
    M = 40;
end

% Don't plot solution if there are too many vectors to 
% see anything.
plot_soln = N <= 256;

if (plot_soln)
    close all
end

% ---------------------------------------------
% Define example domain
% Example 1, 3 : Both Math 427/527
% Example 2 : Math 527
% ---------------------------------------------
example = 1;

if (example == 1)    
    % Domain for uv [ax,bx]x[ay,by]
    u0 = 0;
    u1 = 1;
    v0 = 0;
    v1 = 2*pi;
    
    % Domain for line : t in [t0,t1];
    t0 = 0;
    t1 = 2*pi;
    P = M;
elseif (example == 2)
    % Domain for uv [ax,bx]x[ay,by]
    u0 = 0;
    u1 = 2;
    v0 = -2;
    v1 = 2;
    
    % Domain for line : t in [t0,t1];
    t0 = 0;
    t1 = 2*(2 + 4);
    P = 2*(M + N);
end
domain_uv = [u0,u1,v0,v1];
domain_t = [t0,t1];

% Plot the surface;  return vectors of u and v
[u,v] = plot_surface(domain_uv,N,M);

hold on;

% Plot the curve; return vector of t values
t = plot_curve(domain_t,P);

% Compute F and curl F
F = get_F(t);
curl_F = get_curl_F(u,v);

% Plot the vectors fields
if (plot_soln)
    [x,y,z] = get_r_curve(t);
    plot_vectors(x,y,z,F);
    
    [x,y,z] = get_r_surface(u,v);
    plot_vectors(x,y,z,curl_F);
end

% Compute surface integral
S_surf = surface_integral(u,v,curl_F);

% Compute line integral
S_curve = line_integral(t,F);

fprintf('%20s %24.16f\n','Surface integral',S_surf);
fprintf('%20s %24.16f\n','Line integral',S_curve);
fprintf('%20s %24.4e\n','Difference',abs(S_curve-S_surf));

if (plot_soln)
    daspect([1,1,1]);
    xlabel('x','fontsize',16);
    ylabel('y','fontsize',16);
    zlabel('z','fontsize',16);

    view(3)
    set(gca,'box','on');
    set(gca,'Clipping','off');
    grid on;
    shg
end

end


% -----------------------------------------
% Vector field F and curl F definitions
% -----------------------------------------
function F = get_F(t)

global example

[x,y,z] = get_r_curve(t);

if (example == 1)
    F{1} = exp(y);
    F{2} = exp(z);
    F{3} = exp(x);
elseif (example == 2)
    % TODO : code vector field F
    error('TODO : Code vector field F');
end

end

function curl_F = get_curl_F(u,v)

global example

[x,y,z] = get_r_surface(u,v);

curl_F = cell(1,3);
if (example == 1)
    curl_F{1} = -exp(z);
    curl_F{2} = -exp(x);
    curl_F{3} = -exp(y);
elseif (example == 2)
    error('TODO : Code curl F');
end
    
end

% ------------------------------------------------
% Definitions of the curve r(t) and tangents
% ------------------------------------------------
function [x,y,z] = get_r_curve(t)

global example

if (example == 1)
    x = cos(t);
    y = sin(t);
    z = zeros(size(t));
elseif (example == 2)
    x = zeros(size(t));
    y = zeros(size(t));
    z = zeros(size(t));
    
    tv = [0, 2, 4, 2, 4];  % Edge lengths
    tsum = cumsum(tv);     % Position in t-vector
    
    % Bottom edge
    m = (tsum(1) < t) & (t <= tsum(2));
    tm = t(m) - tsum(1);
    x(m) = tm;
    y(m) = x(m).^2;
    z(m) = -2*ones(size(tm));

    % Right edge
    m = tsum(2) < t & t <= tsum(3);
    tm = t(m) - tsum(2);
    error('TODO : Code r(t) along right edge');
    
    % Top edge
    m = tsum(3) < t & t <= tsum(4);
    tm = t(m) - tsum(3);
    error('TODO : Code r(t) along top edge');
    
    % Left edge
    m = tsum(4) < t & t <= tsum(5);
    tm = t(m) - tsum(4);
    error('TODO : Code r(t) along left edge');
end

end

function [tangents,curve_element] = curve_tangents(t)

global example

% Compute r'(t)
if (example == 1)
    rp1 = -sin(t);
    rp2 = cos(t);
    rp3 = zeros(size(t));
elseif (example == 2)
    rp1 = zeros(size(t));
    rp2 = zeros(size(t));
    rp3 = zeros(size(t));
    
    tv = [0, 2, 4, 2, 4];
    tsum = cumsum(tv);
    
    % Bottom edge
    m = (tsum(1) < t) & (t <= tsum(2));
    tm = t(m) - tsum(1);
    rp1(m) = ones(size(tm));
    rp2(m) = 2*tm;
    rp3(m) = zeros(size(tm));

    % Right edge
    m = tsum(2) < t & t <= tsum(3);
    tm = t(m) - tsum(2);
    error('TODO : Code tangent to curve along right edge (not normalized)');
    
    % Top edge
    m = tsum(3) < t & t <= tsum(4);
    tm = t(m) - tsum(3);
    error('TODO : Code tangent to curve along top edge (not normalized)');
    
    % Left edge
    m = tsum(4) < t & t <= tsum(5);
    tm = t(m) - tsum(4);
    error('TODO : Code tangent to curve along left edge (not normalized)');
end

rp{1} = rp1;
rp{2} = rp2;
rp{3} = rp3;

[tangents,curve_element] = normalize(rp);

end


% ------------------------------------------------------
% Definitions of the surface r(u,v) and surface normals
% ------------------------------------------------------
function [x,y,z] = get_r_surface(u,v)

global example

[um,vm] = meshgrid(u,v);

if (example == 1)
    x = um.*cos(vm);
    y = um.*sin(vm);
    z = zeros(size(vm));
elseif (example == 2)
    error('TODO : Code surface');
end

end

function [n,area_element] = surface_normals(u,v)

global example

[um,vm] = meshgrid(u,v);

% Compute r_u cross r_v
if (example == 1)
    N{1} = zeros(size(um));
    N{2} = zeros(size(um));
    N{3} = ones(size(um));
elseif (example == 2)
    error('TODO : Code surface normal N (not normalized)');
end

[n,area_element] = normalize(N);

end


% -----------------------------------------
% Curve and surface integrals
% -----------------------------------------
function S = line_integral(t,F)

[tangents,curve_element] = curve_tangents(t);

F_dot_dr = curve_element.*(tangents{1}.*F{1} + ...
    tangents{2}.*F{2} + tangents{3}.*F{3});

S = sum(F_dot_dr(:));

end

function S = surface_integral(u,v,curl_F)

[n,area_element] = surface_normals(u,v);

curl_F_dot_n_dA = area_element.*(n{1}.*curl_F{1} + ...
    n{2}.*curl_F{2} + n{3}.*curl_F{3});

S = sum(curl_F_dot_n_dA(:));

end


% ---------------------------------------------------
% Graphics and utilities - Nothing to do below
% ---------------------------------------------------
function plot_vectors(x,y,z,F)

global plot_soln

% Vectors are scaled by 2 for better visualization
s = 2;
if (plot_soln)  
    quiver3(x,y,z,F{1},F{2},F{3},s,'k');
end

end


function [u,v] = plot_surface(domain_uv,N,M)

global plot_soln

ax = domain_uv(1);
bx = domain_uv(2);
ay = domain_uv(3);
by = domain_uv(4);

du = (bx-ax)/N;
dv = (by-ay)/M;

ue = linspace(ax,bx,N+1);
uc = ue(1:end-1) + du/2;

ve = linspace(ay,by,M+1);
vc = ve(1:end-1) + dv/2;

[xep,yep,zep] = get_r_surface(ue,ve);

% Plot the surface
if (plot_soln)
    p = patch(surf2patch(xep,yep,zep));
    set(p,'edgecolor','none');
    set(p,'facecolor',[1,0,0]*0.8);
    set(p,'facealpha',0.8);
end
    
u = uc;
v = vc;

end

function t = plot_curve(domain_t,P)

global plot_soln

t0 = domain_t(1);
t1 = domain_t(2);

dt = (t1-t0)/P;
te = linspace(t0,t1,P+1);
tc = te(1:end-1) + dt/2;
tc1 = tc;
tc1(end+1) = tc(1);
[x,y,z] = get_r_curve(tc1);

if (plot_soln)  
     plot3(x,y,z,'k','linewidth',3);
end

t = tc;

end

function [vunit,v_scale] = normalize(v)

% Compute ||v||
g = sqrt(v{1}.^2 + v{2}.^2 + v{3}.^2);

% Normalize to get unit tangents
vunit = cell(1,3);
for i = 1:3
    vunit{i} = v{i}./g;
end

v_scale = g;

end

