c     this program computes the Bohr model Hbeta level energies,
c     transition energy, and wavelength; it hten computes the same for
c     the Bohr-Sommerfeld model.  a table of the fractional energy
c     corrections for the split and shifted B-S levels are output.

      PROGRAM hw5p1

      implicit none
      double precision pi,e,me,mp,mn,c,Ryd,h,hbar,rbohr,fine
      double precision m_H,m_He,mu_H,mu_He,ratH,ratHe,RydH,RydHe
      double precision a0,aH,aHe

      integer           n,k
      double precision  E4,E2,EHb,wHb
      double precision  dEnk,nn,kk
      dimension         dEnk(5,5)

      double precision  E41,E42,E43,E22,E21
      double precision  EHb4122,EHb4221,EHb4322
      double precision  wHb4122,wHb4221,wHb4322 

c     define physical constants [cgs] always

      pi   = 3.1415926535897932384d0
      e    = 4.803206815d-10
      me   = 9.109389754d-28
      mp   = 1.67262171d-24
      mn   = 1.67492728d-24
      c    = 2.99792458d10
      h    = 6.6260755d-27
      Ryd  = 2.1798741d-11
   
      hbar  = 0.5d0 * h/pi

      rbohr = hbar*hbar/(me * e * e)
      fine  = e*e/(hbar*c)  

c     compute the Rydnerg constant for both the inifnite mass nucleus
c     and the reduced mass; account for reduced mass and nuclear charge
c     for yucks, He is also computed (but not used)

      m_H  = mp
      m_He = 2.0d0*mp + 2.0d0*mn  

      mu_H  = me/(1.0d0+me/m_H)
      mu_He = me/(1.0d0+me/m_He)
      ratH  = mu_H/me
      ratHe = mu_He/me
      RydH  = Ryd*ratH
      RydHe = 4.0d0*Ryd*ratHe

      a0    = rbohr
      aH    = rbohr/ratH
      aHe   = (rbohr/ratHe)/4.0d0

c     write out the reduced mass ratio and Rydberg constants for H

      WRITE(6,*) 'muH/me   = ',ratH
      WRITE(6,*) 'Ryd_inf  = ',Ryd
      WRITE(6,*) 'RydH     = ',RydH

c     compute Bohr Hbeta (infinite mass nucleus)

      WRITE(6,*)
      WRITE(6,*) 'Unreduced Bohr Hbeta'

      E4 = -Ryd/(4.0d0**2)
      E2 = -Ryd/(2.0d0**2)
      EHb = E4 - E2
      wHb = 1.0d8*h*c/Ehb

      WRITE(6,*) 'E_4      = ',E4  ! n=4 energy
      WRITE(6,*) 'E_2      = ',E2  ! n=2 energy
      WRITE(6,*) 'EHb      = ',EHb ! Hbeta trani energy
      WRITE(6,*) 'wHb      = ',wHb ! Hbeta wavelength [anstroms]

c     compute Bohr Hbeta (reduced mass of electron)

      WRITE(6,*)
      WRITE(6,*) 'Reduced (best) Bohr Hbeta'

      E4 = -RydH/(4.0d0**2)
      E2 = -RydH/(2.0d0**2)
      EHb = E4 - E2
      wHb = 1.0d8*h*c/Ehb

      WRITE(6,*) 'E_4      = ',E4  ! n=4 energy
      WRITE(6,*) 'E_2      = ',E2  ! n=2 energy
      WRITE(6,*) 'EHb      = ',EHb ! Hbeta trani energy
      WRITE(6,*) 'wHb      = ',wHb ! Hbeta wavelength [anstroms]

c     compute and communicate the matrix of the energy corrections for
c     the B-S model since it is for H, Z=1 (implicit).  calculation is
c     approximate to order (Z*fine)^2

      WRITE(6,*)
      DO 13 n=1,5
        nn = float(n)
      DO 15 k=1,n
        kk = float(k)
        dEnk(n,k) = (fine**2)*(1.0d0/kk - (3.0d0/(4.0d0*nn)))/nn
 15   CONTINUE
        WRITE(6,600) n,(k,dEnk(n,k),k=1,n)
 13   CONTINUE

c     Bohr-Sommerfeld Hbetas (dk=+/-1 selection rule)
c     n=4,k=1 -> n=2,k=2
c     n=4,k=2 -> n=2,k=1
c     n=4,k=3 -> n=2,k=2

      WRITE(6,*)
      WRITE(6,*) 'Bohr-Sommerfeld Hbetas'

      E41 = E4*(1.0d0+dEnk(4,1))
      E22 = E2*(1.0d0+dEnk(2,2))
      EHb4122 = E41-E22
      wHb4122 = 1.0d8*h*c/Ehb4122

      WRITE(6,*)
      WRITE(6,*) '41 -> 22'
      WRITE(6,*) 'E_41     = ',E41
      WRITE(6,*) 'E_22     = ',E22
      WRITE(6,*) 'EHb_4122 = ',EHb4122
      WRITE(6,*) 'wH_b4122 = ',wHb4122

      E42 = E4*(1.0d0+dEnk(4,2))
      E21 = E2*(1.0d0+dEnk(2,1))
      EHb4221 = E42-E21
      wHb4221 = 1.0d8*h*c/Ehb4221

      WRITE(6,*)
      WRITE(6,*) '42 -> 21'
      WRITE(6,*) 'E_42     = ',E42
      WRITE(6,*) 'E_21     = ',E21
      WRITE(6,*) 'EHb_4221 = ',EHb4221
      WRITE(6,*) 'wH_b4221 = ',wHb4221

      E43 = E4*(1.0d0+dEnk(4,3))
      E22 = E2*(1.0d0+dEnk(2,2))
      EHb4322 = E43-E22
      wHb4322 = 1.0d8*h*c/Ehb4322

      WRITE(6,*)
      WRITE(6,*) '43 -> 22'
      WRITE(6,*) 'E_43     = ',E43
      WRITE(6,*) 'E_22     = ',E22
      WRITE(6,*) 'EHb_4322 = ',EHb4322
      WRITE(6,*) 'wH_b4322 = ',wHb4322


 600  FORMAT(1x,i3,5(i3,1x,1pe13.5))

      STOP
      END