下のソースプログラムをそのままテキストエディタ（メモ帳など）にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(このプログラムは 単純な C 言語なので、このテキストエディタを "filename.c" とセーブしてコンパイルしてほしい。)
ここでは 1 MM = 1 × 10-14 meter という新しい単位を使用している。
このサンプルプログラムでは、実行すると、最初に H2+ の核間距離 ( MM ) を入力する。
その後、電子１の ｙ座標 "r" ( MM )、x 座標 "a" ( MM ) を 入力する。
これら入力値から、このプログラムは H2 の結合エネルギー (eV)、力 F1 と遠心力の比、n2 に作用する 電子と n1 からの力の比 を出力する。
F1 は 電子に作用する 原子核 1 (n1) 方向の力である。
電子に作用する力と その速度から １軌道に含まれるドブロイ波を計算する。
( ここでは ビリアル定理 E = -T = 1/2 V を使用している。)
さらに 原子核 2 の力による x 軸周囲の歳差運動を考慮して このドブロイ波を計算する。
total waves は これらの合計である。
また 原子核 2 からの x 方向の力の影響を考慮した ドブロイ波も計算する (= new total waves )
#include <stdio.h>
#include <math.h>
int main(void)
{
int i,la,lb;
double nnuc,rr,aa,ke,r,nuc,a, rra,rrb,fynb,fxnb,internuc;
double ntwoforce,debrogli,wave,tim,vp,debropre,wavepre,cf,fgen,cosb;
double fnb,fone,goux,gouy,gouxx,gouxy,fonea,diele,fxnbb,gouxxx,foneb,rah;
double newgoux,newx,newradius,newwav,newvelo,newtim,newprewav,twaves;
double poten,ppot, kinetic,velo,binding;
double me=9.1093826e-31;
double pai=3.141592653589793;
double epsi=8.85418781787346e-12;
double h=6.62606896e-34;
double ele=1.60217653e-19;
/* input nuc, x, y coordinate */
printf(" Internuclear distance (MM) of H2 molecule ? ");
scanf("%lf",&nnuc);
printf(" r (MM) ? = y coordinate of electron 1 ");
scanf("%lf",&rr);
printf(" a (MM) ? = x coordinate of electron 1 ");
scanf("%lf", &aa);
printf(" \n");
r = rr * 1.0e-14; nuc = nnuc * 1.0e-14; /* change MM to meter */
ke = (ele*ele)/(4.0*pai*epsi);
for (i=1; i < 10 ;i++) { /* repeat until a=initial a+100 */
a=aa*1.0e-14;
rra=sqrt(a*a+r*r); /* rra =distance between electron and n1 */
rrb=sqrt(r*r+(nuc-a)*(nuc-a)); /* rrb= between electron and n2 */
poten=ke*(-1.0/rra-1.0/rrb+1.0/nuc); /* poten=potential energy (J) */
ppot=poten*6.241509*1.0e18; /* ppot = potential energy (eV) */
kinetic=-0.5*poten; /* kinetic = total kinetic energy of two electrons (J) */
velo=sqrt((2.0*kinetic)/me); /* velo = electron's velocity (m/s) */
binding=-ppot*0.5-13.606; /* binding energy (eV) of H2+ */
fynb=(ke*r)/(rrb*rrb*rrb); /* n2's force acting on electron in -y direction */
fxnb=(ke*(nuc-a))/(rrb*rrb*rrb); /* n2's force acting on electron in x direction */
internuc=ke/(nuc*nuc); /* force between two nuclei */
ntwoforce = fxnb/internuc; /* ratio of forces acting on n2 */
debrogli = h/(me*velo); /* de Broglie wavelength of electron */
wave=(2.0*pai*rra)/debrogli; /* de Broglie's waves around nucleus 1 */
tim=(2.0*pai*rra)/velo; /* tim = time period of electron around n1 */
vp=sqrt((fynb*r)/me); /* vp = electron precession veleocity around x axis */
debropre = h/(me * vp); /* de Broglie wavelength of precession */
wavepre=(vp*tim)/debropre; /* de Broglie's wave number of precession */
cf=(me*velo*velo)/rra; /* cf= centrifugal force around n1 */
fgen=ke/(rra*rra); /* fgen = n1's force acting on the electron */
cosb=(rra*rra+nuc*nuc-rrb*rrb)/(2*rra*nuc); /* cosine b */
fnb=ke/(rrb*rrb);
fnb=(fnb*(rra-nuc*cosb))/rrb; /* fnb= n2's force * cos (theta) */
fone=fnb+fgen; /* total force acting on electron toward n1 (=F1) */
goux=fone*cosb; /* x component of F1 */
gouy=fone*sqrt(1-cosb*cosb); /* y component of F1 */
gouxx=goux-fxnb; /* subtract n2 force from goux */
fonea=sqrt(gouxx*gouxx+gouy*gouy); /* changed F1 */
diele=sqrt((nuc+a)*(nuc+a)+r*r); /* between n2 and electron at position 2 */
fxnbb=(ke*(nuc+a))/(diele*diele*diele); /* x component of force from n2 */
gouxxx=goux+fxnbb; /* add n2 force to goux */
foneb=sqrt(gouxxx*gouxxx+gouy*gouy); /* changed F1 at position 2 */
fonea=(fonea+foneb)/2.0; /* changed average F1 */
rah=fone/cf; /* rah = ratio of F1 to centrifugal force around n1 */
newgoux = sqrt(fonea*fonea-gouy*gouy);
newx=(r*newgoux)/gouy;
newradius=sqrt(r*r+newx*newx); /* radius influened by n2 force in x direction */
newvelo= sqrt((fonea*newradius)/me);
newtim=(2.0*pai*newradius)/newvelo; /* new time period */
newprewav=(vp*newtim*me*vp)/h; /* new de Broglie's wave of precession */
newwav=(2.0*pai*newradius*me*newvelo)/h; /* new de Broglie's waves number */
printf("a:%.1f", aa);
printf(" binding energy (eV): %.4f ", binding);
printf("ratio of F1 to centrifugal: %.4f ", rah);
printf("ratio of forces acting on n2: %.3f \n", ntwoforce);
printf("waves around n1: %.3f ", wave);
printf("+ waves around x axis: %.3f ", wavepre);
la=(int)(1000.0*wave); lb=(int)(1000.0*wavepre);la=la+lb;
twaves=la/1000.0;
printf(" = total waves: %.3f \n ", twaves);
printf("new waves around n1: %.3f ", newwav);
printf("+ new precession wave: %.3f ", newprewav);
la=(int)(1000.0*newwav); lb=(int)(1000.0*newprewav);la=la+lb;
twaves=la/1000.0;
printf(" = new total waves: %.3f \n ", twaves);
printf(" \n");
aa=aa+10;}
return 0;
}