function [] = Fowler_Ellipsograph()
% The shape of the gears is not precise, it creates a decent GIF and a SVG.
%
% 2017-05-30 Jahobr

modul=1;

col.green     = round([0.1 0.7 0.1 ]*255)/255; % pin
col.darkGreen = round([0   0.3 0   ]*255)/255; % frame
col.darkGrey  = round([0.5 0.5 0.65]*255)/255; % wheels
col.grey      = round([0.8 0.8 0.8 ]*255)/255; % wheels
col.grey2     = round([0.7 0.7 0.7 ]*255)/255; % slider
col.E         = round([1   0   0   ]*255)/255; % red
teeth = [60 60];
diameter = modul.*teeth;

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

nFrames = 330; % frames per rotation
scaleReduction = 3; % reduction for nice antialiasing
xSize = 600; % pixel
ySize = 260; % pixel
figHandle = figure(15554461); clf

set(figHandle, 'Units','pixel');
set(figHandle, 'position',[1 1 [xSize ySize]*scaleReduction]); % big start image for antialiasing later [x y width height]
set(figHandle, 'GraphicsSmoothing','on') % requires at least version 2014b
axesHandle = axes;
hold(axesHandle,'on')
xlimValues = [-1.05  2.05]*diameter(1); % ADJUST
ylimValues = [-0.7   0.7] *diameter(1); % ADJUST
set(axesHandle,'position',[-0.05 -0.05 1.1 1.1]); % stretch axis bigger as figure, easy way to get rid of ticks [x y width height]
xlim(xlimValues); ylim(ylimValues); % set axis limits
axis equal; drawnow;

angles_wheels = linspace(0.2*pi, 2.2*pi,nFrames+1); % angles for wheel 
angles_wheels = angles_wheels(1:end-1); % remove last frame, it would be double

leftCrank  = diameter(1)*0.2;
rightCrank = diameter(1)*0.35;
reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate

for iFrame = 1:nFrames
    cla(axesHandle)
    curAngle = angles_wheels(iFrame);
    
    plotBox(0.5*diameter(1), 0,1.4*diameter(1),0.578*diameter(1),0,[1 1 1],col.darkGreen,15) % outer frame

    plotBox(0.5*diameter(1)+leftCrank*cos(-curAngle), 0, 1.05*diameter(1),0.556*diameter(1),0,[0.8 0.8 1],[0 0 0.5],5) % sliging blue box
    
    center_L_wheel = [leftCrank*cos(-curAngle) 0]; % Left wheel
    center_R_wheel = [leftCrank*cos(-curAngle)+mean(diameter)  0]; % Right wheel

    drawCogWheel(axesHandle,center_L_wheel,teeth(1),modul,col.darkGrey,-curAngle); %  left cogwheel
    drawCogWheel(axesHandle,center_R_wheel,teeth(2),modul,col.darkGrey, curAngle-(pi/teeth(2))); % right cogwheel
    
    circlePatch(center_L_wheel(1),center_L_wheel(2),diameter(1)*0.47,col.grey,'k',3) %  grey circle
    circlePatch(center_R_wheel(1),center_R_wheel(2),diameter(2)*0.47,col.grey,'k',3) %  grey circle
    
    plotBox(center_L_wheel(1),center_L_wheel(2),eps,diameter(1)*0.1, curAngle,[1 1 1],'k',3) % missuse box to plot a line
    plotBox(center_R_wheel(1),center_R_wheel(2),eps,diameter(1)*0.1,-curAngle,[1 1 1],'k',3) % missuse box to plot a line
    
    plotBox(center_L_wheel(1),center_L_wheel(2),diameter(1)*0.45,diameter(1)*0.05, curAngle,[1 1 1],'k',3) % plot slit
    plotBox(center_R_wheel(1),center_R_wheel(2),diameter(1)*0.45,diameter(1)*0.05,-curAngle,[1 1 1],'k',3) % plot slit
    
    plotBox(0,leftCrank*sin(curAngle),diameter(1)*0.08,diameter(1)*0.048, curAngle,col.grey2,'k',3) % left crank slider
    circlePatch(0,leftCrank*sin(curAngle),diameter(2)*0.04,col.green,'k',3) % green crank pin

    plotBox(center_R_wheel(1)+cos(curAngle)*rightCrank,center_R_wheel(2)-sin(-curAngle)*rightCrank,diameter(1)*0.08,diameter(1)*0.048,-curAngle,col.grey2,'k',3)  % right pen slider
  
    plot( [1 1]*0.06*diameter(1),[1 -1]*0.57*diameter(1),'-','color',col.darkGreen,'LineWidth',15) % vertical guide
    plot(-[1 1]*0.06*diameter(1),[1 -1]*0.57*diameter(1),'-','color',col.darkGreen,'LineWidth',15) % vertical guide

    ellipse(diameter(1),0,rightCrank+leftCrank,rightCrank,col.E,7)
    plot(center_R_wheel(1)+cos(curAngle)*rightCrank,center_R_wheel(2)-sin(-curAngle)*rightCrank,'.','color',col.E,'MarkerSize',60) % pen marker

   %% save animation
    xlim(xlimValues); ylim(ylimValues); drawnow; % set axis limits
        
    if iFrame == 1 % save SVG
        if ~isempty(which('plot2svg'))
            plot2svg(fullfile(pathstr, [fname '_a.svg']),figHandle) % by Juerg Schwizer, See http://www.zhinst.com/blogs/schwizer/
        else
            disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
        end
    end
    f = getframe(figHandle);
    reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
