REM dif1.bas - jackord@kw.igs.net - BBC BASIC version: richard@rtrussell.co.uk
      REM The Fraunhofer diffraction pattern for light from a helium-neon laser
      REM   striking a circular (d=.05 mm) aperture and falling on a screen 1 m away
      REM   (one pixel representing 1 mm on the screen).
      REM The pattern is calculated directly by superposition of Huygens' wavelets.
      REM The plot is overexposed 1000X to make lower intensity variations visible.
      REM Three algorithms can be selected: A and B are general algorithms,
      REM   C is a specific algorithm for Fraunhofer diffraction
      
      *FLOAT 64
      
      REM Install libraries
      INSTALL @lib$+"WINLIB5"
      
      REM Controls
      push% = FN_button("Run", 4, 4, 40, 20, FN_setproc(PROCplot), 0)
      rad1% = FN_button("Alg A", 52, 4, 58, 20, FN_setproc(PROCq0), 9)
      rad2% = FN_button("Alg B", 114, 4, 58, 20, FN_setproc(PROCq1), 9)
      rad3% = FN_button("Alg C", 176, 4, 58, 20, FN_setproc(PROCq2), 9)
      
      REM Initialize window
      WindowWidth=240                       : REM Pixel Scale 0-240
      WindowHeight=240                      : REM Pixel Scale 0-240
      VDU 23,22,WindowWidth;WindowHeight;8,12,16,128
      SYS "SetWindowText", @hwnd%, "Circular Aperture"
      OFF
      
      ON CLOSE QUIT
      
      k=0
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCq0 : k=0 : ENDPROC
      
      DEF PROCq1 : k=1 : ENDPROC
      
      DEF PROCq2 : k=2 : ENDPROC
      
      DEF PROCplot
      CLS : VDU 4
      DIM c&(169)
      
      la=.633/1000 : n=70                   : REM Helium-neon red
      d=.05: dd=d/(2*n+1): z=1000: zz=z*z
      cm=PI*d*d/dd/dd/4: cm=cm*cm: d4=d*d/4
      tt=TIME
      
      REM Algorithm A
      IF k=0 THEN
        FOR nx=0 TO 168                     : REM Screen X loop
          sr=0: si=0
          FOR j=-n TO n                     : REM Aperture Y loop
            ygsq=j*dd*j*dd
            FOR i=-n TO n                   : REM Aperture X loop
              xg=i*dd
              IF (xg*xg+ygsq)<=d4 THEN
                r=SQR(zz+(nx-xg)*(nx-xg)+ygsq)
                t=2*PI*r/la: a=.5*(1+z/r)*z/r
                sr=sr+a*COS(t): si=si+a*SIN(t)
              ENDIF
            NEXT i
          NEXT j
          cc=INT((sr*sr+si*si)/cm*255000): IF cc>255 THEN c&(nx)=255 ELSE c&(nx)=cc
          PRINT TAB(1,14) "nx = "; nx; "/168";
        NEXT nx
      ENDIF
      
      REM Algorithm B
      IF k=1 THEN
        FOR nx=0 TO 168                     : REM Screen X loop
          r0=SQR(zz+nx*nx): sr=0: si=0
          FOR j=-n TO n                     : REM Aperture Y loop
            ygsq=j*dd*j*dd
            FOR i=-n TO n                   : REM Aperture X loop
              xg=i*dd
              IF (xg*xg+ygsq)<=d4 THEN
                r=SQR(zz+(nx-xg)*(nx-xg)+ygsq)
                t=2*PI*(r-r0)/la: a=.5*(1+z/r)*z/r
                sr=sr+a*COS(t): si=si+a*SIN(t)
              ENDIF
            NEXT i
          NEXT j
          cc=INT((sr*sr+si*si)/cm*255000): IF cc>255 THEN c&(nx)=255 ELSE c&(nx)=cc
          PRINT TAB(1,14) "nx = "; nx; "/168";
        NEXT nx
      ENDIF
      
      REM Algorithm C
      IF k=2 THEN
        FOR nx=0 TO 168                     : REM Screen X loop
          r0=SQR(zz+nx*nx): sr=0: si=0
          FOR j=-n TO n                     : REM Aperture Y loop
            ygsq=j*dd*j*dd
            FOR i=-n TO n                   : REM Aperture X loop
              xg=i*dd
              IF (xg*xg+ygsq)<=d4 THEN
                t=2*PI*nx*xg/r0/la: a=.5*(1+z/r0)*z/r0
                sr=sr+a*COS(t): si=si+a*SIN(t)
              ENDIF
            NEXT i
          NEXT j
          cc=INT((sr*sr+si*si)/cm*255000): IF cc>255 THEN c&(nx)=255 ELSE c&(nx)=cc
          PRINT TAB(1,14); "nx = "; nx; "/168";
        NEXT nx
      ENDIF
      
      ORIGIN 240,240
      GCOL 1
      tc=TIME
      FOR ny=0 TO 120                       : REM Screen Y loop
        FOR nx=0 TO 120                     : REM Screen X loop
          nr=INT(SQR(nx*nx+ny*ny))
          COLOUR 1,c&(nr),0,0
          PLOT 2*nx, 2*ny                   : REM Plot Intensity
          PLOT -2*nx, 2*ny                  : REM assuming symmetry
          PLOT 2*nx, -2*ny
          PLOT -2*nx, -2*ny
        NEXT nx
      NEXT ny
      tp=TIME
      GCOL 15
      VDU 5
      MOVE -230, -203: PRINT "Cal "; (tc-tt)/100; " s     Dis "; (tp-tc)/100; " s";
      ENDPROC