clear all
clf

f = 0;

R = 6.371*10^8; %radius of Earth in centimeters 
theta = 11.5;
inc = 70; %inclination
M = 8*10^25; %magnetic moment of the dipole
Fe = 2*M*cosd(theta)/(R^3)*10^5/sind(inc);
z = 200; %depth to top of dike
d = 10; %radius of dike
A = d^2; %cross sectional area of dike

L = 90;
N = 90;

k = zeros(N+1,N+1);
T = 20*ones(N+1,N+1);
for i = 2:N
    for j = 2:N
            T(i,j) = 1100;
            C = 1*10^5; %Curie constant is material dependent
            mu = 1.257*10^(-6); %permeability of free space
            k(i,j) = C*mu./(T(i,j)-C*mu);%magnetic susceptibility contrast, basalt intruding into sandstone 
    end
end

dx = L/N;
dy = dx;
x=zeros(1,N+1);
y=zeros(1,N+1);
for i=1:N+1
    x(i)=dx*(i-1);
    y(i)=dy*(i-1);
end

Za = zeros(N+1,N+1);
Hax = zeros(N+1,N+1);

for i = 1:N+1
    for j = 1:N+1 
        Za(i,j) = z*k(i,j)*Fe*A/((i-N/2).^2 + (j-N/2).^2 + z.^2).^(3/2);
        Hax(i,j) = i*k(i,j)*Fe*A/((i-N/2).^2 + (j-N/2).^2 + z.^2).^(3/2); %The dike is centered at a position of (50,50)
    end
end


for q = 1:180;
    for j = 1:N+1 
        for i = 1:N+1
            Za(i,j) = (z-q)*k(i,j)*Fe*A/((i-N/2).^2 + (j-N/2).^2 + (z-q).^2).^(3/2);
            Hax(i,j) = i*k(i,j)*Fe*A/((i-N/2).^2 + (j-N/2).^2 + (z-q).^2).^(3/2); %The dike is centered at a position of (50,50)    
        end
    end
    if mod(q,1) == 0
    f = f+1;
    Fa = Za*sind(inc) + Hax*cosd(inc);
    xx=0:L/N:L;
    yy=0:L/N:L;
    [X,Y] = meshgrid(yy,xx);
    surf(X,Y,Fa);
    axis([0 L 0 L 0 6]);
    saveas(gcf, ['magn_',sprintf('%03d',f),'.png'])
    end
end
