﻿ Two electron atom program (Python)

### Python program to calculate two-electron atoms.

If you copy and paste the program source code below into a text editor, you can easily compile and run this.
(This program is Python ( 2.7 version ), so save this text editor as "filename.py" and compile it.)
In this program, first we input the initial x-coordinate r1 (in MM) of electron 1, and the absolute value of the total energy E (in eV) of helium atom.
From the inputted values, this program outputs the y component of electron 1 velocity after a quarter of its orbit, and WN (the number of de Broglie's waves included in one quarter of the orbital).
Here 1 SS = 1 × 10-22 second and 1 MM = 1 × 10-14 meter. Velocity 1 MM/SS = 1 × 108 m/s.
Python is a little slower than C and JAVA, so here we change the definition of 1 SS.
Orbits are a little rough, and calculation value becomes a little different from C and JAVA. ( Though almost same. )
So it's better to run C and JAVA programs, if you can.
( In this program, inputting energy |E| = 79.0054 and r1 = 3074.0 give WN = 0.250000. )
The initial x-coodinate is automatically increased per calculation until +50.

``````
import math
# input initial r1 and |E|
r=float(raw_input("r1 between nucleus and electron 1 (MM) ? "))
E=float(raw_input("total energy |E| of helium atom (eV) ? "))
me=9.1093826e-31;
pai=3.141592653589793; epsi=8.85418781787346e-12
h=6.62606896e-34; ele=1.60217653e-19

Z = 2.0;                  # helium atomic number
nucle = 6.64465650e-27;   # alpha particle
# rm = reduced mass
rm=(2.0*me*nucle)/(2.0*me+nucle); rm=rm*0.5

# repeat until r1 = r1 + 50
for num in range(1,50):
# initial potential energy
poten=-(2.0*Z*ele*ele)/(4.0*pai*epsi*r)+(ele*ele)/(4.0*pai*epsi*2.0*r)
# vya = total kinetic energy
vya=-(E*1.60217646e-19)-poten*1.0e14
if vya > 0:
vyb=math.sqrt(vya/me);  # vyb = initial velocity (m/s)
VY=vyb*1.0e-8;          # change m/s to MM/SS = initial VY
prexx=r; xx = r         # initial x coordinate = r
preyy=0.0               # initial y coordinate = 0
# WN = de Broglie wave number
VX=0.0; WN=0.0          # initial VX ( x velocity ) = 0

while xx > 0:           # repeat until  move 1/4 orbit
xx=prexx+VX; yy=preyy+VY   # coordinate after 1 SS
preVY=VY
vk=VX*VX+VY*VY
leng=math.sqrt(vk)*1.0e-14      # moving distance (m) for 1 SS
wav=h/(rm*math.sqrt(vk)*1.0e8)   # de Broglie wavelength (m)
WN=WN+leng/wav                 #  add de Broglie wave number for 1 SS

ra=math.sqrt(prexx*prexx+preyy*preyy) # distance nucleus and electron
rb=math.sqrt(4.0*prexx*prexx+2.0*preyy*preyy)  # between two electrons
ra=ra*1.0e-14; rb=rb*1.0e-14          # change MM to meter
prexx=prexx*1.0e-14; preyy=preyy*1.0e-14
ac=(ele*ele)/(4.0*pai*epsi*rm)
# add acceleration ( MM/SS^2 )
VX=VX+1.0e-30*ac*prexx*(-Z/(ra*ra*ra)+2.0/(rb*rb*rb))
VY=VY+1.0e-30*ac*preyy*(-Z/(ra*ra*ra)+1.0/(rb*rb*rb))
prexx=xx;preyy=yy
if VY > -0.0001 and VY < 0.0001:    # last VY condition
print "r1: "+str(r)+"   ",
print("VX:%.6f " % VX),
print("VY:%.6f " % VY),
print("preVY:%.6f " % preVY),
print("WN:%.6f " % WN)
r=r+1

``````