clc; close all; clear;
t0=0e-6;
v=input('Enter Ratio v/c between 0 and 0.95: ');
Phi0=input('Enter Magnitude of Potential between 0.5 and 1: ');
D=input('Enter Distant D from Electron > 1e-3[m] and < 10[m]: ');
GeometryFun(v);
ScalarVectorPotentialTogether(v,t0,Phi0);
fprintf(2,'\nPlease be patient. Wait for animation\n');
NewMovingE(v,t0,Phi0);
fprintf(2,'\nPlease be patient. Wait for animation\n');
NewMovingB(v,t0,Phi0);
fprintf(2,'\nPlease be patient. Wait for 3D image\n');
EnergyFlux(v,D);
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function GeometryFun(v);
%% Geometry
ff=figure(1); hold all; grid minor;
movegui(ff,'northwest'); 
set(gca,'xtick',[]);
set(gca,'ytick',[]);
a=[0 0]; 
b=[6 9];
c=[4 0];
plot([a(1) b(1)],[a(2) b(2)],'*b');
line([a(1) b(1)],[a(2) b(2)],'Color','b');
text(b(1)+.2,b(2),'\it(x, y,z)','FontWeight','bold','Fontsize',20);
line([b(1),b(1)],[0,b(2)],'Color','k','LineStyle','--');
%Retarded Position
x3 = [.26 .20];
y3 = [.3 .11];
text(a(1)-0.7,a(2)+.4,'\itP''','Color','b','FontWeight','bold','Fontsize',20);
ar3=annotation('textarrow',x3,y3,'String','RETARDED POSITION');
line([a(1),a(1)],[0,b(2)/2],'Color','k','LineStyle','--');
ar3.Color = 'b';
%Present Position
x0 = [.555 .485 ];
y0 = [.2 .11];
plot([c(1) b(1)],[c(2) b(2)],'*r');
line([c(1) b(1)],[c(2) b(2)],'Color','r');
text(c(1)-0.65,c(2)+.4,'\itP','Color','r','FontWeight','bold','Fontsize',20);
ar=annotation('textarrow',x0,y0,'String','PRESENT POSITION');
ar.Color = 'r';
%x - vt
x1 = [0.625 .485];
y1 = [.08 .08];
ar1=annotation('doublearrow',x1,y1);
ar1.Color = 'k';
ar1.LineWidth = 3;
text(4.3,-.7,'\itx - vt','Color','k','FontWeight','bold','Fontsize',16);
%vt
x1 = [0.131 0.485];
y1 = [.08 .08];
ar1=annotation('arrow',x1,y1);
ar1.Color = 'k';
ar1.LineWidth = 3;
text(1.3,-.7,'\itvt','Color','k','FontWeight','bold','Fontsize',16);
%vt'
x1 = [0.131 0.2];
y1 = [.4 .4];
ar1=annotation('arrow',x1,y1);
ar1.Color = 'k';
ar1.LineWidth = 3;
text(-.95,b(1)/2+1,'\itvt''','Color','k','FontWeight','bold','Fontsize',16);
text(-.95,b(1)/2+2,'\itt'' = t - r''/c','Color','k','FontWeight','bold','Fontsize',16);
%Speed Vector
x2 = [.485 .7];
y2 = [.11 .11];
ar2=annotation('arrow',x2,y2);
ar2.Color = 'r';
ar2.LineWidth = 4;
text(7,.5,['\itV/c = ',num2str(v)],'Color','r','FontWeight','bold','Fontsize',20);
% r-vectors
text(3.8,6.5,'\itr''','Color','b','FontSize',24);
text(4.8,5.5,'\itr','Color','r','FontSize',24);
%Coordinate Label
xlabel('X-direction','Position',[8.5 0],'Fontsize',16,'FontWeight','bold','Color','k');
ylabel('Y-direction','Position',[-1 8],'Fontsize',16,'FontWeight','bold','Color','k');
axis([-1 10 0 10]);
hold off;
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ScalarVectorPotentialTogether(v,t,Phi0)
N=2e2;
[X,Y,Z]=meshgrid(linspace(-1,1,N)*20);
k0=4*pi*8.8541878128e-12;
q=k0; gamma2=(1-v^2);
%% Potential 
Phi=(q/k0)/((X-v*t).^2+gamma2*(Y.^2+Z.^2)).^(1/2);
%%
cc={'k','b','g','y','m','r'}; jc=2;
%%
f2=figure('units','normalized','outerposition',[0.6 0 0.4 1]); 
ax1=subplot(2,1,1); grid on; hold on; set(gcf,'color','w');
%% Charge
    [xs,ys,zs] = sphere;
    if v~=0
        xs = v*xs*0.125; ys = v*ys*0.125; zs = v*zs*0.125;
    else
        xs = xs*0.05; ys = ys*0.05; zs = zs*0.05; 
    end
    hh=surf(xs,ys,zs); set(hh,'FaceColor','k','EdgeColor','none'); 
     h1=patch(isosurface(X,Y,Z,Phi,Phi0)); 
     set(h1,'FaceColor',char(cc(jc)),'FaceAlpha',0.3);
     xd=get(h1,'XData'); yd=get(h1,'YData'); zd=get(h1,'ZData');
     arrow(0,0,0,1.5*max(xd(:)),0,0,6);
     if v~=0; tx=text(1.5*max(xd(:)),0,.2*max(zd(:)),'\bfV');
     else; tx=text(1.5*max(xd(:)),0,.2*max(zd(:)),'\bfX'); end;
     tx.Color=char(cc(6)); tx.FontSize=20;
     arrow(0,0,0,0,1.5*max(yd(:)),0,6);
     ty=text(0,1.5*max(yd(:)),0,'\bfY-axis'); ty.Color=char(cc(6));
     arrow(0,0,0,0,0,1.9*max(zd(:)),6);
     tz=text(0,0,1.9*max(zd(:)),'\bfZ-axis'); tz.Color=char(cc(6));
     arrow(min(xd(:)),0,min(zd(:)),max(xd(:)),0,min(zd(:)),1);
     arrow(max(xd(:)),0,min(zd(:)),min(xd(:)),0,min(zd(:)),1);
     tz=text(0,0,min(zd(:))*1.1,['\bfb = ',num2str(max(xd(:))-min(xd(:)))]);
     tz.Color=char(cc(1));
     arrow(0,min(yd(:)),max(zd(:)),0,max(yd(:)),max(zd(:)),1);
     arrow(0,max(yd(:)),max(zd(:)),0,min(yd(:)),max(zd(:)),1);
     tz=text(-0.3,0,max(zd(:))*1.1,['\bfa = ',num2str(max(yd(:))-min(yd(:)))]);
     tz.Color=char(cc(1));
     tz=text(0,0,min(zd(:))*1.3,['\bfb/a = ',num2str((max(xd(:))-min(xd(:)))/(max(yd(:))-min(yd(:))))]);
     tz.Color=char(cc(1));
     if v~=0;
        tt=text(1.5*max(xd(:)),0,1.4*max(zd(:)),{'\bfScalar Potential Isosurface';'\bfCharge Moving with'; '\bfConstant Speed along X-axis'});
     else;
        tt=text(1.5*max(xd(:)),0,1.4*max(zd(:)),{'\bfScalar Potential Isosurface';'\bfCharge at Rest';'Coulomb Law'});
     end;
     tt.FontSize=12; tt.Color=char(cc(6));
