lambda = 671;
%
% February 2010
% Last updated: March 2012
%
% Paul J. Martin
%
% University of Oregon
% Atom Optics
%
% This script calculates the locations for anamorphic prism pairs
% to be used in the new unibody laser design. Placement will be
% wavelength dependent, but the beam is forced to exit normal to
% the laser's exterior surface at the proper location for fiber 
% coupling.
%
% The bisect algorithm is used three times: to constrain the
% magnification (force a circular beam), the output angle (normal to
% the laser's exterior surface), and the horizontal output location
% (at the fiber coupler).
%
% Currently, supported wavelengths (in nm) are 671, 689, 780, 852, 914, 
% 922, and 994.
%
% WARNING: This code is UGLY! Speed beat out elegance in its
% creation...To make any sense of it, you'll almost certainly
% have to refer to the diagram located in this same directory.
% Calculations of angles and line intersections are not well
% documented, but should be easy enough to reproduce.
%
% DOUBLE WARNING: This code assumes that the grating will be shimmed
% on the far end of the arm (ie, the pivot is set to be the upper
% corner in the assembly drawing), so the diffraction angle needs to  
% be smaller than 45 degrees to make this work without substantial
% changes. If you need to make an adjustment the other way, you'll
% have to fix it accordingly.
%
% To add a new wavelength, the user must alter the code by manually
% inserting (for the new target wavelength) the index of refraction 
% of SF11, the divergence angles of the diode, and the lines per mm
% of the commercially available grating which most nearly reflects
% the zeroth mode at a right angle (again, this should be smaller than
% 90 given how the shimming is assumed to work).
%
% Indices of refraction are for SF11, from http://refractiveindex.info.
% Incident angles are calculated using the grating equation, assuming
% Littrow configuration ( \lambda = 2dsin\theta ).
% The desired magnification is determined from the divergence angles 
% of the laser diodes.
%

if lambda == 780				% /lambda = 780 nm  (Rb) 
    n = 1.76552;				% index of refraction in SF11
    glines = 1800;				% closest grooves/mm for desired geometry
    divpar  = 8.7;				% sharp diode divergences (784nm nominal)
    divperp = 16;
    
elseif lambda == 689			% /lambda = 689 nm  (Sr clock)
    n = 1.77273;
    glines = 2000;
    divpar  = 8.5;				% sacher diode (690 nmnominal)
    divperp = 19;
    
elseif lambda == 994			% /lambda = 994 nm  (Sr repump)
    n = 1.75608;
    glines = 1350;
    divpar  = 10;		% sacher diode (1030nm nominal)
    divperp = 26;
    
elseif lambda == 922            % /lambda = 922 nm  (Sr trap)
    n = 1.75872;
    glines = 1500;
    divpar  = 10;		% sacher diode (920nm nominal)
    divperp = 24;
 
elseif lambda == 914            % /lambda = 914 nm  (Sr magic)
    n = 1.75904;
    glines = 1500;
    divpar  = 10;		% sacher diode (920nm nominal)
    divperp = 24;
    
elseif lambda == 852            % /lambda = 852 nm  (Cs)
    n = 1.76183;
    glines = 1596;
    divpar  = 8.5;
    divperp = 30;
    
elseif lambda == 671            % /lambda = 671 nm  (Li)
    n = 1.77444;
    glines = 2000;
    divpar  = 8;
    divperp = 22;
    
else
	fprintf('Current wavelength not supported.\n\n');
    fprintf('Edit the first line of the script so that lambda takes one of');
    fprintf('\nthe following values:\n');
	fprintf('780 nm\n689 nm\n994 nm\n922 nm\n914 nm\n');
	fprintf('\nTo add a new wavelength, edit this script accordingly.\n\n');
	exit;
end

% angle (referenced to direction of the beam exiting the laser) after
% reflection from the grating)
aincident = pi/2 - 2*asin(lambda*glines/(2e6));

% magnification required to have beam be circular after prisms
% since divergence angles are full-angle FWHM measures, divide by two)
desiredmag = tan(divpar/2*pi/180)/tan(divperp/2*pi/180);

% Refer to diagram p. 15.9 in Melles Griot catalog: a1 and a2 are the
% angles between the back side of prisms and normal to beam
% However, in this version of the calculation, we go from
% right to left in the diagram, with a1 now on the right.

