If you copy and paste the program source code below into a text editor, you can easily compile and run this.
(This class file name is MathMethod, so save this text editor as "MathMethod.java", and compile it.)
In this program, first we input the initial x-coordinate r1 (in MM) of electron 1, and the absolute value of the total energy E (in eV) of helium atom.
From the inputted values, this program computes electrons' orbits by calculating Coulomb force among two electrons and a nucleus, de Broglie wavelength at short time intervals.
Then, this program outputs the y component of electron 1 velocity (= which becomes zero when an orbit is symmetrical around a nucleus ) after the electron has moved a quarter of its orbits, and WN (= the number of de Broglie's waves included in one quarter of the orbital, which should be 0.25000 when one orbit is 1 × de Broglie wavelength, which case uniquely gives the predicted helium energy by this real helium model's program ).
Here 1 SS = 1 × 10-23 second and 1 MM = 1 × 10-14 meter.
The initial x-coodinate is automatically increased per calculation until +100 (= these different initial x-coodinates 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 velocity after electron-1 moves a quarter of its orbit ).
import java.util.Scanner;
class MathMethod {
public static void main(String[] args) {
// input r1 and |E|
Scanner stdIn=new Scanner(System.in);
System.out.println("r1 between nucleus and electron 1 (MM)？ ");
double r=stdIn.nextDouble();
System.out.println("total energy |E| of helium atom (eV) ? ");
double E=stdIn.nextDouble();
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; // Helium atomic number
double nucle = 6.64465650e-27; // He alpha particle
double rm=(2.0*me*nucle)/(2.0*me+nucle); rm=rm*0.5; // reduced mass
for (int i=1;i < 100;i++) { // repeat until r1=initial r1+100
// poten = potential energy
double 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
double vya=-(E*1.60217646e-19)-poten*1.0e14;
if (vya > 0) {
// vyb=electron initial velocity (m/sec)
double vyb=Math.sqrt(vya/me);
double VY=vyb*1.0e-9; // change m/sec to MM/SS
double prexx=r; double VX=0.0; double WN=0.0; double preyy=0.0;
double yy,vk,preVY,preWN,midWN,leng,wav; double xx=0.0;
do {
xx=prexx+VX; yy=preyy+VY; //electron 1 position after 1SS
preVY=VY;preWN=WN ;
vk=VX*VX+VY*VY;
leng=Math.sqrt(vk)*1.0e-14; // moving length (m) for 1 SS
wav=h/(rm*Math.sqrt(vk)*1.0e9); // de Broglie wavelength (m)
WN=WN+leng/wav; // add de Broglie wavelength
//calculation of VX,VY from Coulomb force
double ra=Math.sqrt(prexx*prexx+preyy*preyy); // between nucleus and electron
double rb=Math.sqrt(4.0*prexx*prexx+2.0*preyy*preyy); // between two electrons
// change MM to meter
ra=ra*1.0e-14; rb=rb*1.0e-14;
prexx=prexx*1.0e-14; preyy=preyy*1.0e-14;
double ac=(ele*ele)/(4.0*pai*epsi*rm);
// 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); // electron has moved one quater of an orbit?
if (VY > -0.0001 && VY < 0.0001) { // last VY condition
System.out.print("r1: "+r+" ");
System.out.printf("VX:%.6f", VX);
System.out.printf("VY:%.6f", VY);
System.out.printf("preVY:%.6f", preVY);
midWN=(preWN+WN)/2.0; System.out.printf("midWN:%.6f\n", midWN);
}
} r=r+1;
}}}