% Computes magnetic field variation across a dipole induced in
% the earth

% SI units used.  Both Earth and induced fields in nT

clear all
close all



d=4;     % Depth of centre of sphere (m)
B=50000; % Strength of Earth's field (nT)
k=0.015; % Susceptibility of sphere (SI)
a =2;    % Radius of sphere (m)
vol = 4*pi*a*a*a/3;
mu_0 = 4*pi*1e-7;
DM=vol*k*B/mu_0   % Magnetic moment of sphere (ignores demagnetizing effect)

%Define vectors for geometry plot

  ground_x = [-20,20]
  ground_z = [0,0]
  b_bg =[B,B]


% Start loop to vary inclination and depth
n=360;
M=moviein(n);
for it=1:n
   
  inc = it*2*pi/n;
  b_earth_z =  -B*cos(inc);
  b_earth_x =  -B*sin(inc);
  
  % Define vector for plotting inclination 
  
  % M in sphere
  bplotx1 = [0.5*a*sin(inc),-0.5*a*sin(inc)];
  bplotz1 = [-d+0.5*a*cos(inc),-d-0.5*a*cos(inc)];
  
  % Earth's field at surface
  bplotx2 = [10+1.5*sin(inc),10-1.5*sin(inc)];
  bplotz2 = [8+1.5*cos(inc),8-1.5*cos(inc)];

  npts =100;
  for ipts = 1:npts+1
   sphere_x(ipts) = -a*cos(2*pi*(ipts-1)/npts);
   sphere_z(ipts) = -d-a*sin(2*pi*(ipts-1)/npts);   
  end
  
  % Define points where H is computed on z=0 surface
  ix=0;
  for xdum = -20:0.1:20.0
    ix=ix+1;
    x(ix) =xdum;
  end
  nx=ix;

  for ix =1:nx
    r=sqrt(d*d+x(ix)*x(ix));
    theta = atan(x(ix)/d);
    % Dipole fields
    b_r(ix)=-(mu_0*2*DM*cos(inc-theta))/(r^3);
    b_theta(ix)=-(mu_0*DM*sin(inc-theta))/(r^3);
    
    b_z(ix)= b_r(ix)*cos(theta)+b_theta(ix)*sin(theta);
    b_x(ix)= b_r(ix)*sin(theta)-b_theta(ix)*cos(theta);
 
    % Sum of dipole and Earth's fields  
    b_x_total(ix)=b_x(ix)+b_earth_x;
    b_z_total(ix)=b_z(ix)+b_earth_z;
    b_total(ix) = sqrt(b_x_total(ix)*b_x_total(ix)+b_z_total(ix)*b_z_total(ix));
  end
  
  % PLot magnetic field response
  %subplot (2,1,1)
  plot(x,b_total)
  hold on
  plot(ground_x,b_bg,':')
  xlabel ('x(m)')
  ylabel ('B_T (nT)')
  title ('Total magnetic field over a sphere :  k=0.015')
  axis([-20,20,49000,51000])
  
  %Plot model geometry and inclination
  % In movie, subplots don't work, so rescale
  
  vcent = 49700;
  vscale =75.;
  plot(sphere_x,vcent+vscale*sphere_z)
  plot(bplotx1,vcent+vscale*bplotz1)
  plot(bplotx2,vcent+vscale*bplotz2)  
  plot(bplotx1(2),vcent+vscale*bplotz1(2),'.')
  plot(bplotx2(2),vcent+vscale*bplotz2(2),'.')
  hold off
  M(:,it) = getframe;  
end

movie(M,50)