' mw4.bas - jackord@kw.igs.net - revised 13 Mar 08 - Liberty Basic v4.02
' wave incident (n=1) on a film (n=5/3) on a substrate (n=2-.2i)
' the film thickness is 3/4 of the wavelength in the film

' Initialize Window
  nomainwin
  button#1, "Component", [p1], UL, 20, 10, 80, 20
  button#1, "Resultant", [p2], UL, 120, 10, 80, 20
  WindowWidth=498     ' Scale 0-480
  WindowHeight=278    ' Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Wave Incident on a Film-covered Metal" for graphics_nsb as #1
  #1 "trapclose [quit]"

[waitHere]
  wait

[p1]
  kk=1: goto [motion]
[p2]
  kk=2: goto [motion]

[motion]
  #1 "cls ; down ; color black"
  #1 "line 0 120 480 120 ; line 240 0 240 240 ; line 294 0 294 240"
  #1 "rule XOR"
  dim yi(81): dim yr(81): dim ys(81)
  dim ya(19): dim yb(19): dim yc(19): dim yt(63)
  dim zi(81): dim zr(81): dim zs(81)
  dim za(19): dim zb(19): dim zc(19): dim zt(63)
  pi=4*atn(1): dx=pi/20: dt=dx/2: jd=18
  a1=60: n=5/3: n2=2: k2=.2                         ' Refractive indices
  rp=n2-n*n: ip=0-k2                                ' Complex arithmetic
  num=(rp*rp+ip*ip)^.5
  phasen=atn(ip/rp)+pi
  rp=n2+n*n
  den=(rp*rp+ip*ip)^.5
  phased=atn(ip/rp)
  r1=num/den: phase1=phasen-phased
  r1r=r1*cos(phase1): r1i=r1*sin(phase1)
  num=((1+r1r)*(1+r1r)+r1i*r1i)^.5
  phasen=atn(r1i/(1+r1r))                            '      (ugh!)
  den=(n2*n2+k2*k2)^.5
  phased=atn(0-k2/n2)
  t=n*num/den
  phase2=pi/2+phasen-phased
  a2=t/2/n*((n+n2)*(n+n2)+k2*k2)^.5
  phase3=phase2+atn(0-k2/(n+n2))
  r2=t/2/n*((n-n2)*(n-n2)+k2*k2)^.5
  phase4=phase2+atn(k2/(n-n2))+pi
  r1=a1*r1: a2=a1*a2: r2=a1*r2: t=a1*t              ' Amplitudes (finally)
  for i=0 to 320
    tt=i*dt
    for j=0 to 80
      yi(j)=120-int(a1*sin(tt+(80-j)*dx))
      yr(j)=120-int(r1*sin(tt+(j-80)*dx+phase1))
      ys(j)=yi(j)+yr(j)-120
    next j
    for j=0 to jd
      ya(j)=120-int(a2*sin(tt-n*(j-jd)*dx+phase3))
      yb(j)=120-int(r2*sin(tt+n*(j-jd)*dx+phase4))
      yc(j)=ya(j)+yb(j)-120
    next j
    for j=0 to 62
      yt(j)=120-int(t*exp(0-k2*j*dx)*sin(tt-n2*j*dx+phase2))
    next j
    if kk=1 then
      #1 "color darkpink"
      for j=1 to 80
        if i>0 then #1 "line "; 3*j-3; " "; zi(j-1); " "; 3*j; " "; zi(j)
        #1 "line "; 3*j-3; " "; yi(j-1); " "; 3*j; " "; yi(j)
      next j
      #1 "color red"
      for j=1 to 80
        if i>0 then #1 "line "; 3*j-3; " "; zr(j-1); " "; 3*j; " "; zr(j)
        #1 "line "; 3*j-3; " "; yr(j-1); " "; 3*j; " "; yr(j)
      next j
      #1 "color yellow"
      for j=1 to jd
        if i>0 then #1 "line "; 240+3*j-3; " "; za(j-1); " "; 240+3*j; " "; za(j)
        #1 "line "; 240+3*j-3; " "; ya(j-1); " "; 240+3*j; " "; ya(j)
      next j
      #1 "color green"
      for j=1 to jd
        if i>0 then #1 "line "; 240+3*j-3; " "; zb(j-1); " "; 240+3*j; " "; zb(j)
        #1 "line "; 240+3*j-3; " "; yb(j-1); " "; 240+3*j; " "; yb(j)
      next j
    else
      #1 "color blue"
      for j=1 to 80
        if i>0 then #1 "line "; 3*j-3; " "; zs(j-1); " "; 3*j; " "; zs(j)
        #1 "line "; 3*j-3; " "; ys(j-1); " "; 3*j; " "; ys(j)
      next j
      for j=1 to jd
        if i>0 then #1 "line "; 240+3*j-3; " "; zc(j-1); " "; 240+3*j; " "; zc(j)
        #1 "line "; 240+3*j-3; " "; yc(j-1); " "; 240+3*j; " "; yc(j)
      next j
    end if
    #1 "color blue"
    for j=1 to 62
      if i>0 then #1 "line "; 294+3*j-3; " "; zt(j-1); " "; 294+3*j; " "; zt(j)
      #1 "line "; 294+3*j-3; " "; yt(j-1); " "; 294+3*j; " "; yt(j)
    next j
    for j=0 to 80
      zi(j)=yi(j): zr(j)=yr(j): zs(j)=ys(j)
    next j
    for j=0 to jd
      za(j)=ya(j): zb(j)=yb(j): zc(j)=yc(j)
    next j
    for j=0 to 62
      zt(j)=yt(j)
    next j
  next i
goto [waitHere]

[quit]
  close #1
end