Equation of Motion of Proton (Coulomb Potential)

## A function of the solution of the path
## equation of a proton under the inverse
## square atraction field of an electron
## that takes in the initial seperation distance
## or the position r0 wrt center of the electron
## as an argument and returns the total time
## ttotal it takes for it to reach within 1.0m
## of the electron and plots the velocity v
## vs time graph. usage ttotal=coulomb(r0).
function coulomb(r0)
## constants and initializations
k=9e+9; ## coulomb field constant [Nm^2/C]
q=1.6e-19; ## electronic charge [C]
mp=1.7e-27; ## proton mass [m]
dt=1e-4;; ## increment in time [sec]
r=r0; ## initial seperation[m]
t=0; ## initial time [s]
v0=0; ## initial velocity [m/s]
v=v0; ## 1st entry of v array [m/s]
n=1; ## initialization of loop index
## since we have already n=0 in argument r0.
while(r(n)>0.1);
dr=v(n)*dt; ## increment that is decrease in r
## since v(n) will be negative below.
r=[r;r(n)+dr]; ## decreases r in each step and
## accumulates the results in r array.
dv=-k*q*q*dt/(mp*r(n)*r(n)); ## increment in v to make v more negative.
v=[v;v(n)+dv]; ## increases the magnitude of negative v in each step and
## accumulates the results in v array.
t=[t;t(n)+dt]; ## increases the time in each step and
## accumulates the results in t array
## for both time axis and ttotal.
n++; ## increase n by 1 in each step
endwhile
ttotal=t(n) ## print the last entry of t array gives total time.
plot(t,v,';;')
xlabel('time(sec)'); ## Ola yuppi! I've learned eventually
ylabel('velocity(m/sec)'); ## to label the axes :) happy end!
endfunction