I want to plot an ellipse over the emmitance plot.The ellipse should be the minimum area ellipse that covers the maximum number of points.All the different points in figure 1 need to be covered with a minum area ellipse.Can anyone help?

clc;
close all;
clear all;
%------------------------------------------

%Vacc:-50V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

a = csvread('50V.csv',1,0);
I=a(:,1);%Beam current at z=5mm%
x=a(:,2);
y=a(:,3);
z=a(:,4);
Vx=a(:,5);
Vy=a(:,6);
Vz=a(:,7);

r=sqrt(x.^2+y.^2);
th=(Vx./Vz);
thh=(Vy./Vz);
figure(1)
subplot(4,4,1),scatter_kde(x,th,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('50V');


Ib=[628];
It=[2696];
Ig=[It-Ib];
E=sqrt(mean((x-mean(x)).^2.*(th-mean(th)).^2)-mean((x-mean(x)).*(th-mean(th)).^2));
%------------------------------------------
%Vacc:-100V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

b = csvread('100V.csv',1,0);
I1=b(:,1);%Beam current at z=5mm%
x1=b(:,2);
y1=b(:,3);
z1=b(:,4);
Vx1=b(:,5);
Vy1=b(:,6);
Vz1=b(:,7);

r1=sqrt(x1.^2+y1.^2);
th1=(Vx1./Vz1);

figure(1)
subplot(4,4,2),scatter_kde(x1,th1,'.','LineWidth',2,'MarkerSize',1);

hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('100V');

Ib1=[1028];
It1=[2948];
Ig1=[It1-Ib1];
E1=sqrt(mean((x1-mean(x1)).^2.*(th1-mean(th1)).^2)-mean((x1-mean(x1)).*(th1-mean(th1)).^2));

%------------------------------------------
%Vacc:-150V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

c = csvread('150V1.csv',1,0);
I2=c(:,1);%Beam current at z=5mm%
x2=c(:,2);
y2=c(:,3);
z2=c(:,4);
Vx2=c(:,5);
Vy2=c(:,6);
Vz2=c(:,7);

r2=sqrt(x2.^2+y2.^2);
th2=(Vx2./Vz2);

figure(1)
subplot(4,4,3),scatter_kde(x2,th2,'.','LineWidth',2,'MarkerSize',1);

hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('150V');

Ib2=[1760];
It2=[3092];
Ig2=[It2-Ib2];

E2=sqrt(mean((x2-mean(x2)).^2.*(th2-mean(th2)).^2)-mean((x2-mean(x2)).*(th2-mean(th2)).^2));
%------------------------------------------
%Vacc:-1000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

d = csvread('1000V.csv',1,0);
I3=d(:,1);%Beam current at z=5mm%
x3=d(:,2);
y3=d(:,3);
z3=d(:,4);
Vx3=d(:,5);
Vy3=d(:,6);
Vz3=d(:,7);

r3=sqrt(x3.^2+y3.^2);
th3=(Vx3./Vz3);

figure(1)
subplot(4,4,4),scatter_kde(x3,th3,'.','LineWidth',2,'MarkerSize',2);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('1000V');
Ig3=[0];

E3=sqrt(mean((x3-mean(x3)).^2.*(th3-mean(th3)).^2)-mean((x3-mean(x3)).*(th3-mean(th3)).^2));
%------------------------------------------
%Vacc:-1500V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

e = csvread('1500V.csv',1,0);
I4=e(:,1);%Beam current at z=5mm%
x4=e(:,2);
y4=e(:,3);
z4=e(:,4);
Vx4=e(:,5);
Vy4=e(:,6);
Vz4=e(:,7);

r4=sqrt(x4.^2+y4.^2);
th4=(Vx4./Vz4);

figure(1)
subplot(4,4,5),scatter_kde(x4,th4,'.','LineWidth',2,'MarkerSize',2);

hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('Trace Space');
title('1500V');
Ig4=[0];

E4=sqrt(mean((x4-mean(x4)).^2.*(th4-mean(th4)).^2)-mean((x4-mean(x4)).*(th4-mean(th4)).^2));
%------------------------------------------
%Vacc:-2000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

f = csvread('2000V.csv',1,0);
I5=f(:,1);%Beam current at z=5mm%
x5=f(:,2);
y5=f(:,3);
z5=f(:,4);
Vx5=f(:,5);
Vy5=f(:,6);
Vz5=f(:,7);

r5=sqrt(x5.^2+y5.^2);
th5=(Vx5./Vz5);

figure(1)
subplot(4,4,6),scatter_kde(x5,th5,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('2000V');
Ig5=[0];

E5=sqrt(mean((x5-mean(x5)).^2.*(th5-mean(th5)).^2)-mean((x5-mean(x5)).*(th5-mean(th5)).^2));

%------------------------------------------
%Vacc:-2500V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

g = csvread('2500V.csv',1,0);
I6=g(:,1);%Beam current at z=5mm%
x6=g(:,2);
y6=g(:,3);
z6=g(:,4);
Vx6=g(:,5);
Vy6=g(:,6);
Vz6=g(:,7);

r6=sqrt(x6.^2+y6.^2);
th6=(Vx6./Vz6);

figure(1)
subplot(4,4,7),scatter_kde(x6,th6,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('2500V');
Ig6=[0];

E6=sqrt(mean((x6-mean(x6)).^2.*(th6-mean(th6)).^2)-mean((x6-mean(x6)).*(th6-mean(th6)).^2));

%Vacc:-3500V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

g1 = csvread('3500V.csv',1,0);
I61=g1(:,1);%Beam current at z=5mm%
x61=g1(:,2);
y61=g1(:,3);
z61=g1(:,4);
Vx61=g1(:,5);
Vy61=g1(:,6);
Vz61=g1(:,7);

r61=sqrt(x61.^2+y61.^2);
th61=(Vx61./Vz61);

figure(1)
subplot(4,4,8),scatter_kde(x61,th61,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('3500V');

Ib61=[4080];
It61=[4116];
Ig61=[It61-Ib61];
E61=sqrt(mean((x61-mean(x61)).^2.*(th61-mean(th61)).^2)-mean((x61-mean(x61)).*(th61-mean(th61)).^2));

%------------------------------------------
%Vacc:-4000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

g2 = csvread('4000V.csv',1,0);
I62=g2(:,1);%Beam current at z=5mm%
x62=g2(:,2);
y62=g2(:,3);
z62=g2(:,4);
Vx62=g2(:,5);
Vy62=g2(:,6);
Vz62=g2(:,7);

r62=sqrt(x62.^2+y62.^2);
th62=(Vx62./Vz62);

figure(1)
subplot(4,4,9),scatter_kde(x62,th62,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('4000V');

Ib62=[4152];
It62=[6612];
Ig62=[It62-Ib62];
E62=sqrt(mean((x62-mean(x62)).^2.*(th62-mean(th62)).^2)-mean((x62-mean(x62)).*(th62-mean(th62)).^2));
%------------------------------------------
%Vacc:-5000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

h = csvread('5000V.csv',1,0);
I7=h(:,1);%Beam current at z=5mm%
x7=h(:,2);
y7=h(:,3);
z7=h(:,4);
Vx7=h(:,5);
Vy7=h(:,6);
Vz7=h(:,7);

r7=sqrt(x7.^2+y7.^2);
th7=(Vx7./Vz7);

figure(1)
subplot(4,4,10),scatter_kde(x7,th7,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('5000V');

Ib7=[4216];
It7=[8620];
Ig7=[It7-Ib7];
E7=sqrt(mean((x7-mean(x7)).^2.*(th7-mean(th7)).^2)-mean((x7-mean(x7)).*(th7-mean(th7)).^2));
%------------------------------------------
%Vacc:-6000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

i = csvread('6000V.csv',1,0);
I8=i(:,1);%Beam current at z=5mm%
x8=i(:,2);
y8=i(:,3);
z8=i(:,4);
Vx8=i(:,5);
Vy8=i(:,6);
Vz8=i(:,7);

r8=sqrt(x8.^2+y8.^2);
th8=(Vx8./Vz8);

figure(1)
subplot(4,4,11),scatter_kde(x8,th8,'.','LineWidth',2,'MarkerSize',0.5);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('6000V');

D6=std(th6)*180./3.14;%Divergence in degree%

Ib8=[4292];
It8=[9324];
Ig8=[It8-Ib8];
E8=sqrt(mean((x8-mean(x8)).^2.*(th8-mean(th8)).^2)-mean((x8-mean(x8)).*(th8-mean(th8)).^2));
%------------------------------------------
%Vacc:-8000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

j = csvread('8000V.csv',1,0);
I9=j(:,1);%Beam current at z=5mm%
x9=j(:,2);
y9=j(:,3);
z9=j(:,4);
Vx9=j(:,5);
Vy9=j(:,6);
Vz9=j(:,7);

r9=sqrt(x9.^2+y9.^2);
th9=(Vx9./Vz9);

figure(1)
subplot(4,4,12),scatter_kde(x9,th9,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('8000V');


Ib9=[4352];
It9=[10412];
Ig9=[It9-Ib9];

E9=sqrt(mean((x9-mean(x9)).^2.*(th9-mean(th9)).^2)-mean((x9-mean(x9)).*(th9-mean(th9)).^2));
%------------------------------------------
%Vacc:-10000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

k = csvread('10000V.csv',1,0);
I10=k(:,1);%Beam current at z=5mm%
x10=k(:,2);
y10=k(:,3);
z10=k(:,4);
Vx10=k(:,5);
Vy10=k(:,6);
Vz10=k(:,7);

r10=sqrt(x10.^2+y10.^2);
th10=(Vx10./Vz10);

figure(1)
subplot(4,4,13),scatter_kde(x10,th10,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('10kV');

Ib10=[4392];
It10=[10952];
Ig10=[It10-Ib10];

E10=sqrt(mean((x10-mean(x10)).^2.*(th10-mean(th10)).^2)-mean((x10-mean(x10)).*(th10-mean(th10)).^2));
%------------------------------------------
%Vacc:-15000V Vp:25V Vsc:0V Vcol:0V   Bohmdensity:2E-5/cm2  Te:50000K Meniscus:2.16V%

l = csvread('15000V.csv',1,0);
I11=l(:,1);%Beam current at z=5mm%
x11=l(:,2);
y11=l(:,3);
z11=l(:,4);
Vx11=l(:,5);
Vy11=l(:,6);
Vz11=l(:,7);

r11=sqrt(x11.^2+y11.^2);
th11=(Vx11./Vz11);

figure(1)
subplot(4,4,14),scatter_kde(x11,th11,'.','LineWidth',2,'MarkerSize',1);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('15kV');

Ib11=[4504];
It11=[12096];
Ig11=[It11-Ib11];

E11=sqrt(mean((x11-mean(x11)).^2.*(th11-mean(th11)).^2)-mean((x11-mean(x11)).*(th11-mean(th11)).^2));
%-----------------------------------------------------------------
hv=[50,100,150,1000,1500,2000,2500,3500,4000,5000,6000,8000,10000,15000];
Ig=[Ig/It,Ig1/It1,Ig2/It2,Ig3,Ig4,Ig5,Ig6,Ig61/It61,Ig62/It62,Ig7/It7,Ig8/It8,Ig9/It9,Ig10/It10,Ig11/It11]*100;
figure(2)
plot(hv,Ig,'>k','LineWidth',2,'MarkerSize',12);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
ylabel('Grid-Total current ratio');
xlabel('HV');
title('Perveance');
axis([0 20000 0 100])


%-----------------------------------------------------------------
hv=[50,100,150,1000,1500,2000,2500,3500,4000,5000,6000,8000,10000,15000];
Emmitance=[E,E1,E2,E3,E4,E5,E6,E61,E62,E7,E8,E9,E10,E11];
figure(3)
plot(hv,Emmitance,'Ob','LineWidth',2,'MarkerSize',12);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
ylabel('rms Emmitance(mm-rad)');
xlabel('HV');
title('Rms-emmitance');
axis([0 20000 0 0.2])

figure(52)
fit_ellipse_zero_R1(x5,th5);
hold on
grid on
set (gca,'fontsize',16, 'fontweight', 'b')
xlabel('Radial distance(x-mm)');
ylabel('Theta');
title('Trace Space');