pang = 30*pi/180;   % angle between prism surfaces
mag  = zeros(1,8);  % vector to store beam width at each surface

% use bisection to find a1 so we get desired magnification
% (a2 is tuned to adjust the exiting angle)
a1deglo = -10;
a1deghi = 10;
a1hi = a1deghi*pi/180;
a1lo = a1deglo*pi/180;


while ( abs(a1hi-a1lo) > 1e-6 ),
    
    a1 = (a1hi+a1lo)/2;
    theta1 = asin( sin(a1 + aincident)/n );
    theta2 = asin( n*sin(pang-theta1) );
    
    % bisect to find a2 such that output angle is zero
    a2deglo = 0;
    a2deghi = 40;
    a2hi = a2deghi*pi/180;
    a2lo = a2deglo*pi/180;
    outang=1;
    
    while ( abs(a2hi-a2lo)>3e-6 ),
        a2 = (a2hi+a2lo)/2;
        phi1 = asin( sin(theta2+a1-pang-a2)/n );
        phi2 = asin( n*sin(pang-phi1) );
        outang = -phi2+a2+pang;
        if (outang<0),
            a2hi=a2;
        else
            a2lo=a2;
        end
    end
    
    % Magnification follows from the geometry. Desired magnification 
    % depends on the transverse and longitudinal beam divergences of the 
    % laser diode.
    
    mag(1) = 1/cos(a1+aincident);
    mag(2) = mag(1)*cos(theta1);
    mag(3) = mag(2)/cos(theta1-pang);
    mag(4) = mag(3)*cos(theta2);
    mag(5) = mag(4)/cos(a1+theta2-a2-pang);
    mag(6) = mag(5)*cos(phi1);
    mag(7) = mag(6)/cos(phi1-pang);
    mag(8) = mag(7)*cos(phi2);
    
    if (mag(8) > desiredmag),
        a1hi=a1;
    else
        a1lo=a1;
    end

end

if abs(a1+aincident) < 3*pi/180 || abs(a1-a2+theta2-pang) < 3*pi/180,
    fprintf('WARNING: Danger of back reflection!!!\n');
    fprintf('First interface: %4.2f degrees from normal\n',(a1+aincident)*180/pi); 
    fprintf('Second interface: %4.2f degrees from normal\n',(a1-a2+theta2-pang)*180/pi);
end

% Output reflections angles--regardless of back reflection danger
ref1 = abs(a1+aincident)*180/pi
ref2 = abs(a1-a2+theta2-pang)*180/pi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Calculate the coordinates, using (0,0) as the center of the beam,
% originating from the center of the collimating lens. The coordinates will
% give the location of the center and either end (at half intensity) of the
% beam as it intersects the surfaces of the grating and prisms.
%
% I use the notation t,c,b to mean top, center, and bottom, respectively.
% 
% All lengths are in inches. Dimensions for the prisms are 0.433 x 0.315.
% Dimensions for the grating are 0.492 x 0.126.

% coordinates of the vertices of the grating
g = zeros (5,2);       

g(1,1) = 0.596;        % location of upper right hand corner (pivot)
g(1,2) = 0.263;		   % (here is where we assume an acute total reflection angle)
g(2,1) = g(1,1) + 0.492*sin(pi/4-aincident/2);
g(2,2) = g(1,2) - 0.492*cos(pi/4-aincident/2);
g(3,1) = g(2,1) - 0.126*sin(pi/4+aincident/2);
g(3,2) = g(2,2) - 0.126*cos(pi/4+aincident/2);
g(4,1) = g(1,1) - 0.126*sin(pi/4+aincident/2);
g(4,2) = g(1,2) - 0.126*cos(pi/4+aincident/2);
g(5,1) = g(1,1);
g(5,2) = g(1,2);

% coordinates of the vertices of the prisms
p1 = zeros (5,2);
p2 = zeros (5,2);