end
    
map = createImMap(reducedRGBimage,32,[1 1 1;struct2map(col)]); % 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 '_a.gif']),'DelayTime',1/25,'LoopCount',inf) % save gif
disp([fname '_a.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit

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 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,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,colFa,colEd,linw)
% x  coordinates of the center
% y  coordinates of the center
% r  is the radius of the circle
% colFa  face color  [r g b]
% colEd  edge color  [r g b]
% linw  line width
angleOffPoints = linspace(0,2*pi,300);
xc = x + r*cos(angleOffPoints);
yc = y + r*sin(angleOffPoints);
patch(xc,yc,colFa,'linewidth',linw,'EdgeColor',colEd); %


function plotBox(x,y,wi,hi,rot,colFa,colEd,linw)
% x  coordinates of the center
% y  coordinates of the center
% wi  half of width
% hi  half of height
% rot  dotation in [rad]
% colFa  face color  [r g b]
% colEd  edge color  [r g b]
% linw  line width

% joint in the middle of an edge to get nice corners
xs = [-wi -wi  wi  wi -wi -wi]; % x slit
ys = [ 0  -hi -hi  hi  hi   0]; % y slit
rotM = [cos(-rot) -sin(-rot); sin(-rot) cos(-rot)];      
vecTemp = rotM*[xs; ys]; % rotate slit
xs = vecTemp(1,:)+x;
ys = vecTemp(2,:)+y;
patch(xs,ys,colFa,'EdgeColor',colEd,'LineWidth',linw); %


function ellipse(x,y,a,b,col,linw)
% x coordinates of the center
% y coordinates of the center
% a radius1
% b radius2
% col  color [r g b]
% linw  line width
angleOffPoints = linspace(0,2.001*pi,300);
xe = x + a*cos(angleOffPoints);
ye = y + b*sin(angleOffPoints);
plot(xe,ye,'-','linewidth',linw,'color',col);


function map = struct2map(RGB)
% RGB: struct of depth 1 with [r g b] in each field
fNames = fieldnames(RGB);
nNames = numel(fNames);
map = NaN(nNames,3); % allocate
for iName = 1:nNames
    map(iName,:) = RGB.(fNames{iName}); % 
end


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];