clc

clear all

%% Hard Sphere Collision Routine by Michael N. Macrossan

%function [vx,vy,vz] = HS(vx,vy,vz,plist,badd,N_in);
%global sigdtvol c_in TN_in nsteps gmax % sigdtvol = sigma*dt/vol
%debug = 0; if debug; dbstop if error; end
%Nc = length(badd);
%TN_in = TN_in + N_in; nsteps = nsteps + 1; Nbar = TN_in/nsteps;
%Nsel = floor(0.5*N_in.*Nbar.*gmax*sigdtvol + rand(1,Nc)); % vectorised?
%% note sigma*dt/vol constant = sigdtvol
%for cell = 1:Nc
%N = N_in(cell);
%if N>1 % only for at least 2 particles in cell
%%Nsel = floor(0.5*N*Nbar(cell)*gmax(cell)*sigma*dt/vol + rand); %
%%fractional collision done randomly
%if debug & Nsel(cell)<1
%hs: cell,Nsel ,cell,Nsel(cell)
%N,sigma*dt/vol ,N,sigdtvol
%return
%end
%b0 = badd(cell);
%gm = 0.0; % detects the largest value at each time step
%Ncol = 0;
%if (debug); sumg = 0.0; end
%for sel = 1:Nsel(cell)
%%’ sel = ’,sel
%k1 = ceil(rand*N); % first collision particle
%k2 = k1; while k2 == k1; k2 = ceil(rand*N); end
%p1 = plist(b0 + k1); p2 = plist(b0 + k2);
%gx = vx(p1) - vx(p2); % collision velocity
%gy = vy(p1) - vy(p2); % = rel. velocity
%gz = vz(p1) - vz(p2);
%g = sqrt(gx^2 + gy^2 + gz^2);
%if debug; sumg = sumg + g; end
%gm = max([gm,g]);
%if g/gmax(cell) > rand % accept this pair
%c_in(cell) = c_in(cell) + 1.0;
%if debug; Ncol = Ncol + 1; end
%[vx(p1),vy(p1),vz(p1),vx(p2),vy(p2),vz(p2)] = ...
%newvel(vx(p1),vy(p1),vx(p1),vx(p2),vy(p2),vz(p2),g)
%end % collision
%end % selection
%gmax(cell) = max( [1.1*gm, 0.95*gmax(cell)] );
%% use the latest largest value found, with margin for error.
%% if gmax is decreasing, use smaller value 0.95 factor,
%% but what if low gm found by "chance"?
%end % N > 1
%end % cells
%return


%% Constant quantities

N = 10000; %Number of particles
dt = 1*10^(-5); %time steps (s)
i = 1; % unit vector in x dir.
j = 1; % unit vector in y dir.
k = 1; % unit vector in z dir.
M = 2.21 *10^(-27) % mass of ceasium atom in kg
w = 9.2*10^(7) % angular oscillation frequency of ceasium (Hz)
%B = ; % Space boundary

% Cell Parameters
h = 1*10^(-6);
x = h*i;
y = h*j;
z = h*k;
Vc = h^3; %volume of space
cell = x*y*z;

% Particle Information
r_init = linspace(0,N,1); %position of each particle
v_init = linspace(0,N,1); %velocity of each particle

for i = 0:N
r_new(i) = r_init(i) + v_init(i)*dt

%% POTENTIAL

U = 0.5*M*(w_x^2 * x^2 + w_y^2 * y^2 + w_z^2 * z^2 )

for i = 1:N

v_new(i) = (1/M)*U*dt



%% Constant quantities

N = 10000; %Number of particles
dt = 1*10^(-5); %time steps (s)
i = 1; % unit vector in x dir.
j = 1; % unit vector in y dir.
k = 1; % unit vector in z dir.
M = 2.21 *10^(-27) % mass of ceasium atom in kg
w = 9.2*10^(7) % angular oscillation frequency of ceasium (Hz)
%B = ; % Space boundary

% Cell Parameters
h = 1*10^(-6);
x = h*i;
y = h*j;
z = h*k;
Vc = h^3; %volume of space
cell = x*y*z;


r_init = linspace(0,N,1); %position of each particle
v_init = linspace(0,N,1); %velocity of each particle

for i = 1:N
r_init(i) =
v_init(i) =

for i = 1:N
r_new(i) = r_init(i) + v_init(i)*dt



%% POTENTIAL

U = 0.5*M*(w_x^2 * x^2 + w_y^2 * y^2 + w_z^2 * z^2 )

for i = 1:N

v_new(i) = (1/M)*U*dt