function My_KModeVariationPlot(yy, k, ax)
% KModeVariationPlot: k-th mode of variation plot for FPCA/FPCder
%
% Inputs:
%   yy : FPCA or FPCder object
%   k  : optional, positive integer(s) <= no_opt, default = 1:no_opt
%   ax : optional, axes handle for plotting. If missing, creates new figure for each k.

no_opt = getVal(yy,'no_opt');

if nargin < 2 || isempty(k)
    k = 1:no_opt;
elseif any(k <= 0)
    error('k must be a positive integer!');
elseif any(k > no_opt)
    error('Any k cannot be larger than no_opt = %d', no_opt);
end

lambda = getVal(yy,'lambda');
mu = getVal(yy,'mu');
if iscell(mu)
    mu = mu{1};
end

out1 = getVal(yy,'out1');
phi = getVal(yy,'phi');
if iscell(phi)
    phi = phi{1};
end

alpha = [-2 -1 0 1 2];  % scaling for mode of variation

for i = 1:length(k)
    if nargin < 3 || isempty(ax)
        figure;
        ax_i = gca;
    else
        ax_i = ax;
    end

    cmat = cool(length(alpha));

    axes(ax_i);  % ensure we plot in the given axes
    hold(ax_i,'on');

    for j = 1:length(alpha)
        tmp = mu + alpha(j) * sqrt(lambda(k(i))) * phi(:,k(i))';
        plot(ax_i, out1, tmp, 'Color', cmat(j,:), 'LineWidth', 1.5);
    end

    xlabel(ax_i,'t');
    ylabel(ax_i,['\mu(t) + \alpha \phi_' num2str(k(i)) '(t)']);
    title(ax_i,['K-th mode of variation (k = ' num2str(k(i)) ')']);
    legend(ax_i, {'\alpha=-2','\alpha=-1','\alpha=0','\alpha=1','\alpha=2'}, 'Location','Best');
    grid(ax_i,'on');
    hold(ax_i,'off');
end

end