' qm5lib.bas - jackord@kw.igs.net - revised 18 Feb 08 - Liberty Basic v4.02
' quantum harmonic oscillator
'   Pulse A from the quantum string placed between the ground state
'     turning points of the quantum oscillator
'   direct solution of the time-dependent Schrodinger Equation or
'     expansion in terms of stationary states (a 54-term series)
' differs from Java version: does not use expanded range for the
'   direct solution and hence has boundary reflections
' differs from BBC Basic version: does not show classical mass-spring oscillation

' Initialize Window
  nomainwin
  button#1, "Pulse A", [p1], UL, 5, 5, 60, 20
  button#1, "Motion", [p2], UL, 70, 5, 60, 20
  checkbox#1.c1, "Series", [waitHere], [waitHere], 142, 5, 60, 20
  textbox#1.tx, 215, 5, 80, 20
  textbox#1.ex, 300, 5, 80, 20
  WindowWidth=402     ' Pixel Scale 0-384
  WindowHeight=298    ' Pixel Scale 0-260
  UpperLeftX=10: UpperLeftY=10
  open "Quantum Oscillator" for graphics_nsb as #1
  #1 "trapclose [quit]"
  #1.tx "Frame/320"
  #1.ex "E/gs"
  #1 "place 160 60 ; down"
  #1 "\Busy"
  n=24: n2=n*n: np=384: nc= np/2: nf=320: ns=54: sf=4000: sy=500
  c=160: pi=4*atn(1): dt=n2*pi/c/nf: du=pi/c/nf/n2: w=2*pi/nf
  redim y(ns+1, np+1): redim xt(121): redim yt(121)
  redim yr(np+1): redim yi(np+1)
  redim ar(ns+1): redim ai(ns+1)
  ss=1: dx=0.1
  for k=1 to ns                                        ' Initialize y(k, i)...
    e=2*k-1: z=1+ss: y(k, nc)=z: v=1-ss
    s=z*z/2: a=0-e*z/n2: v=v+a*dx/2
    x=0: jj=0: ii=0: tst=0: tp=n*sqr(e)
    while tst=0
      jj=jj+1: x=x+dx: z=z+v*dx: a=(x*x/n2-e)*z/n2: v=v+a*dx
      if x>tp and v*z>0 then z=0: tst=1
      s=s+z*z
      if jj=10 then
        jj=0: ii=ii+1: if ii<=nc then y(k, nc+ii)=z
      end if
    wend
    s=sqr(2*s*dx)
    for i=0 to nc                                      ' ...and normalize them
      y(k, nc+i)=y(k, nc+i)/s: y(k, nc-i)=ss*y(k, nc+i)
    next i
    ss=0-ss
  next k
  xt(60)=nc: yt(60)=260
  for i=1 to 60                                        ' Potential function
    yt(60+i)=260-4*i: yt(60-i)=yt(60+i)
    xt(60+i)=nc+int(n*sqr(4*i/15)+.5): xt(60-i)=np-xt(60+i)
  next i
  for i=nc-n to nc+n                                   ' Init Pulse B
    phi=(i-(nc-n))*pi/n: yr(i)=(1-cos(phi))/(3*n)^.5
    yi(i)=0-yr(i)*sin(phi): yr(i)=yr(i)*cos(phi)
  next i
  enf=0                                                ' Init enf
  for i=nc-n to nc+n
    s1=yr(i)*(0-n2*(yr(i-1)+yr(i+1)-2*yr(i))+(i-nc)*(i-nc)/n2*yr(i))
    s2=yi(i)*(0-n2*(yi(i-1)+yi(i+1)-2*yi(i))+(i-nc)*(i-nc)/n2*yi(i))
    enf=enf+s1+s2
  next i
  ens=0
  for k=1 to ns-1 step 2                               ' Init ar(k)
    ar(k)=0
    for i=nc-n to nc+n
      ar(k)=ar(k)+yr(i)*y(k, i)
    next i
    ens=ens+ar(k)*ar(k)*(2*k-1)
  next k
  for k=2 to ns step 2                                 ' Init ai(k)
    ai(k)=0
    for i=nc-n to nc+n
      ai(k)=ai(k)+yi(i)*y(k, i)
    next i
    ens=ens+ai(k)*ai(k)*(2*k-1)                        ' Init ens
  next k
  #1 "cls ; place 160 60"
  #1 "\Ready"

[waitHere]
  wait

[p1]
  kk=1: jt=0: tt=time$("ms"): goto [plot]

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

