REM df6bbc.bas - jackord@kw.igs.net - 17 Dec 2010 - BBC Basic 5.50b
      REM Transmitted and Reflected X-ray diffraction patterns from a face-centred
      REM   cubic crystal structure.
      REM Huygens' wavelets are superimposed from an nxnxn close-packed rhombohedral
      REM   array where n is scanned from 2 to 256.
      REM The incident beam is normal to a close-packed plane, and the wavelength
      REM   is set to the cube edge divided by the root of 3 to enhance diffraction
      REM   from (200) and (220) planes.
      
      *FLOAT 64
      
      REM Install libraries
      INSTALL @lib$+"WINLIB5"
      
      REM Controls
      b1% = FN_button("Tramsmit", 5, 5, 70, 20, FN_setproc(PROCq1), 0)
      b2% = FN_button("Reflect", 5, 30, 70, 20, FN_setproc(PROCq2), 0)
      
      REM Initialize window
      WindowWidth=320                                             : REM Pixel Scale 0-320
      WindowHeight=240                                            : REM Pixel Scale 0-240
      VDU 23,22,WindowWidth;WindowHeight;8,12,16,128
      *FONT Arial,10
      OFF: ON CLOSE QUIT
      ORIGIN 400, 240
      a=3.615E-7: la=a/SQR(3): d=a/SQR(2): z=30: zz=z*z: s$="FCC"
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCq1: di=1
      
      DEF PROCq2: di=-1
      
      CLS : VDU 5
      * REFRESH OFF
      FOR n=1 TO 8                                                :REM Scan loop
        FOR ny=-120 TO 120                                        :REM Screen Y loop
          FOR nx=-120 TO 120                                      :REM Screen X loop
            r0=SQR(zz+nx*nx+ny*ny)
            b=2*PI/r0/la
            sr=1: si=0: ni=1: sn=1: dx=d
            FOR i=1 TO n                                          :REM Array X loop
              t=b*ni*dx*nx: rr=COS(t): ri=SIN(t)
              qr=sr+rr*sr-ri*si: si=si+rr*si+ri*sr
              sr=qr: ni=2*ni: sn=2*sn
            NEXT i
            ni=1: dx=d/2: dy=d*SQR(3)/2
            FOR i=1 TO n                                          :REM Array Y loop
              t=b*ni*(dx*nx+dy*ny): rr=COS(t): ri=SIN(t)
              qr=sr+rr*sr-ri*si: si=si+rr*si+ri*sr
              sr=qr: ni=2*ni: sn=2*sn
            NEXT i
            ni=1: dx=d/2: dy=d/2/SQR(3): dz=d/SQR(1.5)
            FOR i=1 TO n                                          :REM Array Z loop
              t=b*ni*(dx*nx+dy*ny+dz*(di*z-r0)): rr=COS(t): ri=SIN(t)
              qr=sr+rr*sr-ri*si: si=si+rr*si+ri*sr
              sr=qr: ni=2*ni: sn=2*sn
            NEXT i
            c=2550*(sr*sr+si*si)/sn/sn: IF c>255 THEN c=255
            GCOL 1: COLOUR 1, 0, c, 0: PLOT 2*nx, 2*ny            :REM Plot Intensity
          NEXT nx
        NEXT ny
        GCOL 15: MOVE -230, -210: PRINT s$; "    n="; 2^n
        * REFRESH
      NEXT n
      * REFRESH ON
      ENDPROC