' df2lib.bas - jackord@kw.igs.net - revised 14 Feb 10 - Liberty Basic v4.03
' The diffraction pattern produced when light from a helium-neon laser
'   strikes a circular (d=.05 mm) aperture and falls on a screen 1000 mm
'   away (125 mm if "z/8" is chosen). One pixel represents 1 mm on the screen.
' The plot is "overexposed" to make lower intensity variations visible.
' The numerical calculation uses either a square grid with n=41 points along d
'   (21 if "n/2" is chosen), or a random array of the same number of points.

' Initialize Window
  nomainwin
  button#1, "Grid", [q1], UL, 5, 5, 50, 20
  button#1, "Rand", [q2], UL, 5, 30, 50, 20
  checkbox#1.c1, "n/2", [waitHere], [waitHere], 5, 55, 50, 20
  checkbox#1.c2, "z/8", [waitHere], [waitHere], 5, 80, 50, 20
  WindowWidth=318     ' Pixel Scale 0-300
  WindowHeight=278    ' Pixel Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Circular Aperture" for graphics_nsb as #1
  #1 "trapclose [quit]"
  la=.633/1000: pi=4*atn(1): d=.05: x0=180: y0=120

[waitHere]
  wait

[q1]
  k=1
goto [plot]

[q2]
  k=2
goto [plot]

[plot]
  #1.c1 "value?": input#1.c1, z$
  if z$="set" then n=10 else n=20
  #1.c2 "value?": input#1.c2, z$
  if z$="set" then z=125 else z=1000
  dd=d/(2*n+1): zz=z*z
  nm=int((2*n+1)*(2*n+1)*pi/4): sc=255000/nm/nm
  dim x(nm+1), y(nm+1)
  #1 "cls ; down"
  tt=time$("ms")
  if k=1 then
    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
        for j=0-n to n                                     ' Aperture Y loop
          yg=j*dd
          for i=0-n to n                                   ' Aperture X loop
            xg=i*dd
            if (xg*xg+yg*yg)<=d*d/4 then
              t=2*pi*(nx*xg+ny*yg)/r0/la: a=.5*(1+z/r0)*z/r0
              sr=sr+a*cos(t): si=si+a*sin(t)
            end if
          next i
        next j
        c=int((sr*sr+si*si)*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 symmetry
        if ny>0 then #1 "set "; x0+nx; " "; y0-ny: #1 "set "; x0-nx; " "; y0-ny
      next nx
      scan
    next ny
  end if
  if k=2 then
    redim x(nm+1): redim y(nm+1): ns=0
    while ns<nm
      xg=d*(rnd(1)-.5): yg=d*(rnd(1)-.5)
      if (xg*xg+yg*yg)<=d*d/4 then ns=ns+1: x(ns)=xg: y(ns)=yg
    wend
    for ny=0-119 to 119                                    ' Screen Y loop
      for nx=0-119 to 119                                  ' Screen X loop
        r0=sqr(zz+nx*nx+ny*ny): sr=0: si=0
        for i=1 to nm                                      ' Array loop
          t=2*pi*(nx*x(i)+ny*y(i))/r0/la: a=.5*(1+z/r0)*z/r0
          sr=sr+a*cos(t): si=si+a*sin(t)
        next i
        c=int((sr*sr+si*si)*sc): if c>255 then c=255
        #1 "color "; c; " 0 0"
        #1 "set "; x0+nx; " "; y0-ny                       ' Plot intensity
      next nx
      scan
    next ny
  end if
  #1 "place 5 235 ; color black"
  #1 "\T="; int((time$("ms")-tt)/1000); " s"
goto [waitHere]

[quit]
  close #1
end