' h7alib.bas - jackord@kw.igs.net - revised 10 Sep 02 - Liberty Basic v3.01
'   Schrodinger equation for the hydrogen atom:
'     bisection algorithm for finding eigenvalues
'     Feynman algorithm for the radial wavefunction
'     optional radial probability density plot

' Initialize Window
  nomainwin
  radiobutton#1.r1, "Eigen", [eigen], [waitHere], 425, 25, 65, 20
  radiobutton#1.r2, "Pdens", [pdens], [waitHere], 490, 25, 65, 20
  button#1, "1s", [b1s], UL, 425, 117, 30, 20
  button#1, "2s", [b2s], UL, 425, 94, 30, 20
  button#1, "2p", [b2p], UL, 458, 94, 30, 20
  button#1, "3s", [b3s], UL, 425, 71, 30, 20
  button#1, "3p", [b3p], UL, 458, 71, 30, 20
  button#1, "3d", [b3d], UL, 491, 71, 30, 20
  button#1, "4s", [b4s], UL, 425, 48, 30, 20
  button#1, "4p", [b4p], UL, 458, 48, 30, 20
  button#1, "4d", [b4d], UL, 491, 48, 30, 20
  button#1, "4f", [b4f], UL, 524, 48, 30, 20
  WindowWidth=570     ' pixel scale 0-560
  WindowHeight=309    ' pixel scale 0-280
  UpperLeftX=10: UpperLeftY=100
  open "Hatom Wavefunctions and Eigenvalues" for graphics_nsb as #1
  #1 "trapclose [quit]"
  #1.r1 "Set": opt=1
  gosub [grid]

[waitHere]
  wait

[eigen]
  opt=1: goto [waitHere]

[pdens]
  opt=2: goto [waitHere]

[b1s]
  n=1: L=0: ymax=1.1: zmax=.368: goto [plot]

[b2s]
  n=2: L=0: ymax=1.1: zmax=.618: goto [plot]

[b2p]
  n=2: L=1: ymax=1: zmax=2.165: goto [plot]

[b3s]
  n=3: L=0: ymax=1.1: zmax=.828: goto [plot]

[b3p]
  n=3: L=1: ymax=1: zmax=2.64: goto [plot]

[b3d]
  n=3: L=2: ymax=5: zmax=36.3: goto [plot]

[b4s]
  n=4: L=0: ymax=1.1: zmax=1.1: goto [plot]

[b4p]
  n=4: L=1: ymax=1: zmax=3.15: goto [plot]

[b4d]
  n=4: L=2: ymax=5: zmax=36.5: goto [plot]

[b4f]
  n=4: L=3: ymax=100: zmax=1200: goto [plot]

[plot]                                                 ' Main Program
  gosub [grid]
  #1 "color blue"
  en=0-1/n/n: de=.00001: zmax=1.05*zmax
  if opt=1 then                                        ' Find Eigenvalue
    a=en*1.1: e=a: gosub [fey]: sa=s
    b=en/1.1: e=b: gosub [fey]: sb=s
    while (b-a)>de                                     ' Bisection Algorithm
      e=(a+b)/2: gosub [fey]: se=s
      if sa=se then a=e else b=e
    wend
    #1 "color red"
    e=(a+b)/2: gosub [fey]
  else                                                 ' Plot Prob Density
    e=en: gosub [fey]
  end if
  #1 "color white ; place 425 2 ; boxfilled 554 22 ; color black"
  #1 "\\E = "; using("##.######", e)
goto [waitHere]

[fey]                                                  ' Feynman Algorithm
  r=0: z=0: u=1: v=0-u/(L+1): dr=n*n*.002
  dv=0-e*u*dr: v=v+dv/2
  if opt=1 then                                        ' Plot R
    if L=0 then y1=int(140*(1-1/ymax)) else y1=140
    #1 "place 0 "; y1
    while z>0-ymax and z<ymax
      r=r+dr: u=u+v*dr: z=u*r^L
      #1 "goto "; int(140*r/n/n); " "; int(140*(1-z/ymax))
      dv=0-(e*u+2/r*(u+(L+1)*(v+dv/2)))*dr: v=v+dv
    wend
  else                                                 ' Plot (rR)^2
    #1 "place 0 140"
    while r<4*n*n
      r=r+dr: u=u+v*dr: z=u*r^(L+1)
      #1 "goto "; int(140*r/n/n); " "; int(140*(1-z*z/zmax/zmax))
      dv=0-(e*u+2/r*(u+(L+1)*(v+dv/2)))*dr: v=v+dv
    wend
  end if
  if u>0 then s=1 else s=0
return

[grid]
  #1 "cls": #1 "down": #1 "color black"
  for i=1 to 3
    x=140*i: #1 "line "; x; " 0 "; x; " 280"
  next i
  y=140: #1 "line 0 "; y; " 560 "; y
  #1 "place 425 2": #1 "\\Energy"
return

[quit]
  close #1
end