function HarmonicDriveAni() % source code for drawing a HarmonicDrive % this is by no means a "simulation". It is a hack job that produces a gif % % 2016-12-05 Jahobr (reworked 2017-09-16)     nTeethOutGear =  42; nTeethFlex = nTeethOutGear-2; modul = 0.1; % modul   colEdge = [0   0   0  ]; % Edge color colFlex = [1   0.2 0.2]; % FlexSpline color colWave = [0.1 0.7 0.1]; % WaveGen color colGear = [0.2 0.2 1  ]; % static OuterGear color   nFrames = 100; frameAngles = linspace(0,-pi,nFrames+1); % rotate clockwise frameAngles = frameAngles(1:end-1); % delete redundant frame   [pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location   figHandle = figure(15674454); clf axesHandle = axes; hold(axesHandle,'on') axis equal xlim([-3 3]) ylim([-3 3]) set(figHandle, 'Units','pixel'); set(figHandle, 'position',[1 1 700 700]); % [x y width height]  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]  set(figHandle,'GraphicsSmoothing','on') % requires at least version 2014b   for iFrame = 1:nFrames          angleWaveGen = frameAngles(iFrame);     angleFlexTeeth = angleWaveGen*(nTeethFlex-nTeethOutGear)/nTeethFlex; %  angle of the flexspline          cla(axesHandle);          %% %%%%%%%%    draw OuterGear (static)   %%%%%%%%     %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          effectiveDiameter = modul*nTeethOutGear;     toothTipDiameter = effectiveDiameter-1.4*modul;     toothBottomDiameter = effectiveDiameter+1.6*modul;          angleBetweenTeeth = 2*pi/nTeethOutGear; % angle between 2 teeth     angleOffPoints = (0:angleBetweenTeeth/8:(2*pi));          %% outerEdge     maxDiameter = toothBottomDiameter*1.2; % definition of outer line     maxXY = samplesEllipse(maxDiameter,maxDiameter,500);     patch(maxXY(:,1),maxXY(:,2),colGear,'EdgeColor',colEdge,'LineWidth',0.5) % full outer disc          %% inner teeth     radiusOffPoints = angleOffPoints; % init          radiusOffPoints(1:8:end) = toothBottomDiameter/2; % middle bottom     radiusOffPoints(2:8:end) = toothBottomDiameter/2; % left bottom     radiusOffPoints(3:8:end) = effectiveDiameter/2; % rising edge     radiusOffPoints(4:8:end) = toothTipDiameter/2; % right top     radiusOffPoints(5:8:end) = toothTipDiameter/2; % middle top     radiusOffPoints(6:8:end) = toothTipDiameter/2; % left top     radiusOffPoints(7:8:end) = effectiveDiameter/2; % falling edge     radiusOffPoints(8:8:end) = toothBottomDiameter/2; % right bottom          [X,Y] = pol2cart(angleOffPoints,radiusOffPoints);          patch(X,Y,[1 1 1],'EdgeColor',colEdge,'LineWidth',0.5) % overlay white area for inner teeth               %% %%%%%%%%%     draw Flexspline        %%%%%%%%%     %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          % % deform estimation based on tooth distance (using the circumferences); could be automated!     % U1 = 42*pi % Circumference of OuterGear     % U1 =     %   131.9469     %     % U2 = pi*sqrt(2*((42/2)^2+(0.9022*42/2)^2))  * 42/40 % Circumference of Flexspline * 42/40     % U2 =     %   131.9435     deform = 0.9022;          deformedDiameter = effectiveDiameter*deform; % scale down, but teeth must still have the same distance          rootEffectiveDia = effectiveDiameter-1.6*modul; % fixed offset     rootDeformedDia  = deformedDiameter-1.6*modul;  % fixed offset          topEffectiveDia = effectiveDiameter+1.4*modul; % fixed offset     topDeformedDia  = deformedDiameter+1.4*modul;  % fixed offset          % % an equidistant sampled ellipse is needed, to keep the tooth distance constant all the way around     offsetOnCircumference = (-angleWaveGen+angleFlexTeeth)/2/pi; % compensation + own_rotation  ,  normalization to "circumference"     equiEffeXY = equidistantSamplesEllipse(effectiveDiameter,deformedDiameter,nTeethFlex*8, offsetOnCircumference); % points on effective diameter     equiRootXY = equidistantSamplesEllipse(rootEffectiveDia, rootDeformedDia, nTeethFlex*8, offsetOnCircumference); % points with inwards offset     equiOutXY  = equidistantSamplesEllipse(topEffectiveDia,  topDeformedDia,  nTeethFlex*8, offsetOnCircumference); % points with outwards offset          toothXY = equiEffeXY; % intit          toothXY(1:8:end,:) = equiOutXY(1:8:end,:); % middle top        I######I     toothXY(2:8:end,:) = equiOutXY(2:8:end,:); % left top          I######+     % toothXY(3:8:end) init did it                                 I####/     toothXY(4:8:end,:) = equiRootXY(4:8:end,:); % right bottom     I##+     toothXY(5:8:end,:) = equiRootXY(5:8:end,:); % middle bottom    I##I     toothXY(6:8:end,:) = equiRootXY(6:8:end,:); % left bottom      I##+     % toothXY(7:8:end) init did it                                 I####\     toothXY(8:8:end,:) = equiOutXY(8:8:end,:); % right top         I######+          [toothXY] = rotateCordiantes(toothXY,angleWaveGen);          patch(toothXY(:,1),toothXY(:,2),colFlex,'EdgeColor',colEdge,'LineWidth',0.5) %draw flexspline with teeth               %% hole     holeEffectiveDia = effectiveDiameter-5*modul; % fixed inwards offset     holeDeformedDia  = deformedDiameter-5*modul;  % fixed inwards offset          holePathXY = samplesEllipse(holeEffectiveDia,holeDeformedDia,500);     holePathXY = rotateCordiantes(holePathXY,angleWaveGen);     patch(holePathXY(:,1),holePathXY(:,2),[1 1 1],'EdgeColor',colEdge,'LineWidth',0.5) % draw hole of deformed ring               %% %%%%%%%%%   draw wave generator      %%%%%%%%%     %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          waveEffectiveDia = holeEffectiveDia; % touch flex spline     waveDeformedDia  = holeDeformedDia-5*modul; % extra air gap to spline, to make it more obvious          wavePathXY = samplesEllipse(waveEffectiveDia,waveDeformedDia,500);     [wavePathXY] = rotateCordiantes(wavePathXY,angleWaveGen);     patch(wavePathXY(:,1),wavePathXY(:,2),colWave,'EdgeColor',colEdge,'LineWidth',0.5) % draw wave generator          %% central shaft     shaftPathXY = samplesEllipse(effectiveDiameter/2.5,effectiveDiameter/2.5,500);     plot(axesHandle,shaftPathXY(:,1),shaftPathXY(:,2),'LineWidth',0.8,'color',colEdge); % draw central shaft outline               %% %%%%%%%%%   save animation     %%%%%%%%%     %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          drawnow;          f = getframe(figHandle);     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                  if ~isempty(which('plot2svg'))             plot2svg(fullfile(pathstr, [fname '_Frame1.svg']),figHandle) % by Juerg Schwizer         else             disp('plot2svg.m not available; see http://www.zhinst.com/blogs/schwizer/');         end     end          imtemp = rgb2ind(f.cdata,map,'nodither');     im(:,:,1,iFrame) = imtemp;      end imwrite(im,map,fullfile(pathstr, [fname '.gif']),'DelayTime',1/30,'LoopCount',inf) % save gif disp([fname '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit   %%% equidistantSamplesEllipse test code % figure(455467);clf;hold on; %  % equidistantXY = equidistantSamplesEllipse(1.5,0.5,40,0.1); % plot(equidistantXY(:,1),equidistantXY(:,2),'bx-') %  %  % equidistantXY = equidistantSamplesEllipse(2,1,40,1); % plot(equidistantXY(:,1),equidistantXY(:,2),'bx-') %  % equidistantXY = equidistantSamplesEllipse(3,2,40,0.5); % plot(equidistantXY(:,1),equidistantXY(:,2),'bx-') %  % equidistantXY = equidistantSamplesEllipse(4,3,40,0); % plot(equidistantXY(:,1),equidistantXY(:,2),'bx-') % pathXY = samplesEllipse(4,3,41); % plot(pathXY(1:end-1,1),pathXY(1:end-1,2),'ro-') %  % plot([4 -4]/2,[0 0],'-k')     function equidistantXY = equidistantSamplesEllipse(diameterH,diameterV,nPoints,offset) % Inputs: %   diameterH  horizontal diameter %   diameterV  vertical diameter %   nPoints    number of resampled points %   offsetFraction between 0 and 1 in circumference of ellipse   pathXY = samplesEllipse(diameterH,diameterV,1000); % create ellipse stepLengths = sqrt(sum(diff(pathXY,[],1).^2,2)); % distance between the points stepLengths = [0; stepLengths]; % add the starting point cumulativeLen = cumsum(stepLengths); % cumulative sum circumference = cumulativeLen(end); finalStepLocs = linspace(0,1, nPoints+1)+offset; % equidistant distribution finalStepLocs = finalStepLocs(1:end-1); % remove redundant point finalStepLocs = mod(finalStepLocs,1)*circumference; % unwrap and scale to circumference equidistantXY = interp1(cumulativeLen, pathXY, finalStepLocs);     function pathXY = samplesEllipse(diameterH,diameterV,nPoints) % point of ellipse; points start on the right, counterclockwise  % first and last points are the same % % Inputs: %   diameterH  horizontal diameter %   diameterV  vertical diameter %   nPoints    number of points   p = linspace(0,2*pi,nPoints)'; pathXY = [cos(p)*diameterH/2 sin(p)*diameterV/2]; % create ellipse     function [xy] = rotateCordiantes(xy,anglee) % [x1 y1; x2 y2; x3 y3; ...] coordinates to rotate % anglee angle of rotation in [rad] rotM = [cos(anglee) -sin(anglee); sin(anglee) cos(anglee)]; xy = (rotM*xy')';