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

leftBar = 1;
xLeftBearing = 0;
yLeftBearing = 0;
centerBar = 2.5;
rightBar = 1.5;
xRightBearing = 2.5;
yRightBearing = -0.5;

RGB.bkgd = [1   1   1  ]; % white background
RGB.edge = [0   0   0  ]; % Edge color
RGB.bars = [0.3 0.3 0.3]; % grey
RGB.icro = [1   1   0.1]; % yellow % Instant center of rotation
RGB.fixC = [0.2 0.2 1  ]; % blue   % fixed Centrode
RGB.movC = [0.1 0.7 0.1]; % green  % moving Centrode
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.12*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 movemet of left bar

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

%% allocate memory
xRightJoint           = Inf(1,nPos); yRightJoint           = Inf(1,nPos);
xInstCentRot          = Inf(1,nPos); yInstCentRot          = Inf(1,nPos);
angleCenterBar        = Inf(1,nPos);
xInstCentRot_MovCoord = Inf(1,nPos); yInstCentRot_MovCoord = Inf(1,nPos);
    
%% calculate geometric values
xLeftJoint = xLeftBearing+cos(anglesLeft)*leftBar; % left bar end point
yLeftJoint = yLeftBearing-sin(anglesLeft)*leftBar; % left bar end point

for iPos = 1:nPos % calculate geometric values for all positions
    
    [xout,yout] = circcirc(xLeftJoint(iPos),yLeftJoint(iPos),centerBar,xRightBearing,yRightBearing,rightBar); % Intersections of circles in Cartesian plane
    if iPos==1
        xRightJoint(iPos) = xout(1); % select inital configuration
        yRightJoint(iPos) = yout(1); % select inital configuration
    else % predict based on last move, % next neighbour based on prediction
        if iPos==2 % next neighbour
            xPred = xRightJoint(iPos-1); % 1 old point is not enough for a prediction; just use last point
            yPred = yRightJoint(iPos-1); % 1 old point is not enough for a prediction; just use last point
        else % make prediction
            xPred = 2*xRightJoint(iPos-1)-xRightJoint(iPos-2); % last point + last delta;(iPos-1)-(iPos-2)
            yPred = 2*yRightJoint(iPos-1)-yRightJoint(iPos-2); % last point + last delta;(iPos-1)-(iPos-2)
        end
        dist1 = norm([xPred-xout(1) yPred-yout(1)]); % distance to prediction
        dist2 = norm([xPred-xout(2) yPred-yout(2)]); % distance to prediction
        if ~isnan(xout(1)) % normal case; calculation worked
            if dist1 < dist2 % choose neighbour
                xRightJoint(iPos) = xout(1);
                yRightJoint(iPos) = yout(1);
            else
                xRightJoint(iPos) = xout(2);
                yRightJoint(iPos) = yout(2);
            end
            
        else % if error; this happens on extreme endpoints of movements
            xRightJoint(iPos) = xPred; % failsafe hack; not nice but it does the job
            yRightJoint(iPos) = yPred; % failsafe hack; not nice but it does the job
        end
    end

    % calculate intersection to dermine instant center of rotation
    [xInstCentRot(iPos),yInstCentRot(iPos)] = intersectionOfLines(...
        xLeftBearing,     yLeftBearing,...     % Point1 Line1
        xLeftJoint(iPos), yLeftJoint(iPos),... % Point2 Line1
        xRightBearing,    yRightBearing,...    % Point1 Line2
        xRightJoint(iPos),yRightJoint(iPos));  % Point2 Line2
        
    %% instant center of rotation in center bar coordinates
    angleCenterBar(iPos) = atan2((yRightJoint(iPos) - yLeftJoint(iPos)) , (xRightJoint(iPos) - xLeftJoint(iPos)));
    rotM = [cos(-angleCenterBar(iPos)) -sin(-angleCenterBar(iPos)); sin(-angleCenterBar(iPos)) cos(-angleCenterBar(iPos))];
    vecTemp = rotM*[xInstCentRot(iPos)-xLeftJoint(iPos);  yInstCentRot(iPos)-yLeftJoint(iPos)];
    xInstCentRot_MovCoord(iPos) = vecTemp(1);
    yInstCentRot_MovCoord(iPos) = vecTemp(2);