h=camlight('left'); set(gca,'visible','off');
daspect([1,1,1]);axis tight; if v==0; view(135,0); end;
%% Vector Potential
    Ax=v*Phi;
    ax2=subplot(2,1,2); grid on; hold on;
    set(gcf,'color','w');
if v==0;   
    A1='Charge at Rest. Vector Potential Ax = v*'; A2 = char(966);
    A3=' = 0.'; A4='  There is no Magnetic Fields';
    fprintf(2,['\n',A1,A2,A3,A4,'\n']);
    title({'\bf',[A1,A2,A3,A4]});
else;    
    hh2=surf(xs,ys,zs);
    set(hh2,'FaceColor','k','EdgeColor','none'); 
    h2 = copyobj(h1,ax2);
     xd=get(h2,'XData'); yd=get(h2,'YData'); zd=get(h2,'ZData');
     xmax=max(xd(:)); ymax=max(yd(:)); zmax=max(xd(:)); 
     arrow(0,0,0,1.5*xmax,0,0,3);
     tx=text(2*xmax,0,.2*zmax,'\bf\itAx');
     tx.Color=char(cc(3)); tx.FontSize=20;
     arrow(min(xd(:)),0,min(zd(:)),max(xd(:)),0,min(zd(:)),1);
     arrow(max(xd(:)),0,min(zd(:)),min(xd(:)),0,min(zd(:)),1);
     tz=text(0,0,min(zd(:))*1.1,['\bfb = ',num2str(max(xd(:))-min(xd(:)))]);
     tz.Color=char(cc(1));
     arrow(0,min(yd(:)),max(zd(:)),0,max(yd(:)),max(zd(:)),1);
     arrow(0,max(yd(:)),max(zd(:)),0,min(yd(:)),max(zd(:)),1);
     tz=text(-0.3*v,0,max(zd(:))*1.1,['\bfa = ',num2str(max(yd(:))-min(yd(:)))]);
     tz.Color=char(cc(1));
     tz=text(0,0,min(zd(:))*1.3,['\bfb/a = ',num2str((max(xd(:))-min(xd(:)))/(max(yd(:))-min(yd(:))))]);
     tz.Color=char(cc(1));
     P=20;
     for jj=0:P-1;
         arrow(-1.5*xmax*v,ymax*sin(2*pi*jj/P)/2,zmax*cos(2*pi*jj/P)/2,1.5*xmax*v,ymax*sin(2*pi*jj/P)/2,zmax*cos(2*pi*jj/P)/2,3);
         arrow(-1.5*xmax*v/2,ymax*sin(2*pi*jj/P)/1.5,zmax*cos(2*pi*jj/P)/1.5,1.5*xmax*v/2,ymax*sin(2*pi*jj/P)/1.5,zmax*cos(2*pi*jj/P)/1.5,3);
         arrow(-1.5*xmax*v/3,1.2*ymax*sin(2*pi*jj/P),1.2*zmax*cos(2*pi*jj/P),1.5*xmax*v/3,1.2*ymax*sin(2*pi*jj/P),1.2*zmax*cos(2*pi*jj/P),3);
     end;   
     arrow(0,0,0,0,1.5*ymax,0,6);
     ty=text(0,1.5*ymax,0,'\bfY-axis'); ty.Color=char(cc(6));
     arrow(0,0,0,0,0,2.4*zmax,6);
     tz=text(0,0,2.4*zmax,'\bfZ-axis'); tz.Color=char(cc(6));
     tt=text(1.5*xmax,0,2*zmax,{'\bfVector Potential Isosurface';'\bfCharge Moving with';'\bfConstant Speed along X-axis'});
     tt.FontSize=12; tt.Color=char(cc(6));
    camlight('left'); set(gca,'visible','off');
    daspect([1,1,1]);axis tight;
    linkprop([ax1,ax2],{'CameraPosition','CameraUpVector'});
    for j1=1:46;
        az=90+(j1-1)*2;
        view(az,0); pause(1);
    end;
