clear
figure(1)
clf
a=input(' Give value of a: ');
nmax=input(' Give number of iterations per initial point: ');
np=input(' Give peppering density of initial points: ');
twopi=2*pi;
delta=twopi/np;
x=0:delta:np*delta;
y=x;
[xx,yy]=meshgrid(x,y);
plot(xx,yy,'.','markersize',2)
hold on
for i=1:nmax
  xxx=mod(xx+yy,twopi);
  yy=mod(yy+a*sin(xx+yy),twopi);
  xx=xxx;
  plot(xx,yy,'.','markersize',2)
end
axis([0 twopi 0 twopi]);
axis square