下のソースプログラムをそのままテキストエディタ（メモ帳など）にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(このプログラムは 単純な C 言語なので、このテキストエディタを "filename.c" とセーブしてコンパイルしてほしい。)
このプログラムでは、最初に 換算質量の条件である 1 から 3 までの数を入力する。
"Normal" の条件 (= 1 ) では、重心が原点以外にあるとき 換算質量を用いる。
"NOT reduced mass" の条件 (= 2 ) では、必ず 通常の電子質量を用いる。
"All reduced mass" の条件 (= 3 ) では、重心が原点にあるときでさえ 必ず 換算質量を用いる。
次に 電子１の初期のｘ座標 r1 (MM) と ヘリウム の全エネルギー E の絶対値 (eV) を入力する。
ここでは、便宜のため、新しい単位を使っている。( 1 MM = 10-14 meter, 1 SS = 10-23 second )
それらを入力すると、電子が、1/4 周した後の電子１の速度 ( VX, VY in MM/SS ) と、軌道 1/4 周に含まれるド・ブロイ波の数 ( midWN ) が画面に表示される。
最初の x座標は 計算ごとに自動的に +30 まで増加していく。
#include <stdio.h>
#include <math.h>
int main(void)
{
int i;
double ree,r,E,rm;
double vya,vyb,poten,VX,VY,prexx,preyy,WN,ra,rb;
double xx,yy,vk,wav,preVY,preWN,midWN,leng,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 nucle = 6.64465650e-27;
double Z = 2.0;
/* input reduced mass condition */
printf("Input 1 (= Normal), 2 (= NOT reduced mass) or 3 (= All reduced mass) ? ");
scanf("%lf",&ree);
if (ree==1 || ree==2 || ree==3) {
/* input r1 and |E| */
printf("r1 between nucleus and electron 1 (MM)？ ");
scanf("%lf",&r);
printf("total energy |E| of helium atom (eV) ? ");
scanf("%lf", &E);
printf(" \n");
rm=(2.0*me*nucle)/(2.0*me+nucle); rm=rm*0.5; /* rm =reduced mass */
for (i=1; i < 30 ;i++) { /* repeat until r1=initial r1+100 */
/* poten = potential energy */
poten=-(2.0*Z*ele*ele)/(4.0*pai*epsi*r)+(ele*ele)/(4.0*pai*epsi*2.0*r);
/* vya= total E-potential energy */
vya=-(E*1.60217646e-19)-poten*1.0e14;
if (vya > 0) {
/* vyb=electron initial velocity (m/sec) */
if (ree==3) {vyb=sqrt(vya/rm);} /* All reducled mass -- ree=3 */
else {vyb=sqrt(vya/me);} /* initial states = usual electron mass */
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 */
preVY=VY;preWN=WN ;
vk=VX*VX+VY*VY;
leng=sqrt(vk)*1.0e-14; /* moving length (m) for 1 SS */
/* wav = de Broglie wavelength (m) */
wav=h/(rm*sqrt(vk)*1.0e9);
if (ree==2) { wav = h/(me*sqrt(vk)*1.0e9);} /* Not using reduced mass */
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 */
rb=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);
if (ree==2) {ac=(ele*ele)/(4.0*pai*epsi*me);} /* Not using reduced mass */
/* acceleration (MM/SS^2) */
VX=VX+1.0e-32*ac*prexx*(-Z/(ra*ra*ra)+2.0/(rb*rb*rb));
VY=VY+1.0e-32*ac*preyy*(-Z/(ra*ra*ra)+1.0/(rb*rb*rb));
prexx=xx;preyy=yy;
} while (xx >= 0); /* repeat above unitl electron 1 arrive at y axis */
if (VY > -0.0001 && VY < 0.0001) { /* last VY condition */
printf("r1= %.2f ", r );
printf("VX= %.6f ", VX);
printf("VY= %.6f ", VY);
printf("preVY= %.6f ", preVY);
midWN=(preWN+WN)/2.0; printf("midWN= %.6f\n", midWN);
}
} r=r+1;
}} return 0;
}