function whitegrating
% Diffraction grating geometry, showing overlap of orders +1 and +2.
% 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 polychromatic 
% light beam. The colored lines are diffracted beams of 7 wavelengths defined in 
% lines 81-89. (Second-order diffraction lines are drawn shorter so that the 
% labels will be clearer). The line labeled "Zero order" is the direction of 
% the zeroth-order diffraction, at the angle of specular reflection from the 
% grating surface.  Press K with the figure window active to see the
% commands.  The labels will be clearer if the figure window is enlarged or
% zoomed to full-screen.
% Tom O'Haver, toh@umd.edu, November 2011

global alphar d RulingDensity

% User-modifiable parameters:
alphar = .75; % initial value of angle of incidence
RulingDensity=600;        % inital value of grating Ruling density, lines/mm,
d=1000000/RulingDensity;  % initial value of groove spacing, in mn

plotangles(alphar,d)
whitebg('k')

% 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 d RulingDensity
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,d);
    case 29
        % angle up when left arrow pressed.    
        alphar=alphar+alphar/20;
         if alphar>1.55,alphar=1.55;end
        plotangles(alphar,d);
    case {97,30}
        % 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,d);
    case {122,31}
        % 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,d);
    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('Ruling density, lines/mm......A/Z or up and down cursor arrows')
    otherwise  
       UnassignedKey=double(key)
       disp('Press k to print out list of keyboard commands')
   end % switch
end % if

function plotangles(alphar,d)
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. 
lambda1=575;  % Wavelength 1, plotted as yellow line
lambda2=500;  % Wavelength 2, plotted as green line
lambda3=400;  % Wavelength 3, plotted as blue line
lambda4=287;  % Wavelength 4, plotted as dotted blue line
lambda5=200;  % Wavelength 5, plotted as dotted blue line
lambda6=700;  % Wavelength 6, plotted as magenta line
lambda7=800;  % Wavelength 7, plotted as dotted magenta line
lambda8=350;  % Wavelength 8, plotted as dotted blue line
lambda9=250;  % Wavelength 9, plotted as dotted blue line
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 each diffracted beam and the X and Y coordinates for
% the endpoint of the line if the angle is 90 degrees or less.
m=1;          % Diffraction Order 
% Wavelength 1, plotted as yellow line
angle1 = asin(m*lambda1/d-sin(alphar));if imag(angle1)==0;x1=-cos(r-angle1);y1=sin(r-angle1);else x1=0;y1=0;end;
% Wavelength 2, plotted as green line
angle2 = asin(m*lambda2/d-sin(alphar));if imag(angle2)==0;x2=-cos(r-angle2);y2=sin(r-angle2);else x2=0;y2=0;end;
% Wavelength 3, plotted as blue line
angle3 = asin(m*lambda3/d-sin(alphar));if imag(angle3)==0;x3=-cos(r-angle3);y3=sin(r-angle3);else x3=0;y3=0;end;
% Wavelength 4, plotted as dotted blue line
angle4 = asin(m*lambda4/d-sin(alphar));if imag(angle4)==0;x4=-cos(r-angle4);y4=sin(r-angle4);else x4=0;y4=0;end;
% Wavelength 5, plotted as dotted blue line
angle5 = asin(m*lambda5/d-sin(alphar));if imag(angle5)==0;x5=-cos(r-angle5);y5=sin(r-angle5);else x5=0;y5=0;end;
% Wavelength 6, plotted as magenta line
angle11 = asin(m*lambda6/d-sin(alphar));if imag(angle11)==0;x11=-cos(r-angle11);y11=sin(r-angle11);else x11=0;y11=0;end;
% Wavelength 7, plotted as dotted magenta line
angle12 = asin(m*lambda7/d-sin(alphar));if imag(angle12)==0;x12=-cos(r-angle12);y12=sin(r-angle12);else x12=0;y12=0;end;

m=2;          % Second Diffraction Order 
angle6 = asin(m*lambda1/d-sin(alphar));if imag(angle6)==0;x6=-.8.*cos(r-angle6);y6=.8.*sin(r-angle6);else x6=0;y6=0;end;
angle7 = asin(m*lambda2/d-sin(alphar));if imag(angle7)==0;x7=-.8.*cos(r-angle7);y7=.8.*sin(r-angle7);else x7=0;y7=0;end;
angle8 = asin(m*lambda3/d-sin(alphar));if imag(angle8)==0;x8=-.8.*cos(r-angle8);y8=.8.*sin(r-angle8);else x8=0;y8=0;end;
angle9 = asin(m*lambda4/d-sin(alphar));if imag(angle9)==0;x9=-.8.*cos(r-angle9);y9=.8.*sin(r-angle9);else x9=0;y9=0;end;
angle10 = asin(m*lambda5/d-sin(alphar));if imag(angle10)==0;x10=-.8.*cos(r-angle10);y10=.8.*sin(r-angle10);else x10=0;y10=0;end;
angle13 = asin(m*lambda6/d-sin(alphar));if imag(angle13)==0;x13=-.8.*cos(r-angle13);y13=.8.*sin(r-angle13);else x13=0;y13=0;end;
angle14 = asin(m*lambda8/d-sin(alphar));if imag(angle14)==0;x14=-.8.*cos(r-angle14);y14=.8.*sin(r-angle14);else x14=0;y14=0;end;
angle15 = asin(m*lambda9/d-sin(alphar));if imag(angle15)==0;x15=-.8.*cos(r-angle15);y15=.8.*sin(r-angle15);else x15=0;y15=0;end;