p1(1,2) = 1.18;        % vert. distance from grating  to right angle of p1
p1(1,1) = (g(3,1)+g(4,1))/2- p1(1,2)*tan(aincident) + .433/2;
p1(4,2) = p1(1,2) + 0.433*sin(a1);
p1(4,1) = p1(1,1) - 0.433*cos(a1);
p1(2,2) = p1(1,2) + 0.315*cos(a1);
p1(2,1) = p1(1,1) + 0.315*sin(a1);
p1(3,2) = p1(2,2) - 0.5*sin(pang-a1);
p1(3,1) = p1(2,1) - 0.5*cos(pang-a1);
p1(5,2) = p1(1,2);
p1(5,1) = p1(1,1);

p1(:,2) = -1.*p1(:,2);


% Before the exact coordinates of the second prism are found, we need to start
% tracing the beam, so that we know about where it exits the laser.
% To find the coordinates, I'll use the solve function to find where two
% lines (of the form y-y0 = m(x-x0) ) intersect.

c = zeros(7,2);   % coordinates for the beam center
m = zeros(10,1);  % slope for the lines which will intersect

% at grating
m(1) = 0;
m(2) = ( g(4,2)-g(3,2) )./( g(4,1)-g(3,1) );
c(2,1) = 1/(m(1)-m(2))*( m(1)*c(1,1) - m(2)*g(3,1) + g(3,2) - c(1,2) );
c(2,2) = m(1)*( c(2,1) - c(1,1) ) + c(1,2); 

% entering 1st prism
m(3) = tan(pi/2-aincident);
m(4) = ( p1(1,2)-p1(4,2) )./( p1(1,1)-p1(4,1) );
c(3,1) = 1/(m(3)-m(4))*( m(3)*c(2,1) - m(4)*p1(1,1) + p1(1,2) - c(2,2) );
c(3,2) = m(3)*( c(3,1) - c(2,1) ) + c(2,2); 

% exiting 1st prism
m(5) = tan(pi/2+a1-theta1);
m(6) = ( p1(2,2)-p1(3,2) )./( p1(2,1)-p1(3,1) );
c(4,1) = 1/(m(5)-m(6))*( m(5)*c(3,1) - m(6)*p1(2,1) + p1(2,2) - c(3,2) );
c(4,2) = m(5)*( c(4,1) - c(3,1) ) + c(3,2); 

% Need to locate the second prism before we can finish the labeling of the
% coordinates. The location of this prism depends on the parameter d, which
% we will adjust to make sure that the beam exits the laser properly.

dmin = 0.1;
dmax = 2;

% Yet another bisection function...

while ( abs(c(7,1)-0.83) > 1e-6 ),
    
    d = (dmin + dmax) /2;
    
    p2(1,2) = c(4,2) - d*sin( pang - a1 - theta2 + pi/2 ) + .433/2*sin(a2);
    p2(1,1) = c(4,1) + d*cos( pang - a1 - theta2 + pi/2 ) + .433/2*cos(a2);
    p2(2,2) = p2(1,2) - 0.433*sin(a2);
    p2(2,1) = p2(1,1) - 0.433*cos(a2);
    p2(3,2) = p2(2,2) - 0.315*cos(a2);
    p2(3,1) = p2(2,1) + 0.315*sin(a2);
    p2(4,2) = p2(3,2) + 0.5*sin(pang+a2);
    p2(4,1) = p2(3,1) + 0.5*cos(pang+a2);
    p2(5,2) = p2(1,2);
    p2(5,1) = p2(1,1);
    
    % entering 2nd prism
    m(7) = tan( a1 + theta2 - pang - pi/2 );
    m(8) = ( p2(1,2)-p2(2,2) )./( p2(1,1)-p2(2,1) );
    c(5,1) = 1/(m(7)-m(8))*( m(7)*c(4,1) - m(8)*p2(1,1) + p2(1,2) - c(4,2) );
    c(5,2) = m(7)*( c(5,1) - c(4,1) ) + c(4,2); 
    
    % exiting 2nd prism
    m(9) = tan( a2 + pi/2 + phi1 );
    m(10) = ( p2(3,2)-p2(4,2) )./( p2(3,1)-p2(4,1) );
    c(6,1) = 1/(m(9)-m(10))*( m(9)*c(5,1) - m(10)*p2(3,1) + p2(3,2) - c(5,2) );
    c(6,2) = m(9)*( c(6,1) - c(5,1) ) + c(5,2); 
    
    % exiting the laser
    c(7,1) = c(6,1);
    c(7,2) = c(6,2) - 1;
    
    if ( c(7,1) > 0.83 ),
        dmax = d;
    else
        dmin = d;
    end
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% To demark the entire width of the beam, just repeat the above
% calculations for the 'center' of the beam, but for the 'top' and
% 'bottom'. There is probably a cleaner way to do this, but since I already
% have the code from above, I'll just cludge it together...ugly but efficient!

