' dif6.bas - jackord@kw.igs.net - 16 Oct 06 - Liberty Basic v4.02
' The intensity graph for a circular aperture of diameter d or a slit
'   width 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.

' Initialize Window
  nomainwin
  button#1, "Circle", [q1], UL, 2, 25, 45, 20
  button#1, "Slit", [q2], UL, 2, 48, 45, 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 "Circle and Slit" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim la(8)
  dim c(121)
  #1.c1 "selectindex 1"
  #1 "cls ; down ; color black ; backcolor white"
  #1 "place 72 0 ; goto 72 240 ; place 168 0 ; goto 168 240"

[waitHere]
  wait

[q1]                                                       ' Select Circle
  kk=1
goto [plot]

[q2]                                                       ' Select Slit
  kk=2
goto [plot]

[plot]
  #1.c1 "selectionindex?": input#1.c1, ind                 ' Wavelength index
  #1 "cls ; down ; color black ; backcolor white"
  #1 "place 72 0 ; goto 72 240 ; place 168 0 ; goto 168 240"
  pi=4*atn(1): n=70: d=96: dd=d/(2*n+1): z=10*d: zz=z*z: cm=0: x0=120
  redim la(8)
  redim c(121)
    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
  for nx=0 to 120                                          ' Screen X loop
    r0=sqr(zz+nx*nx): sr=0: si=0                           ' Algorithm B
    if kk=1 then
      for j=0-n to n                                       ' Circle Y loop
        yg=j*dd
        for i=0-n to n                                     ' Circle 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
    else                                                   ' Slit X loop
      for j=0-n to n
        xg=j*dd
        r=sqr(zz+(nx-xg)*(nx-xg))
        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)
      next j
    end if
    c(nx)=int(sr*sr+si*si): if c(nx)>cm then cm=c(nx)
    if kk=1 then #1 "place 172 20": #1 "\n="; nx; "/120"
    scan
  next nx
  #1 "color red"
  for nx=0 to 120                                          ' Plot Intensity
    x2=240-int(c(nx)/cm*240)
    if (nx<48 or nx>48) and x2<240 then
      #1 "place "; x0+nx; " "; x2: #1 "goto "; x0+nx; " "; 240
      #1 "place "; x0-nx; " "; x2: #1 "goto "; x0-nx; " "; 240
    end if
  next nx
goto [waitHere]

[quit]
  close #1
end