% protonTrajectory.m
% Solution to problem 2 of Computation control 2. 2017/02/13
% J.M.Soler, Feb.2017

% Set parameters
clear all
m = 1.672e-27;       % proton mass (kg)
q = 1.602e-19;       % proton charge (C)
r0 =  [1;2;3]*1e8;   % initial position (m)
v0 = -[1;2;3]*1e7;   % initial velocity (m/s)
tmax = 30;           % trajectory duration (s)
dt = 0.1;            % trajectory time interval (s)

% Integrate trajectory
myForce = @(r,v,t)earthMagneticForce(q,r,v);
r = rungeKutta( myForce, m, r0, v0, dt, tmax );
r = squeeze(r);      % eliminate particle index

% Plot trajectory (and a dot at Earth's position)
plot3(r(1,:),r(2,:),r(3,:),0,0,0,'o')
xlabel('x (m)')
ylabel('y (m)')
zlabel('z (m)')
axis equal
grid on
