' fosdlib.bas - jackord@kw.igs.net - revised 6 Sep 02 - Liberty Basic v3.01
' the driven damped oscillator (initially at rest at origin)
'   - sine or square wave driving force
'   - f/fo equals 1/6, 1/2, 1, 2
'   - Q equals 1/3, 1, 3, inf
'   - sine wave analytic solution superposition option (except f=fo, Q=inf)

' Initialize Window
  nomainwin
  button#1, "Sine", [sin], UL, 5, 5, 55, 20
  button#1, "Square", [sqr], UL, 65, 5, 55, 20
  f$(0)="f=1/6": f$(1)="f=1/2": f$(2)="f=1": f$(3)="f=2"
  q$(0)="Q=1/3": q$(1)="Q=1": q$(2)="Q=3": q$(3)="Q=inf"
  combobox#1.fs, f$(, [fset], 5, 30, 70, 130
  combobox#1.qs, q$(, [qset], 80, 30, 70, 130
  checkbox#1.af, "Anal", [aset], [aclr], 125, 5, 55, 20
  WindowWidth=490     ' pixel scale 0-480
  WindowHeight=409    ' pixel scale 0-380
  UpperLeftX=50: UpperLeftY=50
  open "Driven Damped Oscillator" for graphics_nsb as #1
  #1 "trapclose [quit]"
  wo=atn(1)/10: w2=wo*wo: c=45: d=w2*c: dt=1: n=480
  fs(0)=240: fs(1)=80: fs(2)=40: fs(3)=20
  bs(0)=3: bs(1)=1: bs(2)=1/3: bs(3)=0
  i=2: print#1.fs, "selectindex "; i: nf=fs(i-1)
  #1.qs "selectindex "; i: rq=bs(i-1): b=wo*rq
  pi=4*atn(1)

[waitHere]
  wait

[fset]                                                 ' set frequency
  #1.fs "selectionindex?"
  input#1.fs, index
  nf=fs(index-1)
goto [waitHere]

[qset]                                                 ' set Q
  #1.qs "selectionindex?"
  input#1.qs, index
  rq=bs(index-1): b=wo*rq
goto [waitHere]

[aset]                                                 ' superimpose analytic
  flg=1
goto [waitHere]

[aclr]                                                 ' no analytic
  flg=0
goto [waitHere]

[sin]                                                  ' Sine wave
  #1 "cls ; down ; color black ; line 0 200 480 200"
  w=40*wo/nf
  #1 "color red ; place 0 200"
  for i=1 to n                                         ' Driving force
    x=200-c*sin(i*w*dt)
    #1 "goto "; i; " "; x
  next i
  x=0: v=0: dv=0                                       ' Feynman algorithm
  #1 "color blue ; place 0 200"                        ' no Half-Step (a=0)
  for i=1 to n
    x=x+v*dt: dv=(d*sin(w*i*dt)-b*(v+dv/2)-w2*x)*dt: v=v+dv
    #1 "goto "; i; " "; 200-x
  next i
  if flg=1 and (b>0 or nf>40 or nf<40) then            ' Analytic solution
    #1 "color darkpink ; place 0 200"
    r=w/wo: s=1/r
    if nf=40 then
      a=pi/2
    else
      a=atn(rq/(s-r)): if s<r then a=a+pi
    end if
    cc=c*s/((s-r)*(s-r)+rq*rq)^.5
   if rq>2 then                                        ' ...overdamped...
      wa=0-wo*(rq/2+(rq*rq/4-1)^.5)
      wb=0-wo*(rq/2-(rq*rq/4-1)^.5)
      a2=cc*(wa*sin(a)+w*cos(a))/(wa-wb)
      a1=cc*sin(a)-a2
      for i=1 to n
        y=200-(a1*exp(wa*i*dt)+a2*exp(wb*i*dt)+cc*sin(w*i*dt-a))
        #1 "goto "; i; " "; y
      next i
    else                                               ' ...underdamped
      wa=wo*(1-rq*rq/4)^.5
      a1=cc*sin(a): a2=(a1*b/2-w*cc*cos(a))/wa
      for i=1 to n
        y=200-(exp(0-b*i*dt/2)*(a1*cos(wa*i*dt)+a2*sin(wa*i*dt))+cc*sin(w*i*dt-a))
        #1 "goto "; i; " "; y
      next i
    end if
  end if

goto [waitHere]

[sqr]                                                  ' Square wave
  #1 "cls ; down ; color black ; line 0 200 480 200"
  #1 "color red ; line 0 200 0 155 ; line 480 245 480 200" ' Driving force
  for i=1 to n/nf-1
    #1 "line "; i*nf; " 155 "; i*nf; " 245"
  next i
  for i=0 to n/2/nf-1
    x1=i*2*nf: #1 "line "; x1; " 155 "; x1+nf; " 155"
    #1 "line "; x1+nf; " 245 "; x1+2*nf; " 245"
  next i
  #1 "color blue ; place 0 200"
  x=0: v=0: dv=d*dt: v=v+dv/2                          ' Feynman algorithm
  for i=0 to n/nf-1
    for j=1 to nf
      x=x+v*dt
      dv=(d-b*(v+dv/2)-w2*x)*dt: v=v+dv
      #1 "goto "; j+i*nf; " "; 200-x
    next j
    v=v-d*dt: dv=dv-d*dt                               ' Discontinuity corrn
    d=0-d
  next i
goto [waitHere]

[quit]
  close #1
end