function [] = Rose_Curve_animation_with_Gears()
% The shape of the gears is not precise, it creates a decent GIF and a SVG.
% The image is rendered horizontally, but saved rotated vertical. 
% This allows the rendering of bigger base images, 
% and therefor better (self made) antialiasing. 
% The plot option "GraphicsSmoothing" made this partly redundant.
%
% After rendering the first frame gets a display time of 0.4 sec and 
% the last frame of 2 sec (using "Jasc Animation Shop").
%   
% 2017-07-12 Jahobr (update 2019-04-14 Jahobr)

RGB.black      = [0   0   0  ];
RGB.white      = [1   1   1  ];
RGB.red        = [1   0   0  ];
RGB.blue       = [0   0   1  ];
RGB.darkBlue   = RGB.blue/1.5;   
RGB.brightGrey = [0.8 0.8 0.8]; % gear
RGB.grey       = [0.5 0.5 0.5]; % gear
RGB.darkGrey   = [0.3 0.3 0.3]; % coordinate system

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


versionList = [... % [n1 d1;  n2 d2;  ...]
    1 1; ...
    1 2; ...
    1 3; ...
    1 4; ...
    1 5; ...
    1 6; ...
    1 7; ...
    1 8; ...
    1 9; ...
    1 10;...    

    2 1; ...
%   2 2; ... % see [1 1];
    2 3; ...
%   2 4; ... % see [1 2];
    2 5; ...
%   2 6; ... % see [1 3];
    2 7; ...
%   2 8; ... % see [1 4];
    2 9; ...
%   2 10;... % see [1 5];

    3 1; ...
    3 2; ...
%   3 3; ... % see [1 1];
    3 4; ...
    3 5; ...
%   3 6; ... % see [1 2];
    3 7; ...
    3 8; ...
%   3 9; ... % see [1 3];
    3 10;...    

    4 1; ...
%   4 2; ... % see [2 1];
    4 3; ...
%   4 4; ... % see [1 1];
    4 5; ...
%   4 6; ... % see [2 3];
    4 7; ...
%   4 8; ... % see [1 2];
    4 9; ...
%   4 10;... % see [2 5];

    5 1; ...
    5 2; ...
    5 3; ...
    5 4; ...
%   5 5; ... % see [1 1];
    5 6; ...
    5 7; ...
    5 8; ...
    5 9; ...
%   5 10;... % see [1 2];
    
    6 1; ...
%   6 2; ... % see [3 1];
%   6 3; ... % see [2 1];
%   6 4; ... % see [3 2];
    6 5; ...
%   6 6; ... % see [1 1];
    6 7; ...
%   6 8; ... % see [3 4];
%   6 9; ... % see [2 3];
%   6 10;... % see [3 5];
    ];

module = 1; % gear size

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

figHandle = figure(15554461);
clf
set(figHandle,'Units','pixel');
set(figHandle,'MenuBar','none',  'ToolBar','none'); % free real estate for a maximally large image
set(figHandle,'Color',RGB.white); % white background
axesHandle = axes;
hold(axesHandle,'on')
axis equal
axis off % invisible axes (no ticks)
drawnow;

