function make_grating_plots

  num = 1000;
  start_x = -10.0;
  end_x = 10.0;
  dx = (end_x - start_x) / num;

  // create a nice sine wave
  
  for i = 1 : num
    xarray(i) = start_x + dx*i;
    yarray(i) = sin(xarray(i));
    yinten(i) = yarray(i)*yarray(i);
  end

  if (0 == 1)
    plot(xarray, yarray);
    pause;
  end

  // how far should we offset each wave, in radians
  phase_off = 0.3;

  for i = 1 : num
    yarray2(i) = sin(xarray(i) + phase_off);
    ysum2(i) = yarray(i)+yarray2(i);
    yinten2(i) = ysum2(i) * ysum2(i);
  end

  if (0 == 1)
  clf;
  plot(xarray, yarray, 'b', xarray, yarray2, 'r+', xarray, ysum2, 'k');
  pause;

  clf;
  plot(xarray, yarray, 'b', xarray, yarray2, 'r+', xarray, yinten2, 'k');
  pause;
  end


  // now three waves
  for i = 1 : num
    yarray3(i) = sin(xarray(i) + 2*phase_off);
    ysum3(i) = yarray(i)+yarray2(i)+yarray3(i);
    yinten3(i) = ysum3(i) * ysum3(i);
  end

  if (0 == 1)
  clf;
  plot(xarray, yarray, 'b', xarray, yarray2, 'r+', ...
              xarray, yarray3, 'r', xarray, ysum3, 'k');
  pause;

  clf;
  plot(xarray, yarray, 'b', xarray, yarray2, 'r+', ...
              xarray, yarray3, 'r', xarray, yinten3, 'k');
  pause;
  end


  // now four waves
  for i = 1 : num
    yarray4(i) = sin(xarray(i) + 3*phase_off);
    ysum4(i) = yarray(i)+yarray2(i)+yarray3(i)+yarray4(i);
    yinten4(i) = ysum4(i) * ysum4(i);
  end



  // show normalized version of intensities for 2, 3, 4 waves
  if (0 == 1)
  clf;
  plot(xarray, yinten2/4, 'b', ...
              xarray, yinten3/9, 'r', xarray, yinten4/16, 'k');
  pause;
  end


  // now ten waves
  for i = 1 : num
    yarray5(i) = sin(xarray(i) + 4*phase_off);
    yarray6(i) = sin(xarray(i) + 5*phase_off);
    yarray7(i) = sin(xarray(i) + 6*phase_off);
    yarray8(i) = sin(xarray(i) + 7*phase_off);
    yarray9(i) = sin(xarray(i) + 8*phase_off);
    yarray10(i) = sin(xarray(i) + 9*phase_off);
    ysum10(i) = yarray(i)+yarray2(i)+yarray3(i)+yarray4(i)+yarray5(i) + ...
                    yarray6(i)+yarray7(i)+yarray8(i)+yarray9(i)+yarray10(i);
    yinten10(i) = ysum10(i) * ysum10(i);
  end


  // show normalized version of intensities for 2, 3, 4, 10 waves
  clf;
  plot(xarray, yinten2/4, 'b', ...
              xarray, yinten3/9, 'r', xarray, yinten4/16, 'k', ...
              xarray, yinten10/100, 'k+');
  pause;


  // now compute intensities for waves which have no offset
  for i = 1 : num
    zarray(i) = sin(xarray(i));
    zarray2(i) = 2*zarray(i);
    zarray3(i) = 3*zarray(i);
    zarray4(i) = 4*zarray(i);
    zarray10(i) = 10*zarray(i);

    zinten2(i) = zarray2(i)*zarray2(i);
    zinten3(i) = zarray3(i)*zarray3(i);
    zinten4(i) = zarray4(i)*zarray4(i);
    zinten10(i) = zarray10(i)*zarray10(i);
  end

  // show UN-normalized version of intensities for 2, 3, 4, 10 waves
  clf;
  plot(xarray, yinten2, 'b', xarray, zinten2, 'r');
  pause;

  clf;
  plot(xarray, yinten3, 'b', xarray, zinten3, 'r');
  pause;


  clf;
  plot(xarray, yinten4, 'b', xarray, zinten4, 'r');
  pause;

  clf;
  plot(xarray, yinten10, 'b', xarray, zinten10, 'r');
  pause;








endfunction

