' dif4.bas - jackord@kw.igs.net - 31 Jan 07 - Liberty Basic v4.02
' The diffraction pattern for light from a helium-neon laser
'   striking a rectangular (w=.05 mm, h=.033 mm) or equilateral triangular
'   aperture and falling on a screen 1 m away
'   (one pixel representing 1 mm on the screen).
' The plot is "overexposed" to make lower intensity variations visible.
' When "AlgD" is selected, Algorithm D replaces Algorithm C for rectangle.
' When "la/6" is selected, the wavelength is divided by 6.

' Initialize Window
  nomainwin
  button#1, "Rect", [q1], UL, 2, 2, 50, 20
  button#1, "Tri", [q2], UL, 2, 24, 50, 20
  checkbox#1.cha, "AlgD", [chas], [chac], 2, 46, 50, 20
  checkbox#1.chb, "l/6", [chbs], [chbc], 2, 68, 50, 20
  WindowWidth=250     ' Pixel Scale 0-240
  WindowHeight=269    ' Pixel Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Rectangle and Triangle" for graphics_nsb as #1
  #1 "trapclose [quit]"
  #1.chb, "font symbol 12 bold"
  la1=.633/1000: la2=la1/6: pi=4*atn(1): n=20: nc=13: la=la1: algd=0
  d=.05: dd=d/(2*n+1): x0=120: y0=120: z=1000: zz=z*z

[waitHere]
  wait

[chas]
  algd=1
goto [waitHere]

[chac]
  algd=0
goto [waitHere]

[chbs]
  la=la1/6
goto [waitHere]

[chbc]
  la=la1
goto [waitHere]

[q1]
  k=1
goto [plot]

[q2]
  k=2
goto [plot]

[plot]
  #1 "cls ; down"
  tt=time$("ms")
  if k=1 then                                              ' Rectangle
    nm=(2*n+1)*(2*nc+1): sc=255000/nm/nm
    for ny=0 to 119                                        ' Screen Y loop
      for nx=0 to 119                                      ' Screen X loop
        r0=sqr(zz+nx*nx+ny*ny): sr=0: si=0
        a=.5*(1+z/r0)*z/r0: b=2*pi*dd/r0/la
        if algd=1 then                                     ' Algorithm D
          for i=0-n to n                                   ' Aperture X loop
            t=b*(i*nx-nc*ny)
            sr=sr+cos(t): si=si+sin(t)
          next i
          qr=sr: qi=si
          for i=1 to 2*nc                                  ' Aperture Y loop
            rr=cos(i*b*ny): ri=sin(i*b*ny)
            qr=qr+rr*sr-ri*si: qi=qi+rr*si+ri*sr
          next i
          c=int((qr*qr+qi*qi)*a*a*sc)
        else                                               ' Algorithm C
          for j=0-nc to nc                                 ' Aperture Y loop
            yg=j*dd
            for i=0-n to n                                 ' Aperture X loop
              xg=i*dd
              t=b*(i*nx+j*ny)
              sr=sr+cos(t): si=si*sin(t)
            next i
          next j
          c=int((sr*sr+si*si)*a*a*sc)
        end if
        if c>255 then c=255
        #1 "color "; c; " 0 0"
        #1 "set "; x0+nx; " "; y0+ny                       ' Plot intensity
        #1 "set "; x0-nx; " "; y0+ny                       '   assuming symmetry
        #1 "set "; x0+nx; " "; y0-ny
        #1 "set "; x0-nx; " "; y0-ny
      next nx
      scan
    next ny
  end if
  if k=2 then
    nm=(2*n+1)*(n+1): sc=255000/nm/nm
    for ny=0-119 to 119                                    ' Screen Y loop
      for nx=0 to 119                                      ' Screen X loop
        r0=sqr(zz+nx*nx+ny*ny): sr=0: si=0: nr=0           ' Algorithm C
        a=.5*(1+z/r0)*z/r0: b=2*pi*dd/r0/la
        by=b*sqr(3)/2
        for j=0-int(2*n/3) to int(4*n/3)                   ' Aperture Y loop
          for i=0-n to n-nr
            t=b*nx*(i+.5*nr)+by*ny*j
            sr=sr+cos(t): si=si+sin(t)
          next i
          nr=nr+1
        next j
        c=int((sr*sr+si*si)*a*a*sc): if c>255 then c=255
        #1 "color "; c; " 0 0"
        #1 "set "; x0+nx; " "; y0-ny                       ' Plot intensity
        #1 "set "; x0-nx; " "; y0-ny                       '   assuming x symmetry
      next nx
      scan
    next ny
  end if
  #1 "place 2 235 ; color white ; backcolor black"
  #1 "\Time "; int((time$("ms")-tt)/1000); " s"
  if k=1 then
    #1 "place "; x0-n; " "; y0-nc
    #1 "box "; x0+n; " "; y0+nc
  else
    #1 "place "; x0-n; " "; y0+12
    #1 "goto "; x0+n; " "; y0+12
    #1 "goto "; x0; " "; y0-24
    #1 "goto "; x0-n; " "; y0+12
  end if
goto [waitHere]

[quit]
  close #1
end