clear all
clf

f = 0;

d = 10;
L = 9*d;
n = 60;
N = 9*n;
nHours = 12*30*24;

K = 21*10^(-5);
dt = 60; 
dx = d/n;
eta = dt*K/dx^2;

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); %Earth's Field Intensity
z = 20; %depth to top of dike

C = 1*10^5; %Curie constant is material dependent
mu = 1.257*10^(-6); %permeability of free space

k = zeros(1,N+1); %magnetic susceptibility
T = 20*ones(1,N+1); %Temperature

for i = (4*N/9 + 1):(5*N/9 + 1)
    T(i) = 1100;
    k(i) = C*mu./(T(i)-C*mu);
end
for i = (2*N/9 + 1):(3*N/9 + 1)
    T(i) = 1100;
    k(i) = C*mu./(T(i)-C*mu);
end
for i = (7*N/9 + 1):(8*N/9 + 1)
    T(i) = 600;
    k(i) = C*mu./(T(i)-C*mu);
end

x = zeros(1,N+1);
t = zeros(1,N+1);
 
for i = 1:N+1
    x(i) = (i-1)*dx;
end

Za1 = zeros(1,N+1);
Hax1 = zeros(1,N+1);
Za2 = zeros(1,N+1);
Hax2 = zeros(1,N+1);
Za3 = zeros(1,N+1);
Hax3 = zeros(1,N+1);

Numit = nHours;
    for iter = 1:Numit
        Tnew = T;
        knew = k;
        for i = 2:N
            Tnew(i) = eta*(T(i+1) - 2*T(i) + T(i-1))+T(i);
            if (i >=(4*N/9+1)&& i<=(5*N/9 + 1))
                knew(i) = C*mu./(Tnew(i)-C*mu);
                for j = 1:N+1
                Za1(j) = z*knew(i)*Fe*10/((j-N/2).^2 + z.^2).^(3/2);
                Hax1(j) = j*knew(i)*Fe*10/((j-N/2).^2 + z.^2).^(3/2);
                end
            end
            if (i >=(2*N/9+1)&& i<=(3*N/9 + 1))
                knew(i) = C*mu./(Tnew(i)-C*mu);
                for j = 1:N+1
                Za2(j) = z*knew(i)*Fe*10/((j-10/9*N/4).^2 + z.^2).^(3/2);
                Hax2(j) = j*knew(i)*Fe*10/((j-10/9*N/4).^2 + z.^2).^(3/2);
                end
            end
            if (i >=(7*N/9+1)&& i<=(8*N/9 + 1))
                knew(i) = C*mu./(Tnew(i)-C*mu);
                for j = 1:N+1
                Za3(j) = z*knew(i)*Fe*10/((j-3*10/9*N/4).^2 + z.^2).^(3/2);
                Hax3(j) = j*knew(i)*Fe*10/((j-3*10/9*N/4).^2 + z.^2).^(3/2);
                end
            end
        end
        T = Tnew;
        k = knew;
        Fa = (Za1+Za2+Za3)*sind(inc) + (Hax1+Hax2+Hax3)*cosd(inc); %induced field intensity
        if mod(iter+49,50) == 0
        f = f+1;
        subplot(3,1,1)
        plot(x,T,'-b')
        axis([0 90 0 1500])
        title('Temperature','fontsize', 20)
        xlabel('position (m)')
        ylabel('Temperature (C)')
        pause(.1)
        subplot(3,1,2)
        plot(x,k,'-g')
        axis([0 90 0 10*10^(-4)])
        title('Magnetic Susceptibility','fontsize', 20)
        xlabel('position (m)')
        ylabel('Magnetic Susceptibility (emu)')
        pause(.1)
        subplot(3,1,3)
        plot(x,Fa,'-r')
        axis([0 90 0 10])
        title('Field Intensity','fontsize', 20)
        xlabel('position (m)')
        ylabel('Field Intesity (Nanotesla)')
        saveas(gcf, ['mag_',sprintf('%03d',f),'.png'])
        pause(.1)
        end
    end
    
    




    