for versionNr = 1:size(versionList,1)
    curVers = ['n' num2str(versionList(versionNr,1)) '_d'  num2str(versionList(versionNr,2))]; % 'n1_d1'
    
    n = versionList(versionNr,1);
    d = versionList(versionNr,2);
    
    k = n/d;
    teeth = [60 NaN]; % use "highly composite number"
    teeth(2) = teeth(1)/k;

    if teeth(2)~= round(teeth(2))
        error(curVers)
    end
    
    diameter = teeth.*module;
    
    % % continuous rotation
    % nFrames = 180*max(d,n); % 180 frames per rotation
    % angles_L_wheel = linspace(0, d *2*pi,nFrames); % angles for wheel responsible for Up-Down
    % angles_L_wheel = angles_L_wheel(1:end-1); % remove last frame, it would be double
    
    % stop acc rot deacc stop
    accFrames = 10; % frames for acceleration (first frame will be 0 last at full speed, so practicall it is accFrames-2)
    nFrames = 180*max(d,n)+2*accFrames; % 180 frames per rotation
    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_L_wheel = cumsum(speed)/sum(speed) *d*2*pi; % create position, normalize, scale
    angles_L_wheel(end) = 0; % repalce "d*2*pi" with 0 to avoid rounding errors resulting in minimally crooked lines
    
    center_L_wheel = [0 0]; % wheel for left movement
    center_R_wheel = [mean(diameter) 0]; % wheel for right movement
        
    curveAmplitude = 0.43*diameter(1);
    
    if strcmp(curVers,'n1_d1')
        xLimits = [-0.55*diameter(1)  0.55*diameter(1)+mean(diameter)];  % ADJUST
    else
        xLimits = [-0.55*diameter(1)  1.12*curveAmplitude+mean(diameter)];  % ADJUST
    end
    yLimits = [-0.55  0.55]*diameter(1); % ADJUST
    
    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)
    imageAspectRatio = xRange/yRange;
    MegaPixelTarget = 100*10^6; % Category:Animated GIF files exceeding the 100 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
    scaleReduction = min(...% repeat as often as possible for nice antialiasing
        floor(screenSize(4)/ySize), floor(screenSize(3)/xSize)); 
    if scaleReduction == 0;   error('"MegaPixelTarget" not possible; use smaller target or bigger monitor'); end % check
    
    figPos = [1 1 xSize*scaleReduction ySize*scaleReduction]; % big start image for antialiasing later [x y width height]
    set(figHandle, 'Position', figPos); 
    if ~all(get(figHandle, 'Position') == figPos);   error('figure Position could not be set');   end % check

    liSc = ySize*scaleReduction/600; % LineWidth scale; LineWidth is absolut, a bigger images needs thicker lines to keep them in proportion
    
    axis equal;
    setXYlim(axesHandle,xLimits,yLimits); % reset limits and drawnow

    angles_R_wheel = -angles_L_wheel*k;
    
    if n==4
        angles_R_wheelToothAllign = angles_R_wheel + pi/teeth(2); % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
    else
        angles_R_wheelToothAllign = angles_R_wheel;   % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
    end
    
    yCurvePoints  = curveAmplitude*cos(k*angles_L_wheel).*cos(angles_L_wheel); % rhodonea curve
    xCurvePoints  = curveAmplitude*cos(k*angles_L_wheel).*sin(angles_L_wheel); % rhodonea curve
    
    reducedRGBimage = uint8(ones(xSize,ySize,3,nFrames)); % allocate
    
    for iFrame = 1:nFrames
        cla(axesHandle)
        
        drawSpurWheel(center_L_wheel,teeth(1),module,RGB.brightGrey,2*liSc,RGB.black, angles_L_wheel(iFrame)+pi/teeth(1)); % upper left cogwheel  (fast)
        drawSpurWheel(center_R_wheel,teeth(2),module,RGB.grey      ,2*liSc,RGB.black, angles_R_wheelToothAllign(iFrame)); % lower right cogwheel  (slow)
        
        %     if curveAmplitude >= diameter(2)/2 % if the crank reaches the white area
        %         xCrank = [-module*2   0        module*2       module                  -module];
        %         yCrank = [  0      -module*2    0     curveAmplitude+module  curveAmplitude+module];
        %         [xCrank,yCrank] = rotateCordiantes(xCrank,yCrank,angles_R_wheel(iFrame));
        %         patch(xCrank+center_R_wheel(1),yCrank+center_R_wheel(2),[0.5 0.5 0.5],'EdgeColor',[0 0 0],'LineWidth',2) % raw the crank trapezoid
        %     end
        
        xCoord = [-1 0 0]*curveAmplitude*1.08+center_L_wheel(1);
        yCoord = [ 0 0 1]*curveAmplitude*1.08+center_L_wheel(1);
        [xCoord,yCoord] = rotateCordiantes(xCoord,yCoord,angles_L_wheel(iFrame));
        plot(xCoord,yCoord,'-','Color',RGB.black,'LineWidth',4*liSc) % coordinate system lines
        
        xArrow = [0 -0.8 0.8]*module+center_L_wheel(1);
        yArrow = [1.17 1.07 1.07]*curveAmplitude+center_L_wheel(2);
        [xArrow,yArrow] = rotateCordiantes(xArrow,yArrow,angles_L_wheel(iFrame));
        patch(xArrow,yArrow,RGB.black);
        [xArrow,yArrow] = rotateCordiantes(xArrow,yArrow,pi/2);
        patch(xArrow,yArrow,RGB.black);
        
        xLine = [-1 1]*(diameter(2)/2-module*1.5);
        yLine = [ 0 0]*(diameter(2)/2-module*1.5);
        [xLine,yLine] = rotateCordiantes(xLine,yLine,angles_R_wheel(iFrame));
        plot(xLine+center_R_wheel(1),yLine+center_R_wheel(2),'-','Color',RGB.darkGrey,'LineWidth',5*liSc) % coordinate system lines
        [xLine,yLine] = rotateCordiantes(xLine,yLine,pi/2);
        plot(xLine+center_R_wheel(1),yLine+center_R_wheel(2),'-','Color',RGB.darkGrey,'LineWidth',5*liSc) % coordinate system lines
        
        x_R_Crank = cos(angles_R_wheel(iFrame)+pi/2)*curveAmplitude+center_R_wheel(1); % Up-Down
        y_R_Crank = sin(angles_R_wheel(iFrame)+pi/2)*curveAmplitude+center_R_wheel(2); % Up-Down
        
        plot([center_L_wheel(1) center_L_wheel(1)],[curveAmplitude -curveAmplitude],':','Color',RGB.darkBlue,'LineWidth',6*liSc) % base line
        
        [x,y] = rotateCordiantes(xCurvePoints,yCurvePoints,angles_L_wheel(iFrame));
        plot(x,y,  '-','Color',RGB.red ,'LineWidth',3*liSc) % rose curve
        plot(x(1:max(iFrame,1)),y(1:max(iFrame,1)),  '-','Color',RGB.red ,'LineWidth',7*liSc) % rose curve
        
        plot([center_R_wheel(1) x_R_Crank],[center_R_wheel(2) y_R_Crank],'.-','Color',RGB.darkBlue,'LineWidth',7*liSc) % base line

        plot([0 x_R_Crank],[y_R_Crank y_R_Crank],':','Color',RGB.blue,'LineWidth',6*liSc) % base line
        
        plot(x_R_Crank,y_R_Crank,'.','Color',RGB.blue ,'MarkerSize',60*liSc,'LineWidth',5*liSc) % color marker crank and curve
        plot(0,y_R_Crank,'.','Color',RGB.red  ,'MarkerSize',60*liSc,'LineWidth',7*liSc) % color marker crank and curve
        plot(0,y_R_Crank,'o','Color',RGB.blue ,'MarkerSize',20*liSc,'LineWidth',5*liSc) % color marker crank and curve
        
        %% save animation
        setXYlim(axesHandle,xLimits,yLimits); % reset limits and drawnow
        %  pause(0.01)
        %  return

        f = getframe(figHandle);
        imtemp = imReduceSize(f.cdata,scaleReduction); % allows subpixel lines
        reducedRGBimage(:,:,:,iFrame) = rot90(imtemp,3); % rotate image

%         if iFrame == nFrames % save 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
    
    startMap = cell2mat(struct2cell(RGB)); % struct2colormap; % list of map colors that are not allowed to be changed
    map = createImMap(reducedRGBimage,32,startMap);  % full colormap
    im = uint8(ones(xSize,ySize,1,nFrames)); % allocate
    
    for iFrame = 1:nFrames
        im(:,:,1,iFrame) = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'nodither'); % rgb to colormap image
    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 100 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)
plot(X,Y,'.','MarkerSize',2*linW,'Color',linC); % extra dots make corners look smoother

function [x,y] = rotateCordiantes(x,y,anglee)
% x:      coordinates vertical
% y:      coordinates horizontal
% 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 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 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/scaleReduction ,'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];