' qm1lib.bas - jackord@kw.igs.net - 18 Feb 08 - Liberty Basic v4.02
' quantum rectilinear pluck and pulse wavefunctions
'   and a standing or travelling localized sinusoidal wavefunction
'   in a square well with time dependence calculated either directly
'   from the Schrodinger Equation or from FFSS expansion
' the FFSS expansion runs much faster than the direct solution

' Initialize Window
  nomainwin
  button#1, "Pluck", [p1], UL, 4, 5, 60, 20
  button#1, "Pulse", [p2], UL, 68, 5, 60, 20
  button#1, "PulseS", [p3], UL, 132, 5, 60, 20
  button#1, "PulseA", [p4], UL, 196, 5, 60, 20
  button#1, "Motion", [p5], UL, 260, 5, 60, 20
  checkbox#1.c1, "FFSS", [waitHere], [waitHere], 324, 5, 54, 20
  textbox#1.tx, 5, 30, 80, 20
  textbox#1.ex, 5, 55, 80, 20
  WindowWidth=402     ' Pixel Scale 0-384
  WindowHeight=298    ' Pixel Scale 0-260
  UpperLeftX=10: UpperLeftY=10
  open "Quantum String" for graphics_nsb as #1
  #1 "trapclose [quit] ; color blue; rule xor ; down"
  #1.tx "Frame/320"
  #1.ex "E/gs"
  n=192: nm=64: nf=320: c=600: pi=4*atn(1): dt=n*n/4/pi/c/nf

[waitHere]
  wait

[p1]
  kk=1: jt=0: goto [plot]

[p2]
  kk=2: jt=0: goto [plot]

[p3]
  kk=3: jt=0: goto [plot]

[p4]
  kk=4: jt=0: goto [plot]

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

[plot]
  if jt=0 then
    redim yr(n+1): redim yrr(n+1): redim dyr(n+1)
    redim yi(n+1): redim yii(n+1): redim dyi(n+1)
    redim z(n+1): redim yy(n+1): redim b(nm+1): redim d(nm+1): redim w(nm+1)
    #1 "cls ; color blue"
    if kk=1 then                                       ' Init Pluck
      for i=0 to n/2
        yr(i)=30*i/n: yr(n-i)=yr(i)
      next i
      energy=12/pi/pi
    end if
    if kk=2 then                                       ' Init Pulse
      for i=3*n/8 to n/2
        yr(i)=4*30/n*(i-3*n/8): yr(n-i)=yr(i)
      next i
      energy=16*12/pi/pi
    end if
    if kk=3 or kk=4 then                               ' Init Pulse A
      for i=3*n/8 to 5*n/8
        phi=(i-3*n/8)*pi*8/n: yr(i)=7.5*(1-cos(phi))
        if kk=4 then                                   ' Init Pulse B
          yi(i)=0-yr(i)*sin(phi): yr(i)=yr(i)*cos(phi)
        end if
      next i
      energy=21.333: if kk=4 then energy=64+energy
    end if
    #1.c1 "value?": input#1.c1, z$
    if z$="set" then                                   ' Init FFSS
      s=0: se=0
      for i=1 to nm-1
        for j=1 to n-1
          sn=sin(pi*i*j/n)
          b(i)=b(i)+yr(j)*2/n*sn: d(i)=d(i)+yi(j)*2/n*sn
        next j
        w(i)=i*i*2*pi/8/nf
        s=s+b(i)*b(i)+d(i)*d(i): se=se+(b(i)*b(i)+d(i)*d(i))*i*i
      next i
      energy=se/s
    else
      for i=1 to n-1                                   ' Init look ahead
        dyr(i)=dt*(yi(i-1)+yi(i+1)-2*yi(i))
        dyi(i)=0-dt*(yr(i-1)+yr(i+1)-2*yr(i))
      next i
    end if
    for i=0 to n: yy(i)=260: z(i)=yy(i): next i
    #1.ex using("##.###",energy)
  end if
  if z$="reset" then #1.c1 "reset"
  tt=time$("ms")
  while jt=0 or plt=1                                  ' Main loop
    #1.tx jt; "/320"
    for i=1 to n-1                                     ' plot prob dens
      yy(i)=260-int(yr(i)*yr(i)+yi(i)*yi(i)+.5)
      if yy(i)<z(i) then #1 "line "; 2*i; " "; yy(i); " "; 2*i; " "; z(i)
      if z(i)<yy(i) then #1 "line "; 2*i; " "; z(i); " "; 2*i; " "; yy(i)
      z(i)=yy(i)
    next i
    jt=jt+1
    if kk<4 then                                       ' Pluck/Pulse/Pulse A
      if z$="set" then                                 ' FFSS
        for i=1 to n-1
          yr(i)=0: yi(i)=0
          for k=1 to nm-1 step 2
            sn=sin(pi*i*k/n)
            yr(i)=yr(i)+sn*b(k)*cos(w(k)*jt)
            yi(i)=yi(i)+sn*b(k)*sin(w(k)*jt)
          next k
        next i
      else
        for jj=1 to c                                  ' Accuracy loop
          for i=1 to n/2+1                             ' Look ahead dy/2
            yrr(i)=yr(i)+dyr(i)/2: yii(i)=yi(i)+dyi(i)/2
          next i
          for i=1 to n/2                               ' Schrodinger Equation
            dyr(i)=dt*(yii(i-1)+yii(i+1)-2*yii(i)): dyr(n-i)=dyr(i)
            yr(i)=yr(i)+dyr(i): yr(n-i)=yr(i)
            dyi(i)=0-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i)): dyi(n-i)=dyi(i)
            yi(i)=yi(i)+dyi(i): yi(n-i)=yi(i)
          next i
        next jj
      end if
      for i=1 to n/2
        yy(i)=260-int(yr(i)*yr(i)+yi(i)*yi(i)+.5): yy(n-i)=yy(i)
      next i
    else                                               ' Pulse B
      if z$="set" then                                 ' FFSS
        for i=1 to n-1
          yr(i)=0: yi(i)=0
          for k=1 to nm-1 step 2
            sn=sin(pi*i*k/n)
            yr(i)=yr(i)+sn*b(k)*cos(w(k)*jt)
            yi(i)=yi(i)+sn*b(k)*sin(w(k)*jt)
          next k
          for k=2 to nm-2 step 2
            sn=sin(pi*i*k/n)
            yr(i)=yr(i)-sn*d(k)*sin(w(k)*jt)
            yi(i)=yi(i)+sn*d(k)*cos(w(k)*jt)
          next k
        next i
      else
        for jj=1 to c                                  ' Accuracy loop
          for i=1 to n-1                               ' Project ahead dy/2
            yrr(i)=yr(i)+dyr(i)/2: yii(i)=yi(i)+dyi(i)/2
          next i
          for i=1 to n-1                               ' Schrodinger Equation
            dyr(i)=dt*(yii(i-1)+yii(i+1)-2*yii(i))
            yr(i)=yr(i)+dyr(i)
            dyi(i)=0-dt*(yrr(i-1)+yrr(i+1)-2*yrr(i))
            yi(i)=yi(i)+dyi(i)
          next i
        next jj
      end if
    end if
    scan                                               ' Exit if bored
    if jt=nf+1 then
      plt=0: kk=0
      #1 "color black ; place 280 40"
      #1 "\Time "; (time$("ms")-tt)/1000; " s"
    end if
    #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end