% EXTRA DEMO:   POSITION OF A BICYCLE STEM VS TIME
%
% Program:      curves_of_unity.m
% Written by:   
% For:          EE 311
%               Demo Computer Lab 1
% Date:         6-19-2001
% Purpose:      This program takes the binomial expansion of
%               (cos(theta)+j*sin(theta))^n and plots the sum of 
%               the outer terms. For example, if n = 3, the binomial
%               expansion would be: 
%               (cos(theta))^3 + 3*j*(cos(theta))^2*sin(theta) ...
%               - 3*cos(theta)*(sin(theta))^2 - j*(sin(theta))^3
%
%               From this, there would be plots generated of:
%               z(1) = (cos(theta))^3 - j*(sin(theta))^3, and
%               (2) = 3*j*(cos(theta))^2*sin(theta)- 3*cos(theta)*(sin(theta))^2
%
%               Finally, there would be both curves and their sum, 
%               the unit circle, are plotted on the last graph.
%
% Notes:        For best results, n should be an odd number less than 10.
%
% Functions:    This program requires the following functions:
% Needed:       remove_spaces.m, title_string_left.m, and title_string_right

clear all, close all
j=sqrt(-1);
theta=0:2*pi/200:2*pi;

n=input('Enter a number for your "Roots Of Unity" plots: ');

% Create the polynomial coefficients for the binomial expansion of (x+1)^n
poly_coeff=abs(poly(eye(n)));
num_terms=length(poly_coeff);       % Number of term in binomial expansion


% Create title string,ts, for proper labeling
% Done with an outside functions
tsl=title_string_left(n,poly_coeff);    % Left side of equations
tsr=title_string_right(n,poly_coeff);   % Right side of equations

% Calculate number of graphs to make and make them
num_graphs=floor(num_terms/2);

k=num_terms;                        % k is the index for cos(theta)
                                    % i is the index for sin(theta)
for i=1:k
    z(i,:)=poly_coeff(i)*(cos(theta)).^(k-1).*(j*sin(theta)).^(i-1);
    k=k-1;
end

% PLOTS
% Outside terms first
% Calculate number of graphs to plot first
num_graphs=floor(num_terms/2);

k=num_terms;                        % k is the index for cos(theta)
                                    % i is the index for sin(theta)
for i=1:num_graphs
    z_curve(i,:)=z(i,:)+z(k,:);
    figure
    hold on
    grid
    axis square
    z_current=z_curve(i,:);
    plot(real(z_current),imag(z_current))
    s=remove_spaces([tsl(i,:) tsr(k,:)]);
    s=['z' num2str(i) ' = ' s ];
    title(s)
    k=k-1;
end

% Deal with odd number of terms in the binomial expansion. 
% This occurs when n is even  
if mod(num_terms,2) == 1
    num_graphs=num_graphs+1;        % Correct for final output
    i=i+1;                          % Correct current index
    z_curve=z(i,:);
    figure
    hold on
    grid
    axis square
    z_current=z_curve;
    plot(real(z_current),imag(z_current))
    s=remove_spaces([tsl(i,:)]);
    s=['z' num2str(i) ' = ' s];
    title(s)
end

% Plot the first set of curves and their sum (unit circle) on one graph
figure
hold on
grid
axis square

% Plot the indiviual curves of unity in blue
k=num_terms;
for i=1:num_graphs   
    z_current=z(i,:)+z(k,:);
    plot(real(z_current),imag(z_current),'b')
    k=k-1;
end

% Plot the sum of all the graphs (unit cirlce) in red
circle=sum(z);
plot(circle,'r')
title('The curves of unity and their sum, the unit circle')

% Makes the final sum look like a unit circle when n is an even number
if mod(num_terms,2)==1
    axis equal
end

% Plot miscellaneous curves of unity
for i=1:num_terms
    for k=i:num_terms
        if (k~=i) & ((num_terms-k+1)~=i)
            z_current=z(i,:)+z(k,:);
            % Eliminate graphs of straight lines
            if (sum(real(z_current))~=0) & (sum(imag(z_current))~=0)
                figure
                hold on
                grid
                axis square
                plot(real(z_current),imag(z_current))
                s=remove_spaces([tsl(i,:) tsr(k,:)]);
                title(s)
            end
        end
    end
end

% Plot the sums of 2 different curves of unity
for i=1:(num_graphs-1)
    for k=(i+1):num_graphs
        z_current=z_curve(i,:)+z_curve(k,:);
        figure
        hold on
        grid
        axis square
        plot(real(z_current),imag(z_current))
        s=['z' num2str(i) '+' 'z' num2str(k)];
        title(s)
    end
end