' df1lib.bas - jackord@kw.igs.net - revised 28 Jan 10 - Liberty Basic v4.03
' The Fresnel pattern for a circular aperture of diameter d viewed on a
'   screen 10*d from the aperture for various values of the wavelength.
'  The pattern is calculated directly by superposition of Huygens' wavelets.
' The plot is overexposed 2X.

' Initialize Window
  nomainwin
  button#1, "Run", [q1], UL, 2, 2, 43, 20
  WindowWidth=498     ' Pixel Scale 0-480
  WindowHeight=278    ' Pixel Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  dim a$(6)
  a$(0)="9.6": a$(1)="4.8": a$(2)="2.4"
  a$(3)="1.2": a$(4)="0.8": a$(5)="0.6"
  combobox#1.c1, a$(, [waitHere], 2, 26, 43, 100
  open "Circular Fresnel Patterns" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim c(169)
  dim la(6)
  #1.c1 "selectindex 1"
  #1 "cls ; down ; color black ; backcolor white"
  #1 "line 240 0 240 240 ; line 312 0 312 240 ; line 408 0 408 240"
  #1 "place 120 120 ; circle 48"

[waitHere]
  wait

[q1]
goto [plot]

[plot]
  #1.c1 "selectionindex?": input#1.c1, ind                 ' Wavelength index
  #1 "cls ; down ; color black ; backcolor white"
  #1 "line 240 0 240 240 ; line 312 0 312 240 ; line 408 0 408 240"
  #1 "place 120 120 ; circle 48"
  pi=4*atn(1): n=48: d=96: dd=d/(2*n+1): z=10*d: zz=z*z: cm=0
  x0=120: y0=120
  redim c(169)
  redim la(6)
    la(0)=d/10: la(1)=d/20: la(2)=d/40
    la(3)=d/80: la(4)=d/120: la(5)=d/160
  for nx=0 to 168                                          ' Screen X loop
    r0=sqr(zz+nx*nx): sr=0: si=0
    for j=0 to n                                           ' Aperture Y loop
      ygsq=j*j*dd*dd: na=int(sqr(n*n-j*j))
      for i=0-na to na                                     ' Aperture X loop
        xg=i*dd
        r=sqr(zz+(nx-xg)*(nx-xg)+ygsq)
        t=2*pi*(r-r0)/la(ind-1): a=(1+z/r)*z/r: if j=0 then a=a/2
        sr=sr+a*cos(t): si=si+a*sin(t)
      next i
    next j
    c(nx)=sr*sr+si*si: if c(nx)>cm then cm=c(nx)
  next nx
  for ny=0 to 120                                          ' Screen Y loop
    for nx=0 to 120                                        ' Screen X loop
      nr=int(sqr(nx*nx+ny*ny))
      cc=510*c(nr)/cm: if cc>255 then cc=255
      #1 "color "; cc; " 0 0"
      #1 "set "; x0+nx; " "; y0+ny                         ' Plot Pattern
      #1 "set "; x0-nx; " "; y0+ny
      #1 "set "; x0+nx; " "; y0-ny
      #1 "set "; x0-nx; " "; y0-ny
    next nx
  next ny
  #1 "color white ; place 120 120 ; circle 48"
  #1 "color red": x0=360
  for nx=0 to 120                                          ' Plot Intensity
    x2=240-int(c(nx)/cm*240)
    if nx<>48 and x2<240 then
      #1 "line "; x0+nx; " "; x2; " "; x0+nx; " 240"
      #1 "line "; x0-nx; " "; x2; " "; x0-nx; " 240"
    end if
  next nx
  goto [waitHere]

[quit]
  close #1
end