Plot elipses with a foci based on equations to plot flowers

I am trying to practice with matlab with various equations out there. Does anyone know how to plot this picture (see image below) based on these equations? I have coded the equations (see below picture), but have no idea how to plot it.

I have coded the equations A(k), B(k), C(k), D(k), E(k) but I do not know where to go from here.

Any help would be appreciated.

close all; clc; clear;

sin48 =@(k) sin((48*pi*k)./1000000);
cos48 =@(k) cos((48*pi*k)./1000000);
sin8  =@(k) sin((8*pi*k)./1000000);
cos8  =@(k) cos((8*pi*k)./1000000);
sin24 =@(k) sin((24*pi*k)./1000000);
cos24 =@(k) cos((24*pi*k)./1000000);
sin72 =@(k) sin((72*pi*k)./1000000);
cos72 =@(k) cos((72*pi*k)./1000000);
cos648=@(k) cos((648*pi*k)./1000000);
cos64 =@(k) cos((64*pi*k)./1000000);

Ek = 0.1 + 0.25*sin48(x).^2 + 0.2*cos18(x).^18 *...
    (1 + (1/3)*cos24(x).^8  + 0.1*cos24(x)^20) + (1/6)*sin8(x).^16 +...
    (8/20)*sin8(x).*(1 - 0.4*cos72(x).^4) + 0.0625*sin(x).^18.*...
    sin24(x).^20.*sin72(x).^30 - (0.1*sin24(x).^6 + (1/6)*sin24(x).^30)...
    (1-sin8(x));

Dk = (pi/2) + (1686*pi*x/1000000);

Ck = 0.003 + 0.097*cos648(x).^40;

Bk = -cos64(x).*Ek - 0.1*cos64(x);

Ak = 1.2*sin64(x)*Ek;

[Sam Chak]

Hmm... Didn't you have few examples posted on the last few days?

Can you put the scatter() and colormap(jet) functions in the code?

[Ahmed Mohamed Mansoor]

Yes. However, there are a few issues. These are ellipses with equations involving imaginary numbers.

Answer 1

[Davide Masiello]

The code below does the job drawing the figure, but not the coloring.

I can't run it online (takes too long).

clear,clc
figure(1)
hold on
for k = 1:1000000
plotEllipse(k)
end

function plotEllipse(k)
e       = 98/100;
sin48   = sin((48*pi*k)/1000000);
sin8    = sin((8*pi*k)/1000000);
cos8    = cos((8*pi*k)./1000000);
sin24   = sin((24*pi*k)/1000000);
cos24   = cos((24*pi*k)/1000000);
sin72   = sin((72*pi*k)/1000000);
cos72   = cos((72*pi*k)/1000000);
cos648  = cos((648*pi*k)/1000000);
sin64   = sin((64*pi*k)/1000000);
cos64   = cos((64*pi*k)/1000000);

Ek      = 0.1 + 0.25*sin48^2 + 0.2*cos8^18 *...
            (1 + (1/3)*cos24^8  + 0.1*cos24^20) + (1/6)*sin8^16 +...
            (8/20)*sin8*(1 - 0.4*cos72^4) + 0.0625*sin8^18.*...
            sin24^20*sin72^30 - (0.1*sin24^6 + (1/6)*sin24^30)*...
            (1-sin8);

Dk = (pi/2) + (1686*pi*k/1000000);
Ck = 0.003 + 0.097*cos648^40;
Bk = -cos64*Ek - 0.1*cos64;
Ak = 1.2*sin64*Ek;

F1      = Ak + 1i*Bk + Ck*exp(1i*Dk) ;                                        % Focus 1
F2      = Ak + 1i*Bk - Ck*exp(1i*Dk) ;                                        % Focus 2

centre  = mean([F1,F2]);                                                    % ellipse centre
angle   = atan(abs(imag(F1)-imag(F2))/abs(real(F1)-real(F2)));              % angle between ellipse horizontal axis and x-axis
c       = sqrt((real(F1)-real(F2))^2+(imag(F1)-imag(F2))^2)/2;              % distance from centre to foci
a       = c/e;                                                              % ellipse horizontal axis
b       = a^2-c^2;                                                          % ellipse vertical axis

t       = linspace(0,2*pi,100);                                             % parametric coordinate
x0      = a*cos(t);                                                         % x-coordinates of ellipse if centered with reference axis and horizontal
y0      = b*sin(t);                                                         % y-coordinates of ellipse if centered with reference axis and horizontal
x       = x0*cos(angle)-y0*sin(angle)+real(centre);                         % Absolute x-coordinates (rotation+translation)
y       = y0*cos(angle)+x0*sin(angle)+imag(centre);                         % Absolute y-coordinates (rotation+translation)

plot(x,y)

end

[Davide Masiello]

