function B = magneticField_JMS(r)
% Simple model of the magnetic field of a tokamak
% Input:
%   r(3,1) : cartesian coordinates at which the field is calculated (m)
% Output:
%   B(3,1) : magnetic field at given points (T)
% JMS, Feb.2023

% Set parameters
mu0 = 4*pi*1e-7;          % magnetic permitivity of vacuum (N/A^2)
It = 5e6;                 % total toroidal current (A)
M = [0;0;1.2e7];            % total magnetic moment of central solenoid (A*m^2)

% Find magnetic field
rmod = norm(r);           % distance to torus center
rdir = r/rmod;            % unit vector along r
R = norm(r(1:2));         % distance to torus axis
B = mu0*It/(2*pi*R) * cross(rdir,[0;0;1]) + ...       % toroidal field
    mu0/4/pi * (3*dot(M,rdir)*rdir - M)/rmod^3;       % poloidal field

end
