function Cardan_Gear_linear_movement()
% Source code for drawing Cardan Gears that produces linear movement.
% The shape of the gears is not precise, it creates a decent GIF and a SVG.
%
% 2018-05-28 Jahobr

[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location
 
modul = 16;
 
carrierCol   = round([0.1 0.7 0.1].*255)./255; % green
sunCol       = round([0.95 0.65 0  ].*255)./255; % yellow (obviously)
palnetColOut = round([1 0.2 0.2].*255)./255; % red
palnetColIn  = round([0.2 0.2 1  ].*255)./255; % blue   (obviously)
lineCol      = round([0.6 0.2 0.8].*255)./255; % violet
 
for currentCase = 1:2;

    figHandle = figure(17674755); clf
    set(figHandle, 'Units','pixel');
    set(figHandle, 'Color',[1 1 1]);
    set(figHandle, 'GraphicsSmoothing','on') % requires at least version 2014b

    nFrames = 452;

    switch currentCase
        case 1 % Stationary_Sun
            saveName = [fname '_packed'];
            teethSun  = 32; % must be divisible by 2
            teethPlanIn = 12; % Idler, tooth number does not matter
            scaleReduction = 3; % reduction for nice antialiasing
            ySize = 270; % pixel
        case 2 % Stationary_Sun
            saveName = [fname '_linedUp'];
            teethSun  = 32; % must be divisible by 2
            teethPlanIn = 20; % Idler, tooth number does not matter
            scaleReduction = 3; % reduction for nice antialiasing
            ySize = 250; % pixel
    end

    teethPlanOut = teethSun/2; % necessary condition

    radiusSun  = modul.*teethSun./2; % effective radius
    radiusPlanIn = modul.*teethPlanIn./2; % effective radius
    radiusPlanOut = modul.*teethPlanOut./2; % effective radius

    distInOut = radiusPlanIn+radiusPlanOut; % distance of centers
    distSunIn = radiusSun+radiusPlanIn; % distance of centers

    angleCarrier = -linspace(0,pi*2,nFrames+1); % define gear position in frames
    angleCarrier = angleCarrier(1:end-1); % remove last frame, it would be double
    angleCarrier = circshift(angleCarrier,[0,round(nFrames/9)]);

    anglePlanOut = -angleCarrier; % gear ratio
    anglePlanIn = angleCarrier.*( teethSun/teethPlanIn+1 ); % gear ratio
    angleSun = zeros(size(angleCarrier));

    switch currentCase
        case 1 % Stationary_Sun
            distSunOut = radiusSun+radiusPlanOut+modul*2.5;
            anglePlanOut = anglePlanOut + (pi/teethPlanOut); % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            anglePlanIn = anglePlanIn + (pi/teethPlanIn)*1.3; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            angleSun = angleSun+(pi/teethSun)*1.5 ; % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT

        case 2 % Stationary_Sun
            distSunOut = radiusSun+radiusPlanIn*2+radiusPlanOut;
            anglePlanOut = anglePlanOut + (pi/teethPlanOut); % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            % anglePlanIn = anglePlanIn + (pi/teethPlanIn); % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
            angleSun = angleSun+(pi/teethSun); % ALLIGNMENT; THIS MAY NEED MANUAL ADJUSTMENT
    end
    angPlanInOffset = acos((distInOut^2-distSunIn^2-distSunOut^2) / (-2*distSunIn*distSunOut)); % triangle calculation

    xlimValues = [-1.07 1.07]*(distSunOut*2);
    ylimValues = [-1.07 1.07]*(distSunOut+radiusPlanOut);

    xRange = xlimValues(2)-xlimValues(1);
    yRange = ylimValues(2)-ylimValues(1);

    xSize = round(ySize/10/yRange*xRange)*10; % pixel
 
    set(figHandle, 'position',[1 1 xSize*scaleReduction ySize*scaleReduction]); % big start image for antialiasing later [x y width height]
    axesHandle = axes; hold(axesHandle,'on');
    set(axesHandle,'position',[0 0 1 1]); % stretch axis as big as figure, [x y width height]
    axis equal
    axis off % invisible axes (no ticks)
    xlim(xlimValues); ylim(ylimValues);
    drawnow;

    reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate
 
    x_y_Crank = [-modul*1.5 -modul*1.5 distSunOut+modul distSunOut+modul;... % x crank outline
                 -modul*1.5  modul*1.5            modul           -modul];   % y

    for iFrame = 1:nFrames

        cla(axesHandle) % fresh frame
 
        %% sun
        drawCogWheel(axesHandle,[0 0],teethSun ,modul,sunCol,angleSun(iFrame));

        %% planets
        [Xin,Yin] = pol2cart(angPlanInOffset+angleCarrier(iFrame) ,distSunIn);
        drawCogWheel(axesHandle,[Xin,Yin],teethPlanIn,modul,palnetColIn,anglePlanIn(iFrame)); % inner planetary gear (Idler)

        [Xout,Yout] = pol2cart(angleCarrier(iFrame) ,distSunOut);
        drawCogWheel(axesHandle,[Xout,Yout],teethPlanOut,modul,palnetColOut,anglePlanOut(iFrame)); % outer planetary gear

        %% carrier
        carrierPoints =0:0.1:2*pi-angleCarrier(iFrame); % circle points
        r = min(radiusPlanIn,radiusPlanOut)*0.5; % radius of carrier frame
        xSunCar =  cos(carrierPoints)*r;         ySunCar =  sin(carrierPoints)*r; % points on left circle of carrier; (start on the right; counterClock)
        xOutCar = -cos(carrierPoints)*r+Xout;    yOutCar = -sin(carrierPoints)*r+Yout; % points on right circle of carrier; (start on the left; counterClock)
        xInCar  =  sin(carrierPoints)*r+Xin;     yInCar  = -cos(carrierPoints)*r+Yin; % points on top circle of carrier; (start on the bottom; counterClock)
        xCarrier = [xSunCar xOutCar xInCar];     yCarrier = [ySunCar yOutCar yInCar]; % assemble carrier
        valPoint = convhull(xCarrier, yCarrier); % Convex hull
        xCarrier = xCarrier(valPoint);           yCarrier = yCarrier(valPoint); % remove concave carrier points
        plot(xCarrier,yCarrier,'color',[0 0 0],'LineWidth',13) % black outline
        plot(xCarrier,yCarrier,'color',carrierCol,'LineWidth',10) % carrierCol "filling"

        circlePatch(0,0,r*1.3,carrierCol,2); % sun carrier
        circlePatch(0,0,r*0.8,sunCol,2);     % sun axle

        circlePatch(Xin,Yin,r*1.3,carrierCol,2); % inner planet carrier
        circlePatch(Xin,Yin,r*0.8,palnetColIn,2);    % inner planet axle

        circlePatch(Xout,Yout,r*1.3,carrierCol,2); % outer planet carrier

        rotM = [cos(-angleCarrier(iFrame)) -sin(-angleCarrier(iFrame)); -sin(+angleCarrier(iFrame)) cos(-angleCarrier(iFrame))]; % rotation matrix
        vecTemp = rotM*x_y_Crank; % rotate crank
        patch(vecTemp(1,:)+Xout,vecTemp(2,:)+Yout,palnetColOut,'EdgeColor',[0 0 0],'LineWidth',2) % raw the crank trapezoid

        circlePatch(Xout,Yout,r*0.8,palnetColOut,2); % outer planet axle

        plot([-2 2]*distSunOut,[0 0],'Linewidth',8,'Color',lineCol); % plot the line
        circlePatch(cos(angleCarrier(iFrame))*2*distSunOut,0,modul*2,palnetColOut,2); % outer planet carrier
        circlePatch(cos(angleCarrier(iFrame))*2*distSunOut,0,modul*1.2,lineCol,2); % outer planet carrier

        %% save animation
        f = getframe(figHandle);
        reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing

        if iFrame == 1 % SVG
            if ~isempty(which('plot2svg'))
                plot2svg(fullfile(pathstr, [saveName '_horizontal.svg']),figHandle) % by Juerg Schwizer
            else
                disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
            end
        end
    end

    map = createImMap(reducedRGBimage,16,[0 0 0;1 1 1;carrierCol;sunCol;palnetColOut;palnetColIn;lineCol]); % colormap

    im = uint8(ones(ySize,xSize,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, [saveName '_horizontal.gif']),'DelayTime',1/25,'LoopCount',inf) % save gif
    disp([saveName '_horizontal.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit
    imwrite(fliplr(permute(im,[2 1 3 4])),map,fullfile(pathstr, [saveName '_vertical.gif']),'DelayTime',1/25,'LoopCount',inf) % save gif
    disp([saveName '_vertical.gif  has '   num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit
end

function drawCogWheel(axesHandle,center,toothNumber,modul,colFilling,startOffset)
% DRAWTOOTHEDWHEEL - draw a simple Toothed Wheel
%
%  Input:
%    axesHandle:
%    center:       [x y]
%    toothNumber:  scalar
%    modul:        scalar tooth "size"
%    colFilling:   color of filling [r g b]
%    startOffset:  start rotation (scalar)[rad]
 
effectiveRadius = modul*toothNumber/2; % effective effectiveRadius
 
outsideRadius =     effectiveRadius+1*  modul; %                +---+             +---+
upperRisingRadius = effectiveRadius+0.5*modul; %               /     \           /     \
% effective Radius                             %              /       \         /       \
lowerRisingRadius = effectiveRadius-0.5*modul; %             I I       I I
rootRadius =        effectiveRadius-1.1*modul; %     + - - - +         + - - - +         +
 
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 tooth
angleOffPoints(12:16:end) = angleOffPoints(12:16:end) - 1/toothNumber^1.7; % hack to create slender tooth
 
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,colFilling,'EdgeColor',[0 0 0],'LineWidth',2)
% plot(axesHandle,X,Y,'-x','linewidth',2,'color',[0 0 0]);
 
% %% effective Radius
% [X,Y] = pol2cart(angleOffPoints,effectiveRadius);
% X = X+center(1); % center offset
% Y = Y+center(2); % center offset
% plot(axesHandle,X,Y,'-.','color',[0 0 0]);

function circlePatch(x,y,r,col,linW)
% x coordinates of the center
% y coordinates of the center
% r is the radius of the circle
% col patch color
% linW LineWidth
angleOffPoints = linspace(0,2.001*pi,200);
xc = x + r*cos(angleOffPoints);
yc = y + r*sin(angleOffPoints);
patch(xc,yc,col,'EdgeColor',[0 0 0],'LineWidth',linW);

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:
%  imNew:  unit8([imRows/redSize x imColumns/redSize x nChannel x nStack])
%
% 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 colores; 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 colores
    [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];