' dif1lib.bas - jackord@kw.igs.net - revised 27 Apr 07 - Liberty Basic v4.02
' The Fraunhofer diffraction pattern for light from a helium-neon laser
'   striking a circular (d=.05 mm) aperture and falling on a screen 1 m away
'   (one pixel representing 1 mm on the screen).
' The pattern is calculated directly by superposition of Huygens' wavelets.
' The plot is overexposed 1000X to make lower intensity variations visible.
' Three algorithms can be selected: A and B are general algorithms,
'   C is a specific algorithm for Fraunhofer diffraction

' Initialize Window
  nomainwin
  button#1, "Alg A", [q1], UL, 2, 2, 40, 20
  button#1, "Alg B", [q2], UL, 44, 2, 40, 20
  button#1, "Alg C", [q3], UL, 86, 2, 40, 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]"
  #1 "font arial 10"
  dim c(169)

[waitHere]
  wait

[q1]
  kk=1: goto [plot]

[q2]
  kk=2: goto [plot]

[q3]
  kk=3: goto [plot]

[plot]
  #1 "cls ; down ; color black ; backcolor white"
  redim c(169)
  la=.633/1000: pi=4*atn(1): n=70                          ' Helium-neon red
  d=.05: dd=d/(2*n+1): x0=120: y0=120: z=1000: zz=z*z
  cm=pi*d*d/dd/dd/4: cm=cm*cm
  tt=time$("ms")
  if kk=1 then                                             ' Algorithm A
    for nx=0 to 168                                        ' Screen X loop
      sr=0: si=0
      for j=0-n to n                                       ' Aperture Y loop
        ygsq=j*dd*j*dd
        for i=0-n to n                                     ' Aperture X loop
          xg=i*dd
          if (xg*xg+ygsq)<=d*d/4 then
            r=sqr(zz+(nx-xg)*(nx-xg)+ygsq)
            t=2*pi*r/la: a=.5*(1+z/r)*z/r
            sr=sr+a*cos(t): si=si+a*sin(t)
          end if
        next i
        scan
      next j
      c(nx)=int((sr*sr+si*si)/cm*255000): if c(nx)>255 then c(nx)=255
      #1 "place 140 15": #1 "\nx = "; nx; "/168"
    next nx
  end if
  if kk=2 then                                             ' Algorithm B
    for nx=0 to 168                                        ' Screen X loop
      r0=sqr(zz+nx*nx): sr=0: si=0
      for j=0-n to n                                       ' Aperture Y loop
        ygsq=j*dd*j*dd
        for i=0-n to n                                     ' Aperture X loop
          xg=i*dd
          if (xg*xg+ygsq)<=d*d/4 then
            r=sqr(zz+(nx-xg)*(nx-xg)+ygsq)
            t=2*pi*(r-r0)/la: a=.5*(1+z/r)*z/r
            sr=sr+a*cos(t): si=si+a*sin(t)
          end if
        next i
        scan
      next j
      c(nx)=int((sr*sr+si*si)/cm*255000): if c(nx)>255 then c(nx)=255
      #1 "place 140 15": #1 "\nx = "; nx; "/168"
    next nx
  end if
  if kk=3 then                                             ' Algorithm C
    for nx=0 to 168                                        ' Screen X loop
      r0=sqr(zz+nx*nx): sr=0: si=0
      for j=0-n to n                                       ' Aperture Y loop
        ygsq=j*dd*j*dd
        for i=0-n to n                                     ' Aperture X loop
          xg=i*dd
          if (xg*xg+ygsq)<=d*d/4 then
            t=2*pi*nx*xg/r0/la: a=.5*(1+z/r0)*z/r0
            sr=sr+a*cos(t): si=si+a*sin(t)
          end if
        next i
        scan
      next j
      c(nx)=int((sr*sr+si*si)/cm*255000): if c(nx)>255 then c(nx)=255
      #1 "place 140 15": #1 "\nx = "; nx; "/168"
    next nx
  end if
  tc=time$("ms")
  for ny=0 to 119                                          ' Screen Y loop
    for nx=0 to 119                                        ' Screen X loop
      nr=int(sqr(nx*nx+ny*ny))
      #1 "color "; c(nr); " 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
  tp=time$("ms")
  #1 "place 2 235 ; color white ; backcolor black"
  #1 "\Calculate  "; (tc-tt)/1000; " s    Display  "; (tp-tc)/1000; " s"
goto [waitHere]

[quit]
  close #1
end