little example up to k = 10000

clear,clc
figure(1)
hold on
for k = 1:10000
plotEllipse(k)
end

function plotEllipse(k)
e       = 98/100;
sin48   = sin((48*pi*k)/1000000);
sin8    = sin((8*pi*k)/1000000);
cos8    = cos((8*pi*k)./1000000);
sin24   = sin((24*pi*k)/1000000);
cos24   = cos((24*pi*k)/1000000);
sin72   = sin((72*pi*k)/1000000);
cos72   = cos((72*pi*k)/1000000);
cos648  = cos((648*pi*k)/1000000);
sin64   = sin((64*pi*k)/1000000);
cos64   = cos((64*pi*k)/1000000);
Ek      = 0.1 + 0.25*sin48^2 + 0.2*cos8^18 *...
            (1 + (1/3)*cos24^8  + 0.1*cos24^20) + (1/6)*sin8^16 +...
            (8/20)*sin8*(1 - 0.4*cos72^4) + 0.0625*sin8^18.*...
            sin24^20*sin72^30 - (0.1*sin24^6 + (1/6)*sin24^30)*...
            (1-sin8);
Dk = (pi/2) + (1686*pi*k/1000000);
Ck = 0.003 + 0.097*cos648^40;
Bk = -cos64*Ek - 0.1*cos64;
Ak = 1.2*sin64*Ek;
F1      = Ak + 1i*Bk + Ck*exp(1i*Dk) ;                                        % Focus 1
F2      = Ak + 1i*Bk - Ck*exp(1i*Dk) ;                                        % Focus 2
centre  = mean([F1,F2]);                                                    % ellipse centre
angle   = atan(abs(imag(F1)-imag(F2))/abs(real(F1)-real(F2)));              % angle between ellipse horizontal axis and x-axis
c       = sqrt((real(F1)-real(F2))^2+(imag(F1)-imag(F2))^2)/2;              % distance from centre to foci
a       = c/e;                                                              % ellipse horizontal axis
b       = a^2-c^2;                                                          % ellipse vertical axis
t       = linspace(0,2*pi,100);                                             % parametric coordinate
x0      = a*cos(t);                                                         % x-coordinates of ellipse if centered with reference axis and horizontal
y0      = b*sin(t);                                                         % y-coordinates of ellipse if centered with reference axis and horizontal
x       = x0*cos(angle)-y0*sin(angle)+real(centre);                         % Absolute x-coordinates (rotation+translation)
y       = y0*cos(angle)+x0*sin(angle)+imag(centre);                         % Absolute y-coordinates (rotation+translation)
plot(x,y)
end

[Davide Masiello]

This is faster

clear,clc

tic

n       = 50;                                   % number of points for each ellipse
kmax    = 10000;
x       = zeros(n,kmax);
y       = zeros(n,kmax);

for k = 1:kmax
    [x(:,k),y(:,k)] = ellipseCoordinates(k,n);
end

plot(x,y)

toc

function [x,y] = ellipseCoordinates(k,n)
e       = 98/100;
sin48   = sin((48*pi*k)/1000000);
sin8    = sin((8*pi*k)/1000000);
cos8    = cos((8*pi*k)/1000000);
sin24   = sin((24*pi*k)/1000000);
cos24   = cos((24*pi*k)/1000000);
sin72   = sin((72*pi*k)/1000000);
cos72   = cos((72*pi*k)/1000000);
cos648  = cos((648*pi*k)/1000000);
sin64   = sin((64*pi*k)/1000000);
cos64   = cos((64*pi*k)/1000000);

Ek      = 0.1 + 0.25*sin48^2 + 0.2*cos8^18 *...
            (1 + (1/3)*cos24^8  + 0.1*cos24^20) + (1/6)*sin8^16 +...
            (8/20)*sin8*(1 - 0.4*cos72^4) + 0.0625*sin8^18.*...
            sin24^20*sin72^30 - (0.1*sin24^6 + (1/6)*sin24^30)*...
            (1-sin8);

Dk = (pi/2) + (1686*pi*k/1000000);
Ck = 0.003 + 0.097*cos648^40;
Bk = -cos64*Ek - 0.1*cos64;
Ak = 1.2*sin64*Ek;

F1      = Ak + 1i*Bk + Ck*exp(1i*Dk) ;                                        % Focus 1
F2      = Ak + 1i*Bk - Ck*exp(1i*Dk) ;                                        % Focus 2

centre  = mean([F1,F2]);                                                    % ellipse centre
angle   = atan(abs(imag(F1)-imag(F2))/abs(real(F1)-real(F2)));              % angle between ellipse horizontal axis and x-axis
c       = sqrt((real(F1)-real(F2))^2+(imag(F1)-imag(F2))^2)/2;              % distance from centre to foci
a       = c/e;                                                              % ellipse horizontal axis
b       = a^2-c^2;                                                          % ellipse vertical axis
    
