' qm4lib.bas - jackord@kw.igs.net - 18 Feb 08 - Liberty Basic v4.02
' superposition of the two lowest energy states of a double well (an
'   infinite square well with a central barrier)

' Initialize Window
  nomainwin
  button#1, "P Dens", [b1], UL, 5, 5, 60, 20
  button#1, "Yr Yi", [b2], UL, 70, 5, 55, 20
  button#1, "Motion", [b3], UL, 130, 5, 60, 20
  WindowWidth=418     ' pixel scale 0-400
  WindowHeight=338    ' pixel scale 0-300
  UpperLeftX=10: UpperLeftY=100
  open "Double Well Superposition" for graphics_nsb as #1   ' no slide bars
  #1 "trapclose [quit]"
  n=400: pi=4*atn(1): c=pi/400: c=c*c
  ea=3.85297: eb=4.07503
  gosub [grid]

[waitHere]
  wait

[b1]
  kk=1: jt=0: plt=0: goto [plot]

[b2]
  kk=2: jt=0: plt=0: goto [plot]

[b3]
  if kk>0 then
    plt=1: goto [plot]
  else
    goto [waitHere]
  end if

[plot]
  if jt=0 then
    gosub [grid]
    redim ya(n+1): redim yb(n+1)
    redim z(n+1): redim zi(n+1): redim yy(n+1): redim yz(n+1)
    e=ea: v=0.9: sa=0
    for i=1 to n-1
      ya(i)=ya(i-1)+v: sa=sa+ya(i)*ya(i)
      if i>191 and i<209 then a=(200-e)*c else a=0-e*c
      v=v+a*ya(i)
    next i
    e=eb: v=0.9: sb=0
    for i=1 to n-1
      yb(i)=yb(i-1)+v: sb=sb+yb(i)*yb(i)
      if i>191 and i<209 then a=(200-e)*c else a=0-e*c
      v=v+a*yb(i)
    next i
    sb=sqr(sa/sb)
    for i=1 to n-1: yb(i)=sb*yb(i): next i
    for i=0 to n: yy(i)=300/kk: yz(i)=yy(i): z(i)=yy(i): zi(i)=yy(i): next i
    if kk=1 then nf=320 else nf=960
    a=2*pi/(eb-ea)/nf
  end if
  while jt=0 or plt=1
    #1 "place 290 20 ; color black"
    #1 "\Frame "; jt; "/"; nf
    #1 "color blue"
    for i=1 to n
      yr=ya(i)*cos(jt*ea*a)+yb(i)*cos(jt*eb*a)
      yi=ya(i)*sin(jt*ea*a)+yb(i)*sin(jt*eb*a)
      if kk=1 then                                     ' Prob density
        yy(i)=300-int(.019*(yr*yr+yi*yi)+.5)
        if yy(i)<z(i) then #1 "line "; i; " "; yy(i); " "; i; " "; z(i)
        if z(i)<yy(i) then #1 "line "; i; " "; z(i); " "; i; " "; yy(i)
      else                                             ' Real Comp
        yy(i)=150-int(yr+.5)
        if jt>0 then #1 "line "; i-1; " "; z(i-1); " "; i; " "; z(i)
        #1 "line "; i-1; " "; yy(i-1); " "; i; " "; yy(i)
        #1 "color red"
        yz(i)=150-int(yi+.5)
        if jt>0 then #1 "line "; i-1; " "; zi(i-1); " "; i; " "; zi(i)
        #1 "line "; i-1; " "; yz(i-1); " "; i; " "; yz(i)
        zi(i-1)=yz(i-1)
        #1 "color blue"
      end if
      z(i-1)=yy(i-1)
      #1 "discard"
    next i
    scan
    jt=jt+1: if jt=nf+1 then plt=0: kk=0
  wend
goto [waitHere]

[grid]
  #1 "cls ; rule xor ; down"
  #1 "color green ; backcolor green ; place 192 0 ; boxfilled 208 300"
  #1 "color black ; backcolor white"
return

[quit]
  close #1
end