function spincircles(options)
% Draws hierarchical spinning circles.
% 
% create: 6/4/2014
% modified: 6/30/2014
% copyright @2014 author: Dongxi Zheng
%
% description:
%    This function draws a central symetric pattern. Such pattern features
%    a circle in the center, spinning a number of smaller circles (at the
%    next level) at equal angle intervals, and each of the smaller circles
%    spinning even smaller circles (at an even lower level) in a similar
%    manner and so on so forth.
%
%    inputs:
%       options - a structure containing specifications of the pattern. (optional)
%                 .level - the number of levels of circles.
%                 .numb_branches - the number of circles span by each
%                                  non-leaf circle.
%                 .rt_r - the radius of the root circle.
%                 .r_ratio - the ratio of the radius of level i+1 to that
%                            of level i.
%                 .closeness - the closeness between the circles.
%                 .spinsteps - the number of steps of spinning
%                 .spininterval - the time between each two spins, for gif
%                                 animation
%                 .spinspeeds - a vector containing the numbers of degrees
%                               for each level of circles to spin in each
%                               spin step. The ith element corresponds to
%                               the ith level. If the specified number of
%                               levels is larger than the length of this
%                               vector, the missing will be set 0
%                 .output - the name (or path) of the output gif file
narginchk(0,1);

% set default options
DEFAULTOPTIONS = struct('level', 3, 'numb_branches', 6, ...
                       'rt_r', 8, 'r_ratio', 0.6, 'closeness', 1.2, ...
                       'spinsteps', 60, 'spininterval',0.1, ...
                       'spinspeeds', [1,1,1,1,1], ...
                       'output', 'spincircles.gif');

% read user input option values; if not given, apply the default.
if nargin < 1
    options = DEFAULTOPTIONS;
else
    if isempty(options)
        options = DEFAULTOPTIONS;
    else
        pfields = fieldnames(DEFAULTOPTIONS);
        nFields = size(pfields,1);
        for i=1:nFields
            fldname = pfields{i};
            % if the input params has the needed field
            if isfield(options, fldname)
                % if the field is empty
                if isempty(options.(fldname))
                    options.(fldname) = DEFAULTOPTIONS.(fldname);
                end
            % if the input params does not have the needed field
            else
                options.(fldname) = DEFAULTOPTIONS.(fldname);
            end
        end
    end
end

% Prepare the spinning configuration
spinsteps = options.spinsteps;
spinspeeds = options.spinspeeds;
if length(spinspeeds) < options.level
    tmp = zeros(1, options.level);
    tmp(1,1:length(spinspeeds)) = spinspeeds;
    spinspeeds = tmp;
end

% Calculate the center-to-center distance between each two level of circles.
dists = calculateBranchDistances(options.level, options.rt_r, ...
                                 options.r_ratio, options.closeness);
% Calculate the radius of each level of circles.
radii = options.rt_r*(options.r_ratio.^(0:options.level-1));

% Set the output image size
window_limits = [-1, 1]*(sum(dists)+radii(end));
h = figure; set(h, 'Visible','off');

for i = 1:spinsteps
    % Calculate the rotation
    angles_rotated = spinspeeds*(i-1)/180*pi;
    
    % Start drawing
    clf; % clear the previous drawings
    axis equal
    xlim(window_limits);
    ylim(window_limits);
    hold on
    circle(0, 0, options.rt_r); % Draw the root circle.
    % Recursively draw the higher levels of circles.
    branchout(0, 0, 0, 1, ...
              radii, dists, angles_rotated, ...
              options.level-1, options.numb_branches);

    % Finish drawing
    hold off
    
    % Added the current drawing to the output gif
    im = frame2im(getframe(h));
    [A, map] = rgb2ind(im, 256);
    if i == 1
        imwrite(A, map, options.output, 'gif', ...
                'LoopCount', Inf, 'DelayTime', options.spininterval);
    else
        imwrite(A, map, options.output, 'gif', ...
                'WriteMode', 'append', 'DelayTime', options.spininterval);
    end
end

end

function branchout(cx_cur, cy_cur, angle_cur, level,...
                   radii, dists,  angles_rotated, ...
                   maxlevel, numb_branches)
% Draws the recursive higher levels of circles and spinning lines of the
% current circle.
%
%    inputs:
%       cx_cur - the center x coordinate of the current circle
%       cy_cur - the center y coordinate of the current circle
%       angle_cur - the orientation of the current circle
%       level - the level of the current circle
%       radii - an array containing the radius value of each level of circles
%       dists - an array containing the center-to-center distances of each
%               two successive levels of circles
%       angles_rotated - a vector indicating for each level of circles how
%                        much they have be spun since the beginning of
%                        time
%       maxlevel - the maximum level of circles that have branches
%       numb_branches - the number of next level circles

    % Get the radius and branch out distance for the current level and the
    % radius of the next level
    r_cur = radii(level);
    r_next = radii(level+1);
    dist_cur = dists(level);
    
    % Calculate the baseline angles of the branch circles
    angle_next = (angle_cur + angles_rotated(level) + pi/2) - ...
                    ((0:numb_branches-1)*2*pi/numb_branches);
    
    % Calculate the centers of the branch circles
    cx_next = cx_cur + dist_cur*cos(angle_next);
    cy_next = cy_cur + dist_cur*sin(angle_next);
    
    % Draw branch circles and swing lines
    angle_delt = acos(r_cur/dist_cur);
    for i = 1:numb_branches
        circle(cx_next(i), cy_next(i), r_next);
        angle_off = angle_next(i) + angle_delt;
        end1x = cx_cur + r_cur*cos(angle_off);
        end1y = cy_cur + r_cur*sin(angle_off);
        plot([end1x, cx_next(i)],[end1y, cy_next(i)]);
        
        % Recursively draw next level of branches if the maximum level is
        % not reached yet.
        if (level<maxlevel)
            branchout(cx_next(i), cy_next(i), angle_next(i), level+1, ...
                      radii, dists, angles_rotated, maxlevel, numb_branches);
        end
    end
end

function circle(cx, cy, radius)
% Draws a circle.
%
%    inputs:
%       cx - the x coordinate of the center of the circle
%       cy - the y coordinate of the center of the circle
%       radius - the radius of the circle

    ang = 0:0.01:2*pi;
    xp = radius*cos(ang);
    yp = radius*sin(ang);
    plot(cx+xp, cy+yp);
end

function branchDists = ...
        calculateBranchDistances(level, rt_r, r_ratio, closeness)
% Calculates the center-to-center distances between each successive pair of
% circle levels.
%
%    inputs:
%       level - the total number of levels
%       rt_r - the radius of the root circle
%       r_ratio - the ratio of the radius of level i+1 to that of level i
%       closeness - the closeness between the circles

    radii = rt_r*(r_ratio.^(0:level-1));
    branchDists = zeros(1,level-1);
    branchDists(level-1) = closeness*(sum(radii(level-1:level)));
    % The calculation is bottom-up to ensure that the circles don't
    % intersect when the closeness is larger than 1
    for i = level-2:-1:1
        branchDists(i) = (radii(i)+sum(branchDists)+radii(level))*closeness;
    end
end