function [] = Oil_Damper()
% Programmed in Matlab R2017a.
% Animation of a twin and mono tube oil damper.
% Before anything is drawn the algorithm loops the frames
%    to determine the particle pathes.
% Several versions are created.
%    .gif, .mp4 and .webm
%
% 2020-08-07 Jahobr - Licensing: CC0 1.0 Universal Public Domain Dedication


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

RGB.edge       = [0    0    0  ];
RGB.white      = [1    1    1  ];
RGB.brightGrey = [0.8  0.8  0.8]; % housing
RGB.blueGrey   = [0.75 0.75 0.9]; % moving metal
RGB.oil        = [0.9  0.9  0  ];

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

screenSize = get(groot,'Screensize')-[0 0 5 20]; % [1 1 width height] (minus tolerance for figure borders)

for mode = {'Twin_Tube_Damper','Mono_Tube_Damper'} %
    
    %% y-axis-dimensions (Twin is default)
    %                     _______________________________________________________________
    wall = 0.5;          %                        wall width                             ¦
    outerRad = 4.5;      %  ,-------------------------outerRadius--------------------,   ¦
    %                    %  ¦    ____________________________________________________¦   ¦
    %                    %  ¦   ¦                         wall width                     ¦
    innerRad = 3.0;      %  ¦   ¦  ,---------,====,--innerRadius-(pistonRadius)--,       ¦
    %                    %  ¦   ¦__¦         ¦____¦                              ¦       ¦
    holeY = [2.2  1.5];  %  ¦   :__:  holeY  :____:                              ¦       ¦
    %                    %  ¦   ¦  ¦         ¦    ¦                              ¦       ¦
    rodRad = 1.0;        %  ¦   ¦  ¦         ¦    '--------rodRadius-----------------------------,
    %                    %  ¦   ¦  ¦         ¦                                                   ¦
    %% x-axis-dimensions    + centerline - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
    %                    %  ¦   ¦  0-Position¦                                                   ¦
    rodLenght = 20;      %  ¦   ¦  ¦         ¦    ,<----------------rodLenght------------------->'
    %                    %  ¦   ¦__¦         ¦    ¦                              ¦       ¦
    pistonWidth = 3;     %  ¦   :__:       pistonWidth                           ¦       ¦
    l_inner = 0.0;       %  ¦   ¦  ¦l_inner  ¦    ¦                              ¦       ¦
    r_inner = 20;        %  ¦   ¦  '----------====------------------------r_inner'       ¦
    %                    %  ¦   ¦____________________________________________________,   ¦
    l_outer = -2.5;      %  ¦left_outer                                              ¦   ¦
    r_outer = r_inner+3*wall; %  ---------------------------------------------r_outer'   ¦
    %                    %_______________________________________________________________¦
    
    xLimits = [-4.1 38.2]; % ADJUST
    yLimits = [-5.8  5.8]; % ADJUST

    switch mode{1} %
        case 'Mono_Tube_Damper'
            l_inner = -4.0;
            xLimits = [-5.5 38.2]; % ADJUST
            yLimits = [-4.7  4.7]; % ADJUST
    end
    
    nParticlesOil = 250; % number of moving particles (oil "bubbles")
    nParticlesGas = 50;  % number of moving particles (gas)

    xRange = xLimits(2)-xLimits(1);
    yRange = yLimits(2)-yLimits(1);
    imageAspectRatio = xRange/yRange;

    for form = {'gif'} % {'mp4','gif'} % file format

        switch form{1} %
            case 'gif'
                
                nFrames = 520;
                pistonPosList = cos(linspace(2*pi,0,nFrames)); % [1 -1 1] one loop
                
                delayTime = ones(1,nFrames)*1/30; % delay per frame
                %             delayTime(1) = 0.5;  % keep first frame longer (did not look good)
                %             delayTime(end) = 0.5; % keep last frame longer (did not look good)
                %             delayTime(round(nFrames/2)) = 1; % keep other extreme position frame longer
                
                MegaPixelTarget = 100*10^6; % Category:Animated GIF files exceeding the 100 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
                scaleReduction = min(...% repeat as often as possible for nice antialiasing
                    floor(screenSize(4)/ySize), floor(screenSize(3)/xSize));
                
                if scaleReduction == 0;   error('"MegaPixelTarget" not possible; use smaller target or bigger monitor'); end % check
                
            case 'mp4'
                repetitions = 3;
                nFrames = 520*repetitions;
                %nFrames = 30;
                pistonPosList = cos(linspace(2*pi*repetitions,0,nFrames)); % [1 -1 1] one loop
                
                xSize = screenSize(3); %max size
                ySize = xSize/imageAspectRatio; % gif width
                xSize = floor(xSize/16)*16; ySize = floor(ySize/16)*16; % full pixels, multiple-of-16 boundaries (MPEG splits the video into 16x16 squares called macroblocks)
                scaleReduction = 1; % not usefull for mp4 or webm (gif only)
            otherwise
                error('unknown mode')
        end
        
        oilRed = ceil(xSize*scaleReduction/500); % size reduction factor to speed things up
        
        figPos = [1 1 xSize*scaleReduction ySize*scaleReduction]; % big start image for antialiasing later [x y width height]

        liSc = ySize*scaleReduction/150; % LineWidth scale; LineWidth is absolute, a bigger images needs thicker lines to keep them in proportion
        
        yRange = xRange * ySize/xSize;  % xSize and ySize were rounded; we have the adjust the yRange to stay in the same aspect ratio
        yLimits = [-yRange/2 yRange/2]; % xSize and ySize were rounded; we have the adjust the yRange to stay in the same aspect ratio
        
        figHandle = figure(117854461); clf
        axesHandleOil = axes; hold(axesHandleOil,'on')
        set(figHandle, 'Position', figPos);
        if ~all(get(figHandle, 'Position') == figPos);   error('figure Position could not be set');   end % check
        set(figHandle,'Color',RGB.white); % white background
        set(axesHandleOil,'Position',[0 0 1 1]); % stretch axis as big as figure, [x y width height]
        xlim(xLimits); ylim(yLimits); % set axis limits
        axis off % invisible axes (no ticks)
        axis equal; drawnow;
        
        firstIter = true;

        pistonPosList = (pistonPosList(:)+1)*0.5; % normalize range to [0 1]
        fullPistonPosTravel = r_inner-pistonWidth; % left end is defined as 0
        
        pistonPosList = (pistonPosList*fullPistonPosTravel*0.75) + (fullPistonPosTravel*0.125); % 75% travel
        % pistonPosList = pistonPosList*fullPistonPosTravel; % 100% travel

        ctrs = NaN(nParticlesOil,2,nFrames); % allocate paths of clusters(numberOFCluster, x-y-coordinate, frame)
        
        oilPos = 18-pistonPosList;  % Twin_Tube_Damper
        pistonPos = [pistonPosList pistonPosList+pistonWidth]; % Mono & Twin_Tube_Damper
        
        AirOilSeparatorPos = 5+((pistonPosList-fullPistonPosTravel)/3);      % Mono_Tube_Damper
        AirOilSeparatorPos = [AirOilSeparatorPos-wall*4 AirOilSeparatorPos]; % Mono_Tube_Damper
        
        %% particle path loop
        %               r g b
        OilColStatic = [0 0 0]; % black,   static
        OilColOuter  = [0 0 1]; % blue,    outer chamber moves left slow
        OilColHolesB = [0 1 0]; % green,   holes base moves right fast
        OilColHolesP = [1 0 0]; % red,     holes piston moves left fast
        OilColBase   = [1 1 0]; % yellow,  outer moves left slow
        OilColBaseUp = [1 0 1]; % magenta, vertical movement
        OilColBaseDo = [0 1 1]; % cyan,    vertical movement
        
        for iFrame = 1:nFrames
            
            figure(figHandle); cla
                        
            plotBox([pistonPos(iFrame,2) r_inner], [ innerRad  rodRad], OilColStatic, 'none', liSc) % inner oil chamber above the rod
            plotBox([pistonPos(iFrame,2) r_inner], [-rodRad -innerRad], OilColStatic, 'none', liSc) % inner oil chamber below the rod
            
            if strcmp(mode{1},'Twin_Tube_Damper')
                plotBox([l_inner pistonPos(iFrame,1)], [-1 1]*innerRad, OilColBase, 'none', liSc) % inner oil chamber, in front of the piston
            else % Mono_Tube_Damper
                plotBox([AirOilSeparatorPos(iFrame,2) pistonPos(iFrame,1)], [-1 1]*innerRad, OilColBase, 'none', liSc) % inner oil chamber, in front of the piston
            end
            
            plotBox(pistonPos(iFrame,:)+[-wall wall]*2, +holeY, OilColHolesP, 'none', liSc) % upper hole, through the piston
            plotBox(pistonPos(iFrame,:)+[-wall wall]*2, -holeY, OilColHolesP, 'none', liSc) % lower hole, through the piston           
            
            if strcmp(mode{1},'Twin_Tube_Damper')
                plotBox([l_outer oilPos(iFrame)], [innerRad+wall outerRad],OilColOuter, 'none' ,liSc) % upper outer oil chamber
                plotBox([l_outer oilPos(iFrame)],-[innerRad+wall outerRad],OilColOuter, 'none', liSc) % lower outer oil chamber
                
                plotBox([l_outer l_inner-2*wall], [-1 1]*(innerRad+wall),     OilColStatic, 'none', liSc) % outer oil chamber bottom, static oil
                plotBox([l_outer l_inner-2*wall], [-outerRad+wall -holeY(2)], OilColBaseUp, 'none', liSc) % outer oil chamber bottom, oil up
                plotBox([l_outer l_inner-2*wall], [+outerRad-wall +holeY(2)], OilColBaseDo, 'none', liSc) % outer oil chamber bottom, oil down
                
                plotBox([l_outer l_inner+wall*2],+holeY, OilColHolesB, 'none', liSc) % upper hole, connecting inner and outer volume
                plotBox([l_outer l_inner+wall*2],-holeY, OilColHolesB, 'none', liSc) % lower hole, connecting inner and outer volume
            end
            
            xlim(xLimits); ylim(yLimits); % set axis limits
            axis equal; drawnow;
            f = getframe(figHandle);
            OilImage = f.cdata;
            
            isOilMap = ~and(and(OilImage(:,:,1)==uint8(255),OilImage(:,:,2)==uint8(255)),OilImage(:,:,3)==uint8(255));
            isOilMap = imerode(isOilMap,strel('rectangle',[1 1].*round(ySize*scaleReduction/35))); % paricles away from border
            
            OilImage = randSubsampling(OilImage, oilRed); % reduce size
            isOilMap = randSubsampling(isOilMap, oilRed); % reduce size
            
            figure(23446); clf
            image(isOilMap*255)
            
            [row,col] = find(isOilMap);

            optsFirst = statset('Display', 'final', 'UseParallel', true, 'MaxIter', 3000);
            optsIter  = statset('Display', 'final', 'UseParallel', true, 'MaxIter', 1);
            if firstIter
                [idx,ctrs(:,:,iFrame)] = kmeans([row,col],nParticlesOil,'Distance','sqEuclidean','start','uniform',...
                    'Options',optsFirst,'replicates',40);
                firstIter = false;
            else
                
                disp('nextIter')
                movSpeed = pistonPosList(iFrame)-pistonPosList(iFrame-1);
                movSpeedPx = movSpeed/xRange*xSize/oilRed*scaleReduction;
                ctrs(:,:,iFrame) = ctrs(:,:,iFrame-1); % use previous state
                
                for iCluster = 1:nParticlesOil
                    colorParticle = OilImage(round(ctrs(iCluster,1,iFrame)),round(ctrs(iCluster,2,iFrame)),:);
                    colorParticle = colorParticle(:)';
                    if(all(colorParticle==OilColStatic*255)) % black static
                        % do nothing
                    elseif(all(colorParticle==OilColOuter*255))  % = [0 0 1]; blue, outer chamber moves left slow
                        ctrs(iCluster,2,iFrame) = ctrs(iCluster,2,iFrame)-movSpeedPx;
                    elseif(all(colorParticle==OilColHolesB*255)) % = [0 1 0]; green, holes base moves right fast
                        ctrs(iCluster,2,iFrame) = ctrs(iCluster,2,iFrame)+1.5*movSpeedPx;
                    elseif(all(colorParticle==OilColHolesP*255)) % = [1 0 0]; red, holes piston moves left fast
                        ctrs(iCluster,2,iFrame) = ctrs(iCluster,2,iFrame)-3.5*movSpeedPx; % 3 3.2
                    elseif(all(colorParticle==OilColBase*255))   % = [1 1 0]; yellow, outer moves left slow
                        ctrs(iCluster,2,iFrame) = ctrs(iCluster,2,iFrame)+0.4*movSpeedPx; % 0.3333
                        
                    elseif(all(colorParticle==OilColBaseUp*255)) % = [1 0 1]; magenta, vertical movement down
                        ctrs(iCluster,1,iFrame) = ctrs(iCluster,1,iFrame)-movSpeedPx*0.8;
                    elseif(all(colorParticle==OilColBaseDo*255)) % = [0 1 1]; cyan, vertical movement up
                        ctrs(iCluster,1,iFrame) = ctrs(iCluster,1,iFrame)+movSpeedPx*0.8;
                        
                    end
                end

                [idx,ctrs(:,:,iFrame)] = kmeans([row,col],nParticlesOil,'Distance','sqEuclidean','start',ctrs(:,:,iFrame),...
                    'Options',optsIter);
            end
        end
        
        %% filter paths
        
        % figure(2376574); clf; hold on
        % xlim([1 xSize]/oilRed); ylim([1 ySize]/oilRed); % set axis limits
        wind = 4;
        for iCluster = 1:nParticlesOil % for all cluster-enters (particles)
            %     cla
            %     plot(permute(ctrs(iCluster,2,:),[3,2,1]),permute(ctrs(iCluster,1,:),[3,2,1]),'b-');
            ctrs(iCluster,1,:) = filtfilt(ones(1,wind)/wind,1,ctrs(iCluster,1,:)); % filter particle path x-coordinate
            ctrs(iCluster,2,:) = filtfilt(ones(1,wind)/wind,1,ctrs(iCluster,2,:)); % filter particle path y-coordinate
            %     plot(permute(ctrs(iCluster,2,:),[3,2,1]),permute(ctrs(iCluster,1,:),[3,2,1]),'m-');
            %     pause()
        end

        switch form{1} %
            case 'gif'
                
                reducedRGBimage = ones([ySize,xSize,3,nFrames],'uint8'); % allocate
                
            case 'mp4'
                vid = VideoWriter(mode{1},'MPEG-4'); % Prepare the new file
                vid.FrameRate = 30;
                vid.Quality = 100; % quality in [0% 100%]
                open(vid);
        end
        
        gasColMap = cool(100);
        minOilPos = min(oilPos);
        maxOilPos = max(oilPos);        

        if strcmp(mode{1},'Twin_Tube_Damper')
            gasOffset_x = 0.75*((rand(nParticlesGas/2,2).*wall)-wall/2);  % gas positions distribution.
            gasOffset_x(1,:)   = 0; % avoid the leftmost particle get over the line
            gasOffset_x(end,:) = 0; % avoid the rightmost particle get over the line
            gasOffset_x = gasOffset_x(:);
            gasOffset_y = (rand(nParticlesGas,1).*wall*0.2)+0.09; % gas positions distribution.
            gasOffset_y = gasOffset_y .* (2*rem( (1:nParticlesGas)', 2) - 1); % alternate 1 -1 1 -1 offset
            colorShading = [ones(nParticlesGas/2,1)*0.5; ones(nParticlesGas/2,1)*0.7];

            gasPos_x = linspace(0, 1, (nParticlesGas)+1)';
            gasPos_x = [gasPos_x(2:2:end); gasPos_x(2:2:end)];
            gasPos_y = mean(outerRad-wall, innerRad)* ones(nParticlesGas/2,1);
            gasPos_y_scale = [gasPos_y; -gasPos_y];
        else
            XY_gasPos = net(haltonset(2,'Skip',1),nParticlesGas*1.2); % generate quasi random numbers, 'Skip',1 to remove the point at 0,0 
            clusterIndex = clusterdata(XY_gasPos,nParticlesGas); % XY_gasPos contains still some very close neighbours (touching particles)
            [~,representativeOfCluster] = unique(clusterIndex); % get only one member of each cluster
            gasPos_x = XY_gasPos(representativeOfCluster,1); % extract nicely distributed quasi-random random numbers
            gasPos_y = XY_gasPos(representativeOfCluster,2); % extract nicely distributed quasi-random random numbers
            colorShading = interp1([0 1],[0.7 0.5],gasPos_y);
            gasPos_y_scale = (gasPos_y-0.5)*2*(innerRad-wall/2); % scale to full radius
            gasOffset_x = (rand(nParticlesGas,1)-0.5).*wall*0.3; % 
            gasOffset_y = (rand(nParticlesGas,1)-0.5).*wall*0.3; % 
        end
        
        %%  plot loop
        for iFrame = 1:nFrames
            
            figure(figHandle); cla
            
            gasTempPercent = round((oilPos(iFrame)-minOilPos)/(maxOilPos-minOilPos)*99)+1;
            
            if strcmp(mode{1},'Twin_Tube_Damper')
                plotBox([l_outer-2*wall l_outer], [-1 1]*(outerRad+wall), RGB.brightGrey, RGB.edge ,liSc) % outer body left, bottom
                
                plotCylinder([l_outer oilPos(iFrame)], [-1 1]*outerRad, RGB.oil,   RGB.edge, liSc, false) % outer body inner void (oil)
                plotCylinder([oilPos(iFrame) r_outer], [-1 1]*outerRad, gasColMap(gasTempPercent,:), RGB.edge, liSc, false) % outer body inner void (air)
                
                xLid = [...
                    r_outer r_outer+2*wall;... box1
                    r_inner r_inner+3*wall]; %  box2, both boxes will be unified into one object
                yLid =  [...
                    [-1 1]*(outerRad+wall);... box1
                    [-1 1]*innerRad]; %  box2, both boxes will be unified into one object
                plotBox(xLid, yLid, RGB.brightGrey, RGB.edge ,liSc) % outer body right, lid
                
                plotBox([l_outer-wall r_outer+wall],  [outerRad outerRad+wall], RGB.brightGrey, RGB.edge ,liSc) % outer body upper wall
                plotBox([l_outer-wall r_outer+wall], -[outerRad outerRad+wall], RGB.brightGrey, RGB.edge ,liSc) % outer body lower wall
                
                plotCylinder([l_inner r_inner], [-1 1]*innerRad,RGB.oil ,RGB.edge, liSc, false) % inner body oil filling
            else % Mono_Tube_Damper
                plotCylinder([AirOilSeparatorPos(iFrame,2) r_inner], [-1 1]*innerRad, RGB.oil ,RGB.edge, liSc, false) % inner body oil filling
                plotBox(AirOilSeparatorPos(iFrame,:),[-1 1]*innerRad, RGB.blueGrey, RGB.edge, liSc) % 2nd piston between oil and air
                plotCylinder([l_inner AirOilSeparatorPos(iFrame,1)], [-1 1]*innerRad, gasColMap(gasTempPercent,:) ,RGB.edge, liSc, false) % inner body air filling
                
                plotBox(AirOilSeparatorPos(iFrame,:)+[wall -wall], [0 -wall/2]+innerRad, [0.1 0.1 0.1], RGB.edge, liSc) % flat piston seal
                plotBox(AirOilSeparatorPos(iFrame,:)+[wall -wall], [0 +wall/2]-innerRad, [0.1 0.1 0.1], RGB.edge, liSc) % flat piston seal

                plotBox([r_inner r_outer+3*wall], [-1 1]*(innerRad+wall), RGB.brightGrey, RGB.edge ,liSc) % outer body right, lid
                
            end
            
            plotBox([l_inner-2*wall l_inner], [-1 1]*(innerRad+wall),RGB.brightGrey, RGB.edge ,liSc) % inner body left, bottom
            
            plotBox([l_inner-wall r_inner+3*wall],  [innerRad innerRad+wall], RGB.brightGrey, RGB.edge ,liSc) % inner body upper wall
            plotBox([l_inner-wall r_inner+3*wall], -[innerRad innerRad+wall], RGB.brightGrey, RGB.edge ,liSc) % inner body upper wall
            
            if strcmp(mode{1},'Twin_Tube_Damper')
                plotCylinder([l_inner-2*wall l_inner], +holeY, RGB.oil ,RGB.edge ,liSc, false) % hole in inner lid
                plotCylinder([l_inner-2*wall l_inner], -holeY, RGB.oil ,RGB.edge ,liSc, false) % hole in inner lid
            end
            
            plotSeal([1 2]*wall+r_outer, [0  wall]+rodRad, RGB.brightGrey, RGB.edge, liSc, '\') % plot upper wiper seal
            plotSeal([1 2]*wall+r_outer, [0 -wall]-rodRad, RGB.brightGrey, RGB.edge, liSc, '/') % plot lower wiper seal
            
            plotSeal([1 2]*wall+r_inner, [0  wall]+rodRad, RGB.brightGrey, RGB.edge, liSc, 'o') % plot o-Ring
            plotSeal([1 2]*wall+r_inner, [0 -wall]-rodRad, RGB.brightGrey, RGB.edge, liSc, 'o') % plot o-Ring
            
            plotCylinder(pistonPos(iFrame,2)+[0 rodLenght],[-1 1]*rodRad, RGB.blueGrey, RGB.edge, liSc, true) % rod
            plotBox(pistonPos(iFrame,:),[-1 1]*innerRad, RGB.blueGrey, RGB.edge, liSc) % piston
            
            plotCylinder(pistonPos(iFrame,:),+holeY, RGB.oil, RGB.edge, liSc, false) % hole in piston
            plotCylinder(pistonPos(iFrame,:),-holeY, RGB.oil, RGB.edge, liSc, false) % hole in piston
            
            plotBox(pistonPos(iFrame,:)+[wall -wall], [0 -wall/2]+innerRad, [0.1 0.1 0.1], RGB.edge, liSc) % flat piston seal
            plotBox(pistonPos(iFrame,:)+[wall -wall], [0 +wall/2]-innerRad, [0.1 0.1 0.1], RGB.edge, liSc) % flat piston seal
            
            figure(figHandle)
            
            scatter(... % plot oil particles
                (ctrs(:,2,iFrame)*oilRed/scaleReduction-1)*(xRange/xSize)+xLimits(1),... % x-pos
                (ctrs(:,1,iFrame)*oilRed/scaleReduction-1)*(yRange/ySize)+yLimits(1),... % y-pos
                (liSc^2)*10, 'o',...  % size and shape
                'LineWidth',  liSc,...
                'MarkerEdgeColor', RGB.oil.*0.6,...
                'MarkerFaceColor', RGB.oil.*0.9,...
                'MarkerEdgeAlpha', 0.5,...
                'MarkerFaceAlpha', 0.5)
            
            if strcmp(mode{1},'Twin_Tube_Damper')
                gasPos_x_scale = interp1([0 1],[oilPos(iFrame)+wall/4, r_outer-wall/4],gasPos_x);
            else % Mono_Tube_Damper
                gasPos_x_scale = interp1([0 1],[l_inner+wall/2, AirOilSeparatorPos(iFrame,1)-wall/2],gasPos_x);
            end
            
            scatter(... % plot gas particles
                gasPos_x_scale + gasOffset_x*gasTempPercent/100 + randn(nParticlesGas,1)*0.05*gasTempPercent/100,... % x-pos
                gasPos_y_scale + gasOffset_y*gasTempPercent/100 + randn(nParticlesGas,1)*0.05*gasTempPercent/100,... % y-pos
                200 + (pistonPos(iFrame,2)^2 * (liSc^2)*0.8),...
                colorShading * gasColMap(gasTempPercent,:),...
                '.')
            
            %% save animation
            xlim(xLimits); ylim(yLimits); % set axis limits
            axis equal; drawnow;
            f = getframe(figHandle);
            
            switch form{1} % Init frame(s)
                case {'mp4'}
                    writeVideo(vid, rot90(f.cdata)); % save video frames
                case {'gif'}
                    reducedRGBimage(:,:,:,iFrame) = imReduceSize(f.cdata,scaleReduction); % the size reduction: adds antialiasing
            end
            
            if iFrame == nFrames
                plot((col*oilRed/scaleReduction-1)*(xRange/xSize)+xLimits(1), (row*oilRed/scaleReduction-1)*(yRange/ySize)+yLimits(1),'.b') % check points for clustering (debug)
            end
            
        end % plot loop
        
        switch form{1} %
            case 'gif'
                
                colormapImage = permute(reducedRGBimage(:,:,:,:),[1 2 4 3]); % step1; make unified image
                colormapImage = reshape(colormapImage,[ySize,xSize*(nFrames) 3 1]); % step2; make unified image
                
                map = createImMap(colormapImage,256,[RGB.edge;RGB.white;RGB.brightGrey;RGB.blueGrey;RGB.oil;min(colorShading)*gasColMap(end,:)]); % colormap
                
                for iFrame = 1:nFrames
                    im = rgb2ind(reducedRGBimage(:,:,:,iFrame),map,'dither');
                    
                    if iFrame == 1
                        imwrite(rot90(im),map,fullfile(pathstr, [mode{1} '.gif']), 'LoopCount',Inf,'DelayTime',delayTime(iFrame)); % individual timings
                    else
                        imwrite(rot90(im),map,fullfile(pathstr, [mode{1} '.gif']), 'WriteMode','append','DelayTime',delayTime(iFrame)); % individual timings
                    end
                end
                
                disp([mode{1} '.gif  has ' num2str(numel(reducedRGBimage)/3/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 100 MP limit
                
            case 'mp4'
                close(vid ); disp(['MP4 completed: ' mode{1} '.mp4']);
                disp(['Conversion of ' mode{1} ' from MP4 to WebM (Matlab currently does not support webm)']);
                if exist('ffmpeg.exe','file') == 2 % check if it is a valid path to an existing file
                    try
                        dosCommand = ['!ffmpeg -i ' mode{1} '.mp4 '...   % Command and source
                            '-c:v libvpx-vp9 -crf 30 -b:v 0 -deadline best '... % constant (best) quality; see: https://trac.ffmpeg.org/wiki/Encode/VP9
                            '-metadata title="' mode{1} '" '... % see: https://wiki.multimedia.cx/index.php/FFmpeg_Metadata
                            '-metadata author="Jahobr" -metadata copyright="CC0 1.0 Universal Public Domain Dedication" '...
                            '-metadata license="CC0 1.0 Universal Public Domain Dedication" '...
                            mode{1} '.webm'];                            % conversion target
                        eval(dosCommand); % run conversion from mp4 to webm
                        % delete(fullfile(pathstr,[mode{1} '.mp4'])); % delete mp4 version after conversion (to save disc space)
                    catch ME
                        disp(ME)
                    end
                else
                    disp(['"ffmpeg.exe" not available in path "' pathstr '" (or other problem); see https://ffmpeg.zeranoe.com/builds/']);
                end
        end
    end
end
  

function disc(x,y,r,colFa)
% x coordinates of the center
% y coordinates of the center
% r is the radius of the circle
% colFa face color   [r g b]

angleOffPoints = linspace(0,2*pi,300);
xc = x + r*sin(angleOffPoints);
yc = y + r*cos(angleOffPoints);

patch(xc,yc,'k','LineWidth',1,'EdgeColor','none','FaceColor',colFa);


function plotSeal(x,y,colFa,colEd,linw,type)
% x [left right]
% y [top bottom]
% colFa  face color  [r g b]
% colEd  edge color  [r g b]
% linw  line width

x = sort(x); y = sort(y); % just nor make sure
plotBox(x, y, colFa, colEd, linw)
switch type
    case {'o','O'}
        disc(mean(x), mean(y), (y(2)-y(1))/2.3, [0.1 0.1 0.1])
    case '/'
        plot(x,           y, 'Color', colEd, 'LineWidth', linw*1.5)
    case '\'
        plot([x(2) x(1)], y, 'Color', colEd, 'LineWidth', linw*1.5)
end


function plotBox(x,y,colFa,colEd,linw)
% x [left right]
% y [top bottom]
%  Or 
% x [left right;... box1
%    left right] %  box2, both boxes will be unified int one object
% y [left right;... box1
%    left right] %  box2, both boxes will be unified int one object
% colFa  face color  [r g b]
% colEd  edge color  [r g b]
% linw  line width

xs = [ ]; % x
ys = [ ]; % y

for iBox = 1:size(x,1)
    xAdd = [x(iBox,2)       x(iBox,2) x(iBox,1) x(iBox,1) x(iBox,2)       x(iBox,2)]; % x
    yAdd = [mean(y(iBox,:)) y(iBox,1) y(iBox,1) y(iBox,2) y(iBox,2) mean(y(iBox,:))]; % y joint in the middle of an edge to get nice corners
    [xs,ys] = polybool('union', xs, ys, xAdd, yAdd);
end

if strcmp(colFa,'none')
   plot(xs,ys,'Color',colEd,'LineWidth',linw) %
else
   patch(xs,ys,colFa,'EdgeColor',colEd,'LineWidth',linw) %
end


function plotCylinder(x,y,colFa,colEd,linw,material_or_hole)
% x [left right]
% y [top bottom]
% colFa  face color  [r g b]
% colEd  edge color  [r g b]
% linw  line width
% material_or_hole:    1 Material, Highlight in top
%                      0 Hole, Highlight in bottom
x = sort(x);
y = sort(y);
if material_or_hole
   y = flip(y);
end
y_center = mean(y);
y_highlight = mean([y(1) y_center]);

vertex_matrix = [ % clockwise                                8----1
    x(1)  y(1);        % 1 right top                         ¦    ¦
    x(1)  y_highlight; % 2 right upper highlight             7----2
    x(1)  y_center;    % 3 right center                      ¦    ¦
    x(1)  y(2);        % 4 right bottom                      6----3
    x(2)  y(2);        % 5 left  bottom                      ¦    ¦
    x(2)  y_center;    % 6 left  center                      ¦    ¦
    x(2)  y_highlight; % 7 left  upper highlight             ¦    ¦
    x(2)  y(1)];       % 8 left  top                         5----4


vertex_color = [ % color matrix
    colFa;           % 1 top 
    1-(1-colFa)*0.6; % 2 upper highlight (brighter)
    colFa;           % 3 center 
    colFa*0.8];      % 4 bottom (darker)

vertex_color = [vertex_color; flipud(vertex_color)]; % colors 5 to 8 in inverse order

faces_matrix = [ % top to bottom
    1 2 7;  % top to highlight,    upper right triangle
    1 7 8;  % top to highlight,    lower left triangle
    2 3 6;  % highlight to center, upper right triangle
    2 6 7;  % highlight to center, lower left triangle
    3 4 5;  % center to bottom,    upper right triangle
    3 5 6]; % center to bottom,    lower left triangle

patch('Vertices',vertex_matrix,'Faces',faces_matrix,...
      'FaceVertexCData',vertex_color,'FaceColor','interp',...
      'Edgecolor','none');
plotBox(x,y,'none',colEd,linw)


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


function im = randSubsampling(im, imRed)
% Input:
%  Subsampling image. Colors are not mixed, no new colors are created.
%  im:      image, [imRows x imColumns x nChannel] (unit8)
%                      imRows, imColumns: must be divisible by redSize
%                      nChannel: usually 3 (RGB) or 1 (grey)
%  redSize: 2 = half the size (quarter of pixels)
%           3 = third the size (ninth of pixels)
%           ... and so on
%  Example: imRed=3 results in use of one      010 000 000 000 100 000 000
%           random specimen chosen from        000 000 001 010 000 000 001
%           the 3x3-origal pixel block         000 100 000 000 000 100 000   

[nRow,nCol,nChannel] = size(im);

index = reshape(1:nRow*nCol, nRow, nCol); % build index matrix using linear indexing
index = index(1:imRed:end,1:imRed:end);   % reduce size (this uses effectively the upper left corner of each sub-box)
[nRowTarget, nColTarget] = size(index);   % new size

randRowOffset =  randi([0 imRed-1],nRowTarget,nColTarget); % rand row-offset within pixel block
indexCheckLimRow = nRowTarget*imRed+randRowOffset(end,:) > nRow; % find creation of index out of original image size
randRowOffset(end,indexCheckLimRow) = 0; % avoid creation of index out of original image size

randColOffset =  randi([0 imRed-1],nRowTarget,nColTarget); % rand col-offset within pixel block
indexCheckLimCol = nColTarget*imRed+randColOffset(:,end) > nCol; % find creation of index out of original image size
randColOffset(indexCheckLimCol,end) = 0; % avoid creation of index out of original image size

index = index + randRowOffset + randColOffset*nRow; % create linear indexing that chooses random pixel in each pixel block

if nChannel == 3
    index = repmat(index,1,1,3);        % for each color channel
    index(:,:,2) = index(:,:,2) + nRow*nCol;   % for each color channel
    index(:,:,3) = index(:,:,3) + nRow*nCol*2; % for each color channel
end

im = im(index);