function [] = Archimedean_Spiral_animation_with_Gears()
% Source code for drawing gears that produces archimedean spiral.
% The shape of the gears is not precise, it creates a decent GIF and a SVG.
%
% 2019-02-17 Jahobr

RGB.bkgd       = [1   1   1  ]; % white background
RGB.black      = [0   0   0  ]; %
RGB.red        = [1   0.2 0.2]; %
RGB.blue       = [0.2 0.2 1  ]; %
RGB.brightGrey = [0.8 0.8 0.8]; % gear
RGB.grey       = [0.5 0.5 0.5]; % rack
RGB.darkGrey   = [0.3 0.3 0.3]; % marks
RGB.green      = [0.1 0.7 0.1]; %
RGB.violet     = [0.6 0.2 0.8]; % 
RGB.cyan       = [0   0.8 0.8]; %

RGB = structfun(@(q)round(q*255)/255, RGB, 'UniformOutput',false); % round to values that are nicely uint8 compatible

[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location

teeth = 16; % use "highly composite number"
module = 1/8; % gear size

for versionNr = 1:3

    if versionNr == 1
        curVers = 'positive_tracing';
        nFrames = 355;
        accFrames = 10; % frames for acceleration (first frame will be 0 last at full speed, so practicall it is accFrames-2)
        speed = [linspace(0,1,accFrames) ones(1,nFrames+1 -2*accFrames) linspace(1,0,accFrames)];
        speed = speed(1:end-1); % last speed is 0, this does nothing in cumsum; (compensated by +1 frames in center)
        angles = cumsum(speed)/sum(speed) *4*pi;
        xLimits = [-13.2 13.2];  % ADJUST
        yLimits = [-11.6  8.6];  % ADJUST
    elseif versionNr == 2
        curVers = 'positive_and_negative';
        nFrames = 155;
        accFrames = 10; % frames for acceleration (first frame will be 0 last at full speed, so practicall it is accFrames-2)
        speed = [linspace(0,1,accFrames) ones(1,nFrames+1 -2*accFrames) linspace(1,0,accFrames)];
        speed = speed(1:end-1); % last speed is 0, this does nothing in cumsum; (compensated by +1 frames in center)
        angles = cumsum(speed)/sum(speed) *2*pi;
        xLimits = [-13.2 13.2];  % ADJUST
        yLimits = [-10.1 10.1];  % ADJUST
    elseif versionNr == 3
        curVers = 'in_and_out_tracing';
        nFrames = 355;
        accFrames = 10; % frames for acceleration (first frame will be 0 last at full speed, so practicall it is accFrames-2)
        speed = [linspace(0,1,accFrames) ones(1,nFrames+1 -2*accFrames) linspace(1,0,accFrames)];
        speed = speed(1:end-1); % last speed is 0, this does nothing in cumsum; (compensated by +1 frames in center)
        angles = cumsum(speed)/sum(speed) *4*pi;
        xLimits = [-13.2 13.2];  % ADJUST
        yLimits = [-10.1 10.1];  % ADJUST
    end
    
    xRange = xLimits(2)-xLimits(1);
    yRange = yLimits(2)-yLimits(1);
    
    screenSize = get(groot,'Screensize')-[0 0 5 20]; % [1 1 width height] (minus tolerance for figure borders)
    screenAspectRatio = screenSize(3)/screenSize(4); % width/height
    imageAspectRatio = xRange/yRange;
    MegaPixelTarget = 51*10^6; % Category:Animated GIF files exceeding the 50 MP limit
    pxPerImage = MegaPixelTarget/nFrames; % pixel per gif frame
    ySize = sqrt(pxPerImage/imageAspectRatio); % gif height
    xSize = ySize*imageAspectRatio; % gif width
    xSize = floor(xSize); ySize = floor(ySize); % full pixels
    if imageAspectRatio > screenAspectRatio % width will be the problem
        scaleReduction = floor(screenSize(3)/xSize); % repeat as often as possible
    else % height will be the problem
        scaleReduction = floor(screenSize(4)/ySize); % repeat as often as possible
    end
    
    figHandle = figure(25554461); clf
    set(figHandle,'Units','pixel');
    set(figHandle,'Position',[1 1 xSize*scaleReduction ySize*scaleReduction]); % big start image for antialiasing later [x y width height]
    set(figHandle,'MenuBar','none',  'ToolBar','none'); % free real estate for a maximally large image
    set(figHandle,'Color',RGB.bkgd); % white background
    axesHandle = axes;
    hold(axesHandle,'on')
    axis off % invisible axes (no ticks)
    axis equal; 
    setXYlim(axesHandle,xLimits,yLimits); % set limits and drawnow;
    
    liSc = ySize*scaleReduction/800; % LineWidth scale; LineWidth is absolut, a bigger images needs thicker lines to keep them in proportion

    reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate
    
    for iFrame = 1:nFrames
        cla(axesHandle)

        if versionNr == 1 % 'positive_tracing'
            % a moving dotted line jitters, and looks bad; so it is drawn static an the unused part is covered up
            plot([0              xLimits(2)+5],[0 0],':','Color',RGB.blue,'LineWidth', 6*liSc) % horizontal base line 
            plot([angles(iFrame) xLimits(2)+5],[0 0],'-','Color',RGB.bkgd,'LineWidth',10*liSc) % cover up unseen part
        elseif versionNr == 2 % 'positive_and_negative'
            plot(xLimits+[-5 5],[0 0],':','Color',RGB.blue,'LineWidth',6*liSc) % horizontal base line
        elseif versionNr == 3 % 'in_and_out_tracing'
            % a moving dotted line jitters, and looks bad; so it is drawn static an the unused part is covered up
            plot([ 0              xLimits(2)+5],[0 0],':','Color',RGB.blue,  'LineWidth', 6*liSc) % horizontal base line
            plot([ 0              xLimits(1)-5],[0 0],':','Color',RGB.violet,'LineWidth', 6*liSc) % horizontal base line
            plot([ angles(iFrame) xLimits(2)+5],[0 0],'-','Color',RGB.bkgd,  'LineWidth',10*liSc) % cover up unseen part
            plot([-angles(iFrame) xLimits(1)-5],[0 0],'-','Color',RGB.bkgd,  'LineWidth',10*liSc) % cover up unseen part
        end

        drawSpurWheel([0 0],teeth,module,RGB.brightGrey,2*liSc,RGB.black,-angles(iFrame) ); % cogwheel
        sideOffset = angles(iFrame)*teeth*module/2;
        drawRack([0 teeth*module/2],150,module,RGB.grey,2*liSc,RGB.black,sideOffset+ pi/16,1);
        
        if versionNr == 3
            drawRack([0 -teeth*module/2],150,module,RGB.grey,2*liSc,RGB.black,-sideOffset+ pi/16,0);
        end
        
        for iRotMark = (-4:4)*teeth*module*pi
            plot([1 1]* (sideOffset+iRotMark), [10 17]*module,'-','Color',RGB.darkGrey,'LineWidth',4*liSc) % vertical marker on upper rack
            if versionNr ==3
                plot([1 1]* (-sideOffset+iRotMark), [10 17]*-module,'-','Color',RGB.darkGrey,'LineWidth',4*liSc) % vertical marker on lower rack
            end
        end
        
        xCoord = [4 0 0];
        yCoord = [0 0 4];
        
        rotAngle = -angles(iFrame);
        if iFrame == nFrames
            rotAngle = 0; % fix slightly angled lines
        end
        [xCoord,yCoord] = rotateCordiantes(xCoord,yCoord,rotAngle);
        plot(xCoord,yCoord,'-','Color','k','LineWidth',4*liSc) % coordinate system lines
        
        xArrow = [4.5 4.0  4.0];
        yArrow = [0  -0.15 0.15];
        
        [xArrow,yArrow] = rotateCordiantes(xArrow,yArrow,rotAngle);
        patch(xArrow,yArrow,'k');
        [xArrow,yArrow] = rotateCordiantes(xArrow,yArrow,pi/2);
        patch(xArrow,yArrow,'k');
        
        curvePoints = linspace(0,8*pi,500);
        [x,y] = pol2cart(curvePoints-angles(iFrame),curvePoints);

        if versionNr == 1 % 'positive_tracing'
            
            plot(x,y,  '-','Color',RGB.red ,'LineWidth',3*liSc) % curve
            
            [x,y] = pol2cart(angles(1:iFrame)-angles(iFrame),angles(1:iFrame));
            plot( x,y,  '-','Color',RGB.red ,'LineWidth',7*liSc) % trace curve
            
            plot(sideOffset*[1 1],[0 8*module],'-','Color',RGB.blue,'LineWidth',6*liSc) % verical base line

            plot(sideOffset,0,'.','Color',RGB.red,  'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            plot(sideOffset,0,'o','Color',RGB.blue ,'MarkerSize',20*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            
        elseif versionNr == 2 % 'positive_and_negative'

            [xa,ya] = pol2cart(-curvePoints-angles(iFrame),-curvePoints);
            plot(xa,ya,  '-','Color',RGB.cyan ,'LineWidth',7*liSc) % anti-curve
            plot( x, y,  '-','Color',RGB.red  ,'LineWidth',7*liSc) % curve

            for iRotMark = (-4:4)*teeth*module*pi
                plot(sideOffset*[1 1]+iRotMark,[0 8*module],'-','Color',RGB.blue,'LineWidth',6*liSc) % verical base line
            end
            
            for iRotMark = (-4:4)*teeth*module*pi
                if sideOffset+iRotMark >= 0
                    plot(sideOffset+iRotMark,0,'.','Color',RGB.red ,'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
                else
                    plot(sideOffset+iRotMark,0,'.','Color',RGB.cyan ,'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
                end
                plot(sideOffset+iRotMark,0,'o','Color',RGB.blue ,'MarkerSize',20*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            end
            if iFrame == nFrames
                plot(0,0,'.','Color',RGB.cyan ,'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
                plot(0,0,'o','Color',RGB.blue ,'MarkerSize',20*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            end
            
        elseif versionNr == 3 % 'in_and_out_tracing'
            
            plot(-x,-y,'-','Color',RGB.green,'LineWidth',3*liSc) % anti curve
            plot( x, y,'-','Color',RGB.red,  'LineWidth',3*liSc) % curve

            [x,y] = pol2cart(angles(1:iFrame)-angles(iFrame),angles(1:iFrame));
            plot(-x,-y,'-','Color',RGB.green ,'LineWidth',7*liSc) % trace anti curve
            plot( x, y,'-','Color',RGB.red,   'LineWidth',7*liSc) % trace curve
            
            plot( sideOffset*[1 1], [0 8*module],'-','Color',RGB.blue,  'LineWidth',6*liSc) % verical base line
            plot(-sideOffset*[1 1],-[0 8*module],'-','Color',RGB.violet,'LineWidth',6*liSc) % verical base line
            
            plot(-sideOffset,0,'.','Color',RGB.green,  'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            plot(-sideOffset,0,'o','Color',RGB.violet ,'MarkerSize',20*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            
            plot(sideOffset,0,'.','Color',RGB.red,  'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker on trace point of curve
            plot(sideOffset,0,'o','Color',RGB.blue ,'MarkerSize',20*liSc,'LineWidth',5*liSc) % color marker on trace point of curve

        end
        

        %% save animation
        setXYlim(axesHandle,xLimits,yLimits); % reset limits and drawnow
        
        f = getframe(figHandle);
        reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
        
        if iFrame == nFrames % SVG
            if ~isempty(which('plot2svg'))
                plot2svg(fullfile(pathstr, [fname '_' curVers '.svg']),figHandle) % by Juerg Schwizer
            else
                disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
            end
        end
    end
    
    RGBcurrent = RGB; % copy colormap; to be pruned to the used colors
    if versionNr == 1
        RGBcurrent = rmfield(RGBcurrent,{'green','violet','cyan'}); % remove unused colors
        nColors = 32;
    elseif versionNr == 2
        RGBcurrent = rmfield(RGBcurrent,{'green','violet'}); % remove unused colors
        nColors = 64;
    elseif versionNr == 3
        RGBcurrent = rmfield(RGBcurrent,{'cyan'}); % remove unused colors
        nColors = 64;
    end

    startMap = cell2mat(struct2cell(RGBcurrent)); % struct2colormap; % list of map colors that are not allowed to be changed
    map = createImMap(reducedRGBimage,nColors,startMap);  % full colormap
    
    im = uint8(ones(ySize,xSize,1,nFrames)); % allocate
    for iFrame = 1:nFrames
        im(:,:,1,iFrame) = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'nodither');
    end
    
    imwrite(im,map,fullfile(pathstr, [fname '_' curVers '.gif']),'DelayTime',1/25,'LoopCount',inf) % save gif
    disp([fname '_' curVers '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit
end

function drawSpurWheel(center,toothNumber,module,fillC,linW,linC,startOffset)
% DRAWSPURWHEEL - draw a simple Toothed Wheel
%    center:       [x y]
%    toothNumber:  scalar
%    module:       scalar tooth "size"
%    fillC:        color of filling [r g b]
%    linW:         LineWidth
%    linC:         LineColor
%    startOffset:  start rotation (scalar)[rad]

effectiveRadius = module*toothNumber/2; % effective Radius

outsideRadius =     effectiveRadius+1*  module; %                +---+             +---+
upperRisingRadius = effectiveRadius+0.5*module; %               /     \           /     \
% effective Radius                              %              /       \         /       \
lowerRisingRadius = effectiveRadius-0.5*module; %             I         I       I         I
rootRadius =        effectiveRadius-1.1*module; %     + - - - +         + - - - +         +

angleBetweenTeeth = 2*pi/toothNumber; % angle between 2 teeth
angleOffPoints = (0:angleBetweenTeeth/16:(2*pi));
angleOffPoints = angleOffPoints+startOffset; % apply rotation offset

angleOffPoints(7:16:end) =  angleOffPoints(7:16:end)  + 1/toothNumber^1.2; % hack to create smaller tooth tip
angleOffPoints(11:16:end) = angleOffPoints(11:16:end) - 1/toothNumber^1.2; % hack to create smaller tooth tip

angleOffPoints(8:16:end)  = (angleOffPoints(7:16:end) +  angleOffPoints(9:16:end))/2; % shift the neighbouring tip point in accordingly
angleOffPoints(10:16:end) = (angleOffPoints(11:16:end) + angleOffPoints(9:16:end))/2; % shift the neighbouring tip point in accordingly

angleOffPoints(6:16:end) =  angleOffPoints(6:16:end)  + 1/toothNumber^1.7; % hack to create slender upperRisingRadius
angleOffPoints(12:16:end) = angleOffPoints(12:16:end) - 1/toothNumber^1.7; % hack to create slender upperRisingRadius

radiusOffPoints = angleOffPoints; % allocate with correct site

radiusOffPoints( 1:16:end) = rootRadius;        % center bottom         I
radiusOffPoints( 2:16:end) = rootRadius;        % left bottom           I
radiusOffPoints( 3:16:end) = rootRadius;        % left bottom corner    +
radiusOffPoints( 4:16:end) = lowerRisingRadius; % lower rising bottom      \
radiusOffPoints( 5:16:end) = effectiveRadius;   % rising edge                 \
radiusOffPoints( 6:16:end) = upperRisingRadius; % upper rising edge              \
radiusOffPoints( 7:16:end) = outsideRadius;     % right top  corner                 +
radiusOffPoints( 8:16:end) = outsideRadius;     % right top                         I
radiusOffPoints( 9:16:end) = outsideRadius;     % center top                        I
radiusOffPoints(10:16:end) = outsideRadius;     % left top                          I
radiusOffPoints(11:16:end) = outsideRadius;     % left top  corner                  +
radiusOffPoints(12:16:end) = upperRisingRadius; % upper falling edge             /
radiusOffPoints(13:16:end) = effectiveRadius;   % falling edge                /
radiusOffPoints(14:16:end) = lowerRisingRadius; % lower falling edge       /
radiusOffPoints(15:16:end) = rootRadius;        % right bottom  corner  +
radiusOffPoints(16:16:end) = rootRadius;        % right bottom          I

[X,Y] = pol2cart(angleOffPoints,radiusOffPoints);
X = X+center(1); % center offset
Y = Y+center(2); % center offset
patch(X,Y,fillC,'EdgeColor',linC,'LineWidth',linW)


function drawRack(center,toothNumber,module,fillC,linW,linC,startOffset,top)
%    center:       [x y]
%    toothNumber:  scalar
%    module:       scalar tooth "size"
%    fillC:        color of filling [r g b]
%    linW:         LineWidth
%    linC:         LineColor
%    startOffset:  initial shift
%    top:          1=top 0=bottom
x = (0:toothNumber*4-1)*pi*module/4;
x = x-mean(x)+center(1)+startOffset;

y = zeros(size(x));

y(1:4:end) = y(1:4:end)+1.1*module; %    +###I bottom
y(2:4:end) = y(2:4:end)-1  *module; % +######I tip
y(3:4:end) = y(3:4:end)-1  *module; % +######I tip
y(4:4:end) = y(4:4:end)+1.1*module; %    +###I bottom

x(1:4:end) = x(1:4:end)-0.14*module; % bottom smaller
x(2:4:end) = x(2:4:end)+0.14*module; % tip smaller
x(3:4:end) = x(3:4:end)-0.14*module; % tip smaller
x(4:4:end) = x(4:4:end)+0.14*module; % bottom smaller

x = [x(1) x x(end)];
y = [y(1)+9*module y y(end)+9*module];

if ~top
    y = -y; % flip
end
y = y+center(2);
patch(x,y,fillC,'EdgeColor',linC,'LineWidth',linW);


function setXYlim(axesHandle,xLimits,yLimits)
% set limits; practically the axis overhangs the figure all around, to
% hide rendering error at line-ends.
% Input:
%   axesHandle:        
%   xLimits, yLimits:  [min max]
overh = 0.05; % 5% overhang all around; 10% bigger in x and y
xlim([+xLimits(1)*(1+overh)-xLimits(2)*overh  -xLimits(1)*overh+xLimits(2)*(1+overh)])
ylim([+yLimits(1)*(1+overh)-yLimits(2)*overh  -yLimits(1)*overh+yLimits(2)*(1+overh)])
set(axesHandle,'Position',[-overh -overh  1+2*overh 1+2*overh]); % stretch axis as bigger as figure, [x y width height]
drawnow;


function [x,y] = rotateCordiantes(x,y,anglee)
% x coordinates of the center
% y coordinates of the center
% anglee angle of rotation in [rad]
rotM = [cos(anglee) -sin(anglee); sin(anglee) cos(anglee)];
x_y = rotM*[x(:)';y(:)'];
x = x_y(1,:);
y = x_y(2,:);


function im = imReduceSize(im,redSize)
% Input:
%  im:      image, [imRows x imColumns x nChannel x nStack] (unit8)
%                      imRows, imColumns: must be divisible by redSize
%                      nChannel: usually 3 (RGB) or 1 (grey)
%                      nStack:   number of stacked images
%                                usually 1; >1 for animations
%  redSize: 2 = half the size (quarter of pixels)
%           3 = third the size (ninth of pixels)
%           ... and so on
% Output:
%  im:     [imRows/redSize x imColumns/redSize x nChannel x nStack] (unit8)
%
% an alternative is: imNew = imresize(im,1/reduceImage,'bilinear');
%        BUT 'bicubic' & 'bilinear'  produces fuzzy lines
%        IMHO this function produces nicer results as "imresize"
 
[nRow,nCol,nChannel,nStack] = size(im);

if redSize==1;  return;  end % nothing to do
if redSize~=round(abs(redSize));             error('"redSize" must be a positive integer');  end
if rem(nRow,redSize)~=0;     error('number of pixel-rows must be a multiple of "redSize"');  end
if rem(nCol,redSize)~=0;  error('number of pixel-columns must be a multiple of "redSize"');  end

nRowNew = nRow/redSize;
nColNew = nCol/redSize;

im = double(im).^2; % brightness rescaling from "linear to the human eye" to the "physics domain"; see youtube: /watch?v=LKnqECcg6Gw
im = reshape(im, nRow, redSize, nColNew*nChannel*nStack); % packets of width redSize, as columns next to each other
im = sum(im,2); % sum in all rows. Size of result: [nRow, 1, nColNew*nChannel]
im = permute(im, [3,1,2,4]); % move singleton-dimension-2 to dimension-3; transpose image. Size of result: [nColNew*nChannel, nRow, 1]
im = reshape(im, nColNew*nChannel*nStack, redSize, nRowNew); % packets of width redSize, as columns next to each other
im = sum(im,2); % sum in all rows. Size of result: [nColNew*nChannel, 1, nRowNew]
im = permute(im, [3,1,2,4]); % move singleton-dimension-2 to dimension-3; transpose image back. Size of result: [nRowNew, nColNew*nChannel, 1]
im = reshape(im, nRowNew, nColNew, nChannel, nStack); % putting all channels (rgb) back behind each other in the third dimension
im = uint8(sqrt(im./redSize^2)); % mean; re-normalize brightness: "scale linear to the human eye"; back in uint8


function map = createImMap(imRGB,nCol,startMap)
% createImMap creates a color-map including predefined colors.
% "rgb2ind" creates a map but there is no option to predefine some colors,
%         and it does not handle stacked images.
% Input:
%   imRGB:     image, [imRows x imColumns x 3(RGB) x nStack] (unit8)
%   nCol:      total number of colors the map should have, [integer]
%   startMap:  predefined colors; colormap format, [p x 3] (double)

imRGB = permute(imRGB,[1 2 4 3]); % step1; make unified column-image (handling possible nStack)
imRGBcolumn = reshape(imRGB,[],1,3,1); % step2; make unified column-image

fullMap = double(permute(imRGBcolumn,[1 3 2]))./255; % "column image" to color map 
[fullMap,~,imMapColumn] = unique(fullMap,'rows'); % find all unique colors; create indexed colormap-image
% "cmunique" could be used but is buggy and inconvenient because the output changes between "uint8" and "double"

nColFul = size(fullMap,1);
nColStart = size(startMap,1);
disp(['Number of colors: ' num2str(nColFul) ' (including ' num2str(nColStart) ' self defined)']);

if nCol<=nColStart;  error('Not enough colors');        end
if nCol>nColFul;   warning('More colors than needed');  end

isPreDefCol = false(size(imMapColumn)); % init
 
for iCol = 1:nColStart
    diff = sum(abs(fullMap-repmat(startMap(iCol,:),nColFul,1)),2); % difference between a predefined and all colors
    [mDiff,index] = min(diff); % find matching (or most similar) color
    if mDiff>0.05 % color handling is not precise
        warning(['Predefined color ' num2str(iCol) ' does not appear in image'])
        continue
    end
    isThisPreDefCol = imMapColumn==index; % find all pixel with predefined color
    disp([num2str(sum(isThisPreDefCol(:))) ' pixel have predefined color ' num2str(iCol)]);
    isPreDefCol = or(isPreDefCol,isThisPreDefCol); % combine with overall list
end
[~,mapAdditional] = rgb2ind(imRGBcolumn(~isPreDefCol,:,:),nCol-nColStart,'nodither'); % create map of remaining colors
map = [startMap;mapAdditional];