' dif5.bas - jackord@kw.igs.net - 16 Oct 06 - Liberty Basic v4.02
' 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, 25, 35, 20
  WindowWidth=250     ' Pixel Scale 0-240
  WindowHeight=269    ' Pixel Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  dim a$(8): a$(0)="la=d/5": a$(1)="la=d/10": a$(2)="la=d/20": a$(3)="la=d/40"
  a$(4)="la=d/80": a$(5)="la=d/120": a$(6)="la=d/160": a$(7)="la=d/200"
  combobox#1.c1, a$(, [waitHere], 2, 2, 80, 120
  open "Circular Fresnel Patterns" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim c(169)
  dim la(8)
  #1.c1 "selectindex 1"
  #1 "cls ; down ; color black ; backcolor white"
  #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 "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(8)
    la(0)=d/5: la(1)=d/10: la(2)=d/20: la(3)=d/40
    la(4)=d/80: la(5)=d/120: la(6)=d/160: la(7)=d/200
  tt=time$("ms")
  for nx=0 to 168                                          ' Screen X loop
    r0=sqr(zz+nx*nx): sr=0: si=0                           ' Algorithm B
    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
          r=sqr(zz+(nx-xg)*(nx-xg)+yg*yg)
          t=2*pi*(r-r0)/la(ind-1): a=.5*(1+z/r)*z/r
          sr=sr+a*cos(t): si=si+a*sin(t)
        end if
      next i
    next j
    c(nx)=int(sr*sr+si*si): if c(nx)>cm then cm=c(nx)
    #1 "place 2 235": #1 "\nx = "; nx; "/168"
    scan
  next nx
  for nx=0 to 168
    c(nx)=int(510*c(nx)/cm): if c(nx)>255 then c(nx)=255   ' Overexpose 2X
  next nx
  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
      #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; " ms  Display "; tp-tc; " ms"
  #1 "place 120 120 ; circle 48"
goto [waitHere]

[quit]
  close #1
end