end

%% filter values out of visilbe range; this cleans up the plot from "unwrapping" lines
xInstCentRotPlot = xInstCentRot;
xInstCentRotPlot(xInstCentRot> 50) = NaN;
xInstCentRotPlot(xInstCentRot<-50) = NaN;

yInstCentRotPlot = yInstCentRot;
yInstCentRotPlot(yInstCentRot> 50) = NaN;
yInstCentRotPlot(yInstCentRot<-50) = NaN;

xInstCentRot_MovCoordPlot = xInstCentRot_MovCoord;
xInstCentRot_MovCoordPlot(xInstCentRot_MovCoord> 50) = NaN;
xInstCentRot_MovCoordPlot(xInstCentRot_MovCoord<-50) = NaN;

yInstCentRot_MovCoordPlot = yInstCentRot_MovCoord;
yInstCentRot_MovCoordPlot(yInstCentRot_MovCoord> 50) = NaN;
yInstCentRot_MovCoordPlot(yInstCentRot_MovCoord<-50) = NaN;

%% 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

%% plot loop
for currentCase = 1:5
    
    switch currentCase
        case 1 %
            saveName = [fname '_speedVector'];
        case 2 %
            saveName = [fname '_InstantCenterRotation']; % Momentanpol
        case 3 %
            saveName = [fname '_fixedCentrode']; % Rastpolbahn
        case 4 %
            saveName = [fname '_movingCentrode']; % Gangpolbahn
        case 5 %
            saveName = [fname '_Centrodes'];
    end
    
    xLimits = [-1.4 5.0];
    yLimits = [-2.7 2.7];

    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)
       
    if currentCase == 1
        xLimits = [-1.4 4.0]; % reduced range
        yLimits = [-1.4 1.4]; % reduced range
        xRange = xLimits(2)-xLimits(1);
        yRange = yLimits(2)-yLimits(1);
        ySize = round(xSize/xRange*yRange); % ony adapt image hight
    end
    
    reducedRGBimage = uint8(ones(ySize,xSize,3,nFrames)); % allocate
    
    iFrame = 0;
    for iPos = 2:intermediatePoints:nPos % leave out first frame, it would be double
        
        iFrame = iFrame+1;
        cla(axesHandle) % fresh frame

        sizze = 0.3;
        bearing(xLeftBearing ,yLeftBearing ,sizze,RGB.edge)
        bearing(xRightBearing,yRightBearing,sizze,RGB.edge)
        
        if currentCase~=1
            plot([xLeftJoint(iPos)  xInstCentRot(iPos)], [yLeftJoint(iPos)  yInstCentRot(iPos)],'--','LineWidth',2.5,'Color',RGB.bars) % line to intersection
            plot([xRightJoint(iPos) xInstCentRot(iPos)], [yRightJoint(iPos) yInstCentRot(iPos)],'--','LineWidth',2.5,'Color',RGB.bars) % line to intersection
        end
        
        plot(xLeftJoint,yLeftJoint,'--','LineWidth',2.5,'Color',RGB.RofM) % range of movment
        
        [~,xMinInd] = min(xRightJoint); % find ranges to avoid doubled lines
        [~,xMaxInd] = max(xRightJoint); % find ranges to avoid doubled lines
        rangeInd = sort([xMinInd xMaxInd]);
        plot(xRightJoint(rangeInd(1):rangeInd(2)),yRightJoint(rangeInd(1):rangeInd(2)),'--','LineWidth',2.5,'Color',RGB.RofM) % range of movment
        
        xBeams = [xLeftBearing xLeftJoint(iPos) xRightJoint(iPos) xRightBearing];
        yBeams = [yLeftBearing yLeftJoint(iPos) yRightJoint(iPos) yRightBearing];
        plot(xBeams,yBeams,'-','MarkerSize',15,'LineWidth',8,'Color',RGB.bars); % all bars
        plot(xBeams(2:3),yBeams(2:3),'-','MarkerSize',15,'LineWidth',12,'Color',RGB.edge); % main bar
        
        if or( currentCase==3 , currentCase==5 ) % fixed Centrode
            plot(xInstCentRotPlot,yInstCentRotPlot,'-','LineWidth',4,'Color',RGB.fixC);
        end
        
        if or( currentCase==4 , currentCase==5 ) % moving Centrode
            rotM = [cos(angleCenterBar(iPos)) -sin(angleCenterBar(iPos)); sin(angleCenterBar(iPos)) cos(angleCenterBar(iPos))];
            vecTemp = rotM*[xInstCentRot_MovCoordPlot; yInstCentRot_MovCoordPlot];
            plot(vecTemp(1,:)+xLeftJoint(iPos),vecTemp(2,:)+yLeftJoint(iPos),'-','LineWidth',4,'Color',RGB.movC);
        end
        
        plot(xBeams,yBeams,'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
        normalVecRight = [xRightJoint(iPos),yRightJoint(iPos)] - [xRightJoint(iPos-1),yRightJoint(iPos-1)]; % lazy way to approximate the speed vector
        
        speedScale = intermediatePoints*35;
        
        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 normalVecRight(1)]*speedScale,yRightJoint(iPos)+[0 normalVecRight(2)]*speedScale,    'Color',RGB.sVec,'LineWidth',4) % speed Vector
        plot(xRightJoint(iPos)+   normalVecRight(1)*speedScale, yRightJoint(iPos)+   normalVecRight(2) *speedScale,'.','Color',RGB.sVec,'LineWidth',4,'MarkerSize',20) % speed Vector end marker
        
        if currentCase~=1
            plot(xInstCentRotPlot(iPos),yInstCentRotPlot(iPos),'o','LineWidth',2,'MarkerEdgeColor',RGB.edge,'MarkerFaceColor',RGB.icro,'MarkerSize',14); % Instant center of rotation
            plot(xInstCentRotPlot(iPos),yInstCentRotPlot(iPos),'.','MarkerSize',5,'Color',RGB.edge); % Instant center of rotation
        end
        
        %% resize figure
        axis equal;
        set(axesHandle,'Position',[-0.01 -0.01 1.02 1.02]); % stretch axis bigger as figure, to have nice line ends [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
    
    switch currentCase
            case 1 %
                map = createImMap(reducedRGBimage,32,[RGB.bkgd;RGB.edge;RGB.sVec;RGB.bars;RGB.RofM]); % only speed Vector
            case 2 %
                map = createImMap(reducedRGBimage,32,[RGB.bkgd;RGB.edge;RGB.sVec;RGB.bars;RGB.RofM;RGB.icro]); % Instant Center Rotation
            case 3 % 
                map = createImMap(reducedRGBimage,64,[RGB.bkgd;RGB.edge;RGB.sVec;RGB.bars;RGB.RofM;RGB.icro;RGB.fixC]); % fixed Centrode
            case 4 % 
                map = createImMap(reducedRGBimage,64,[RGB.bkgd;RGB.edge;RGB.sVec;RGB.bars;RGB.RofM;RGB.icro;         RGB.movC]); % moving Centrode
            case 5 %
                map = createImMap(reducedRGBimage,64,[RGB.bkgd;RGB.edge;RGB.sVec;RGB.bars;RGB.RofM;RGB.icro;RGB.fixC;RGB.movC]); % both Centrodes
    end
    
    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
end

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 [xs,ys] = intersectionOfLines(x1,y1, x2,y2, x3,y3, x4,y4)
% Point1_Line1 % Point2_Line1 % Point1_Line2 % Point2_Line2
xs = ((x4-x3)*(x2*y1-x1*y2) - (x2-x1)*(x4*y3-x3*y4)) / ((y4-y3)*(x2-x1) - (y2-y1)*(x4-x3));
ys = ((y1-y2)*(x4*y3-x3*y4) - (y3-y4)*(x2*y1-x1*y2)) / ((y4-y3)*(x2-x1) - (y2-y1)*(x4-x3));

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