REM mw4bbc.bas - jackord@kw.igs.net - 22 Apr 07 - BBC BASIC 5.50b
      REM wave incident (n=1) on a film (n=5/3) on a substrate (n=2-.2i)
      REM the film thickness is 3/4 of the wavelength in the film
      
      *FLOAT64
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      button1% = FN_button("Component", 10, 5, 80, 20, FN_setproc(PROCp1), 0)
      button2% = FN_button("Resultant", 110, 5, 80, 20, FN_setproc(PROCp2), 0)
      WindowWidth=480     : REM Scale 0-480
      WindowHeight=240    : REM Scale 0-240
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Film on a Dielectric"
      OFF
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCp1 : kk=1
      DEF PROCp2 : kk=2
      
      ORIGIN 0, 240
      b=PI/120: c=PI/60: jd=54
      a1=120: n=5/3: n2=2: k2=.2                     : REM Refractive indices
      rp=n2-n*n: ip=-k2                              : REM Complex arithmetic
      num=(rp*rp+ip*ip)^.5
      phn=ATN(ip/rp)+PI
      rp=n2+n*n
      den=(rp*rp+ip*ip)^.5
      phd=ATN(ip/rp)
      r1=num/den: ph1=phn-phd
      r1r=r1*COS(ph1): r1i=r1*SIN(ph1)
      num=((1+r1r)*(1+r1r)+r1i*r1i)^.5
      phn=ATN(r1i/(1+r1r))                           : REM      (ugh!)
      den=(n2*n2+k2*k2)^.5
      phd=ATN(-k2/n2)
      t=n*num/den
      ph2=PI/2+phn-phd
      a2=t/2/n*((n+n2)*(n+n2)+k2*k2)^.5
      ph3=ph2+ATN(-k2/(n+n2))
      r2=t/2/n*((n-n2)*(n-n2)+k2*k2)^.5
      ph4=ph2+ATN(k2/(n-n2))+PI
      r1=a1*r1: a2=a1*a2: r2=a1*r2: t=a1*t           : REM Amplitudes (finally)
      tt=TIME: del=1.5
      * REFRESH OFF
      FOR i=0 TO 960
        CLS: GCOL 0
        LINE 0, 0, 960, 0 : LINE 480, -240, 480, 240
        LINE 588, -240, 588, 240
        IF kk=1 THEN
          GCOL 9: MOVE 0, INT(a1*SIN(i*b))           : REM Plot components
          FOR j=1 TO 240
            DRAW 2*j, INT(a1*SIN(i*b-j*c))
          NEXT j
          GCOL 6: MOVE 0, INT(r1*SIN(i*b+ph1))
          FOR j=1 TO 240
            DRAW 2*j, INT(r1*SIN(i*b+j*c+ph1))
          NEXT j
          GCOL 13: MOVE 480, INT(a2*SIN(i*b+jd*n*c+ph3))
          FOR j=1 TO jd
            DRAW 480+2*j, INT(a2*SIN(i*b-(j-jd)*n*c+ph3))
          NEXT j
          GCOL 10: MOVE 480, INT(r2*SIN(i*b-jd*n*c+ph4))
          FOR j=1 TO jd
            DRAW 480+2*j, INT(r2*SIN(i*b+(j-jd)*n*c+ph4))
          NEXT j
        ELSE
          GCOL 12: MOVE 0, INT(a1*SIN(i*b)+r1*SIN(i*b+ph1))
          FOR j=1 TO 240                              : REM Plot resultant
            DRAW 2*j, INT(a1*SIN(i*b-j*c)+r1*SIN(i*b+j*c+ph1))
          NEXT j
          FOR j=1 TO jd
            DRAW 480+2*j, INT(a2*SIN(i*b-(j-jd)*n*c+ph3)+r2*SIN(i*b+(j-jd)*n*c+ph4))
          NEXT j
        ENDIF
        GCOL 12: MOVE 588, INT(t*SIN(i*b+ph2))      : REM Plot transmitted wave
        FOR j=jd+1 TO 240
          DRAW 480+2*j, INT(t*EXP(-k2*(j-jd)*c)*SIN(i*b-(j-jd)*n2*c+ph2))
        NEXT j
        tt=tt+del: WHILE tt>TIME: ENDWHILE
        * REFRESH
      NEXT i
      * REFRESH ON
      ENDPROC