end;
hold off;      
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function NewMovingE(v,t,Phi0)
N=1e2;
[X,Y,Z]=meshgrid(linspace(-1,1,N)*10);
k0=4*pi*8.8541878128e-12;
q=1.602176634e-19; 
%% Potential 
gamma2=(1-v^2);
Phi=(q/k0)/((X-v*t).^2+gamma2*(Y.^2+Z.^2)).^(1/2);
Phi=Phi/max(Phi(:));
[Ex,Ey,Ez]=gradient(-Phi,20/N);
Ex=Ex*(1-v^2);
E=sqrt(Ex.^2+Ey.^2+Ez.^2);
E=E./max(E(:));
FF=log(E.^2);
%
f3=figure('units','normalized','outerposition',[0.3 0.3 0.6 0.6]); 
hold on; grid on; view(-210,30); axis square;
% xlim([min(X(:)) max(X(:))]); ylim([min(X(:)) max(X(:))]); zlim([min(X(:)) max(X(:))]);
title('\bfElectric Field ({\color{red}vector lines}) and Equipotential Lines ({\color{blue}blue})');
if v~=0;
    arrow(0,0,0,1.5*max(X(:)),0,0,6);
    tx=text(1.7*max(X(:)),0,0.1*max(Z(:)),'\bfV');
    tx.FontSize=20; tx.Color='r';
    t1=text(1.7*max(X(:)),0,-.1*max(Z(:)),['\bfv/c = ',num2str(v)]);
    t1.Color='r'; t1.FontSize=12;
