下のソースプログラムをそのままテキストエディタ（メモ帳など）にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(この プログラムは 単純な C 言語なので、このテキストエディタを "ファイル名.c" とセーブしてコンパイルしてほしい。)
このサンプルプログラムでは、実行すると、最初に 原子番号 Z ( 1 = 水素原子、2 = ヘリウムイオン ) と 換算質量の条件 ( 0 = 換算質量、1 = 通常の電子質量 ) を入力するように表示される。
これらの値から 水素様原子の 全エネルギーと、（ボーア）半径の理論値が出力される。
次に、電子 のスタート地点のｘ座標（単位ＭＭ）、と 水素様原子の全エネルギー （ eV ) の絶対値（正の値）を入力するように画面に表示される。
それらを入力すると、このプログラムは、軌道計算後（電子が４分の１周した後）の状態における、電子の速度のｙ成分（ゼロに近いとき、電子軌道は安定する）、y 座標、ＷＮ（4分の1軌道に含まれるド・ブロイ波の数）を画面に表示する。
ここでは、1 MM = 1 × 10-14 meter、 1 SS = 1 × 10-23 second で計算している。
各計算ごとに r1 は +100 まで自動的に増加していく。
ここの方法は 中性ヘリウムのプログラムと同じ方法を用いている。
水素原子 (H) では、 全エネルギー = -13.60569 eV、 ボーア半径 = 5291.77 MM。
換算質量を用いると、 全エネルギー = 13.59829 eV、半径 = 5294.65 MM。
ヘリウムイオン (He+) では、 全エネルギー = -54.42277 eV、半径 = 2645.89 MM。
換算質量を用いると、全エネルギー = -54.41531 eV、 半径 = 2646.25 MM。
これらの値 ( |E| = |全エネルギー|、r1 = 半径 )で、last VY はちょうどゼロになる。一度試してほしい。
#include <stdio.h>
#include <math.h>
/* hydrogen like atoms (ions) */
int main(void)
{
int i;
double Z,whi,toenergy,radius,rm,r,E;
double vya,vyb,poten,VX,VY,prexx,preyy,WN,ra;
double xx,yy,vk,preVY,preWN,midWN,leng,wav,ac;
double me=9.1093826e-31;
double pai=3.141592653589793;
double epsi=8.85418781787346e-12;
double h=6.62606896e-34;
double ele=1.60217653e-19;
double mp=1.67262171e-27;
double alph = 6.64465650e-27;
/* input atomic number Z and reduced mass condition */
printf("Atomic number Z ( H = 1, He+ = 2 )？ ");
scanf("%lf",&Z);
printf("You use reduced mass (= 0 ) or usual electron mass (= 1 ) ?");
scanf("%lf", &whi);
rm=me;
if (Z==1) { rm=(me*mp)/(me+mp); } /* hydrogen reduced mass */
if (Z==2) { rm=(me*alph)/(me+alph);} /* helium ion reduced mass */
if (whi == 1 ) { rm=me; }
/* theoretical values of total energy and (Bohr) radius */
toenergy = -(Z*Z*rm*ele*ele*ele*ele)/(8.0*epsi*epsi*h*h);
radius = (epsi*h*h)/(pai*rm*Z*ele*ele);
toenergy = toenergy * 6.241509e18; /* J to eV */
radius = radius * 1.0e14; /* meter to MM */
printf("Total energy: %.5f ", toenergy);
printf("Radius: %.2f \n", radius);
/* input r1 and |E| */
printf("r1 between nucleus and electron 1 (MM)？ ");
scanf("%lf",&r);
printf("total energy |E| of hydrogen like atoms (eV) ? ");
scanf("%lf", &E);
for (i=1; i < 100 ;i++) { /* repeat until r1=initial r1+100 */
/* poten = potential energy */
poten=-(Z*ele*ele)/(4.0*pai*epsi*r);
/* vya= total E-potential energy */
vya=-(E*1.60217646e-19)-poten*1.0e14;
if (vya > 0) {
/* vyb=electron initial velocity (m/sec) */
vyb=sqrt((2*vya)/rm);
VY=vyb*1.0e-9; /* change m/sec to MM/SS */
prexx=r; VX=0.0; WN=0.0; preyy=0.0;
do {
xx=prexx+VX; yy=preyy+VY; /* electron 1 position after 1SS */
preWN=WN ;
vk=VX*VX+VY*VY;
leng=sqrt(vk)*1.0e-14; /* moving length (m) for 1 SS */
wav=h/(rm*sqrt(vk)*1.0e9); /* de Broglie wavelength (m) */
WN=WN+leng/wav; /* add de Broglie wavelength */
/* calculation of VX,VY from Coulomb force */
ra=sqrt(prexx*prexx+preyy*preyy); /* between nucleus and electron */
ra=ra*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);
/* acceleration (MM/SS^2) */
VX=VX+1.0e-32*ac*prexx*(-Z/(ra*ra*ra));
VY=VY+1.0e-32*ac*preyy*(-Z/(ra*ra*ra));
prexx=xx;preyy=yy;
} while (xx >= 0); /* electron has moved one quater of an orbit? */
if (VY > -0.0001 && VY < 0.0001) { /* last VY condition */
printf("r1= %.2f ", r );
printf("VX= %.6f ", VX);
printf("VY= %.6f ", VY);
printf("last y = %.2f ", yy);
midWN=(preWN+WN)/2.0; printf("midWN= %.6f\n", midWN);
}
} r=r+1;
} return 0;
}