REM fosdb.bas - jackord@kw.igs.net - revised 27 Feb 2011 - BBC Basic 5.92a
      REM the driven damped oscillator (initially at rest at the origin)
      REM   - sine or square wave driving force
      REM   - f/fo equals 1/6, 1/2, 1, 2
      REM   - Q equals 1/3, 1, 3, inf
      REM   - for the sinusoidal force the analytic solution is superimposed after 2 s
      REM       (except for f=fo and Q=inf)
      
      *FLOAT 64
      INSTALL @lib$+ "WINLIB5"
      
      b1% = FN_button("Sine", 5, 5, 55, 25, FN_setproc(PROCsin), 0)
      b2% = FN_button("Square", 65, 5, 55, 25, FN_setproc(PROCsqr), 0)
      c1% = FN_combobox("", 125, 5, 65, 0, 0, 3)
      c2% = FN_combobox("", 195, 5, 65, 0, 0, 3)
      WindowWidth=480: WindowHeight=380
      VDU 23, 22, WindowWidth; WindowHeight; 8, 16, 16, 128
      SYS "SetWindowText", @hwnd%, "The Driven Damped Oscillator"
      SYS "SendMessage", c1%, &143, 0, "f=1/6"
      SYS "SendMessage", c1%, &143, 1, "f=1/2"
      SYS "SendMessage", c1%, &143, 1, "f=1"
      SYS "SendMessage", c1%, &143, 1, "f=2"
      SYS "SendMessage", c1%, &14E, 2, 0
      SYS "SendMessage", c2%, &143, 0, "Q=1/3"
      SYS "SendMessage", c2%, &143, 1, "Q=1"
      SYS "SendMessage", c2%, &143, 1, "Q=3"
      SYS "SendMessage", c2%, &143, 1, "Q=inf"
      SYS "SendMessage", c2%, &14E, 1, 0
      DIM fs(4), bs(4)
      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
      wo=ATN(1)/10: w2=wo*wo: c=90: d=w2*c: dt=1: n=480
      VDU 5: OFF
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCsin
      SYS "SendMessage", c1%, &147, 0, 0 TO fi: nf=fs(fi): w=40*wo/nf
      SYS "SendMessage", c2%, &147, 0, 0 TO qi: rq=bs(qi): b=wo*rq
      CLS: GCOL 0: LINE 0, 360, 960, 360: GCOL 9: MOVE 0, 360
      FOR i=1 TO n
        x=c*SIN(i*w*dt): DRAW 2*i, 360+x
      NEXT i
      GCOL 12: MOVE 0, 360
      x=0: v=0: dv=0                                    REM No half-step since a=0
      FOR i=1 TO n                                      REM Feynman loop
        x=x+v*dt: dv=(d*SIN(w*i*dt)-b*(v+dv/2)-w2*x)*dt: v=v+dv: DRAW 2*i, 360+x
      NEXT i
      IF b>0 OR nf<>40 THEN
        tt=TIME: WHILE tt+200>TIME: ENDWHILE
        GCOL 13: MOVE 0, 360                            REM Analytic solution
        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
        cc=c*s/SQR((s-r)*(s-r)+rq*rq)
        IF rq>2 THEN
          wa=-wo*(rq/2+SQR(rq*rq/4-1))                  REM ...overdamped
          wb=-wo*(rq/2-SQR(rq*rq/4-1))
          a2=cc*(wa*SIN(a)+w*COS(a))/(wa-wb)
          a1=cc*SIN(a)-a2
          FOR i=1 TO n
            y=360+(a1*EXP(wa*i*dt)+a2*EXP(wb*i*dt)+cc*SIN(w*i*dt-a))
            DRAW 2*i, y
          NEXT i
        ELSE
          wa=wo*SQR(1-rq*rq/4)                          REM ...underdamped
          a1=cc*SIN(a): a2=(a1*b/2-w*cc*COS(a))/wa
          FOR i=1 TO n
            y=360+(EXP(-b*i*dt/2)*(a1*COS(wa*i*dt)+a2*SIN(wa*i*dt))+cc*SIN(w*i*dt-a))
            DRAW 2*i, y
          NEXT i
        ENDIF
      ENDIF
      ENDPROC
      
      DEF PROCsqr
      SYS "SendMessage", c1%, &147, 0, 0 TO fi: nf=fs(fi): w=40*wo/nf
      SYS "SendMessage", c2%, &147, 0, 0 TO qi: rq=bs(qi): b=wo*rq
      CLS: GCOL 0: LINE 0, 360, 960, 360: GCOL 9
      FOR i=1 TO n/nf-1
        LINE 2*i*nf, 270, 2*i*nf, 450
      NEXT i
      FOR i=0 TO n/2/nf-1
        x1=4*i*nf: LINE x1, 450, x1+2*nf, 450: LINE x1+2*nf, 270, x1+4*nf,270
      NEXT i
      GCOL 12: MOVE 0, 360
      x=0: v=0: dv=d*dt: v=v+dv/2                       REM Feynman half-step
      FOR i=0 TO n/nf-1                                 REM Feynman loop
        FOR j=1 TO nf
          x=x+v*dt: dv=(d-b*(v+dv/2)-w2*x)*dt: v=v+dv: DRAW 2*nf*i+2*j, 360+x
        NEXT j
        v=v-d*dt: dv=dv-d*dt: d=-d                      REM Discontinuity correction
      NEXT i
      ENDPROC