function ParallelogramLinkage1()
% source code that produces a GIF and a SVG
%
% 2017-04-26 Jahobr (update 2019-02-04 Jahobr)

leftBar = 2; 
xLeftBearing = -1.5;
yLeftBearing = 0;
centerBar = 3; 
rightBar = leftBar; %  parallel by definition
xRightBearing = xLeftBearing+centerBar; % parallel by definition
yRightBearing = yLeftBearing; % parallel by definition

RGB.bkgd = [1   1   1  ]; % white background
RGB.edge = [0   0   0  ]; % Edge color
RGB.bars = [0.3 0.3 0.3]; % grey
RGB.RofM = [1   0   1  ]; % magenta% range of movment
RGB.sVec = [1   0   0  ]; % red    % speed Vector

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


nFrames = 400;
startFrame = round(0.86*nFrames); % get a nice first frame with everything visible
intermediatePoints = 1; % extra points between to smooth lines
nPos = nFrames*intermediatePoints+1; % number of positions that have to be calculated

anglesLeft = linspace(0,2*pi,nPos); % define gear position in frames

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

%% calculate geometric values
xLeftJoint = xLeftBearing+cos(anglesLeft)*leftBar; % left bar end point
yLeftJoint = yLeftBearing-sin(anglesLeft)*leftBar; % left bar end point

xRightJoint = xRightBearing+cos(anglesLeft)*rightBar; % Right bar end point
yRightJoint = yRightBearing-sin(anglesLeft)*rightBar; % left bar end point

% extra bar to enforce parallel motion
offAng = pi/2;
xLeftJointExtra = xLeftBearing+cos(anglesLeft+offAng)*leftBar/2; % left bar end point
yLeftJointExtra = yLeftBearing-sin(anglesLeft+offAng)*leftBar/2; % left bar end point

xRightJointExtra = xRightBearing+cos(anglesLeft+offAng)*rightBar/2; % Right bar end point
yRightJointExtra = yRightBearing-sin(anglesLeft+offAng)*rightBar/2; % left bar end point


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

xLimits = [-4.0 4.0];
yLimits = [-2.5 2.5];

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
scaleReduction = 2; % not auto mode. (Line width is not programmed adaptive! I want all verions to look similar)
reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate

iFrame = 0;

