% tokamak.m
% Simulate a charge trajectory in a tokamak
% JMS, Feb.2023

% Set parameters
clear
Rmin = 2;                % minor torus radius (m)
Rmax = 4;                % major torus radius (m)
r0 = [3.5;0;0.3];        % initial position (m)
v0 = [1.5;0.2;2.0]*1e6;  % initial velocity (m/s) 
% v0 = 1.5e6*randn(3,1);   % initial velocity (m/s) 
dt = 1e-8;               % time step (s)
tmax = 1e-4;             % final time (s)
m = 1.6726e-27;          % proton mass (kg)
q = 1.6022e-19;          % proton charge (C)

% Find proton trajectory
lorenzForce = @(r,v,t)q*cross(v,magneticField_JMS(r));
r = newtonRK( lorenzForce, m, r0, v0, dt, tmax );

% Plot trajectory
figure(11)
dphi = 2*pi/30;
phi = 0:dphi:2*pi;
a = (Rmax+Rmin)/2;
b = (Rmax-Rmin)/2;
plot3(r(1,:),r(2,:),r(3,:),'-','LineWidth',1.5); hold on
plot3(Rmin*cos(phi),Rmin*sin(phi),  zeros(size(phi)),'k')
plot3(Rmax*cos(phi),Rmax*sin(phi),  zeros(size(phi)),'k')
plot3(   a*cos(phi),   a*sin(phi),+b*ones(size(phi)),'k')
plot3(   a*cos(phi),   a*sin(phi),-b*ones(size(phi)),'k')
hold off
axis equal
grid on
