% generating the magnetic anomaly from a magnetized infinite prism 
% of square section, size a, depth d and magnetic moment J making angle th 
% with horizontal (th, positive counter-clockwise)
% J.A. Olive, Dec. 2010
clear all


% buried square 
d=200; % depth (m)
a=40; % size (m)

% magnetization
J=10; % norm of magnetization vector (A/m)
th=90; % angle with horizontal (deg, >0 counter-clockwise)

Jx=J*cosd(th);
Jz=J*sind(th);

% observation profile 
Nobs=1000; % number of observation points
x0=linspace(-5*d,5*d,Nobs);


x=[-a/2 a/2 a/2 -a/2];
z=[d-a/2 d-a/2 d+a/2 d+a/2];

V=zeros(1,Nobs);
H=zeros(1,Nobs);

for u=1:Nobs
    
    xa=x+x0(u);
    za=z;
    [Vt, Ht]=quadmag(xa,za,Jx,Jz);
    V(u)=V(u)+Vt;
    H(u)=H(u)+Ht;
    
end
H=-H;

y0=zeros(1,Nobs);
smp=15; % selecting only a few points of the profile 
% to plot magnetic field vectors

subplot(411)
plot(x0,H)
hold on
grid on
ylabel('H (A/m)','fontsize',15)
set(gca,'PlotBoxAspectRatio',[2 1 1])

subplot(412)
plot(x0,V)
hold on
grid on
ylabel('V (A/m)','fontsize',15)
set(gca,'PlotBoxAspectRatio',[2 1 1])

subplot(413)
plot(x0,sqrt(V.^2+H.^2))
hold on
grid on
ylabel('||B|| (A/m)','fontsize',15)
set(gca,'PlotBoxAspectRatio',[2 1 1])

subplot(414)
quiver(x0(1:smp:end),y0(1:smp:end),H(1:smp:end),V(1:smp:end))
hold on
plot(x0,y0,'k')
plot(x,-z,'k.')
axis equal
quiver(0,-d,Jx,Jz,20,'r')
xlim([-5*d 5*d])
xlabel('x (m)','fontsize',15)
ylabel('z (m)','fontsize',15)
set(gca,'PlotBoxAspectRatio',[2 1 1])