end;
text(max(X(:)),max(Y(:)),max(Z(:)),{['max(Ex) = ',num2str(max(Ex(:))),'[V/m]'];...
    ['max(Ey) = ',num2str(max(Ey(:))),'[V/m]'];...
    ['max(Ez) = ',num2str(max(Ez(:))),'[V/m]']});
%
h1=contourslice(X,Y,Z,FF,[],0,[],60); 
hs1=streamslice(X,Y,Z,Ex,Ey,Ez,[],0,[],2); set(hs1,'Color','r');
h2 = contourslice(X,Y,Z,FF,0,[],[],60); 
hs2=streamslice(X,Y,Z,Ex,Ey,Ez,0,[],[],2); set(hs2,'Color','r');
c=jet; c=c(1:3:end,:); colormap(c); hc=colorbar; shading interp;
hc.TickLabels = num2cell(FF);
ylabel(hc,'\bf Ln(Relative E-field Energy Intensity)','FontSize',20);
zlabel('\bfZ-axis'); ylabel('\bfY-axis'); xlabel('\bfX-axis');
h3 = contourslice(X,Y,Z,FF,[],[],min(Z(:)),60); 
h4 = slice(X,Y,Z,FF,[],[],min(Z(:))); set(h4,'FaceColor','interp','EdgeColor','none');
hs3=streamslice(X,Y,Z,Ex,Ey,Ez,[],[],min(Z(:)),2); 
set(hs3,'Color','r');view(20,21);
vals = linspace(min(Z(:)),max(Z(:)),11);
% daObj=VideoWriter(FileName,'MPEG-4')
% daObj.Quality = 95;
% open(daObj);
for id = vals; 
    delete(h3); delete(hs3); delete(h4); 
    h3=contourslice(X,Y,Z,FF,id,[],[],60); 
    h4 = slice(X,Y,Z,FF,id,[],[]); set(h4,'FaceColor','interp','EdgeColor','none');
    hs3=streamslice(X,Y,Z,Ex,Ey,Ez,id,[],[],2); 
    set(hs3,'Color','r'); view(20,21); axis square;
    drawnow;
%     writeVideo(daObj,getframe(gcf));
    pause(1);
end;
% close(daObj);
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function NewMovingB(v,t,Phi0)
k0=4*pi*8.8541878128e-12;
q=k0;
N=1e2;
[X,Y,Z]=meshgrid(linspace(-1,1,N)*10);
cc={'k','b','g','y','m','r'};
%% Potential 
gamma2=(1-v^2);
Ax=Phi0*v*(q/k0)/((X-v*t).^2+gamma2*(Y.^2+Z.^2)).^(1/2);
[Bx,By,Bz]=curl(X,Y,Z,Ax,0*Ax,0*Ax);
B=sqrt(Bx.^2+By.^2+Bz.^2);
R=sqrt(X.^2+Y.^2+Z.^2);
%%
figure('units','normalized','outerposition',[0.3 0.1 0.6 0.6]); 
grid minor; axis('square','tight'); hold on; view(3);
xlabel('\bfX-axis'); ylabel('\bfY-axis'); zlabel('\bfZ-axis');
axis tight;
%% Vector Potential
vals = linspace(min(X(:)),max(X(:)),21);
if v==0;   
    A1='Charge at Rest. Vector Potential Ax = v*'; A2 = char(966);
    A3=' = 0.'; A4='  There is no Magnetic Fields';
    fprintf(2,['\n',A1,A2,A3,A4,'\n']);
    title({'\bf',[A1,A2,A3,A4]});
else;
B=B/max(B(:));
B=log(B.^2);
%Arrow
arrow(min(X(:)),0,0,max(X(:))+5,0,0,4);
tx=text(max(X(:))+5,0,.1*max(Z(:)),'\bfV');
tx.FontSize=20; tx.Color='y';
t1=text(max(X(:)),max(Y(:)),-.1*max(Z(:)),['\bfv/c = ',num2str(v)]);
t1.Color='r'; t1.FontSize=12;
text(max(X(:)),max(Y(:)),max(Z(:)),{['max(Bx) = ',num2str(max(Bx(:))),'[T]'];...
    ['max(By) = ',num2str(max(By(:))),'[T]'];...
    ['max(Bz) = ',num2str(max(Bz(:))),'[T]']});
