' qm2lib.bas - jackord@kw.igs.net - 18 Feb 08 - Liberty Basic v4.02
' wave pulse on a quantum string
' time dependence calculated either directly from the Schrodinger Equation
'   or from FFSS expansion

' Initialize Window
  nomainwin
  button#1, "P Dens", [p1], UL, 5, 5, 70, 20
  button#1, "R Comp", [p2], UL, 85, 5, 70, 20
  button#1, "Motion", [p3], UL, 165, 5, 70, 20
  checkbox#1.c1, "FFSS", [waitHere], [waitHere], 245, 5, 70, 20
  textbox#1.tx, 5, 30, 80, 20
  WindowWidth=338     ' Pixel Scale 0-320
  WindowHeight=258    ' Pixel Scale 0-220
  UpperLeftX=10: UpperLeftY=10
  open "Quantum Pulse" for graphics_nsb as #1
  #1.tx "Frame/320"
  #1 "trapclose [quit] ; color blue ; rule xor"
  n=160: nm=60: nf=320: nk=4: c=8: jt=0: pi=4*atn(1): dt=n*n/4/pi/c/nf/nk/nk

[waitHere]
  wait

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

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

[p3]
  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 ; down"
    for i=3*n/8 to 5*n/8
      phi=(i-3*n/8)*pi*8/n: yr(i)=(1-cos(phi))/2
      yi(i)=0-yr(i)*sin(phi*nk): yr(i)=yr(i)*cos(phi*nk)
    next i
    #1.c1 "value?": input#1.c1, z$
    if z$="set" then
      for i=1 to nm-1                                   ' Init FFSS
        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/16
      next i
    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)=220/kk: z(i)=yy(i): next i
  end if
  if z$="reset" then #1.c1 "reset"
  while jt=0 or plt=1                                  ' Main loop
    #1.tx jt; "/320"
    for i=1 to n
      if kk=1 then                                     'Prob density
        yy(i)=220-int(160*(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)
      else                                             'Real Comp
        yy(i)=110-int(80*yr(i)+.5)
        if jt>0 then #1 "line "; 2*(i-1); " "; z(i-1); " "; 2*i; " "; z(i)
        #1 "line "; 2*(i-1); " "; yy(i-1); " "; 2*i; " "; yy(i)
      end if
      z(i-1)=yy(i-1)
    next i
    jt=jt+1
    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
    scan                                               ' Exit if bored
    if jt=nf+1 then plt=0: kk=0
    #1 "discard"
  wend
goto [waitHere]

[quit]
  close #1
end