REM nm5bbc.bas - jackord@kw.igs.net - revised 15 Feb 08 - BBC Basic v5.50b
      REM pulse on n-segment string with fixed ends for n=48, 96, 192, 384
      REM time dependence from Newton's Law or FFSS (lowest 64 modes)
      REM   and dispersion relation
      REM symmetry used to speed FFSS calculation
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      b1% = FN_button("Pulse A", 10, 5, 60, 20, FN_setproc(PROCp1), 0)
      b2% = FN_button("Pulse B", 80, 5, 60, 20, FN_setproc(PROCp2), 0)
      b3% = FN_button("Motion", 150, 5, 60, 20, FN_setproc(PROCp3), 0)
      c1% = FN_combobox("", 220, 5, 70, 120, 0, 3)
      c2% = FN_combobox("", 300, 5, 70, 40, 0, 3)
      WindowWidth=392
      WindowHeight=260
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Pulse on an N-Segment String"
      SYS "SendMessage", c1%, &143, 0, "n=48"
      SYS "SendMessage", c1%, &143, 1, "n=96"
      SYS "SendMessage", c1%, &143, 1, "n=192"
      SYS "SendMessage", c1%, &143, 1, "n=384"
      SYS "SendMessage", c1%, &14E, 0, 0
      SYS "SendMessage", c2%, &143, 0, "F=ma"
      SYS "SendMessage", c2%, &143, 1, "FFSS"
      SYS "SendMessage", c2%, &14E, 0, 0
      DIM y(384): DIM v(384): DIM xx%(384): DIM yy%(384)
      DIM b(64): DIM d(64): DIM w(64)
      kk=0: plt=0: jt=0: nf=768
      OFF: VDU 5
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCp1: kk=1: jt=0                         : REM Pulse A
      
      DEF PROCp2: kk=2: jt=0                         : REM Pulse B
      
      DEF PROCp3: IF kk>0 THEN plt=1 ELSE ENDPROC    : REM Motion
      
      IF jt=0 THEN
        SYS "SendMessage", c1%, &147, 0, 0 TO ch1%
        SYS "SendMessage", c2%, &147, 0, 0 TO ch2%
        n=48*2^ch1%
        IF n<64 THEN nm=n ELSE nm=64
        FOR i=0 TO n                                 : REM Initialize
          xx%(i)=8+768/n*i: y(i)=0: v(i)=0
        NEXT i
        dt=n/384
        FOR i=3*n/8 TO 5*n/8
          y(i)=40*(1-COS((i-3*n/8)*PI/n*8))
          v(i)=-40*PI/n*8*SIN((i-3*n/8)*PI/n*8)
          IF kk=2 THEN
            v(i)=v(i)*COS(4*(i-3*n/8)*PI/n*8)+y(i)*4*PI/n*8*SIN(4*(i-3*n/8)*PI/n*8)
            y(i)=y(i)*COS(4*(i-3*n/8)*PI/n*8)
          ENDIF
        NEXT i
        FOR i=0 TO n
          yy%(i)=230+INT(2*y(i))
        NEXT i
        IF ch2%=1 THEN
          FOR i=1 TO nm-1                            : REM FFSS
            b(i)=0: d(i)=0: w(i)=2*SIN(PI*i/2/n)
            FOR j=1 TO n-1
              b(i)=b(i)+2*y(j)*SIN(PI*i*j/n)/n
              d(i)=d(i)+2*v(j)/w(i)*SIN(PI*i*j/n)/n
            NEXT j
          NEXT i
        ELSE
          dt=dt/4                                    : REM F=ma Feynman half-step
          FOR i=1 TO n-1
            v(i)=v(i)+(y(i-1)+y(i+1)-2*y(i))*dt/2
          NEXT i
        ENDIF
      ENDIF
      IF plt=1 THEN * REFRESH OFF
      tt=TIME
      WHILE jt=0 OR plt=1
        CLS: GCOL 0: MOVE xx%(0), 190: DRAW xx%(0), 270: MOVE xx%(n), 190: DRAW xx%(n), 270
        GCOL 12: MOVE 8, 260
        FOR i=0 TO n: yy%(i)=230+INT(2*y(i)+1): NEXT i
        FOR i=0 TO n: DRAW xx%(i), yy%(i): NEXT i    : REM Draw string
        IF n=48 THEN
          FOR i=0 TO n                               : REM Point masses if n=48
            CIRCLE FILL xx%(i)+2, yy%(i)-2, 8
          NEXT i
        ENDIF
        IF plt=1 THEN * REFRESH
        jt=jt+1
        IF ch2%=1 THEN
          FOR i=1 TO n-1                             : REM FFSS
            y(i)=0
            FOR k=1 TO nm-1 STEP 2
              y(i)=y(i)+b(k)*SIN(PI*i*k/n)*COS(w(k)*jt*dt)
            NEXT k
            FOR k=2 TO nm-2 STEP 2
              y(i)=y(i)+d(k)*SIN(PI*i*k/n)*SIN(w(k)*jt*dt)
            NEXT k
          NEXT i
        ELSE
          FOR j=1 TO 4                               : REM F=ma
            FOR i=1 TO n-1
              y(i)=y(i)+v(i)*dt
            NEXT i
            FOR i=1 TO n-1
              v(i)=v(i)+(y(i-1)+y(i+1)-2*y(i))*dt
            NEXT i
          NEXT j
        ENDIF
        WHILE (TIME-tt)<2*jt: WAIT 1: ENDWHILE
        IF jt=nf+1 THEN
          GCOL 0: MOVE 10, 40
          PRINT "Run time = "; (TIME-tt)/100; " sec"
          plt=0: kk=0: * REFRESH ON
        ENDIF
      ENDWHILE
      ENDPROC