% 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

% Set initial inclination (degrees)
inc =90;
rad = 180/pi;
inc = (90-inc)/rad
% inc = 90 is a vertical B field


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;
M=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
irun =1
while irun ~= -1
 
  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 = [2+1.5*cos(inc),2-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*M*cos(inc-theta))/(r^3);
    b_theta(ix)=-(mu_0*M*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,':')
  hold off
  %plot(x,h_z_total,'-')
  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
  subplot (2,1,2)
  plot(ground_x,ground_z)
  axis([-20,20, -10,5])
  hold on
    plot(sphere_x,sphere_z)
    plot(bplotx1,bplotz1)
    plot(bplotx2,bplotz2)  
    plot(bplotx1(2),bplotz1(2),'.')
    plot(bplotx2(2),bplotz2(2),'.')
    xlabel ('x(m)')
    ylabel ('z(m)')
  hold off
  
  
% Change these parameters if desired 
 choice=menu('Options', 'Change i','Change depth','Quit');
  if choice ==1  
   90-inc*rad
   inc=input('Enter inc:')
   inc = (90-inc)/rad;
  end
  if choice ==2
     d
     d=input('Enter depth:')
  end
  if choice > 2
   irun=-1;
   print -dpsc  mag_dipole.ps
  end
end