% rename the center beam, so c is free to use again.
center = c;     

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up initial 'top' point
c(1,1) = 0;
c(1,2) = 0.177*tan( divperp/2 * pi/180 );

% Now, through the gauntlet...
% at grating
m(1) = 0;
m(2) = ( g(4,2)-g(3,2) )./( g(4,1)-g(3,1) );
c(2,1) = 1/(m(1)-m(2))*( m(1)*c(1,1) - m(2)*g(3,1) + g(3,2) - c(1,2) );
c(2,2) = m(1)*( c(2,1) - c(1,1) ) + c(1,2); 

% entering 1st prism
m(3) = tan(pi/2-aincident);
m(4) = ( p1(1,2)-p1(4,2) )./( p1(1,1)-p1(4,1) );
c(3,1) = 1/(m(3)-m(4))*( m(3)*c(2,1) - m(4)*p1(1,1) + p1(1,2) - c(2,2) );
c(3,2) = m(3)*( c(3,1) - c(2,1) ) + c(2,2); 

% exiting 1st prism
m(5) = tan(pi/2+a1-theta1);
m(6) = ( p1(2,2)-p1(3,2) )./( p1(2,1)-p1(3,1) );
c(4,1) = 1/(m(5)-m(6))*( m(5)*c(3,1) - m(6)*p1(2,1) + p1(2,2) - c(3,2) );
c(4,2) = m(5)*( c(4,1) - c(3,1) ) + c(3,2); 

% entering 2nd prism
m(7) = tan( a1 + theta2 - pang - pi/2 );
m(8) = ( p2(1,2)-p2(2,2) )./( p2(1,1)-p2(2,1) );
c(5,1) = 1/(m(7)-m(8))*( m(7)*c(4,1) - m(8)*p2(1,1) + p2(1,2) - c(4,2) );
c(5,2) = m(7)*( c(5,1) - c(4,1) ) + c(4,2); 

% exiting 2nd prism
m(9) = tan( a2 + pi/2 + phi1 );
m(10) = ( p2(3,2)-p2(4,2) )./( p2(3,1)-p2(4,1) );
c(6,1) = 1/(m(9)-m(10))*( m(9)*c(5,1) - m(10)*p2(3,1) + p2(3,2) - c(5,2) );
c(6,2) = m(9)*( c(6,1) - c(5,1) ) + c(5,2); 

% exiting the laser
c(7,1) = c(6,1);
c(7,2) = c(6,2) - 1;

% again, rename so that I can reuse the same code
top = c;    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up initial 'bottom' point
c(1,1) = 0;
c(1,2) = -0.177*tan( divperp/2 * pi/180 ); 

% Now, through the gauntlet...
% at grating
m(1) = 0;
m(2) = ( g(4,2)-g(3,2) )./( g(4,1)-g(3,1) );
c(2,1) = 1/(m(1)-m(2))*( m(1)*c(1,1) - m(2)*g(3,1) + g(3,2) - c(1,2) );
c(2,2) = m(1)*( c(2,1) - c(1,1) ) + c(1,2); 

% entering 1st prism
m(3) = tan(pi/2-aincident);
m(4) = ( p1(1,2)-p1(4,2) )./( p1(1,1)-p1(4,1) );
c(3,1) = 1/(m(3)-m(4))*( m(3)*c(2,1) - m(4)*p1(1,1) + p1(1,2) - c(2,2) );
c(3,2) = m(3)*( c(3,1) - c(2,1) ) + c(2,2); 

% exiting 1st prism
m(5) = tan(pi/2+a1-theta1);
m(6) = ( p1(2,2)-p1(3,2) )./( p1(2,1)-p1(3,1) );
c(4,1) = 1/(m(5)-m(6))*( m(5)*c(3,1) - m(6)*p1(2,1) + p1(2,2) - c(3,2) );
c(4,2) = m(5)*( c(4,1) - c(3,1) ) + c(3,2); 