% Draw lines from 0,0 to the endpoints
plot([0 xi],[0 yi],'w',[0 x0],[0 y0],'w',...
    [0 x1],[0 y1],'y',[0 x2],[0 y2],'g',[0 x3],[0 y3],'c',...
    [0 x4],[0 y4],'b',[0 x5],[0 y5],'b:',[0 x11],[0 y11],'r',[0 x12],[0 y12],'m:',[0 x6],[0 y6],'y',...
    [0 x7],[0 y7],'g',[0 x8],[0 y8],'c',[0 x9],[0 y9],'b',[0 x10],[0 y10],'c:',...
    [0 x13],[0 y13],'r',[0 x14],[0 y14],'b:',[0 x15],[0 y15],'b:',[0 0],[0 1.3],'w:','LineWidth',3)
whitebg('k')

% Add labels and title
degrees=alphar*360/(2*pi);
xlabel(['Angle {\leftarrow \rightarrow} ' num2str(round(10*degrees)/10) ' degrees.   Lines/mm (\uparrow \downarrow  or A/Z) ' num2str(round(10*1000000/d)/10) ])
title('Diffraction grating with white incident beam, illustrating overlap of orders +1 and +2.');
text(xi,yi,'Incident');
text(xi,yi-yi/20,'beam');
text(x0,y0,'Zero');
text(x0,y0-y0/20,'order');
% Label diffracted beams if it is on scale (90 degrees or less)
if imag(angle1)==0;text(x1,y1,[num2str(lambda1) ' nm.1']);end;
if imag(angle1)==0;text(x1,y1-y1/25,'m=1');end;
if imag(angle2)==0;text(x2,y2,[num2str(lambda2) ' nm.']);end;
if imag(angle2)==0;text(x2,y2-y2/25,'m=1');end;
if imag(angle3)==0;text(x3,y3,[num2str(lambda3) ' nm.']);end;
if imag(angle3)==0;text(x3,y3-y3/25,'m=1');end;
if imag(angle4)==0;text(x4,y4,[num2str(lambda4) ' nm.']);end;
if imag(angle4)==0;text(x4,y4-y4/25,'m=1');end;
if imag(angle5)==0;text(x5,y5,[num2str(lambda5) ' nm.']);end;
if imag(angle5)==0;text(x5,y5-y5/25,'m=1');end;
if imag(angle6)==0;text(x6,y6,[num2str(lambda1) ' nm.']);end;
if imag(angle6)==0;text(x6,y6-y6/20,'m=2');end;
if imag(angle7)==0;text(x7,y7,[num2str(lambda2) ' nm.']);end;
if imag(angle7)==0;text(x7,y7-y7/20,'m=2');end;
if imag(angle8)==0;text(x8,y8,[num2str(lambda3) ' nm.']);end;
if imag(angle8)==0;text(x8,y8-y8/20,'m=2');end;
if imag(angle9)==0;text(x9,y9,[num2str(lambda4) ' nm.']);end;
if imag(angle9)==0;text(x9,y9-y9/20,'m=2');end;
if imag(angle10)==0;text(x10,y10,[num2str(lambda5) ' nm.']);end;
if imag(angle10)==0;text(x10,y10-y10/20,'m=2');end;
if imag(angle11)==0;text(x11,y11,[num2str(lambda6) ' nm.']);end;
if imag(angle11)==0;text(x11,y11-y11/20,'m=1');end;
if imag(angle12)==0;text(x12,y12,[num2str(lambda7) ' nm.']);end;
if imag(angle12)==0;text(x12,y12-y12/20,'m=1');end;
if imag(angle13)==0;text(x13,y13,[num2str(lambda6) ' nm.']);end;
if imag(angle13)==0;text(x13,y13-y13/20,'m=2');end;
if imag(angle14)==0;text(x14,y14,[num2str(lambda8) ' nm.']);end;
if imag(angle14)==0;text(x14,y14-y14/20,'m=2');end;
if imag(angle15)==0;text(x15,y15,[num2str(lambda9) ' nm.']);end;
if imag(angle15)==0;text(x15,y15-y15/20,'m=2');end;
axis([-1 1 0 1.5]);

text(-.5,1.45,['Wavelength diffracted at grating normal:']);
text(-.5,1.4,['  (first order)= '  num2str(round(d*sin(alphar))) ' nm']);
text(-.5,1.35,['  (second order)= '  num2str(round(d/2*sin(alphar))) ' nm']);
