import math from datetime import datetime def run_ptv_calculation(): # --- PHYSICAL CONSTANTS (CODATA June 2019) --- alpha1 = 0.0072973525693 alpha2 = alpha1 * alpha1 Rydfreq = 3289841960.2508 Me = 9.1093837015e-31 clight = 299792458 M = 0.000544617021488 pi = 3.141592653589793 h = 6.62607015e-34 # NIST experimental data: (label, NIST0, NIST1, n, L, Sp) states = { 1: ("1S1/2", 3288087922.4160, 3288086502.0102, 1, 0, 0.5), 2: ("2S1/2", 822025577.0922, 822025399.5354, 2, 0, 0.5), 3: ("3S1/2", 365343617.9043, 365343565.2949, 3, 0, 0.5), 4: ("4S1/2", 205505309.9525, 205505287.7579, 4, 0, 0.5), 5: ("5S1/2", 131523180.9882, 131523169.6246, 5, 0, 0.5), 6: ("6S1/2", 91335431.6017, 91335425.0256, 6, 0, 0.5), 7: ("2P1/2", 822026546.1359, 822026486.9664, 2, 1, -0.5), 8: ("3P1/2", 365343906.4710, 365343888.9392, 3, 1, -0.5), 9: ("4P1/2", 205505431.9299, 205505424.5337, 4, 1, -0.5), 10: ("5P1/2", 131523243.5005, 131523239.7136, 5, 1, -0.5), 11: ("6P1/2", 91335467.7972, 91335465.6058, 6, 1, -0.5), 12: ("7P1/2", 67103543.4423, 67103542.0628, 7, 1, -0.5), 13: ("2P3/2", 822015547.4967, 822015523.8451, 2, 1, 0.5), 14: ("3P3/2", 365340647.6119, 365340640.6040, 3, 1, 0.5), 15: ("4P3/2", 205504057.1004, 205504054.1439, 4, 1, 0.5), 16: ("5P3/2", 131522539.5886, 131522538.0749, 5, 1, 0.5), 17: ("6P3/2", 91335060.4413, 91335059.5653, 6, 1, 0.5), 18: ("7P3/2", 67103286.9157, 67103286.3640, 7, 2, -0.5), 19: ("3D3/2", 365340651.1931, 365340646.9865, 3, 2, -0.5), 20: ("4D3/2", 205504058.6495, 205504056.8748, 4, 2, -0.5), 21: ("5D3/2", 131522540.3924, 131522539.4837, 5, 2, -0.5), 22: ("6D3/2", 91335060.9101, 91335060.3843, 6, 2, -0.5), 23: ("7D3/2", 67103287.2124, 67103286.8813, 7, 2, -0.5), 24: ("8D3/2", 51375939.6283, 51375939.4065, 8, 2, -0.5), 25: ("3D5/2", 365339566.8037, 365339564.1003, 3, 2, 0.5), 26: ("4D5/2", 205503601.1722, 205503600.0315, 4, 2, 0.5), 27: ("5D5/2", 131522306.1639, 131522305.5800, 5, 2, 0.5), 28: ("6D5/2", 91334925.3612, 91334925.0233, 6, 2, 0.5), 29: ("7D5/2", 67103201.8522, 67103201.6394, 7, 2, 0.5), 30: ("8D5/2", 51375882.4437, 51375882.3011, 8, 2, 0.5), } # emult (gamma) values: [mm=1 (positive), mm=2 (negative)] emult_table = { # nn 1-6 nS1/2 'S': [ 1.006942700, -0.549012], # nn 7-12 nP1/2 'P12': [ 0.55773890261, -0.20481644002], # nn 13-18 nP3/2 'P32': [ 0.48878536995, -0.23361014241], # nn 19-24 nD3/2 'D32': [ 0.3957755, -0.1836735], # nn 25-30 nD5/2 'D52': [ 0.37459705, -0.1950811], } # A and B values for exitsh1 = A * dist^B, keyed by group and emult sign magshift_table = { 'S': {'+': (-28923705.16563, -3.0000421859595), '-': ( 15666540.39306, -2.999976494243)}, 'P12': {'+': (-15981616.6169, -3.00000946384), '-': ( 5870876.7450, -3.000001596923)}, 'P32': {'+': (-14008000.6385, -3.00000976361), '-': ( 6693635.451, -3.000002303498)}, 'D32': {'+': (-11341280.0376, -3.000004673972), '-': ( 5263782.149, -3.000002402484)}, 'D52': {'+': (-10734985.1348, -3.000004778242), '-': ( 5590077.3226, -3.0000023717)}, } def get_group(nn): if 1 <= nn <= 6: return 'S' if 7 <= nn <= 12: return 'P12' if 13 <= nn <= 18: return 'P32' if 19 <= nn <= 24: return 'D32' return 'D52' def get_mult(nn): if 1 <= nn <= 6: return 1 if 7 <= nn <= 12: return 1/3 if 13 <= nn <= 18: return 2/15 if 19 <= nn <= 24: return 2/25 return 25/486 def sommerfeld_fs(na, nr, alpha1, alpha2, Rydfreq, M): """[Sommerfeldfs] subroutine""" redmass1 = 1.0 / (1.0 + M) Somm04 = alpha1 / (math.sqrt(na*na - alpha2) + nr) Somm07 = 1.0 - 1.0 / math.sqrt(1.0 + Somm04**2) SommEnergy = Somm07 * redmass1 * Rydfreq * 2.0 / alpha2 return SommEnergy, redmass1 def coordinates(n, na, nr, alpha1, alpha2, M, Y, Y2): """[coordinates] subroutine — brute-force search for Xstore""" R2 = alpha1 / math.sqrt(1.0 - alpha2) Xstore = 0.0 for acc in range(-1, 15): storemark = None for msearch in range(0, 151): Xinc = msearch / (10.0**acc) X = Xstore + Xinc FX2 = n * math.sqrt(1.0 + alpha2 / Y2) * (1.0 - alpha2) Fint1 = 1.0 - M Fint2 = 1.0 + (X*X + R2*R2) / (Fint1*Fint1) Fint3 = 1.0 / (Fint1 * math.sqrt(Fint2)) Fint4 = R2*R2 * (X*X + Fint1*Fint1) / (Fint1**4 * Fint2**2) Fint5 = Fint4 * Fint4 Fint = Fint3 * (1.0 + (3.0/4.0)*Fint4 + (105.0/64.0)*Fint5) FX = (1.0 / FX2) * Fint Felectron = (0.5 / (Y2 + alpha2)) * (1.0 / math.sqrt(1.0 - alpha2 / (2.0*(Y2 + alpha2)))) Fsearch = FX - Felectron mark = 0 if Fsearch < 0 else 1 if msearch > 0 and storemark is not None and mark != storemark: Xstore = Xstore + (msearch - 1) / (10.0**acc) break storemark = mark return Xstore def mvr_fs2(na, nr, alpha1, alpha2, Rydfreq, M, evel3): """[MVRfs2] subroutine — PTV fine structure energies""" # Without evel3 (D=1 case for VortexEnergy0) evel3a = nr + math.sqrt(na*na - alpha2) evel3b = evel3a * evel3a evel3c = alpha2 / evel3b evel3d = evel3c / (2.0 * (1.0 + evel3c)) evel3e = math.sqrt(1.0 - evel3d) evel3f = evel3b + alpha2 # With evel3 (D multiplier) evel3aa = nr + math.sqrt(na*na - alpha2*evel3) evel3bb = evel3aa * evel3aa evel3cc = alpha2*evel3 / evel3bb evel3dd = evel3cc / (2.0 * (1.0 + evel3cc)) evel3ee = math.sqrt(1.0 - evel3dd) evel3ff = evel3bb + alpha2*evel3 k2 = M / evel3ee redmass3 = 1.0 / (1.0 + k2) VortexEnergy0 = Rydfreq / (evel3e * evel3f) VortexEnergy1 = redmass3 * Rydfreq / (evel3ee * evel3ff) VortexEnergy2 = redmass3 * evel3 * Rydfreq / (evel3ee * evel3ff) return VortexEnergy0, VortexEnergy1, VortexEnergy2, redmass3 def hf_model(mult, n, Rydfreq, alpha1, alpha2, dist, Y, Y2): """[Hfmodel] subroutine""" Hfsummary1 = 2.0 * n**2 * pi * Rydfreq * alpha1**3 Hfsummary2 = (dist**3) * (1.0 - alpha2)**1.5 * Y * math.sqrt(Y2 + alpha2) * math.sqrt(2.0) Hfshift = mult * Hfsummary1 / Hfsummary2 return Hfshift output_file = "PTV.txt" print(f"Writing results to {output_file}...") with open(output_file, "w") as ptv: ptv.write(f"\r\n{datetime.now().strftime('%H:%M:%S %b %d, %Y')}\r\n\r\n") ptv.write("*** A program to calculate hydrogen atom hyperfine levels using the new photonic toroidal vortex (PTV) model ***\r\n") ptv.write(" written in Python based on logic by Dr Barry R. Clarke\r\n\r\n") ptv.write(" All equations and tables refer to the following papers.\r\n\r\n") ptv.write(" PTV1: Clarke, Barry R., A photonic toroidal vortex model of the hydrogen atom fine structure, Quantum Studies:\r\n") ptv.write(" Mathematics and Foundations, 12, article 18 (2025)\r\n") ptv.write(" PTV2: Clarke, Barry R., The hydrogen atom hyperfine structure from a photonic toroidal vortex model,\r\n") ptv.write(" in peer-review (2026)\r\n\r\n") ptv.write(" EXPERIMENTAL DATA\r\n") ptv.write(" Horbatsch, M. and E. A. Hessels, Tabulation of the bound-state energies of atomic hydrogen, Physical Review A 93, 022513 (2016)\r\n\r\n") ptv.write("_____________________________________________________________________________________________________\r\n\r\n") # Outer loop mm = 1 to 2 (two gamma values per state group), inner loop nn = 1 to 30 for mm in range(1, 3): for nn in range(1, 31): label, NIST0, NIST1, n, L, Sp = states[nn] group = get_group(nn) mult = get_mult(nn) # Sommerfeld quantum numbers na = L + Sp + 0.5 nr = n - na Y = nr + math.sqrt(na*na - alpha2) Y2 = Y * Y # Sommerfeld fine structure SommEnergy, redmass1 = sommerfeld_fs(na, nr, alpha1, alpha2, Rydfreq, M) # Coordinate search (independent of emult) Xstore = coordinates(n, na, nr, alpha1, alpha2, M, Y, Y2) r2f2 = Xstore**2 + (1.0 - M)**2 dist = math.sqrt(r2f2) # Select emult (gamma) emult = emult_table[group][mm - 1] # evel3 = D multiplier: (1 + emult*M/n)^2 evel3 = (1.0 + emult * M / n)**2 # PTV fine structure energies VortexEnergy0, VortexEnergy1, VortexEnergy2, redmass3 = \ mvr_fs2(na, nr, alpha1, alpha2, Rydfreq, M, evel3) # Hyperfine shift Hfshift = hf_model(mult, n, Rydfreq, alpha1, alpha2, dist, Y, Y2) # Magnetic potential adjustment sign_key = '+' if emult > 0 else '-' A, B = magshift_table[group][sign_key] exitsh1 = A * (dist**B) # --- OUTPUT (mirrors BASIC [printout1]) --- # Liberty BASIC formatting helpers: # template1$ = '###################.######' => {:27.6f} with space-sign (27 total) # template2$ = '#############.##########' => {:25.10f} with space-sign (25 total) # template3$ = '######.################' => {:24.16f} with space-sign (24 total) # evel3 printed with BASIC default: 8dp, trailing zeros stripped def t1(v): return f'{v: 27.6f}' def t2(v): return f'{v: 25.10f}' def t3(v): return f'{v: 24.16f}' def evel3_fmt(v): return f'{v:.8f}'.rstrip('0').rstrip('.') # Quantum number formatting: int if whole number, else as-is def qn(v): return int(v) if v == int(v) else v ptv.write(f"{label}\r\n") ptv.write("QUANTUM NUMBERS\r\n") ptv.write(f"(1) Quantum Mechanics: n = {qn(n)} L = {qn(L)} Sp = {Sp} Sommerfeld: n(phi) = {qn(na)} n(r) = {qn(nr)}\r\n\r\n") ptv.write("BOUND STATE DISTANCES (as a muliplier of a fixed radius for all states)\r\n") ptv.write(f"(2) Horizontal distance between p-e Sp-2 centers ={t3(Xstore)} PTV1 x-bar Eqs (84)\r\n") ptv.write(f"(3) Shortest distance between p-e Sp-2 centers ={t3(dist)} PTV1 d2-bar Eq (86)\r\n\r\n") ptv.write(f" gamma ={t3(emult)}\r\n\r\n") ptv.write("FINE STRUCTURE (ALL ENERGIES IN MHz)\r\n") ptv.write(f"(4) Sommerfeld f-s energy no adjustments ={t1(SommEnergy/redmass1)}\r\n") ptv.write(f"(5) New PTV f-s energy no adjustments ={t1(VortexEnergy0)} PTV1 Eq (54)\r\n") ptv.write(f"(6) Sommerfeld f-s energy reduced mass ={t1(SommEnergy)}\r\n") ptv.write(f"(7) Relativistic reduced mass multiplier = {t3(redmass3)} PTV2 (1 + Mr)^(-1) Eq (13)\r\n") ptv.write(f"(8) New PTV f-s energy relativ. reduced mass ={t1(VortexEnergy1)} PTV2 (D = 1) Eq (13)\r\n") ptv.write(f"(9) Velocity D multiplier = {evel3_fmt(evel3)} PTV2 Eq (13)\r\n") ptv.write("(10) New PTV f-s energy with relativistic reduced mass \r\n") ptv.write(f" and adjusted velocity ={t1(VortexEnergy2)} PTV2 Eqs (12)/(13)\r\n\r\n") ptv.write("MAGNETIC POTENTIAL ADJUSTMENT (ALL ENERGIES IN MHz)\r\n") ptv.write(f"(11) Average h-f energy ={t1((NIST0 + NIST1)/2)} Horbatsch and Hessels (2016)\r\n") diff_mag = abs(VortexEnergy2 - (NIST0 + NIST1)/2) ptv.write("(12) Difference between PTV f-s with corrections\r\n") ptv.write(f" and average h-f energy (magnitude) ={t1(diff_mag)}\r\n") ptv.write(f"(13) Magnetic potential adjustment ={t1(exitsh1)} PTV2 Eq (22) & Tables 2 and 3\r\n") if emult > 0: # F=1 branch (lower frequency state) partscalc0 = NIST1 / (NIST1 - (VortexEnergy2 - Hfshift + exitsh1)) aa01 = abs(VortexEnergy2 - (NIST0 + NIST1)/2 - exitsh1) ptv.write(f"(14) H-f mid-point error magnitude ={t1(aa01)} PTV2 Eq (24) & Table 2\r\n") ptv.write(f"(15) H-f mid-point PTV ={t1(abs(VortexEnergy2 + exitsh1))}\r\n\r\n") ptv.write("HYPERFINE SHIFT (ALL ENERGIES IN MHz)\r\n") ptv.write(f"(16) Lower h-f energy from experiment ={t1(NIST1)} Horbatsch and Hessels (2016)\r\n") ptv.write(f"(17) Lower h-f energy calculated from PTV ={t1(VortexEnergy2 + exitsh1 - Hfshift)} PTV2 Eqs (12)/(13) plus Eq. (24) minus Eq. (35)\r\n") ptv.write(f"(18) Deviation of PTV from experiment is one part in{t1(abs(partscalc0))}\r\n") ptv.write(f"(19) H-f difference taken from experiment ={t1((NIST0 - NIST1)/2)} Horbatsch and Hessels (2016)\r\n") ptv.write(f"(20) H-f difference calculated from PTV ={t1(Hfshift)} PTV2 Eq. (35)\r\n") ptv.write(f"(21) H-f error magnitude ={t1(abs(Hfshift - (NIST0 - NIST1)/2))}\r\n") ptv.write(f"(22) A/gamma = {t2(A/emult)}\r\n") else: # F=0 branch (higher frequency state) partscalc1 = NIST0 / (NIST0 - (VortexEnergy2 + Hfshift + exitsh1)) bb01 = abs(VortexEnergy2 - (NIST0 + NIST1)/2 + exitsh1) ptv.write(f"(14) H-f mid-point error magnitude ={t1(bb01)}\r\n") ptv.write(f"(15) H-f mid-point PTV ={t1(abs(VortexEnergy2 + exitsh1))}\r\n\r\n") ptv.write("HYPERFINE SHIFT (ALL ENERGIES IN MHz)\r\n") ptv.write(f"(16) Higher h-f energy from experiment ={t1(NIST0)} Horbatsch and Hessels (2016)\r\n") ptv.write(f"(17) Higher h-f energy calculated from PTV ={t1(VortexEnergy2 + exitsh1 + Hfshift)} PTV2 Eqs (12)/(13) plus Eq. (22) minus Eq. (35)\r\n") ptv.write(f"(18) Deviation of PTV from experiment is one part in{t1(abs(partscalc1))}\r\n") ptv.write(f"(19) H-f difference taken from experiment ={t1((NIST0 - NIST1)/2)} Horbatsch and Hessels (2016)\r\n") ptv.write(f"(20) H-f difference calculated from PTV ={t1(Hfshift)} PTV2 Eq. (35)\r\n") ptv.write(f"(21) H-f error magnitude ={t1(abs(Hfshift - (NIST0 - NIST1)/2))}\r\n") ptv.write(f"(22) A/gamma = {t2(A/emult)}\r\n") ptv.write("\r\n_____________________________________________________________________________________________________\r\n\r\n") ptv.write(" *** Free license, Barry R. Clarke (aleteller@barryispuzzled.com), 26 April 2026 ***\r\n") print("Done.") if __name__ == "__main__": run_ptv_calculation()