This is the output of Hbeta.f (with some additional comments)

C the reduced mass and Rydberg constants (computed from fundamental constants)

 muH/me   =   0.999455678993799     
 Ryd_inf  =   2.179874100000000E-011
 RydH     =   2.178687548736497E-011

C Bohr model Hbeta energies and wavelength for unreduced electron mass

 Unreduced Bohr Hbeta
 E_4      =  -1.362421312500000E-012
 E_2      =  -5.449685250000000E-012
 EHb      =   4.087263937500000E-012
 wHb      =    4860.09098363636     

C Bohr model Hbeta level energies, transition energy, and wavelength for 
C reduced electron mass: note the 2 angstrom longer wavelength; the energies 
C show that that the electron is less bound when accounting for non-infinite 
c mass nucleus

 Reduced (best) Bohr Hbeta
 E_4      =  -1.361679717960311E-012
 E_2      =  -5.446718871841243E-012
 EHb      =   4.085039153880932E-012
 wHb      =    4862.73787400883     

C Bohr-Sommerfeld fractional energy shifts and splittings, accounting for reduced 
C mass of electron

   n  k   dE_nk/E_n    k   dE_nk/E_n    k   dE_nk/E_n    k   dE_nk/E_n    k   dE_nk/E_n
   1  1   1.33128E-05
   2  1   1.66411E-05  2   3.32821E-06
   3  1   1.33128E-05  2   4.43761E-06  3   1.47920E-06
   4  1   1.08167E-05  2   4.16026E-06  3   1.94146E-06  4   8.32053E-07
   5  1   9.05273E-06  2   3.72760E-06  3   1.95255E-06  4   1.06503E-06  5   5.32514E-07

C The Bohr-Sommerfeld Hbeta level energies, transition energies, and wavelengths
C using the empirical Sommerfeld selection rule dk=+/-1

 Bohr-Sommerfeld Hbetas

 41 -> 22
 E_41     =  -1.361694446817958E-012
 E_22     =  -5.446736999666040E-012
 EHb_4122 =   4.085042552848081E-012
 wH_b4122 =    4862.73382795886     

 42 -> 21
 E_42     =  -1.361685382905560E-012
 E_21     =  -5.446809510965229E-012
 EHb_4221 =   4.085124128059669E-012
 wH_b4221 =    4862.63672477950     

 43 -> 22
 E_43     =  -1.361682361601427E-012
 E_22     =  -5.446736999666040E-012
 EHb_4322 =   4.085054638064613E-012
 wH_b4322 =    4862.71944205795     

C the fractional % differences from the Bohr model (not computed by code)

wH_b4122 - WHb = -0.004047 (smallest discrepency)
wH_b4221 - Whb = -0.101150 (largest discrepency)
wH_b4322 - Whb = -0.018432

C nk=41 is most eccentric orbit for n=4; nk=22 is a circular orbit for n=2 
C nk=42 is second most eccentircity for n=4, nk=21 is most eccentric orbit for n=2