% 
h1=contourslice(X,Y,Z,B,0,[],[],20); set(h1,'linewidth',1.5);
h2=contourslice(X,Y,Z,B,[],0,[],20); set(h2,'linewidth',2);
h3=contourslice(X,Y,Z,B,[],[],0,20); set(h3,'linewidth',1.5);
%
hs1=streamslice(X,Y,Z,Bx,By,Bz,v*t,[],[],1); set(hs1,'Color','b');
hs2=streamslice(X,Y,Z,Bx,By,Bz,[],0,[],2); set(hs2,'Color','b');
hs3=streamslice(X,Y,Z,Bx,By,Bz,[],[],0,2); set(hs3,'Color','b');
title('\bfB-Field ({\color{blue}vector lines}) and Equipotential Lines ({\color{red}red})');
c=hot; c=c(12:2:end,:); colormap(c);  
hc=colorbar; hc.TickLabels = num2cell(B);
ylabel(hc,'\bf Ln(Relative B-field Energy Intensity)','FontSize',20);
view(20,21);
%
for id = vals;
    delete(hc);
    h=contourslice(X,Y,Z,B,id,[],[],20);
    hs=streamslice(X,Y,Z,Bx,By,Bz,id,[],[],1); set(hs,'Color','b');
    h4 = slice(X,Y,Z,B,id,[],[]); set(h4,'FaceColor','interp','EdgeColor','none');
    hc=colorbar; hc.TickLabels = num2cell(B); 
    ylabel(hc,'\bf Ln(Relative B-field Energy Intensity','FontSize',20);
    pause(.1);
    if id~=vals(end);
        delete(h); delete(hs); delete(hc); delete(h4);
    end;