% entering 2nd prism
m(7) = tan( a1 + theta2 - pang - pi/2 );
m(8) = ( p2(1,2)-p2(2,2) )./( p2(1,1)-p2(2,1) );
c(5,1) = 1/(m(7)-m(8))*( m(7)*c(4,1) - m(8)*p2(1,1) + p2(1,2) - c(4,2) );
c(5,2) = m(7)*( c(5,1) - c(4,1) ) + c(4,2); 

% exiting 2nd prism
m(9) = tan( a2 + pi/2 + phi1 );
m(10) = ( p2(3,2)-p2(4,2) )./( p2(3,1)-p2(4,1) );
c(6,1) = 1/(m(9)-m(10))*( m(9)*c(5,1) - m(10)*p2(3,1) + p2(3,2) - c(5,2) );
c(6,2) = m(9)*( c(6,1) - c(5,1) ) + c(5,2); 

% exiting the laser
c(7,1) = c(6,1);
c(7,2) = c(6,2) - 1;

% might as well rename this one, too
bottom = c;     

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Octave plot of grating, the prisms, and the three beam paths (or, more
% specifically, the coordinates of these)...mostly for debugging and setting
% the boundaries correctly for .eps output.
%
hold on;
axis equal;
% set plot boundaries (in inches)
bdys = [-0.5 1.5 -3 0.5];
axis(bdys);
plot(g(:,1),g(:,2));
plot(p1(:,1),p1(:,2));
plot(p2(:,1),p2(:,2));

plot(center(:,1),center(:,2));
plot(top(:,1),top(:,2));
plot(bottom(:,1),bottom(:,2));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Create .eps of beam path and prism location

% function to plot an nx2 array in the postscript file using lineto
function yout = psplotarray(file, array)
  fprintf(file, "newpath\n");
  fprintf(file, "%f in %f in moveto\n", array(1,1), array(1,2));
  for j=2:length(array);
    fprintf(file, "%f in %f in lineto\n", array(j,1), array(j,2));
  end %for
  fprintf(file, "stroke\n");
end %function


%%% write eps file, to anamorphic.lambda.eps 
psfilename = sprintf("anamorphic.%dnm.ps", lambda);

file = fopen(psfilename, "wt");
fprintf(file, "%%!PS-Adobe-2.0\n");
fprintf(file, "/in {72 mul} bind def  %% conversion to inches\n");
fprintf(file, "%f in  %f in translate\n", -bdys(1), -bdys(3));
fprintf(file, ".2 setlinewidth\n");
fprintf(file, "1 setlinecap\n");

% plot grating, prism 1, prism 2, beam path
psplotarray(file, g);
psplotarray(file, p1);
psplotarray(file, p2);
psplotarray(file, center);
psplotarray(file, top);
psplotarray(file, bottom);

% plot pocket-sized guide (manually insert coordinates...measured in VW)
% arc command in postscript: xcenter, ycenter, radius, 
% starting angle, final angle (ccw)

fprintf(file, "newpath\n");
fprintf(file, "%f in %f in moveto\n",0.155 ,-1.101 );
fprintf(file, "%f in %f in %f in 90 0 arcn\n",1.507, -1.294, 0.193 );
fprintf(file, "%f in %f in %f in 0 -90 arcn\n",1.507, -1.544, 0.193 );
fprintf(file, "%f in %f in lineto\n",1.382, -1.736 );
fprintf(file, "%f in %f in %f in 0 -90 arcn\n",1.194, -2.030, 0.188 );
fprintf(file, "%f in %f in %f in 270 180 arcn\n",0.467, -2.030, 0.188 );
fprintf(file, "%f in %f in lineto\n",0.280, -1.586 );
fprintf(file, "%f in %f in %f in 270 180 arcn\n",0.155, -1.394, 0.193 );
fprintf(file, "%f in %f in %f in 180 90 arcn\n",0.155, -1.294, 0.193 );
fprintf(file, "stroke\n");


fclose(file);
cmd = sprintf("ps2eps -f %s", psfilename);
system(cmd);
unlink(psfilename);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%