' df3lib.bas - jackord@kw.igs.net - 15 Feb 10 - Liberty Basic v4.03
' The diffraction pattern produced when light from a helium-neon laser strikes either
'   (1) an equilateral triangular aperture (side d=.05 mm),
'   (2) a rectangular aperture (.05x.025 mm), or
'   (3) a square aperture (side .05 mm)
'   and falls on a screen 1000 mm away (125 mm if "z/8" is chosen for the triangle.
' One pixel represents 1 mm on the screen.
' The plot is "overexposed" to make lower intensity variations visible.
' The numerical calculation uses
'   (1) a close-packed array for the triangle with 50 points along d (25 if "n/2" is chosen)
'   (2) a square 1024x512 grid for the rectangle
'   (3) a square 1024x1024 grid for the square.
' The rectangle and square use the fast algorithm.

' Initialize Window
  nomainwin
  button#1, "Tri", [tri], UL, 5, 5, 50, 20
  checkbox#1.ch1, "n/2", [waitHere], [waitHere], 5, 30, 50, 20
  checkbox#1.ch2, "z/8", [waitHere], [waitHere], 5, 55, 50, 20
  button#1, "Rec", [rec], UL, 5, 80, 50, 20
  checkbox#1.ch3, "Sqr", [waitHere], [waitHere], 5, 105, 50, 20
  WindowWidth=318     ' Pixel Scale 0-300
  WindowHeight=278    ' Pixel Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Rectangle and Triangle" for graphics_nsb as #1
  #1 "trapclose [quit]"
  la=.633/1000: d=.05: x0=180: y0=120: pi=4*atn(1)

[waitHere]
  wait

[tri]                                                    ' Triangle
  #1.ch1 "value? result$": if result$="set" then n=25 else n=50
  #1.ch2 "value? result$": if result$="set" then z=125 else z=1000
  #1.ch3 "reset"
  #1 "cls ; down"
  tt=time$("ms")
  dd=d/n: zz=z*z: nm=n*(n+1)/2: sc=255000/nm/nm
  for ny=0-120 to 120                                    ' Screen Y loop
    for nx=0 to 120                                      ' 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
      by=b*sqr(3)/2
      for j=0 to n-1                                     ' Aperture Y loop
        for i=1 to n-j
          t=b*nx*(i+.5*j)+by*ny*j
          sr=sr+cos(t): si=si+sin(t)
        next i
      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
  #1 "place 2 235 ; color black"
  #1 "\T="; int((time$("ms")-tt)/1000); " s"
  #1 "color white"
    #1 "place "; x0-20; " "; y0+12
    #1 "goto "; x0+20; " "; y0+12
    #1 "goto "; x0; " "; y0-24
    #1 "goto "; x0-20; " "; y0+12
goto [waitHere]

[rec]                                                    ' Rectangle and Square
  #1.ch1 "reset"
  #1.ch2 "reset"
  #1.ch3 "value? result$": if result$="set" then nc=10 else nc=9
  #1 "cls ; down"
  tt=time$("ms")
  nb=10: n=1024: dd=d/n: z=1000: zz=z*z: sc=255000
  for ny=0 to 120                                        ' Screen Y loop
    for nx=0 to 120                                      ' Screen X loop
      r0=sqr(zz+nx*nx+ny*ny)
      a=.5*(1+z/r0)*z/r0: b=2*pi*dd/r0/la
      sr=1: si=0: ni=1: sn=1
      for i=1 to nb                                      ' Aperture X loop
        t=b*ni*nx
        rr=cos(t): ri=sin(t)
        qr=sr+rr*sr-ri*si: si=si+rr*si+ri*sr
        sr=qr: ni=2*ni: sn=2*sn
      next i
      ni=1
      for i=1 to nc                                      ' Aperture X loop
        t=b*ni*ny
        rr=cos(t): ri=sin(t)
        qr=sr+rr*sr-ri*si: si=si+rr*si+ri*sr
        sr=qr: ni=2*ni: sn=2*sn
      next i
      c=int((sr*sr+si*si)*a*a*sc/sn/sn): 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
      if ny>0 then #1 "set "; x0+nx; " "; y0-ny: #1 "set "; x0-nx; " "; y0-ny
    next nx
    scan
  next ny
  #1 "place 2 235 ; color black"
  #1 "\T="; int((time$("ms")-tt)/1000); " s"
  #1 "color white"
  if nc=9 then
    #1 "place "; x0-20; " "; y0-10: #1 "box "; x0+20; " "; y0+10
  else
    #1 "place "; x0-20; " "; y0-20: #1 "box "; x0+20; " "; y0+20
  end if
goto [waitHere]

[quit]
  close #1
end