C language program to calculate the Helium orbits (manipulating reduced mass).

top page
Back to no helium compound page

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 redeuced mass except when the center of mass is at the origin.
In "Not reduced mass" (= 2 ) condition, we always use the usual mass.
In "All reduced mass" (= 3 ) condition, we always use the reduced mass even when the center of mass is at the origin.

Next 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. From the inputted values, this program outputs the y component of electron 1 velocity in Fig. 13, and WN (the number of de Broglie's waves contained in one quarter of the orbital). Here 1 SS = 1 × 10-23 second.
The initial x-coodinate is automatically increased per calculation until +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;
   }