function ordergrating
% Diffraction grating with monochromatic incident beam, with a single selected 
% order shown. Cross section of the geometry of a diffraction grating, a common
% illustration in textbooks of optics, spectroscopy, and analytical chemistry.
% The grating surface is at the bottom of the diagram, along the x axis.
% The line labeled "Incident beam" is the direction of the incoming light
% beam.  The dotted line is the diffracted (outgoing) beam. The line  
% labeled "Order 0" is the direction of the zeroth-order diffraction, 
% at the angle of specular reflection from the grating surface.   
% Tom O'Haver, toh@umd.edu, November 2011

global alphar wavelength d RulingDensity m

% User-modifiable parameters:
wavelength = 389;  % initial value of wavelength
alphar = 1;   % initial value of angle of incidence
RulingDensity=45;       % initial value of grating ruling density, lines/mm,
d=1000000/RulingDensity;  % initial value of groove spacing, in mn
m=104;  % initial value of order

plotangles(alphar,wavelength,d,m)

% Attaches KeyPress test function to the figure.
set(gcf,'KeyPressFcn',@ReadKey)
uicontrol('Style','text')
% end of outer function
% ----------------------------SUBFUNCTIONS--------------------------------
function ReadKey(obj,eventdata)
global alphar wavelength d RulingDensity m
key=get(gcf,'CurrentCharacter');
if isscalar(key),
  switch double(key),
    case 28
        % angle down when right arrow pressed.
          alphar=alphar-alphar/20;
          if alphar>1.55,alphar=1.55;end
          plotangles(alphar,wavelength,d,m);
    case 29
        % angle up when left arrow pressed.    
        alphar=alphar+alphar/20;
         if alphar>1.55,alphar=1.55;end
        plotangles(alphar,wavelength,d,m);
    case 30
        % Zooms up when up arrow pressed.
        wavelength=wavelength+wavelength/20;
        plotangles(alphar,wavelength,d,m);
    case 31
        % Zooms down when down arrow pressed.
        wavelength=wavelength-wavelength/20;
        if wavelength<0,wavelength=0;end
        plotangles(alphar,wavelength,d,m);
    case 97
        % When 'a' key is pressed, INcreases "RulingDensity" by 1 or 5%
        if RulingDensity==0,RulingDensity=1;end
        if RulingDensity>20,
            RulingDensity=round(RulingDensity+RulingDensity./20);
        else
            RulingDensity=RulingDensity+1;
        end
        d=1000000./RulingDensity;
        plotangles(alphar,wavelength,d,m);
    case 122
        % When 'z' key is pressed, DEcreases "RulingDensity" by 1 or 5%
        if RulingDensity==0,RulingDensity=1;end
        if RulingDensity>20,
            RulingDensity=round(RulingDensity-RulingDensity./20);
        else
            RulingDensity=RulingDensity-1;
        end
        d=1000000./RulingDensity;
        plotangles(alphar,wavelength,d,m);
    case 115
        % When 's' key is pressed, INcreases order (m) by 1 or 5%
        if m>20,
            m=round(m+m./10);
        else
            m=m+1;
        end
        plotangles(alphar,wavelength,d,m);
    case 120
        % When 'x' key is pressed, DEcreases  order (m) by 1 or 5%
        if m>20,
            m=round(m-m./10);
        else
            m=m-1;
        end
        plotangles(alphar,wavelength,d,m);
     case 107
        % When 'K' key is pressed, prints out keyboard commands
        disp('Keyboard commands:')
        disp('Angle of incidence, degrees...left and right cursor arrows')
        disp('Wavelength, nm................up and down cursor arrows')
        disp('Ruling density, lines/mm......A/Z')
        disp('Diffraction order.............S/X')
    otherwise  
       UnassignedKey=double(key)
       disp('Press k to print out list of keyboard commands')
   end % switch
end % if

function plotangles(alphar,lambda,d,m)

r=pi/2; % Constant used below

% Compute the angle of the zero-order beam and the X and Y coordinates for
% the endoint of the incident and zero-order beams. 
xi=-cos(r-alphar);yi=sin(r-alphar); % Incident beam
angle0 = asin(-sin(alphar));x0=-cos(r-angle0);y0=sin(r-angle0); % Zero order

% Compute the angle of the diffracted beam and the X and Y coordinates for
% its endpoint if the angle is on scale (90 degrees or less)
angle1 = asin(m*lambda/d-sin(alphar));if imag(angle1)==0;x1=-cos(r-angle1);y1=sin(r-angle1);else x1=0;y1=0;end;
angle2 = asin(m*(lambda+1)/d-sin(alphar));if imag(angle2)==0;x2=-cos(r-angle2);y2=sin(r-angle2);else x2=0;y2=0;end;

% Draw line from 0,0 to the endpoint
plot([0 xi],[0 yi],'k',[0 x0],[0 y0],'k--',[0 x1],[0 y1],'k:',[0 x2],[0 y2],'k:','LineWidth',2)
whitebg('w')

% Compute and display the angular dispersion
AngularDispersison=asin(m*(lambda+1)/d-sin(alphar))-asin(m*lambda/d-sin(alphar));
text(-.4,1.34,['Angular Dispersion = '  num2str(AngularDispersison*360/(2*pi)) ' degrees/nm']);
text(-.4,1.4,['Diffraction order shown (S/X)= '  num2str(m)]);

% Add labels and title
text(xi,yi,'Incident beam');
if imag(angle0)==0;text(x0/2,y0/2,'Order 0');end;
if imag(angle1)==0;text(x1,y1,['Order ' num2str(m)]);end;
% Note: when the wavelength is adjusted, the order label is replaced with
% the wavelength.
title('Diffraction grating with monochromatic incident beam, single selected order shown.');
degrees=alphar*360/(2*pi);
xlabel(['Angle {\leftarrow \rightarrow} ' num2str(round(10*degrees)/10) ' degrees.   Wavelength {\uparrow \downarrow} '  num2str(round(10*lambda)/10) ' nm.  Lines/mm A/Z ' num2str(round(10*1000000/d)/10) ])
whitebg('w')
h2=gca;axis([-1 1 0 1.5]);