' dif4x.bas - jackord@kw.igs.net - 31 Jan 07 - Liberty Basic v4.02
' The diffraction pattern from an array in the form of a cube
'   or close-packed rhomb
' The 41x41 2-D arrays can be assembled into stacks of 1, 3, 11, or 41
'   layers (41 layers are required to make a complete cube or rhomb)
' Wavelengths are chosen to enhance diffraction from (111) planes in the
'   simple cube and from (200) planes in the rhomb (fcc lattice)

' Initialize Window
  nomainwin
  button#1, "Cube", [cube], UL, 2, 2, 58, 20
  button#1, "Rhomb", [rhomb], UL, 2, 24, 58, 20
  dim a$(4): a$(0)="N=1": a$(1)="N=3": a$(2)="N=11": a$(3)="N=41"
  combobox#1.c1, a$(, [waitHere], 2, 46, 58, 55
  WindowWidth=310     ' Pixel Scale 0-300
  WindowHeight=269    ' Pixel Scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Cube and Rhomb" for graphics_nsb as #1
  #1 "trapclose [quit]"
  #1.c1 "selectindex 1"
  dim nl(4): nl(0)=0: nl(1)=1: nl(2)=5: nl(3)=20
  pi=4*atn(1): n=20: dd=2.55E-7: la=la1: algd=0
  x0=180: y0=120: z=30: zz=z*z

[waitHere]
  wait

[cube]
  #1.c1 "selectionindex?": input #1.c1, ind: np=nl(ind-1)
  #1 "cls ; down"
  tt=time$("ms")
  nm=(2*n+1)*(2*n+1)*(2*np+1): sc=255000/nm/nm
  la=2*dd/3
  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
      b=2*pi*dd/r0/la
      for i=0-n to n                                     ' X loop Algorithm D
        t=b*i*nx
        sr=sr+cos(t): si=si+sin(t)
      next i
      qr=sr: qi=si
      for i=0-n to n                                     ' Y loop Algorithm D
        rr=cos(i*b*ny): ri=sin(i*b*ny)
        qr=qr+rr*sr-ri*si: qi=qi+rr*si+ri*sr
      next i
      if np>0 then
        sr=qr: si=qi
        for i=0-np to np
          if i=0 then i=i+1
          rr=cos(i*b*(z-r0)): ri=sin(i*b*(z-r0))
          qr=qr+rr*sr-ri*si: qi=qi+rr*si+ri*sr
        next i
      end if
      c=int((qr*qr+qi*qi)*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
  #1 "place 2 235"
  #1 "\T = "; int((time$("ms")-tt)/1000); " s"
goto [waitHere]

[rhomb]
  #1.c1 "selectionindex?": input #1.c1, ind: np=nl(ind-1)
  #1 "cls ; down"
  tt=time$("ms")
  nm=(2*n+1)*(2*n+1)*(2*np+1): sc=255000/nm/nm
  la=.8165*dd
  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
      b=2*pi*dd/r0/la
      by=b*sqr(3)/2
      for i=0-n to n                                     ' X loop Algorithm D
        t=b*i*nx
        sr=sr+cos(t): si=si+sin(t)
      next i
      qr=sr: qi=si
      for i=0-n to n                                     ' Y loop Algorithm D
        t=i*b*nx/2+i*by*ny
        rr=cos(t): ri=sin(t)
        qr=qr+rr*sr-ri*si: qi=qi+rr*si+ri*sr
      next i
      if np>0 then
        sr=qr: si=qi
        for i=0-np to np                                 ' Z loop Algorithm D
          if i=0 then i=i+1
          t=i*b*nx/2+i*by*ny/3+i*b*(z-r0)/sqr(1.5)
          rr=cos(t): ri=sin(t)
          qr=qr+rr*sr-ri*si: qi=qi+rr*si+ri*sr
        next i
      end if
      c=int((qr*qr+qi*qi)*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
  #1 "place 2 235"
  #1 "\T = "; int((time$("ms")-tt)/1000); " s"
goto [waitHere]

[quit]
  close #1
end