Create design and algorithm for multiple slits diffraction

[samreid]

In https://github.com/phetsims/wave-interference/issues/7 I implemented slits with spacing and angle. I'm using the same 2D FFT as the other patterns but I'm not sure whether the pattern is correct--it seems odd and would be good to have @ariel-phet take a look. Please test a variety of angles because it angle=0 and angle >0 seem asymmetrical and probably one of them is exhibiting artifacts.

Dev version is here for review: https://www.colorado.edu/physics/phet/dev/html/wave-interference/1.0.0-dev.3/wave-interference_en.html

For example:

[ariel-phet]

@samreid the pattern is a bit difficult to evaluate because we don't really have a diffraction grating yet. A few changes we should make:

1. Instead of "spacing" we should be thinking about a "number of lines"
2. The "slits" should be much thinner. The slits should all be the minimal allowable width. By "allowable" I mean the minimum width we can use for computation (in this case I am guessing 1 pixel wide)
3. As the number of lines changes the lines should ideally always be symmetrically spaced about the center. So for instance, 2 slits would be have equal black regions on the left and right hand side.
4. Ideally the number of lines would go from 2 up to a value where the diffracting object looks like a grating.

Lets start there and see how the pattern is looking

[ariel-phet]

@samreid also once we have the horizontal pattern looking good, we may decide that angle is not necessary and it would be better to focus on being able to change wavelength as the primary learning goal. Angle is interesting, but given the way the computation is done we might get lots of artifacts as you suggest.

[samreid]

I switched to slit size = 1 pixel, centered it and converted to use numberOfLines instead of slit separation. I left the rotation control to help us test for artifacts, but rotation is no longer centering the slits. Let's discuss and decide where to go: https://www.colorado.edu/physics/phet/dev/html/wave-interference/1.0.0-dev.4/wave-interference_en.html

[ariel-phet]

@samreid something is not right...

First of all the diffraction pattern should be horizontal, but it is showing up as vertical.

If shining a "gaussian-esque" laser beam through a vertical diffraction grating, we should see a horizontal pattern something like this:

We might have to play some games to hollywood that

[samreid]

I'll try adding a gaussian distribution over the slits and see what that looks like.

[samreid]

Gaussians don't seem to be helping:

Poor man's branch here in case we want to go back to this:

// Copyright 2017, University of Colorado Boulder

/**
 * @author Sam Reid (PhET Interactive Simulations)
 */
