' dif2.bas - jackord@kw.igs.net - 31 Jan 07 - Liberty Basic v4.02
' The diffraction pattern produced when light from a helium-neon laser
'   strikes a circular (d=.05 mm) aperture and falls on a screen 1 m away
'   (one pixel representing 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, or a random array of n*n points spread over a square of side d.

' Initialize Window
  nomainwin
  button#1, "Grid", [q1], UL, 2, 2, 50, 20
  button#1, "Rand", [q2], UL, 2, 24, 50, 20
  WindowWidth=250     ' Pixel Scale 0-240
  WindowHeight=269    ' 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): n=20                          ' Helium-neon red
  d=.05: dd=d/(2*n+1): x0=120: y0=120: z=1000: zz=z*z
  nm=int((2*n+1)*(2*n+1)*pi/4): sc=255000/nm/nm
  dim x(nm+1), y(nm+1)

[waitHere]
  wait

[q1]
  k=1
goto [plot]

[q2]
  k=2
goto [plot]

[plot]
  #1 "cls ; down"
  tt=time$("ms")
  if k=1 then                                              ' Algorithm C
    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
        #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                 ' Algorithm C
        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 2 235 ; color white ; backcolor black"
  #1 "\Time "; (time$("ms")-tt)/1000; " s"
goto [waitHere]

[quit]
  close #1
end