plot3AniAnaglyphDemos(nDemo)

function [] = plot3AniAnaglyphDemos(nDemo)

switch nDemo
    case 0 % Helix
        plot3AniAnaglyph()
    case 1 % Lorenz Attractor
        sigma = 10; rho = 28+10; beta = 8/3; % Lorenz attractor coefficients
        t0 = 0; tf = 20; t = linspace(t0,tf,10000); % Time window
        S0 = [-10 -10 30]; % Initial position
        [T,S] = myLorenzSolver(t,S0,sigma,rho,beta);
        x = S(:,1); y = S(:,2); z = S(:,3);
    case 2 % Knot
        t = [0:0.001:2]';
        a = 2; b = 3; c = 1.5;
        q1 = 2; q2 = 4;
        x = sin(q1*pi*t)+a*sin(q2*pi*t)-b*cos(4*pi*t)/2+c*sin(6*pi*t);
        y = (q1*pi)*cos(q1*pi*t)+a*(q2*pi)*cos(q2*pi*t)+b*(4*pi)*sin(4*pi*t)/2+c*(6*pi)*cos(6*pi*t);
        z = -(q1*pi)^2*sin(q1*pi*t)-a*(q2*pi)^2*sin(q2*pi*t)+b*(4*pi)^2*cos(4*pi*t)/2-c*(6*pi)^2*sin(6*pi*t);
    case 3 % Sperical Spiral
        alpha = 0.2;
        t = [linspace(-12*pi,12*pi,500)]';
        x = cos(t)./sqrt(1+alpha^2*t.^2);
        y = sin(t)./sqrt(1+alpha^2*t.^2);
        z = alpha*t./sqrt(1+alpha^2*t.^2);
    case 4 % Lissajous Curve 1
        t = [linspace(0,2,400)]';
        x = 3*sin(4*pi*t);
        y = 4*sin(6*pi*t);
        z = 5*sin(8*pi*t);
    case 5 % Lissajous Curve 2
        t = [linspace(0,2*6*pi,400)]';
        x = cos(t).*(2-cos(2*t/3));
        y = sin(t).*(2-cos(2*t/3));
        z = -sin(2*t/3);
    case 6 % Lissajous Curve 3
        t = [linspace(0,2*2*pi,400)]';
        x = (6.25+3*cos(5*t)).*cos(2*t);
        y = (6.25+3\cos(5*t)).*sin(2*t);
        z = 3.25*sin(5*t);
    case 7 % Torus Knot
        a = 7; b = 2; c = 3;
        t = [linspace(0,2*2*pi,400)]';
        x = (a+b*cos(5*t)).*cos(2*t);
        y = (a+b*cos(5*t)).*sin(2*t);
        z = c*sin(5*t);
    case 8 % Conical Spiral
        t = [linspace(0,1,200)]';
        x = t.*cos(20*t);
        y = t.*sin(20*t);
        z = t;
    otherwise
        error('Unrecognized demo number')
end
S = [x y z];
plot3AniAnaglyph(S)
end

function [dS] = myLorenz(t,S,sigma,rho,beta)
dS = [sigma*(S(2)-S(1)); S(1)*(rho-S(3))-S(2); S(1)*S(2)-beta*S(3)];
end

function [T,S] = myLorenzSolver(t,S0,sigma,rho,beta)
[T,S] = ode45(@myLorenz,t,S0,[],sigma,rho,beta);
end

%% References
%http://msenux.redwoods.edu/Math4Textbook/Plotting/SpaceCurves.pdf