下のソースプログラムをそのままテキストエディタ（メモ帳など）にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(この プログラムは Python なので、このテキストエディタを "ファイル名.py" とセーブしてコンパイルしてほしい。)
このサンプルプログラムでは、実行すると、最初に電子１のスタート地点のｘ座標（単位ＭＭ）、とヘリウム原子の全エネルギー（ eV )の絶対値（正の値）を入力するように画面に表示される。
それらを入力すると、このプログラムは、軌道計算後（電子が４分の１周した後）の状態における、電子１の速度のｙ成分（ゼロに近いとき、電子軌道は安定する）と、ＷＮ（4分の1軌道に含まれるド・ブロイ波の数）を画面に表示する。
ここでは、1 MM = 1 × 10-14 meter、 1 SS = 1 × 10-22 second、速度 1 MM/SS = 1 × 108 m/s で計算している。
Python は C や JAVA などと比較して 速度が少し遅いため 1 SS の定義を少し変更した。
そのため 計算軌道が 少し粗いものになったため、計算結果も C や JAVA のものと少し異なる ( ただほとんど同じだが。)
できるなら なるべく C や JAVA のプログラムで実行することをお勧めする。
( このプログラムでは エネルギー |E| = 79.0054、r1 = 3074.0 で WN = 0.250000 の結果がでる。)
最初の x 座標は 自動的に +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