define( function( require ) {
  'use strict';

  // modules
  var Circle = require( 'SCENERY/nodes/Circle' );
  var HBox = require( 'SCENERY/nodes/HBox' );
  var Image = require( 'SCENERY/nodes/Image' );
  var inherit = require( 'PHET_CORE/inherit' );
  var LaserPointerNode = require( 'SCENERY_PHET/LaserPointerNode' );
  var Matrix3 = require( 'DOT/Matrix3' );
  var NumberControl = require( 'SCENERY_PHET/NumberControl' );
  var Panel = require( 'SUN/Panel' );
  var Property = require( 'AXON/Property' );
  var RadioButtonGroup = require( 'SUN/buttons/RadioButtonGroup' );
  var Range = require( 'DOT/Range' );
  var Rectangle = require( 'SCENERY/nodes/Rectangle' );
  var ResetAllButton = require( 'SCENERY_PHET/buttons/ResetAllButton' );
  var ScreenView = require( 'JOIST/ScreenView' );
  var Util = require( 'DOT/Util' );
  var VBox = require( 'SCENERY/nodes/VBox' );
  var waveInterference = require( 'WAVE_INTERFERENCE/waveInterference' );

  // constants
  var width = 256;
  var height = width;

  /**
   * @param {number} x0
   * @param {number} y0
   * @param {number} sigmaX
   * @param {number} sigmaY
   * @param {number} x
   * @param {number} y
   * @returns {number}
   */
  function gaussian( x0, y0, sigmaX, sigmaY, x, y ) {
    var dx = x - x0;
    var dy = y - y0;
    var a = dx * dx / sigmaX / sigmaX;
    var b = dy * dy / sigmaY / sigmaY;
    return Math.pow( Math.E, -(a + b) / 2 );
  }

  /**
   * @param {WaveInterferenceModel} diffractionModel
   * @constructor
   */
  function DiffractionScreenView( diffractionModel ) {

    var self = this;
    ScreenView.call( this );

    this.onProperty = new Property( true );
    var laserPointerNode = new LaserPointerNode( this.onProperty, {
      left: 10, centerY: 50
    } );

    // Reset All button
    var resetAllButton = new ResetAllButton( {
      listener: function() {
        diffractionModel.reset();
      },
      right: this.layoutBounds.maxX - 10,
      bottom: this.layoutBounds.maxY - 10
    } );
    this.addChild( resetAllButton );

    this.squareWidthProperty = new Property( 10 );
    this.squareHeightProperty = new Property( 10 );

    this.numberOfLinesProperty = new Property( 10 );
    this.angleProperty = new Property( 0 );

    this.sigmaXProperty = new Property( 1.5 );
    this.sigmaYProperty = new Property( 1000 );
    this.gaussianMagnitudeProperty = new Property( 250 ); // TODO: is this useful for students?

    this.sceneProperty = new Property( 'rectangle' );
    var toggleButtonsContent = [ {
      value: 'rectangle',
      node: new Rectangle( 0, 0, 20, 20, { fill: 'black' } )
    }, {
      value: 'circle',
      node: new Circle( 10, { fill: 'black' } )
    }, {
      value: 'slits',
      node: new HBox( {
        spacing: 2,
        children: _.range( 1, 8 ).map( function( r ) {
          return new Rectangle( 0, 0, 2, 20, { fill: 'black' } );
        } )
      } )
    } ];
    var radioButtonGroup = new RadioButtonGroup( this.sceneProperty, toggleButtonsContent, {
      left: 10,
      bottom: this.layoutBounds.bottom - 10
    } );

    this.placeholderImage = document.createElement( 'canvas' );
    this.placeholderImage.width = width;
    this.placeholderImage.height = height;

    var context = this.placeholderImage.getContext( '2d' );
    context.fillStyle = 'black';
    context.fillRect( 0, 0, width, height );

    var imageScale = 1.5;
    this.apertureImage = new Image( this.placeholderImage, { scale: imageScale, top: 100, left: 140 } );
    self.addChild( this.apertureImage );

    this.diffractionImage = new Image( this.placeholderImage, {
      right: this.layoutBounds.right - 10,
      scale: imageScale,
      top: 100
    } );
    self.addChild( this.diffractionImage );

    var ICON_SCALE = 0.2;
    this.apertureIcon = new Image( this.placeholderImage, {
      scale: ICON_SCALE,
      centerY: laserPointerNode.centerY,
      centerX: this.apertureImage.centerX,
      matrix: Matrix3.affine( 1, 0.25,
        0, 1,
        0, 0 )
    } );

    this.diffractionIcon = new Image( this.placeholderImage, {
      scale: ICON_SCALE,
      centerY: laserPointerNode.centerY,
      centerX: this.diffractionImage.centerX,
      matrix: Matrix3.affine( 1, 0.25,
        0, 1,
        0, 0 )
    } );
    self.addChild( this.diffractionIcon );

    var updateCanvases = function() {
      self.updateCanvases();
    };
    this.sceneProperty.lazyLink( updateCanvases );

    this.addChild( radioButtonGroup );

    this.squareWidthProperty.lazyLink( updateCanvases );
    this.squareHeightProperty.lazyLink( updateCanvases );
    this.sigmaXProperty.lazyLink( updateCanvases );
    this.sigmaYProperty.lazyLink( updateCanvases );
    this.onProperty.lazyLink( updateCanvases );
    this.numberOfLinesProperty.lazyLink( updateCanvases );
    this.angleProperty.lazyLink( updateCanvases );
    this.gaussianMagnitudeProperty.lazyLink( updateCanvases );
    this.squareControlPanel = new Panel( new VBox( {
      children: [
        new NumberControl( 'width', this.squareWidthProperty, new Range( 2, 100 ), {
          delta: 2 // avoid odd/even artifacts
        } ),
        new NumberControl( 'height', this.squareHeightProperty, new Range( 2, 100 ), {
          delta: 2 // avoid odd/even artifacts
        } ) ]
    } ), {
      centerTop: this.apertureImage.centerBottom.plusXY( 0, 5 )
    } );
    this.addChild( this.squareControlPanel );

    this.gaussianControlPanel = new Panel( new HBox( {
      children: [
        new VBox( {
          children: [
            new NumberControl( 'sigmaX', this.sigmaXProperty, new Range( 0, 200 ) ),
            new NumberControl( 'sigmaY', this.sigmaYProperty, new Range( 0, 200 ) )
          ]
        } ), new NumberControl( 'magnitude', this.gaussianMagnitudeProperty, new Range( 1, 1000 ) ) ],
    } ), {
      leftTop: this.apertureImage.leftBottom.plusXY( 0, 5 )
    } );
    this.addChild( this.gaussianControlPanel );

    this.slitsControlPanel = new Panel( new VBox( {
      children: [
        new NumberControl( 'number of lines', this.numberOfLinesProperty, new Range( 2, 200 ) ),
        new NumberControl( 'angle', this.angleProperty, new Range( 0, Math.PI * 2 ), {
          delta: 0.01
        } )
      ]
    } ), {
      leftTop: this.apertureImage.leftBottom.plusXY( 0, 5 )
    } );
    this.addChild( this.slitsControlPanel );

    this.sceneProperty.link( function( scene ) {
      self.squareControlPanel.visible = scene === 'rectangle';
      self.gaussianControlPanel.visible = scene === 'circle';
      self.slitsControlPanel.visible = scene === 'slits';
    } );

    var beamWidth = 40;
    var incidentBeam = new Rectangle( laserPointerNode.right, laserPointerNode.centerY - beamWidth / 2, this.apertureIcon.centerX - laserPointerNode.right, beamWidth, {
      fill: 'gray',
      opacity: 0.7
    } );

    var transmittedBeam = new Rectangle( this.apertureIcon.centerX, laserPointerNode.centerY - beamWidth / 2, this.diffractionIcon.centerX - this.apertureIcon.centerX, beamWidth, {
      fill: 'gray',
      opacity: 0.7
    } );

    this.onProperty.linkAttribute( incidentBeam, 'visible' );
    this.onProperty.linkAttribute( transmittedBeam, 'visible' );

    this.addChild( transmittedBeam );
    self.addChild( this.apertureIcon );
    this.addChild( incidentBeam );
    this.addChild( laserPointerNode );

    updateCanvases();
  }

  waveInterference.register( 'DiffractionScreenView', DiffractionScreenView );

  return inherit( ScreenView, DiffractionScreenView, {

    //TODO Called by the animation loop. Optional, so if your view has no animation, please delete this.
    // @public
    step: function( dt ) {
      //TODO Handle view animation here.
    },

    updateCanvases: function() {

      // Usage code from JS-Fourier-Image-Analysis/js/main.js
      var dims = [ width, height ]; // will be set later
      var cc = 9e-3; // contrast constant
      var apertureCanvas;
      var diffractionCanvas;
      var diffractionContext;

      var start = +new Date();
      // make each canvas the image's exact size
      apertureCanvas = document.createElement( 'canvas' );
      apertureCanvas.width = dims[ 0 ];
      apertureCanvas.height = dims[ 1 ];
      var apertureContext = apertureCanvas.getContext( '2d' );

      diffractionCanvas = document.createElement( 'canvas' );
      diffractionCanvas.width = dims[ 0 ];
      diffractionCanvas.height = dims[ 1 ];
      diffractionContext = diffractionCanvas.getContext( '2d' );

      apertureContext.fillStyle = 'black';
      apertureContext.fillRect( 0, 0, width, height );

      apertureContext.fillStyle = 'white';

      var i;

      if ( this.sceneProperty.value === 'rectangle' ) {
        var rectWidth = this.squareWidthProperty.value;
        var rectHeight = this.squareHeightProperty.value;
        apertureContext.fillRect( width / 2 - rectWidth / 2, width / 2 - rectHeight / 2, rectWidth, rectHeight );
      }
      else if ( this.sceneProperty.value === 'circle' ) {
        for ( i = 0; i < width; i++ ) {
          for ( var k = 0; k < height; k++ ) {
            var v = Util.clamp( Math.floor( gaussian( width / 2, height / 2, this.sigmaXProperty.value, this.sigmaYProperty.value, i, k ) * this.gaussianMagnitudeProperty.value ), 0, 255 );
            apertureContext.fillStyle = 'rgb(' + v + ',' + v + ',' + v + ')';
            apertureContext.fillRect( i, k, 1, 1 );
          }
        }
      }
      else if ( this.sceneProperty.value === 'slits' ) {

        apertureContext.rotate( this.angleProperty.value );
        var numberOfLines = this.numberOfLinesProperty.value;
        var lineSpacing = width / (numberOfLines + 1);

        for ( var m = 0; m < width; m++ ) {
          for ( var n = 0; n < height; n++ ) {

            var sum = 0;

            // Even number of lines
            if ( numberOfLines % 2 === 0 ) {

              for ( i = 0; i < numberOfLines / 2; i++ ) {

                sum += this.gaussianMagnitudeProperty.value*gaussian( width / 2 + i * lineSpacing + lineSpacing / 2, height / 2, this.sigmaXProperty.value, this.sigmaYProperty.value, m,n )
                sum += this.gaussianMagnitudeProperty.value*gaussian( width / 2 - i * lineSpacing - lineSpacing / 2, height / 2, this.sigmaXProperty.value, this.sigmaYProperty.value, m,n )

                // apertureContext.fillRect( width / 2 + i * lineSpacing + lineSpacing / 2, -1000, slitWidth, 2000 );
                // apertureContext.fillRect( width / 2 - i * lineSpacing - lineSpacing / 2, -1000, slitWidth, 2000 );
              }
            }
            else {

              // odd number of lines
              for ( i = 0; i < numberOfLines / 2; i++ ) {
                sum += this.gaussianMagnitudeProperty.value*gaussian( width / 2 + i * lineSpacing, height / 2, this.sigmaXProperty.value, this.sigmaYProperty.value, m,n );
                // apertureContext.fillRect( width / 2 + i * lineSpacing, -1000, slitWidth, 2000 );
                if ( i !== 0 ) { // only draw central line once
                  // apertureContext.fillRect( width / 2 - i * lineSpacing, -1000, slitWidth, 2000 );
                  sum += this.gaussianMagnitudeProperty.value*gaussian( width / 2 - i * lineSpacing, height / 2, this.sigmaXProperty.value, this.sigmaYProperty.value, m,n );
                }
              }
            }

            // console.log(m,n,sum);
            var v = Util.clamp( Math.floor( sum ), 0, 255 );
            apertureContext.fillStyle = 'rgb(' + v + ',' + v + ',' + v + ')';
            apertureContext.fillRect( m, n, 1, 1 );
          }
        }
      }

      // circle
      // apertureContext.beginPath();
      // apertureContext.arc( 75, 75, 5, 0, 2 * Math.PI );
      // apertureContext.fill();

      // grab the pixels
      var imageData = apertureContext.getImageData( 0, 0, dims[ 0 ], dims[ 1 ] );
      var h_es = []; // the h values
      for ( var ai = 0; ai < imageData.data.length; ai += 4 ) {

        // greyscale, so you only need every 4th value
        h_es.push( imageData.data[ ai ] );
      }

      // initialize the h values
      var h = function( n, m ) {
        if ( arguments.length === 0 ) {
          return h_es;
        }

        var idx = n * dims[ 0 ] + m;
        return h_es[ idx ];
      }; // make it a function so the code matches the math

      var duration = +new Date() - start;
      console.log( 'It took ' + duration + 'ms to draw the image.' );

      start = +new Date();

      // compute the h hat values
      var h_hats = [];
      Fourier.transform( h(), h_hats );
      h_hats = Fourier.shift( h_hats, dims );

      // get the largest magnitude
      var maxMagnitude = 0;
      for ( ai = 0; ai < h_hats.length; ai++ ) {
        var mag = h_hats[ ai ].magnitude();
        if ( mag > maxMagnitude ) {
          maxMagnitude = mag;
        }
      }

      Fourier.filter( h_hats, dims, NaN, NaN );

      // store them in a nice function to match the math
      var $h = function( k, l ) {
        if ( arguments.length === 0 ) {
          return h_hats;
        }

        var idx = k * dims[ 0 ] + l;
        return h_hats[ idx ];
      };

      // draw the pixels
      var currImageData = diffractionContext.getImageData( 0, 0, dims[ 0 ], dims[ 1 ] );
      var logOfMaxMag = Math.log( cc * maxMagnitude + 1 );
      for ( k = 0; k < dims[ 1 ]; k++ ) {
        for ( var l = 0; l < dims[ 0 ]; l++ ) {
          var idxInPixels = 4 * (dims[ 0 ] * k + l);
          currImageData.data[ idxInPixels + 3 ] = 255; // full alpha
          var color = Math.log( cc * $h( k, l ).magnitude() + 1 );
          color = Math.round( 255 * (color / logOfMaxMag) );
          // RGB are the same -> gray
          for ( var c = 0; c < 3; c++ ) { // lol c++
            currImageData.data[ idxInPixels + c ] = color;
          }
        }
      }
      diffractionContext.putImageData( currImageData, 0, 0 );

      this.apertureImage.image = apertureCanvas;
      this.diffractionImage.image = this.onProperty.value ? diffractionCanvas : this.placeholderImage;

      this.apertureIcon.image = apertureCanvas;
      this.diffractionIcon.image = this.onProperty.value ? diffractionCanvas : this.placeholderImage;

      duration = +new Date() - start;
      console.log( 'It took ' + duration + 'ms to compute the FT.' );
    }
  } );
} );

[samreid]

@ariel-phet it is not clear to me why the Fourier Transform is failing for this case. Any ideas? Want to try using dsin(theta)=m*lambda to synthesize the result? By the way, I found this page a useful reference: https://en.wikipedia.org/wiki/Diffraction#Diffraction_grating

[ariel-phet]

@samreid yes, I was thinking just using dsin(theta)=mlambda would be easiest in this case.

My guess is that to make the FFT work you need to play some games thinking about how to make a grating somewhat infinite and working in the illumination. Doesn't seem worth it when we should be able to hollywood fairly well with the diffraction grating equation.

[samreid]

Labeled for design meeting discussion to decide whether to try to compute the pattern through FFT or through analytical equation.

[samreid]

We decided performance should be better on iPad air, whether through a smaller lattice or different algorithm.

[samreid]

Doesn't seem worth it when we should be able to hollywood fairly well with the diffraction grating equation.

@arouinfar would you like to create an initial design and propose an algorithm for this, or should I put something together?

[arouinfar]

The current design for the Diffraction screen does not include a diffraction grating, so I think we can close this issue.