If you copy and paste the program source code below into a text editor, you can easily compile and run this.
(This program is simple C language, so save this text editor as "filename.c", and compile it.)
In this program, we first input the number ( 1-3 ) of reduced mass conditions.
In "Normal" condition (= 1 ), we use the reduced mass except when the center of mass of two electrons is just at the nucleus (= use reduced mass except when at the start where two electrons just start from x-axis which uses ordinary electron mass because the helium nucleus transiently stops only at the start ).
In "Not reduced mass" (= 2 ) condition, we always use the ordinary electron mass.
In "All reduced mass" (= 3 ) condition, we always use the reduced mass even when the helium nucleus transiently stops
In this program, first we input the initial x-coordinate r1 (in MM) of electron 1 (= the electron-1 and electron-2 start from x-axis ), and the absolute value of the total energy E (in eV) of helium atom.
↑ This initial input r1 value increases per calculation of the electron's 1/4 orbit until r1 + 30.
From the inputted values, this program computes electrons' 1/4 orbits (= electron-1 moves from x-axis to y-axis on x-y place, while electron-2 moves from x-axis to z-axis on the perpendicular x-z plane symmetrically ) by rigorously calculating Coulomb forces among two electrons and a nucleus, de Broglie wavelength, updating electrons' positions, velocities at short time intervals of 1SS (= 1 × 10-23 second ).
Then, this program outputs the x (= VX ) and y (= VY ) components of electron-1's last velocity after the electron has moved a quarter of its orbits from x-axis to y-axis (= electron-1 must cross the y-axis perpendicularly with the last VY velocity being zero to form the most stable symmetrical orbit around the nucleus = so we should choose the uniquely-determined r1 which gives the zero last VY in each different helium total energy E ).
This program also outputs WN (= the number of de Broglie's waves included in one quarter of the orbital, which should be 0.250000 when one orbit is 1 × de Broglie wavelength, which case uniquely gives the predicted helium energy by this real helium model's program ).
Here we use new convenient units of time = 1 SS = 1 × 10-23 second, distance = 1 MM = 1 × 10-14, velocity = 1MM/SS = 1 × 109 m/s, acceleration = 1MM/SS^2 = 1032 m/s^2.
The initial x-coordinate is automatically increased per calculation until +30 (= these different initial x-coordinates give different orbits, out of which we can uniquely determine only one initial x-coordinate that gives a symmetrical orbit around a nucleus with zero y-component of electron-1's last velocity (= VY = 0.000000 ) after electron-1 moves a quarter of its orbit ).
#include <stdio.h>
#include <math.h>
int main(void)
{
/* use reduced mass or electron mass */
int i;
double r,E,rm, ree, nucle;
double vya,vyb,poten,VX,VY,prexx,preyy,WN,ra,rb;
double xx,yy,vk,leng,wav,ac;
/* me = electron mass kg, pai = pi,
epsi= electron constant, h = Planck constant,
ele = electron charge, Z = nuclear charge
*/
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 Z = 2.0;
/*
choose "always use reduce mass = 3",
"always use ordinary electron mass = 2"
"Normal = 1 = use reduced mass except when two electrons are just at the opposite sides of the nucleus (= use ordinary elecrton mass only at starting point where electrons-1 and 2 just start from x-axis )"
*/
printf("Input 1 (= Normal), 2 (= NOT reduced mass) or 3 (= All reduced mass) ? ");
scanf("%lf",&ree);
/* input initial electron-1's x-coordinate r1 (MM = 10^-14 meter ) */
printf("r1 between nucleus and electron 1 (MM)？ ");
scanf("%lf",&r);
/* input absolute value of helium total energy (eV) */
printf("absolute value of total energy |E| of helium atom (eV) ? ");
scanf("%lf", &E);
/* nucle = nuclear mass (kg),
rm = reduced mass
*/
nucle = 6.64465650e-27;
rm =(2.0*me*nucle)/(2.0*me+nucle);
rm=rm*0.5;
/* when Not using reduced mass */
if ( ree > 1.5 && ree < 2.5 ) { rm = me ;}
for (i=1; i < 30 ;i++) {
/* repeat until initial r1 increases to r1+30 */
/* poten = initial Coulomb potential energy (J) */
poten=-(2.0*Z*ele*ele)/(4.0*pai*epsi*r*1.0e-14)+(ele*ele)/(4.0*pai*epsi*2.0*r*1.0e-14);
/* vya= total energy E - potential energy (J) */
vya=-(E*1.60217646e-19)-poten;
if (vya > 0) {
/* vyb=electron's initial velocity (m/sec) */
vyb=sqrt(vya/me);
/* when always using reduced mass */
if ( ree > 2.5 ) { vyb=sqrt(vya/rm); }
/* VX, VY = electron-1's x,y-velocity (MM/SS = 10^9 m/s ) */
VY=vyb*1.0e-9; VX = 0.0;
/* prexx, preyy = initial electron-1's x,y-coordinate */
prexx=r; preyy=0.0;
/* WN = number of de Broglie wavelength */
WN=0.0;
do {
/* xx,yy = electron-1's coordinate after 1SS (= 10^-23 second ) */
xx=prexx+VX; yy=preyy+VY;
vk=VX*VX+VY*VY;
/* leng = electron's moving length (= meter ) for 1SS */
leng=sqrt(vk)*1.0e-14;
/* wav = de Broglie wavelength = h/mv (= meter ) */
wav=h/(rm*sqrt(vk)*1.0e9);
/* WN = sum of de Broglie wavelength */
WN=WN+leng/wav;
/* ra = distance between electron and nucleus */
ra=sqrt(prexx*prexx+preyy*preyy);
/* rb = distance between two electrons */
rb=sqrt(4.0*prexx*prexx+2.0*preyy*preyy);
/* change unit of length from MM into meter */
ra=ra*1.0e-14; rb=rb*1.0e-14;
prexx=prexx*1.0e-14; preyy=preyy*1.0e-14;
ac=(ele*ele)/(4.0*pai*epsi*rm);
/* calculating acceleration 1MM/SS^2 = 10^32 m/s^2 */
/* update x-velocity VX of electron-1 from Coulomb force */
VX=VX+1.0e-32*ac*prexx*(-Z/(ra*ra*ra)+2.0/(rb*rb*rb));
/* update y-velocity VY of electron-1 from Coulomb force */
VY=VY+1.0e-32*ac*preyy*(-Z/(ra*ra*ra)+1.0/(rb*rb*rb));
prexx=xx;preyy=yy;
} while (xx >= 0);
/* until electron-1 has moved its 1/4 orbit to reach y-axis */
/* display initial electron-1's x-coordinate r1 */
printf("r1= %.2f ", r );
/* display electron-1's last velocity VX, VY */
printf("VX= %.6f ", VX);
printf("VY= %.6f ", VY);
/* display total de Broglie wavelength WN of 1/4 orbit */
printf("WN= %.6f\n", WN);
} r=r+1;
} return 0;
}