%% plot loop
for iPos = 2:intermediatePoints:nPos % leave out first frame, it would be double
    
    saveName = [fname '_speedVector'];
    
    iFrame = iFrame+1;
    cla(axesHandle) % fresh frame
    
    sizze = 0.4;
    bearing(       xLeftBearing ,yLeftBearing ,sizze,RGB.edge)
    slidingBearing(xRightBearing,yRightBearing,sizze,RGB.edge)
    
    plot(xLeftJoint, yLeftJoint, '--','LineWidth',2.5,'Color',RGB.RofM) % range of movment
    plot(xRightJoint,yRightJoint,'--','LineWidth',2.5,'Color',RGB.RofM) % range of movment
    
    xLeftBeams  = [xLeftBearing xLeftJointExtra(iPos) xLeftJoint(iPos) xLeftBearing];
    yLeftBeams  = [yLeftBearing yLeftJointExtra(iPos) yLeftJoint(iPos) yLeftBearing];
    plot(xLeftBeams,yLeftBeams,'-','MarkerSize',15,'LineWidth',8,'Color',RGB.bars); % left triangle
    xRightBeams = [xRightBearing xRightJoint(iPos) xRightJointExtra(iPos) xRightBearing];
    yRightBeams = [yRightBearing yRightJoint(iPos) yRightJointExtra(iPos) yRightBearing];
    plot(xRightBeams,yRightBeams,'-','MarkerSize',15,'LineWidth',8,'Color',RGB.bars); % right triangle
    plot([xLeftJointExtra(iPos) xRightJointExtra(iPos)],[yLeftJointExtra(iPos) yRightJointExtra(iPos)],'-','MarkerSize',15,'LineWidth',8,'Color',RGB.bars); % extra bar to enforce parallel motion
    plot([xLeftJoint(iPos) xRightJoint(iPos)],[yLeftJoint(iPos) yRightJoint(iPos)],'-','MarkerSize',15,'LineWidth',12,'Color',RGB.edge); % main bar
    
    
    plot([xLeftBeams xRightBeams],[yLeftBeams yRightBeams],'o','LineWidth',3,'MarkerEdgeColor',RGB.edge,'MarkerFaceColor',RGB.bkgd,'MarkerSize',12); % joints
    
    diffVecLeft  = [xLeftJoint(iPos), yLeftJoint(iPos)]  - [xLeftJoint(iPos-1), yLeftJoint(iPos-1)];  % lazy way to approximate the speed vector
    diffVecRight = [xRightJoint(iPos),yRightJoint(iPos)] - [xRightJoint(iPos-1),yRightJoint(iPos-1)]; % lazy way to approximate the speed vector
    
    speedScale = intermediatePoints*30;
    
    line(xLeftJoint(iPos)+[0 diffVecLeft(1)]*speedScale, yLeftJoint(iPos)+[0 diffVecLeft(2)]*speedScale,    'Color',RGB.sVec,'LineWidth',4) % speed Vector
    plot(xLeftJoint(iPos)+   diffVecLeft(1) *speedScale, yLeftJoint(iPos)+   diffVecLeft(2) *speedScale,'.','Color',RGB.sVec,'LineWidth',4,'MarkerSize',20) % speed Vector end marker
    
    line(xRightJoint(iPos)+[0 diffVecRight(1)]*speedScale,yRightJoint(iPos)+[0 diffVecRight(2)]*speedScale,    'Color',RGB.sVec,'LineWidth',4) % speed Vector
    plot(xRightJoint(iPos)+   diffVecRight(1)*speedScale, yRightJoint(iPos)+   diffVecRight(2) *speedScale,'.','Color',RGB.sVec,'LineWidth',4,'MarkerSize',20) % speed Vector end marker
    
    
    %% resize figure
    axis equal;
    set(axesHandle,'Position',[0 0 1 1]); % stretch axis as big as figure, [x y width height]
    set(figHandle, 'Position',[1 1 xSize*scaleReduction ySize*scaleReduction]); % big start image for antialiasing later [x y width height]
    xlim(xLimits); ylim(yLimits);
    drawnow;
    pause(0.005)
    
    % save SVG
    if iFrame == startFrame
        if ~isempty(which('plot2svg'))
            plot2svg(fullfile(pathstr, [saveName '.svg']),figHandle) % by Juerg Schwizer
        else
            disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');
        end
    end
    
    %% save animation
    f = getframe(figHandle);
    reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
end

reducedRGBimage = circshift(reducedRGBimage,1-startFrame,4); % shift animation do get nice start frame

map = createImMap(reducedRGBimage,32,[RGB.bkgd;RGB.edge;RGB.sVec;RGB.bars;RGB.RofM]); % only speed Vector

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 '.gif']),'DelayTime',1/25,'LoopCount',inf) % save gif
disp([saveName '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit


if ispc; dos(['explorer ' pathstr]); end % open folder with files in it
return
%%

function bearing(x,y,sizze,col)
% x coordinates of the center
% y coordinates of the center
% size
plot([0 -0.5 0.5 0]*sizze+x,[0 -0.8660 -0.8660 0]*sizze+y,'k','LineWidth',3,'Color',col); % Triangle %  0.8660 = sqrt(3)*0.5
plot([-0.7 0.7]*sizze+x,[-0.87 -0.87]*sizze+y,'k','LineWidth',3,'Color',col); % base line
for iLine = -0.6:0.2:0.7
    plot(([-0.1 0.1]+iLine)*sizze+x,[-1.07 -0.87]*sizze+y,'k','LineWidth',2,'Color',col); % Hatching
end

function slidingBearing(x,y,sizze,col)
% x coordinates of the center
% y coordinates of the center
% size
plot([0 -0.5 0.5 0]*sizze+x,[0 -0.8660 -0.8660 0]*sizze+y,'k','LineWidth',3,'Color',col); % Triangle %  0.8660 = sqrt(3)*0.5
plot([-0.7 0.7]*sizze+x,[-1.05 -1.05]*sizze+y,'k','LineWidth',3,'Color',col); % base line
for iLine = -0.6:0.2:0.7
    plot(([-0.1 0.1]+iLine)*sizze+x,[-1.25 -1.05]*sizze+y,'k','LineWidth',2,'Color',col); % Hatching
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:
%  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];