function [] = Black_and_white_wiggler()
% source code for drawing the animation
%
% 2017-04-24 Jahobr

fps = 50;
nFrames = 200;

nBands = 30; % black + white
mergeCircles = 8; % circles merged to create a band
nCircles = nBands*mergeCircles; % individual circle patches

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

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

xyLim = [-1.5 1.5];

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.05 -0.05 1.1 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;

%% Radii of the circles, big -> small
radii = linspace(1.2,0.25,nCircles).^4; % 1.2^4 = 2.1 covers the whole image ... 0.25^4 = 0.004 is still visible;
% radii = logspace(log10(2),log10(0.003),40);

%% lag of each circle, big -> small
phaseShift = linspace(1.2,0.2,nCircles).^4;
% phaseShift = logspace(log10(0.5),log10(0.0005),nCircles);
% phaseShift = linspace(1.,0.0000000001,nCircles);

%% scale the movement bigger circles are slowed down, big -> small
movementScale = (1-cos(linspace(0,pi,nCircles)))*0.5;
% movementScale = ones(1,nCircles)
% movementScale = linspace(0.3,1,nCircles);
% movementScale = nthroot(linspace(0.3^3,1,nCircles),3);

angleOffPoints = ((2*pi) : -pi/200 : 0);

xTrajectory = zeros(1,nFrames); % allocate
yTrajectory = zeros(1,nFrames); % allocate

for iFrame = 1:nFrames
    cla(axesHandle) % fresh frame
    col = [0 0 0]; % start black
    
    for iCirc = 1:nCircles % draw stack of circles big to small
        
        nPatch = 1+mergeCircles-rem(iCirc,mergeCircles);
        
        % draw a figure-8
        xc=sin(   phaseShift(iCirc)+k(iFrame)) *movementScale(iCirc);
        yc=sin(2*(phaseShift(iCirc)+k(iFrame)))*movementScale(iCirc);
        
%         % draw a more complex movement
%         xc=sin(2*(phaseShift(iCirc)+k(iFrame)))*movementScale(iCirc);
%         yc=cos(3*(phaseShift(iCirc)+k(iFrame)))*movementScale(iCirc);

        [X(nPatch,:),Y(nPatch,:)] = pol2cart(angleOffPoints,radii(iCirc));

        X(nPatch,:) = X(nPatch,:)+xc; % center offset
        Y(nPatch,:) = Y(nPatch,:)+yc; % center offset        


        if nPatch == mergeCircles
            
            MergeX = X(nPatch,:);
            MergeY = Y(nPatch,:);
            
            for iMerge = 2:mergeCircles % draw stack of circles big to small
                [MergeX,MergeY] = polybool('union', MergeX,MergeY,...
                    X(iMerge,:),Y(iMerge,:));
            end
            patch(MergeX,MergeY,col,'EdgeColor',1-col);
            
            col = 1-col; % flip black and white
        end
    end
    
    xTrajectory(iFrame) = xc; % store trajectory of smallest circle
    yTrajectory(iFrame) = yc; % store trajectory of smallest circle
    
    %% save animation
    xlim(xyLim); ylim(xyLim); % set axis limits
    drawnow % update figure window and execute pending callbacks
    
    f = getframe(figHandle);
    f.cdata = imresize(f.cdata,0.5); % the size reduction: adds antialiasing, edge line thinner
    if iFrame == 1 % create colormap
        map = gray(8); % create own color map %
        im = rgb2ind(f.cdata,map,'nodither'); %
        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) % save gif
disp([fname '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit


%% Helper Figures
% check trajectory of the tip
figure(15634469); clf; hold on
plot(xTrajectory,yTrajectory)
title('trajectory')

figure(13234469); clf; hold on
subplot(1,3,1)
plot(radii)
title('radii')
subplot(1,3,2)
plot(phaseShift)
title('phase shift')
subplot(1,3,3)
plot(movementScale)
title('movement scale')