Skip to main content

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

Comments

Popular posts from this blog

NEWTON-RAPSON METHOD FOR HEAT FLOW

##Constants and initializations
a=5.67E-8; ## Stefan-Boltzman constant[Watt/meter^2Kelvin^4]
e=0.8; ## Rod surface emissivity [Dimensionless]
h=20; ## Heat transfer coefficient of air flow [W/m^2-K]
Tinf=Ts=25; ## Temperature of air and the walls of the close[Celcius]
D=0.1; ## Diameter of the rod[meter]
I2R=100; ## Electric power dissipated in rod (Ohmic Heat)[W]
T=[]; ## Temperature of the rod[*C]
T(1)=25; ## Initial guess of the temperature of the rod[*C]
Q=[]; ## Heat function [W]
Qp=[]; ## First derivative of Q wrt T [W/C*].
for i=1:100
Q(i)=pi*D*(h*(T(i)-Tinf)+e*a*(T(i)^4-Ts^4))-I2R;
Qp(i)=pi*D*(h+4*e*a*T(i)^3);
T(i+1)=T(i)-Q(i)/Qp(i); ## Newton-Rapson Method
endfor
printf('The steady state temperature is %f\n',T(i+1))
save -text HeatFlowTemp.dat
## The plot
t=1:100; ##temperature
for n=1:100
H(n)=pi*D*(h*(t(n)-Tinf)+e*a*(t(n)^4-Ts^4))-I2R;
endfor
plot(t,H)
xlabel('T(Celcius)');
ylabel('Q(Watt)');
legend('Q(T)');
title('Heat flow vs Temperatu…