[plot]
  if jt=0 then
    redim yr(np+1): redim yrr(np+1): redim dyr(np+1)
    redim yi(np+1): redim yii(np+1): redim dyi(np+1)
    redim yz(np+1): redim yy(np+1)
    for i=nc-n to nc+n                                 ' Init Pulse B
      phi=(i-(nc-n))*pi/n: yr(i)=(1-cos(phi))/(3*n)^.5
      yi(i)=0-yr(i)*sin(phi): yr(i)=yr(i)*cos(phi)
    next i
    #1.c1 "value?": input#1.c1, z$
    if z$="set" then energy=ens else energy=enf
      for i=nc-n to nc+n                               ' Init dy's
        dyr(i)=dt*(yi(i-1)+yi(i+1)-2*yi(i))-du*(i-nc)*(i-nc)*yi(i)
        dyi(i)=0-dt*(yr(i-1)+yr(i+1)-2*yr(i))+du*(i-nc)*(i-nc)*yr(i)
      next i
    for i=0 to np: yy(i)=260: yz(i)=yy(i): next i
    #1.ex using("##.###",energy)
    xtp=n*sqr(energy)
    yen=260-int(15*energy+.5)
  end if
  while jt=0 or plt=1                                  ' Main loop
    if plt=0 then
      for jj=1 to 2
        #1 "cls"
        if z$="set" then
          for i=0 to np                                ' Series
            y(0, i)=0
            for k=1 to ns
              if jj=1 then s1=ar(k) else s1=ai(k)
              y(0, i)=y(0, i)+s1*y(k,i)
            next k
          next i
          #1 "color green"
          for k=1 to ns                                ' Components
            #1 "place 0 130"
            for i=1 to np
              if jj=1 then s2=ar(k)*y(k, i) else s2=ai(k)*y(k, i)
              #1 "goto "; i; " "; 130-int(sy*s2+.5)
            next i
          next k
        else
          for i=0 to np
            if jj=1 then s1=yr(i) else s1=yi(i)
            y(0, i)=s1
          next i
        end if
        #1 "color blue ; place 0 130"
        for i=1 to np
          #1 "goto "; i; " "; 130-int(sy*y(0,i)+.5)
        next i
        #1 "place "; 140-10*jj; " 50"
        if jj=1 then
          #1 "\REAL COMPONENT"
        else
          #1 "\IMAGINARY COMPONENT"
        end if
        timer 2000, [cont]
        wait
        [cont]
        timer 0
      next jj
    end if
    if z$="set" then
      for i=0 to np                                    ' Series prob density
        re=0: im=0
        for k=1 to ns
          re=re+y(k, i)*(ar(k)*cos((k-1)*jt*w)-ai(k)*sin((k-1)*jt*w))
          im=im+y(k, i)*(ar(k)*sin((k-1)*jt*w)+ai(k)*cos((k-1)*jt*w))
        next k
        yy(i)=260-int(sf*(re*re+im*im)+.5)
      next i
    else
      for i=0 to np                                    ' Schrodinger prob density
        yy(i)=260                                      ' Avoid underflow error
        if abs(yr(i))>.001 then yy(i)=yy(i)-int(sf*yr(i)*yr(i)+.5)
        if abs(yi(i))>.001 then yy(i)=yy(i)-int(sf*yi(i)*yi(i)+.5)
      next i
    end if
    if jt=0 then
      #1 "cls ; down ; rule over ; color red ; place "; xt(0); " "; yt(0)
      for i=1 to 120                                   ' Plot potential and ..
        #1 "goto "; xt(i); " "; yt(i)
      next i                                           ' .. gs energy and turning points
      #1 "line 168 245 216 245 ; line 216 245 216 259 ; line 168 245 168 259 "
      #1 "color green ; line "; nc-xtp; " "; yen; " "; nc+xtp; " "; yen
      #1 "line "; nc+xtp; " "; yen; " "; nc+xtp; " 260"
      #1 "line "; nc-xtp; " "; yen; " "; nc-xtp; " 260"
      #1 "color blue ; rule xor"
    end if
    #1.tx jt; "/320"
    for i=1 to np-1                                    ' Plot prob dens
      if yy(i)<yz(i) then #1 "line "; i; " "; yy(i); " "; i; " "; yz(i)
      if yz(i)<yy(i) then #1 "line "; i; " "; yz(i); " "; i; " "; yy(i)
      yz(i)=yy(i)
    next i
    if z$="reset" then
      for jj=0 to c
        for i=1 to np-1                                ' Project ahead
          yrr(i)=yr(i)+dyr(i)/2: yii(i)=yi(i)+dyi(i)/2
        next i
        for i=1 to np-1                                ' Schrodinger equation
          dyr(i)=dt*(yii(i-1)+yii(i+1)-2*yii(i))-du*(i-nc)*(i-nc)*yii(i)
          yr(i)=yr(i)+dyr(i)
          dyi(i)=0-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i))+du*(i-nc)*(i-nc)*yrr(i)
          yi(i)=yi(i)+dyi(i)
        next i
      next jj
    end if
    jt=jt+1
    scan                                               ' Exit if bored
    #1 "discard"
    if jt=nf+1 then plt=0: kk=0: #1.tx "Time "; int((time$("ms")-tt)/1000); " s"
  wend
goto [waitHere]

[quit]
  close #1
end