%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                Losningsforslag til Potensial rundt dipol (øv 6 TFY4155)     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; % Nullstille alle variabler og deres dimensjoner
% Velg omraadet -2 < x < 2 og -2 < z < 2 for plotting, med skrittlengde
% 2/50 langs baade x og z, dvs ca 100 x 100 punkter i alt
% z forskjøvet 0.01 for å unngå Ve(0,1/2)=0 og division by zero for Relavvik
[x,z] = meshgrid(-2:2/50:2,-2.01:2/50:2);
% Ve er det eksakte potensialet fra dipolen
Ve = (1./(sqrt(x.^2+(z-1/2).^2))-1./(sqrt(x.^2+(z+1/2).^2)));
% Va er dipolpotensialet til ledende orden i a/r naar r >> a
Va = z.*(x.^2.+z.^2).^(-3/2);
% Relavvik er relativ differanse mellom Ve og Va i prosent
Relavvik = 100*abs((Ve-Va)./(Ve));
figure(1); % Ber om figur
subplot(2,2,1); % Viser opptil fire figurer i 2x2-mønster
% mesh(x,y,z) tegner opp et 3D-plott av z som funksjon av x og y
% obs i vår notasjon: y -> z og z -> Ve:
mesh(x,z,Ve);
% Kommandoen axis([a b c d e f]) setter aksene for 3D-plot slik:
% a < x < b, c < y < d, e < z < f
axis([-2 2 -2 2 -10 10]);
% Kommandoen caxis([zmin zmax]) setter fargeskalaen slik at blaatt
% tilsvarer zmin og rodt tilsvarer zmax
caxis([-10 10]);
xlabel('x/a'); % Tekst paa aksene
ylabel('z/a');
zlabel('Ve');
subplot(2,2,2); % Plotter Va ved siden av Ve
mesh(x,z,Va); 
axis([-2 2 -2 2 -10 10]);
caxis([-10 10]);
xlabel('x/a');
ylabel('z/a');
zlabel('Va');
subplot(2,2,3); % Plotter Relavvik under Ve og Va
mesh(x,z,Relavvik); 
axis([-2 2 -2 2 0 100]); % Plotter Relavvik fra 0% til 100%:
% Fargeskala settes for Relavvik mellom 0 og 100
caxis([0 100]);
xlabel('x/a');
ylabel('z/a');
zlabel('|(Ve - Va)/Ve|*100%');