function [x,B] = mdipole(A, z, Be, inc)
% [x,B] = mdipole(A, z, Be, inc)
% on input:
%    A = dipole moment of dipole, in nT.m^3
%    Be = Earth's field magnitude, nT, always positive
%    inc = Earth field inclination, degrees, positive (down) in the northern hemisphere
%    z = depth of dipole, meters
% on output:
%    x = station position
%    B = magnetic field magnitude
% this version corrected for missing 4-pi (Feb 2024)

x = linspace(-20, 20, 401); % create a line of stations spaced 0.1 m
incr = pi*inc/180; % copy across inclination
r = sqrt(z.^2 + x.^2);
theta = acos(x./r);
% note that you need the '.' in the division (x and r are both vectors)
%figure;plot(x,theta); % you can make sure the angles work by plotting them
% if you do this you will see that (a) theta is always positive and (b) it 
% goes to zero at large positive x and 180 at large negative x.  
% it is equal to Inc + phi in the figure in the dipole notes. Thus...
phi = theta - incr;

% since the azimuth is measured clockwise from the positive pole, the 
% azimuth will be 180 + phi
azim = pi + phi;
%figure;plot(x,azim);
Br = A*2*cos(azim)./(4*pi*r.^3); % radial field from dipole equation
Bt = A*sin(azim)./(4*pi*r.^3);   % azimuthal field from dipole equation
Tr = -Be*cos(phi);   % Earth's field projected onto radial direction
Tt = Be*sin(phi);   % Earth's field projected onto azimuthal direction
B = sqrt((Br+Tr).^2 + (Bt+Tt).^2);
%figure;plot(x,Br,'r');
%figure;plot(x,Bt,'m');
%figure; plot (x, B-50000,'b')