t       = linspace(0,2*pi,n);                                               % parametric coordinate
x0      = a*cos(t);                                                         % x-coordinates of ellipse if centered with reference axis and horizontal
y0      = b*sin(t);                                                         % y-coordinates of ellipse if centered with reference axis and horizontal
x       = x0*cos(angle)-y0*sin(angle)+real(centre);                         % Absolute x-coordinates
y       = y0*cos(angle)+x0*sin(angle)+imag(centre);                         % Absolute y-coordinates

x       = x';
y       = y';

end

[Davide Masiello]

Actually, the whole thing can be done without loop

clear,clc

tic
n       = 50;       
k       = 1:1000000;
e       = 98/100; 

Ek      = 0.1 + 0.25*sin((48*pi*k)/1000000).^2 + 0.2*cos((8*pi*k)/1000000).^18.*...
            (1 + (1/3)*cos((24*pi*k)/1000000).^8  + 0.1*cos((24*pi*k)/1000000).^20) + (1/6)*sin((8*pi*k)/1000000).^16 +...
            (8/20)*sin((8*pi*k)/1000000).*(1 - 0.4*cos((72*pi*k)/1000000).^4) + 0.0625*sin((8*pi*k)/1000000).^18.*...
            sin((24*pi*k)/1000000).^20.*sin((72*pi*k)/1000000).^30 - (0.1*sin((24*pi*k)/1000000).^6 + (1/6)*sin((24*pi*k)/1000000).^30).*...                                                                                                
            (1-sin((8*pi*k)/1000000));

Dk = (pi/2) + (1686*pi*k/1000000);
Ck = 0.003 + 0.097*cos((648*pi*k)/1000000).^40;
Bk = -cos((64*pi*k)/1000000).*Ek - 0.1*cos((64*pi*k)/1000000);
Ak = 1.2*sin((64*pi*k)/1000000).*Ek;

F1      = Ak + 1i*Bk + Ck.*exp(1i*Dk) ;                                        % Focus 1
F2      = Ak + 1i*Bk - Ck.*exp(1i*Dk) ;                                        % Focus 2

centre  = mean([F1;F2],1);                                                    % ellipse centre
angle   = atan(abs(imag(F1)-imag(F2))./abs(real(F1)-real(F2)));              % angle between ellipse horizontal axis and x-axis
c       = sqrt((real(F1)-real(F2)).^2+(imag(F1)-imag(F2)).^2)/2;              % distance from centre to foci
a       = c/e;                                                              % ellipse horizontal axis
b       = a.^2-c.^2;                                                          % ellipse vertical axis
    
t       = linspace(0,2*pi,n)';                                               % parametric coordinate
x0      = a.*cos(t);                                                         % x-coordinates of ellipse if centered with reference axis and horizontal
y0      = b.*sin(t);                                                         % y-coordinates of ellipse if centered with reference axis and horizontal
x       = x0.*cos(angle)-y0.*sin(angle)+real(centre);                         % Absolute x-coordinates
y       = y0.*cos(angle)+x0.*sin(angle)+imag(centre);                         % Absolute y-coordinates

plot(x,y)
toc

[Ahmed Mohamed Mansoor]

This is just amazing. I am very impressed. There is an issue from my side. When I run the code. It overwhelms my RAM and swap memory and matlab then crashes. I tried to limit it to 10000 like you did initialy, and it worked. But when I ran the whole thing, it keeps crashing the software.

Note: I am currently using a MATLAB on ubuntu running on a Dell Precision 5820 Tower with 16 Gigs of Ram and 1.3 TB of disk capacity.

[Davide Masiello]

If you run the code in my last comment without the line

plot(x,y)

then it completes the calculation in less than a second (on my PC).

After that, I tried to execute the plot command and it took around 4 mins to run, but still the picture does not show up, so I guess the problem has to do with graphic memory.

[Ahmed Mohamed Mansoor]

@Davide Masiello I ran it without the plot command and it was able to run in less than a second. The same things happens to me once I start plotting it. Matlab crashes after my RAM gets overwhelmed. Ill probably try to run on this on a computer with a better GPU. Again thank you for the help. I might be back with a few more questions about the code itself.

[Davide Masiello]

@Ahmed Mohamed Mansoor Happy to help, please post the plot if you ever manage to produce it. Also, since you are interested in this cool plots, I suggest you check this out.

https://ch.mathworks.com/matlabcentral/communitycontests/contests/5/entries

[Ahmed Mohamed Mansoor]

@Davide Masiello The coincidence!!!! I was literally looking at them right now and was just mesmerised at the different plots. There's just so much to learn!!

and yes i'll try to get the plot when possible.