Plot of a butterfly using mathematical equations

Hello, I saw these set of equations that plot a butterfly using circles (See Image Below). I was able to code it down and the plot seems to be okay. However, there are a few issues with my code and it would be great if anybody could help.

My questions are:

1) Is the code correct?

2) Is there a faster method of plotting these circles. The current code takes quite a while for it to plot.

3) How can i get a color scheme based on the picture above?

The code that I attempted was as follows:

clear; close all; clc;

sin1  =@(k) sin((pi*k)./40000);
cos1  =@(k) cos((pi*k)./40000);
sin141=@(k) sin((141*pi*k)./40000);
cos141=@(k) cos((141*pi*k)./40000);
sin2  =@(k) sin((2*pi*k)./40000);
cos2  =@(k) cos((2*pi*k)./40000);
cos32 =@(k) cos((32*pi*k)./40000);
sin6  =@(k) sin((6*pi*k)./40000);
sin18 =@(k) sin((18*pi*k)./40000);
cos3  =@(k) cos((3*pi*k)./40000);
cos21 =@(k) cos((21*pi*k)./40000);
cos12 =@(k) cos((12*pi*k)./40000);

x = 1:40000;

Xk = (3/2)*cos141(x).^9 .* (1-0.5*sin1(x)) .* (1 - 0.25*cos2(x).^30 .* ...
    (1+cos32(x).^20)) .* (1 - 0.5*sin2(x).^30 .* sin6(x).^10 .* ...
    (0.5 + 0.5*sin18(x).^20));
Yk = cos2(x) .* cos141(x).^2 .* (1+0.25*(cos1(x).*cos3(x).*cos21(x)).^24);
Rk = 0.025 .* (cos141(x).^14 + sin141(x).^6) .* (1- (cos1(x).*cos3(x).*...
    cos12(x)).^16) + 0.01;

axis([-1.5 1.5 -1 1.5])
for i = 1:length(x)
    viscircles([Xk(i), Yk(i)], Rk(i),'Color',[0.7 0.7 0.7]);
    drawnow
    pause(10^-10000)
end

Your help would be appreciated.

Thank you

Answer 1

[KSSV]

Don't use loop...instead use scatter. Read about it.

clc; clear all ;
clear; close all; clc;

sin1  =@(k) sin((pi*k)./40000);
cos1  =@(k) cos((pi*k)./40000);
sin141=@(k) sin((141*pi*k)./40000);
cos141=@(k) cos((141*pi*k)./40000);
sin2  =@(k) sin((2*pi*k)./40000);
cos2  =@(k) cos((2*pi*k)./40000);
cos32 =@(k) cos((32*pi*k)./40000);
sin6  =@(k) sin((6*pi*k)./40000);
sin18 =@(k) sin((18*pi*k)./40000);
cos3  =@(k) cos((3*pi*k)./40000);
cos21 =@(k) cos((21*pi*k)./40000);
cos12 =@(k) cos((12*pi*k)./40000);

x = 1:40000;

Xk = (3/2)*cos141(x).^9 .* (1-0.5*sin1(x)) .* (1 - 0.25*cos2(x).^30 .* ...
    (1+cos32(x).^20)) .* (1 - 0.5*sin2(x).^30 .* sin6(x).^10 .* ...
    (0.5 + 0.5*sin18(x).^20));
Yk = cos2(x) .* cos141(x).^2 .* (1+0.25*(cos1(x).*cos3(x).*cos21(x)).^24);
Rk = 0.025 .* (cos141(x).^14 + sin141(x).^6) .* (1- (cos1(x).*cos3(x).*...
    cos12(x)).^16) + 0.01;

scatter(Xk,Yk,[],Rk)
colormap(jet)
axis([-1.5 1.5 -1 1.5])

[Ahmed Mohamed Mansoor]

Thank you so much!. This was very helpful.

This is an additional question out of curiosity. If there was a way to do a quick animation using either the 'drawnow' or 'animated line' function on MATLAB (or any other way). Would that be possible?

[KSSV]

REad about comet

[Ahmed Mohamed Mansoor]

Thank you again.

I came up with this. Is there a way to speed this up?

for i =  1:length(x)
    comet(Xk(i),Yk(i),Rk(i))
    colormap(copper)
    axis([-1.5 1.5 -1 1.5])
    hold on
    pause(0.000005)
end