function Rotating_black_and_white_squares()
% source code for drawing the animation
%
% 2017-04-26 Jahobr

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

fps = 50;
nFrames = 200;

rotShift = cos(linspace(0,pi,nFrames+1))*0.5+0.5;
rotShift = rotShift(1:end-1); % 0 and 2pi are the same image

xyLim = [-1 1];

figHandle = figure(15124455);
clf
axesHandle = axes;
hold(axesHandle,'on')
set(figHandle, 'Units','pixel');
set(figHandle, 'position',[1 1 1000 1000]); % big start image for antialiasing later [x y width hight]
set(axesHandle,'position',[0 0 1 1]); % stetch axis bigger as figure, easy way to get rid of ticks [x y width hight]
set(figHandle,'GraphicsSmoothing','on') % requires at least version 2014b
xlim(xyLim); ylim(xyLim); % set axis limits
axis equal; drawnow;

for iFrame = 1:nFrames
    cla(axesHandle) % fresh frame
    
    col = [0 0 0]; % start black
    
    x =[-1 -1 1  1];
    y =[-1  1 1 -1];

    k = rotShift(iFrame);
    
    while norm([x(1) y(1)]) > 0.5/900 % squares smaller than a pixel
        
        patch([x x(1)],[y y(1)],col,'EdgeColor','none');
        
        col = 1-col; % flip black and white

        x = [k*x(1)+(1-k)*x(2)  k*x(2)+(1-k)*x(3)  k*x(3)+(1-k)*x(4)  k*x(4)+(1-k)*x(1)]; % create next square
        y = [k*y(1)+(1-k)*y(2)  k*y(2)+(1-k)*y(3)  k*y(3)+(1-k)*y(4)  k*y(4)+(1-k)*y(1)]; % create next square
        
        linearScale = 0.9;
        constantOffset = 0.02;
        
        cornerRadius = norm([x(1) y(1)]);
        offsetScale =  (cornerRadius-constantOffset)/cornerRadius;
        
        % usedScale = offsetScale;
        % usedScale = linearScale;
        usedScale = mean([offsetScale linearScale]);
        
        x = usedScale*x;
        y = usedScale*y;
        
        if usedScale<0
            break
        end
    end
    
    %% save animation
    xlim(xyLim); ylim(xyLim); % set axis limits
    drawnow % update figure window and execute pending callbacks
    pause(0.01)
    
    f = getframe(figHandle);
    f.cdata = imresize(f.cdata,0.5); % the size reduction: adds antialiasing
    if iFrame== 1
        map=gray(8); % 8 colores % create color map % or use : [im,map] = rgb2ind(f.cdata,8,'nodither'); % 
        im = rgb2ind(f.cdata,map,'nodither'); % create first image
        im(1,1,1,nFrames) = 0; % allocate
    end
    
    imtemp = rgb2ind(f.cdata,map,'nodither');
    im(:,:,1,iFrame) = imtemp;
    
end
imwrite(im,map,fullfile(pathstr, [fname '.gif']),'DelayTime',1/fps,'LoopCount',inf) %
disp([fname '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit

end