function f = devils(k);
% The function calculates the Devil's Staircase
% for a value of the parameter K > 0 real
% Call: devils(K)
% The program draws the staircase graph;
% Original by J Neubauer / Humboldt U Berlin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameter settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
omegamin = 0.0;
omegamax = 1.0;
nomegastep = 1000;
trans = 400; % number of transient points
out = 400; % number of plotted points
x0 = 0.3; % starting values x_0 for iterations for all parameters
omegastep = (omegamax-omegamin)/(nomegastep-1);
omega = (omegamin:omegastep:omegamax);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = x0 * ones( 1, length(omega) ); % row vector of IVs of length r
xstart = x;
for i = 1:trans % transient iterations for all selected parameters
x = x + omega - k .* sin(2.0 * pi * x) / 2.0 / pi; % circle map
end;
xin = x;
for i = 1:out % output iterations
x = x + omega - k .* sin(2.0 * pi * x) / 2.0 / pi; % circle map
w = (x - xin)/i; % winding number
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(gcf), plot(omega,w,'b .','markersize',3);
xlabel('Bifurcation parameter');
ylabel('Winding number');
title('Devil''s staircase');
grid on;
set(gca,'xlim',[omegamin omegamax]);
set(gca,'ylim',[-.10 1.1]);
%
uicontrol('Style','pushbutton','Units','normalized',...
'Position',[.75 .01 .15 .05],'String','CLOSE',...
'BackGroundColor',[.0 .0 .9],'ForeGroundColor','w',...
'Fontsize',12,'Callback','uiresume;');
uiwait;
close;