%
% demonstrate various ways of plotting fields in Matlab
%
% start with a clean slate
close all; clear;
% create a meshgrid of sample points
dx = 0.075; dy = 0.075; dz=0.075;
[x,y,z] = meshgrid(-2:dx:2,-2:dy:2,-2:dz:2);
% calculate the potential function, V(x,y,z)
% in this case two equal and opposite charges separated by 2d
% choose k=q=d=1 to give convenient sized numbers
k=1; q=1; d=1;
V=(k*q./(x.^2+(y-d).^2+z.^2).^(1/2))-(k*q./(x.^2+(y+d).^2+z.^2).^(1/2));
% for a simpler example (single point charge) replace with the following
%V = k./(x.^2+y.^2+z.^2).^(1/2);
% place 3-D plots in figure 1
figure(1);
% show the volume as a solid
slice(x,y,z,V,[-2 1.95],[-2 1.95],[-2 1.95]);shading interp;
colormap hsv;
xlabel('x');
ylabel('y');
zlabel('z');
title('Electric Field of two point charges');
axis square;
colorbar;
rotate3d on;
pause;
% view a single slice
slice(x,y,z,V,[],[],0); shading interp;
xlabel('x'); ylabel('y'); zlabel('z');
title('Slice through z=0');
colorbar;
pause;
% add slices in vertical planes
slice(x,y,z,V,0,1,0);shading interp;
xlabel('x'); ylabel('y'); zlabel('z');
title('Add slices through x=0 and y=1');
colorbar;
pause;
%place 2-D plots in figure 2
figure(2);
x2=[-2:dx:2];
y2=[-2:dy:2];
% plot the slice through z=0
plane=V(1:54,1:54,27);
surf(x2,y2,plane);
xlabel('x'); ylabel('y'); zlabel('z');
title('A mountain plot of the plane z=0');
pause;
% plot the z=0 plane as a contour plot
contour(x2,y2,plane,20);
xlabel('x'); ylabel('y'); zlabel('z');
title('A contour plot of the plane z=0');
pause;
% plot two isosurfaces
figure(1);
isosurface(x,y,z,V,1);
isosurface(x,y,z,V,-1);
xlabel('x'); ylabel('y'); zlabel('z');
title('Isosurfaces V=1 and V=-1');
pause;
% add lighting to show surfaces more clearly
camlight;
lighting gouraud;
pause;
hold off;
% compute the field as the gradient of the potential
[Ex,Ey,Ez] = gradient(V);
% plot the magnitude of the field
M=(Ex.^2+Ey.^2+Ez.^2).^.5;
slice(x,y,z,abs(M),0,1,0);shading interp;
xlabel('x'); ylabel('y'); zlabel('z');
title('Magnitude of the E field');
colormap(jet);
colorbar;
pause
hold off
slice(x,y,z,V,0,1,0);shading interp;
[cx,cy,cz]=meshgrid([-1.5 : .25 : 1.5]);
coneplot(x,y,z,Ex,Ey,Ez,cx,cy,cz,20);
xlabel('x'); ylabel('y'); zlabel('z');
title('Cone plot of the field');
pause;
figure(2);
Px=Ex(1:54,1:54,27);
Py=Ey(1:54,1:54,27);
quiver(Px,Py);
title('quiver plot of plane z=0');
pause;
streamslice(Px,Py);
title('streamline plot of z=0');
pause;
figure(3);
slice(x,y,z,V,0,1,0);shading interp;
alpha(.5);
xlabel('x'); ylabel('y'); zlabel('z');
%title('streamline plot of entire volume');
axis square;
colorbar;
rotate3d on;
%pause;
[sx,sy,sz] = meshgrid(-.075:dx/2:.075,-1.075:dy/2:-.925,-.075:dz/2:.075);
streamline(x,y,z,Ex,Ey,Ez,sx,sy,sz);
[sx,sy,sz] = meshgrid(-.075:dx/2:.075,1.075,-.075:dz/2:.075);
streamline(x,y,z,-Ex,-Ey,-Ez,sx,sy,sz);
xlabel('x'); ylabel('y'); zlabel('z');
title('streamline plot');
pause;
figure(1);
streamslice(x,y,z,Ex,Ey,Ez,0,1,0);
xlabel('x'); ylabel('y'); zlabel('z');
title('streamslice plot on several planes');
%pause;