end;
end;
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function EnergyFlux(v,D)
cc={'k','b','g','y','m','r'};
c=299792458; %Light Speed [m/s]
gamma2=(1-v^2);
eps0=8.8541878128e-12;
k0=4*pi*eps0;
mu0=4*pi*1e-7;
q=1.602176634e-19;
%
N=1.5e2;
[X,Y,Z]=meshgrid(linspace(-1,1,N)*D);
R=sqrt(X.^2+Y.^2+Z.^2);
%E- and B-field. Energy Flux g
    t=0*c;
    p=((X-v*t).^2+(1-v^2)*(Y.^2+Z.^2)).^(3/2);
    Ex=(q/k0)*(1-v^2)*(X-v*t)./p; 
    Ey=(q/k0)*(1-v^2)*Y./p; 
    Ez=(q/k0)*(1-v^2)*Z./p; 
    Bx=0*Ex;
    By=v*(1-v^2 )*Z./p;
    Bz=-v*(1-v^2)*Y./p;
    ux=(Ey.*Bz-Ez.*By)/mu0;
    uy=(Ez.*Bx-Ex.*Bz)/mu0;
    uz=(Ex.*By-Ey.*Bx)/mu0;
    U=sqrt(ux.^2+uy.^2+uz.^2);
    MaxU=max(U(:));
    MinU=min(U(:));
    UL=log10(U); 
    %
    ff=figure('units','normalized','outerposition',[0 0 1 1]); 
    grid minor; axis('square','tight'); hold on; view(-12,21);
    xlabel('\bfX-axis'); ylabel('\bfY-axis'); zlabel('\bfZ-axis');
    title('\bfEM Energy Flux [W/m^2] and its Isosurfaces (green) where |u| = const.')
    axis tight
    %
    h1=contourslice(X,Y,Z,UL,0,[],[],40); set(h1,'linewidth',1.5)
    cm = get(gca,'Colormap'); hc = colorbar; 
    ylabel(hc,'\bflog10(Intensity of Energy Flux due to EM Field) [W/m^2]','FontSize',12)
    cslicecolor2linestyle(h1)
    hh1 = slice(X,Y,Z,UL,0,[],[]); set(hh1,'FaceColor','interp','EdgeColor','none','FaceAlpha',0.5);
    h2=contourslice(X,Y,Z,UL,[],0,[],40); set(h2,'linewidth',2)
    cslicecolor2linestyle(h2)
    hh2 = slice(X,Y,Z,UL,[],0,[]); set(hh2,'FaceColor','interp','EdgeColor','none','FaceAlpha',0.5);
    h3=contourslice(X,Y,Z,UL,[],[],0,40); set(h3,'linewidth',1.5)
    cslicecolor2linestyle(h3)
    hh3 = slice(X,Y,Z,UL,[],[],0); set(hh3,'FaceColor','interp','EdgeColor','none','FaceAlpha',0.5);
    %
    hp=patch(isosurface(X,Y,Z,UL,log10(MinU*1e3/2))); 
    text(max(X(:))*1.2,max(Y(:)),max(Z(:)),['Flux Isovalue = ',num2str((MinU*1e3/2)),' [W/m^2]'])
    set(hp,'FaceColor',char(cc(3)),'EdgeColor','none','FaceAlpha',0.6);
    hp=patch(isosurface(X,Y,Z,UL,log10(MinU*1e3/4)));
    text(max(X(:))*1.2,max(Y(:)),max(Z(:))*.9,['Flux Isovalue = ',num2str((MinU*1e3/4)),' [W/m^2]'])
    set(hp,'FaceColor',char(cc(3)),'EdgeColor','none','FaceAlpha',0.5);
    hp=patch(isosurface(X,Y,Z,UL,log10(MinU*1e3/8)));
    text(max(X(:))*1.2,max(Y(:)),max(Z(:))*.8,['Flux Isovalue = ',num2str((MinU*1e3/8)),' [W/m^2]'])
    set(hp,'FaceColor',char(cc(3)),'EdgeColor','none','FaceAlpha',0.4);
    hp=patch(isosurface(X,Y,Z,UL,log10(MinU*1e3/16)));
    text(max(X(:))*1.2,max(Y(:)),max(Z(:))*.7,['Flux Isovalue = ',num2str((MinU*1e3/16)),' [W/m^2]'])
    set(hp,'FaceColor',char(cc(3)),'EdgeColor','none','FaceAlpha',0.3);
    %
    hs1=streamslice(X,Y,Z,ux,uy,uz,0,[],[],2); set(hs1,'Color','k','LineWidth',2);
    hs2=streamslice(X,Y,Z,ux,uy,uz,[],0,[],2); set(hs2,'Color','k','LineWidth',2);
    hs3=streamslice(X,Y,Z,ux,uy,uz,[],[],0,2); set(hs3,'Color','k','LineWidth',2);
    % 
function cslicecolor2linestyle(h)
% Greg Aloe (2020),https://www.mathworks.com/matlabcentral/fileexchange/1157-cslicecolor2linestyle 
% Find all the unique colors
colors = get(h,'cdata');
colors=unique(cat(1,colors{:}));
colors=colors(~isnan(colors));
% Loop through all the patches returned by CONTOURSLICE, and designate a linestyle for each
% Define the line style (note that this can be changed since the code is written generally)
linespec = {'-','--',':','-.'};
linestyles = repmat(linespec,1,ceil(length(colors)/length(linespec)));
linestyles = {linestyles{1:length(colors)}};
for n=1:length(h)
    % Find the unique color number associated with the handle
    color = find(max(get(h(n),'cdata'))==colors);
    % Convert the color to the associated linestyle
    linestyle = linestyles{color};
    set(h(n),'linestyle',linestyle);
end
end
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = arrow(x0, y0, z0, x1, y1, z1, c)
    X = [x0, x1, nan]; Y = [y0, y1, nan]; Z = [z0, z1, nan];
    u = [x1 - x0; y1 - y0; z1 - z0]; uN = norm(u);
    for i = 1:200;
        rad = cross(u, ones(3, 1) - 2*rand(3, 1)); 
        rad = rad/norm(rad)*uN/20;   
        X = [X, x1, x1 - (x1 - x0)/5 + rad(1), nan];
        Y = [Y, y1, y1 - (y1 - y0)/5 + rad(2), nan];
        Z = [Z, z1, z1 - (z1 - z0)/5 + rad(3), nan];
    end;
    cc={'k','b','g','y','m','r'};
    a = line(X, Y, Z, 'linewidth', 2, 'color', char(cc(c)));
end