function Tusi_couple_vs_Trammel()
% source code for drawing Tusi_couple_vs_Trammel
% produces a GIF and a SVG
%
% 2017-04-10 Jahobr

p = 1;
nFrames = 170;
angleSmall = linspace(0,2*pi,nFrames+1); % define gear position in frames
angleSmall = angleSmall(1:end-1); % remove last frame, it would be double

figHandle = figure(15674455);
clf
axesHandle = axes;
hold(axesHandle,'on')
set(figHandle, 'Units','pixel');


for currentCase = 1:3

    axis equal; drawnow;
    
    for iFrame = 1:nFrames;
        
        boxN = 0.1*p; % narrow side
        boxL = 0.2*p; % long side
        
        currentAngle = angleSmall(iFrame);
        
        cla(axesHandle) % fresh frame
        plot([-p,-boxN/2 -boxN/2],[ boxN/2  boxN/2  p],'-','color',[0.8 0.8 0.8],'linewidth',3) % horizontal line
        plot([ p, boxN/2  boxN/2],[ boxN/2  boxN/2  p],'-','color',[0.8 0.8 0.8],'linewidth',3) % horizontal line
        plot([ p, boxN/2  boxN/2],[-boxN/2 -boxN/2 -p],'-','color',[0.8 0.8 0.8],'linewidth',3) % horizontal line
        plot([-p,-boxN/2 -boxN/2],[-boxN/2 -boxN/2 -p],'-','color',[0.8 0.8 0.8],'linewidth',3) % horizontal line
        
        barCol = [0.5 0.5 0.5]; % color of trammel bar
        
        
        switch currentCase
            case 1 % no ellipses
                set(figHandle, 'position',[1 1 1000 1000]); % big start image for antialiasing later [x y width hight]
                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(p*[-1.3 1.3]);
                ylim(p*[-1.3 1.3]);
                
                plot([p*cos(currentAngle)  0],[0 p*sin(currentAngle)],'.-','color',[0.5 0.5 0.5],'linewidth',35,'MarkerSize',100); % grey trammel bar
            case 2 % ellipses vertical
                set(figHandle, 'position',[1 1 650 1000]); % big start image for antialiasing later [x y width hight]
                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(p*[-1.3 1.3]);
                ylim(p*[-2 2]);
                plot([p*cos(currentAngle) -0.7*p*cos(currentAngle)],[0 1.7*p*sin(currentAngle)],'.-','color',[0.5 0.5 0.5],'linewidth',35,'MarkerSize',100); % grey trammel bar
            case 3 % ellipses horizontal
                set(figHandle, 'position',[1 1 1100 715]); % big start image for antialiasing later [x y width hight]
                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(p*[-2 2]);
                ylim(p*[-1.3 1.3]);
                plot([1.7*p*cos(currentAngle) 0],[-0.7*p*sin(currentAngle) p*sin(currentAngle)],'.-','color',[0.5 0.5 0.5],'linewidth',35,'MarkerSize',100); % grey trammel bar
        end
               
        plot(p*cos(currentAngle),p*sin(currentAngle),'.m','MarkerSize',40) % touching point / Instant centre of rotation
        
        circle(0,0,p) % outer circle
        circle(cos(currentAngle)*0.5,sin(currentAngle)*0.5,p/2)  % inner circle
        
        colHori = [0.1 0.7 0.1]; % green
        rectangle(p*cos(currentAngle),0,boxL,boxN,colHori)
        plot(p*cos(currentAngle),0,'.','MarkerSize',30,'linewidth',4,'color',barCol) % overwrite line with trammel bar color
        plot([p,-p],[0 0],':','MarkerSize',25,'linewidth',4,'color',colHori) % horizontal line
        plot(p*cos(currentAngle),0,'.','MarkerSize',20,'linewidth',4,'color',colHori) % small dot to indicate the precise position
        
        rectangle(0,p*sin(currentAngle),boxN,boxL,[0 0 1])
        plot(0,p*sin(currentAngle),'.','MarkerSize',30,'linewidth',4,'color',barCol) % overwrite line with trammel bar color
        plot([0 0],[p,-p],'b:','MarkerSize',25,'linewidth',4) % vertical line
        plot(0,p*sin(currentAngle),'b.','MarkerSize',20,'linewidth',4) % small dot to indicate the precise position
        
        col1   = [0  0.9 0.9]; % cyan
        ellOneOffset = 0.2;
        col2 = [1 0 0]; % red
        ellTwoOffset = 0.7;
        switch currentCase
            case 1 % no ellipses
                saveName = 'Tusi_couple_vs_Trammel';
                if iFrame == 10
                    plot2svg([saveName '.svg'],figHandle) % by Juerg Schwizer, See http://www.zhinst.com/blogs/schwizer/
                end
            case 2 % ellipses vertical
                saveName = 'Tusi_couple_vs_Trammel_plus_Ellipses';
                ellipse(0,0,ellOneOffset,p+ellOneOffset,col1,4);
                plot(-ellOneOffset*cos(currentAngle),(p+ellOneOffset)*sin(currentAngle),'.','MarkerSize',40,'linewidth',3,'color',col1);
                ellipse(0,0,ellTwoOffset,p+ellTwoOffset,col2,4);
                plot(-ellTwoOffset*cos(currentAngle),(p+ellTwoOffset)*sin(currentAngle),'.','MarkerSize',40,'linewidth',3,'color',col2);
                if iFrame == 12
                    plot2svg([saveName '.svg'],figHandle) % by Juerg Schwizer, See http://www.zhinst.com/blogs/schwizer/
                end
            case 3 % ellipses horizontal
                saveName = 'Tusi_couple_vs_Trammel_plus_Ellipses2';
                ellipse(0,0,p+ellOneOffset,ellOneOffset,col1,4);
                plot((p+ellOneOffset)*cos(currentAngle),-ellOneOffset*sin(currentAngle),'.','MarkerSize',40,'linewidth',3,'color',col1);
                ellipse(0,0,p+ellTwoOffset,ellTwoOffset,col2,4);
                plot((p+ellTwoOffset)*cos(currentAngle),-ellTwoOffset*sin(currentAngle),'.','MarkerSize',40,'linewidth',3,'color',col2);
                if iFrame == 30
                    plot2svg([saveName '.svg'],figHandle) % by Juerg Schwizer, See http://www.zhinst.com/blogs/schwizer/
                end
        end

        %% save animation
        drawnow
        pause(0.05)
        f = getframe(figHandle);
        f.cdata = imresize(f.cdata,0.5); % the size reduction: adds antialiasing
        if iFrame == 1 % create colormap
            [im,map] = rgb2ind(f.cdata,32,'nodither'); % 32 colors % create color map %% THE FIRST FRAME MUST INCLUDE ALL COLORES !!!
            % FIX WHITE, rgb2ind sets white to [0.9961    0.9961    0.9961], which is annoying
            [~,wIndex] = max(sum(map,2)); % find "white"
            map(wIndex,:) = 1; % make it truly white
            im(1,1,1,nFrames) = 0; % allocate
        end
     
        imtemp = rgb2ind(f.cdata,map,'nodither');
        im(:,:,1,iFrame) = imtemp;
        
    end
    imwrite(im,map,[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
return
%%

function circle(x,y,r)
% x coordinates of the center
% y coordinates of the center
% r is the radius of the circle
angleOffPoints = linspace(0,2*pi,300);
xc = x + r*cos(angleOffPoints);
yc = y + r*sin(angleOffPoints);
plot(xc,yc,'k','linewidth',4);

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

function rectangle(x,y,width,height,col)
% x coordinates of the center
% y coordinates of the center
% width
% height
plot([x-width/2 x+width/2 x+width/2 x-width/2 x-width/2],[y-height/2 y-height/2 y+height/2 y+height/2 y-height/2],'-','linewidth',3,'color',col);