function opticstrain(filename)
% Models the paraxial propagation of TEM beams through an optics train specified in a text file
%-------------------------------------------------------------------------
% PROGRAM NAME: opticstrain
% AUTHOR: Andri M. Gretarsson
% UPDATES: AMG -- First complete version 1/8/2010
%
% SYNTAX: opticstrain(filename)
%
% DESCRIPTION:
%
% The Matlab program opticstrain.m plots the width of a Gaussian beam as a 
% function of propagation distance.  The location and type of optics 
% encountered are listed in a control file with four columns: Optic location (z), 
% Type of optic, Optic parameters, Program control commands. The locations z 
% of all optics are specified with respect to the same origin. All lengths 
% should be specified in the same units as the beam parameter 
% (complex curvature) q. For example, if q is specified in millimeters, 
% so must z, and all other length parameters.
%
% See the file "opticstrain help.pdf" for further details.
%
% INPUT VARIABLES:
% filename  =   name of the control file that specifies the optics chain.
%
% OUTPUT VARIABLES: none
% 
% EXAMPLE: opticstrain('demo_optics.xlsx');
%
%-------------------------------------------------------------------------
% SYNTAX: opticstrain(filename)
%-------------------------------------------------------------------------

% Parse the input file and if not .txt file then convert to a text file.
[pathstr, name, ext] = fileparts(filename);
if strcmp(ext,'.xls') || strcmp(ext,'.xlsx')
    disp('Converting MSExcel file to text format.');
    xls2txt(filename);  
else
	if ~strcmp(ext,'.txt')
        error('Filename extension must be ''.xls'', ''.xlsx'', or ''.txt''');
	end
end
controlfile=[filename(1:end-length(ext)),'.txt'];
cf=fopen(controlfile);


% Default values of optic parameters changeable via the control file
mode=[0,0];                             % Hermite-Gaussian, TEM_00
q=29.5262i;                             % Collimated beam (e.g. w=1 mm if lambda=1064 nm)
lambda=1.064e-3;                        % mm (1.064 microns expressed in millimeters)
zobj=0;                                 % location of optic
a=1;                                    % amplitude of beam
n_med=1;                                % index of propagation medium (usually air)
lensdiam=25.4;                          % default diameter of all optics in train


% Default values of program control parameters changeable via the control file
profiledisp=0;

% Values for internal program parameters
num_bounces=150;                        % # bounces in cavity used to calculate cav. behavior
npts=300;                               % # 
nxpts=500;
yoffset=0;                              % y-position of optic centers in the optics train figure
zoffset=0;                              % z-...
linespacing=1.25;                       % Spacing between optics train branches in units of lensdiam
ctemp=[0;0;Inf;Inf;0;0;0;0;1.46;n_med;n_med;0;0]; % defaults for cavity in case length(values)<12 when cavity is called

% Initialize variables
xdom=linspace(-lensdiam/4,lensdiam/4,nxpts);
linefeed=(max(xdom)-min(xdom))*linespacing;
try close(1); end
trainfig=figure(1);                     % open a new optics train figure
fignum=1;
controlkey={};                          % structure stores the control command strings
bd=0;                                   % branch depth
zstore=[]; qstore=[];                   % heap for storing the beam parameters
linenum=0; commandnum=0;                % used in parsing the command file and command strings

cline=fgetl(cf);                        % get the first control file line
while cline ~= -1
    linenum=linenum+1;
    cline=strtrim(cline);               % remove spaces
    if ~isempty(cline) && cline(1)~='%' % ignore empty lines and lines starting with '%'
        commandnum=commandnum+1;

        % PARSE THE CONTROL FILE INPUT LINE
        cline(cline==','|cline==';')=' ';                                   % Replace punctuation with spaces
        oldz=zobj;                                                          % oldz is the position of the last optic
        [zobj n ee nn]=sscanf(cline,'%f'); cline=cline(nn:end);             % zobj is the pos. of the current optic
        [object n ee nn]=sscanf(cline,'%s',1); cline=cline(nn:end);         % get type of optic
        [values n ee nn]=sscanf(cline,'%f'); cline=cline(nn:end);           % get optic parameters
        temp=' ';
        if oldz==zobj, k=length(controlkey); else k=0; controlkey={}; end   % all control keys at same z-location are used simultaneously
        while ~isempty(temp)                                                % Parse control commands into cell string array
            k=k+1;
            [temp n ee nn]=sscanf(cline,'%s',1); cline=cline(nn:end);       % read in next string (spaces are separators) and then remove it from the line
            if ~isempty(temp), controlkey(k,1)={temp}; end                  %#ok<AGROW> %Then store the control key in the control key structure
        end

        % PARSE CONTROL KEY TO GET PROGRAM ACTIONS REQUESTED
        if sum(strcmp(controlkey,'profile')),profiledisp=1;                 % look for the word "profile" in the control key list
        else profiledisp=0; end
        if sum(strcmp(controlkey,'pref')), a=1; end                         % "pref"
        yoffset=yoffset-sum(strcmp(controlkey,'down'))*linefeed;            % "down"
        yoffset=yoffset+sum(strcmp(controlkey,'up'))*linefeed;              % "up"
        if sum(strcmp(controlkey,'pop'))                                    % "pop" ("push" is handled at the end of this loop)
            q=qstore(:,bd);
            oldz=zstore(:,bd);
            bd=bd-1;
            qstore=qstore(:,1:bd);
            zstore=zstore(:,1:bd);
        end
        
        % PERFORM ACTIONS APPROPRIATE TO THE OPTIC SPECIFIED

        % lambda
        if strcmp(object,'lambda'), lambda=values;

        % mode
        elseif strcmp(object,'mode'), mode=values;
        
        %lens diameter set
        elseif strcmp(object,'opticdiam'), lensdiam=values;
        
        % q set
        elseif strcmp(object,'q'), q=values(1)+i*values(2);

        % index set
        elseif strcmp(object,'index')||strcmp(object,'nmedium') , n_med=values;

        % lens
        elseif strcmp(object,'lens'),
            if length(values)>1, lensdiam=values(2); end
            figure(trainfig); hold on;
            if ~isempty(oldz)                                               % when a previous optic was present in figure
                zdom=linspace(0,abs(zobj-oldz),npts);                       % domain over which to plot beam
                beamplot(q,zdom,oldz,yoffset,0,n_med,lambda);               % oldz is the zoffset (z-origin) for the beam plotting.
                plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
                [q,a]=prop(q,free(zobj-oldz),mode,a);                       % update the beam parameters to reach the new optic
            end
            lensplot([lensdiam values(1)/2],5,zobj,yoffset);                % draw the new optic
            xlabel('z'); ylabel('r');
            [q,a]=prop(q,lens(values(1)),mode,a);                           % uptdate the beam param's to propagate through the new optic


        % mirror
        elseif strcmp(object,'mirr'),   [q,a]=prop(q,mirr(values)*free(zobj-oldz),mode,a);
            if length(values)>1, lensdiam=values(2); end
            figure(trainfig); hold on;
            if ~isempty(oldz)
                zdom=linspace(0,abs(zobj-oldz),npts);
                beamplot(q,zdom,oldz,yoffset,0,n_med,lambda);
                plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
                plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
                [q,a]=prop(q,free(zobj-oldz),mode,a);
            end
            lensplot(lensdiam,0,zobj,yoffset);
            xlabel('z'); ylabel('r');
            [q,a]=prop(q,lens(values(1)),mode,a);


        % window
        elseif strcmp(object,'window')
            zdomwin=linspace(0,values(2),npts);
            if ~isempty(oldz)
                zdom=linspace(0,abs(zobj-oldz),npts);
                beamplot(q,zdom,oldz,yoffset,0,n_med,lambda);
                plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
                [q,a]=prop(q,free(zobj-oldz),mode,a);
            end
            [q,a]=prop(q,fdie(n_med,values(1)),mode,a);
            beamplot(q,zdomwin,zobj,yoffset,0,values(1),lambda);
            [q,a]=prop(q,fdie(values(1),n_med)*free(values(2)),mode,a);
            plot([zdomwin(1),zdomwin(end)]+zobj,[0,0]+yoffset,'k:');
            xlabel('z'); ylabel('r');
            plot([zobj zobj],[-lensdiam/2 lensdiam/2]+yoffset,'k-');
            plot([zobj+values(2) zobj+values(2)],[-lensdiam/2 lensdiam/2]+yoffset,'k-');
            plot([zobj zobj+values(2)],[lensdiam/2 lensdiam/2]+yoffset,'k-');
            plot([zobj zobj+values(2)],[-lensdiam/2 -lensdiam/2]+yoffset,'k-');
            zobj=zobj+values(2);

        % beamstop
        elseif strcmp(object,'beamstop')
            zdom=linspace(0,abs(zobj-oldz),npts);
            if ~isempty(oldz)
                beamplot(q,zdom,oldz,yoffset,0,n_med,lambda);
                plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
                [q,a]=prop(q,free(zobj-oldz),mode,a);
            end
            xlabel('z'); ylabel('r');
            plot([zobj zobj],[-lensdiam/2 lensdiam/2]+yoffset,'k-');

        % cavity
        elseif strcmp(object(1:end-1),'cavity')||strcmp(object,'cavity')
            if ~isempty(values)                                             % if no values submitted, use zcav, qrefl, qinside and qtrans from last cavity call
                zcav=zobj;
                cvals=[values;ctemp(length(values)+1:end)];                 % use default values for cavity parameters not specified
                L=round(cvals(1)/lambda)*lambda+cvals(2);                   % round cavity length to an integer number of wavelengths and add dL
                figure(trainfig); hold on;
                if ~isempty(oldz)
                    zdom=linspace(0,abs(zcav-oldz),npts);
                    beamplot(q,zdom,oldz,yoffset,0,n_med,lambda);
                    plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
                end
                lensplot(lensdiam,0,zcav,yoffset);
                lensplot(lensdiam,0,zcav+L,yoffset);
                plot([zcav, zcav+L],[0,0]+yoffset,'k:');
                xlabel('z'); ylabel('r');
                if ~isempty(oldz)
                    [q,a]=prop(q,free(zcav-oldz),mode,a);
                end
                [qrefl,qtrans,qinside,arefl,atrans,ainside]=...
                    cav(q,mode(1),mode(2),lambda,L,cvals(3),cvals(4),...
                    cvals(5),cvals(6),cvals(7),cvals(8),num_bounces,...
                    [cvals(9),cvals(10),cvals(11)],[cvals(12),cvals(13)],a);
                zin=HermiteGaussianE([mode(1),mode(2),...                   % field profile of input beam to the cavity
                    q,lambda,a],xdom);      
                pin=trapz(xdom,abs(zin).^2);
                uno=ones(size(qrefl));
            else
                zobj=zcav;                                                  % if no cavity parameters revert to the z of the last cavity results (qref, qtrans, etc. still in memory)
            end

            % cavityI
            if strcmp(object,'cavityI')
                qinputbeam=prop(q,sdie(cvals(3),cvals(9),cvals(10))*...
                    free(cvals(12))*fdie(cvals(11),cvals(9)),mode,a);
                q=sum(qinside.*ainside)./sum(ainside);
                a=sum(ainside);
                aa=ainside; qq=qinside;                                     % contains params for _each_ traversal in cavity. Used during plotting.
                zobjnew=zcav+cvals(12);
                
            % cavityR
            elseif strcmp(object,'cavityR')
                qinputbeam=q;
                a=sum(arefl);
                aa=arefl; qq=qrefl; zobjnew=zcav;
             
            % cavityI
            else q=sum(qtrans.*atrans)./sum(atrans);
               qinputbeam=prop(q,sdie(cvals(4),cvals(10),cvals(9))*...      % end mirror
                   free(cvals(13))*fdie(cvals(9),cvals(11))*...
                   free(L)*...                                              % cavity traversal
                   sdie(cvals(3),cvals(9),cvals(10))*free(cvals(12))*...    % input mirror
                   fdie(cvals(11),cvals(9)),mode,a);
                a=sum(atrans);
                aa=atrans; qq=qtrans; zobjnew=zcav+L+cvals(12)+cvals(13);
            end
                       
            if profiledisp                                                  % cavity profile display must be handled separately due to use of qq, etc.
                fignum=fignum+1;
                figure(fignum);
                z=sum(nan2zero(HermiteGaussianE([mode(1)*uno,mode(2)*uno,...
                    qq,lambda*uno,aa],xdom)),2);                           
                p=trapz(xdom,abs(z).^2);
                trunclogic=abs(xdom)<2.5*w_(qq(1),lambda); xdomtrunc=xdom(trunclogic);
                ztrunc=z(trunclogic); zintrunc=zin(trunclogic);
                ang=(unwrap(angle([zintrunc,ztrunc])));ang=ang-ang(round(size(ang,1)/2),1);
                ang(:,2)=ang(:,2)-floor(ang(round(size(ang,1)/2),2)/2/pi)*2*pi;                
                [hax,h1,h2]=plotyy(xdom,abs([z,zin]).^2,xdomtrunc,ang);  
                set(hax(2),'xaxislocation','top');  fighandle=get(hax(1),'parent');
                set(h1,'linewidth',2);  set(h2,'linestyle',':');
                axchil=get(hax(1),'children'); 
                set(hax(1),'children',axchil(end:-1:1),'color','none','box','off');
                set(hax(2),'color','w');
                figchil=get(fighandle,'children'); set(fighandle,'children',figchil(end:-1:1));
                xlabel('Radial position'); ylabel('Intensity (power/area)');
                hleg1=legend(hax(1),'Input Intens.','Intensity',2);
                hleg2=legend(hax(2),'Input Phase','Phase',1);
                origfont=get(hax(2),'fontname');
                set(get(hax(2),'Ylabel'),'String','Phase (radians)') 
                ypiticks(hax(2),-250:1:250);
                set(hleg1,'color','w'); set(hleg2,'color','w','fontname',origfont);                
                if length(num2str(zobj,'%0.2f'))>=9,
                    figtext(0.5,1.5,['z = ',num2str(zobjnew,'%0.2e')]);
                else
                    figtext(0.5,1.5,['z = ',num2str(zobjnew,'%0.2f')]);
                end
                figtext(0.5,1,['P/P_{in} = ',num2str(p/pin,'%0.2g')]);
                hold off;
            end
            figure(trainfig); hold on;
            zdom=linspace(0,abs(L),npts);
            h=beamplot(qinside(2:end),zdom,zcav,yoffset,0,n_med,lambda); hold on;
            cc=colormap; set(h,'color',cc(1,:));
            lensplot(lensdiam,0,zcav,yoffset);
            if ~isempty(oldz)
                plot([zdom(1),zdom(end)]+oldz,[0,0]+yoffset,'k:');
            end

            
            zobj=zobjnew;
            
        % optic not recognized    
        else error(['Optic in line ',num2str(linenum),...
            ' of the control file ''',controlfile,...
            ''' was not recognized. Check syntax.']);
        end
        
        % DISPLAY BEAM PROFILE (for all optics other than "cavity")
        if  profiledisp && ~(strcmp(object(1:end-1),'cavity')||strcmp(object,'cavity'))
            fignum=fignum+1; figure(fignum);
            z=nan2zero(HermiteGaussianE([mode(1),mode(2),q,lambda,a],xdom));
            p=trapz(xdom,abs(z).^2);
            trunclogic=abs(xdom)<2.5*w_(q(1),lambda); xdomtrunc=xdom(trunclogic);
            ztrunc=z(trunclogic);
            ang=(unwrap(angle(ztrunc)));    ang=ang-ang(round(size(ang,1)/2),1);
            ang=ang-floor(ang(round(size(ang,1)/2))/2/pi)*2*pi;
            [hax,h1,h2]=plotyy(xdom,abs(z).^2,xdomtrunc,ang);
            set(hax(2),'xaxislocation','top');  fighandle=get(hax(1),'parent');
            set(h1,'linewidth',2); set(h2,'linestyle',':');
            axchil=get(hax(1),'children');
            set(hax(1),'color','none','box','off');
            set(hax(2),'color','w');
            figchil=get(fighandle,'children'); set(fighandle,'children',figchil(end:-1:1));
            xlabel('Radial position'); ylabel('Intensity (power/area)');
            hleg1=legend(hax(1),'Intensity',2);
            hleg2=legend(hax(2),'Phase',1);
            origfont=get(hax(2),'fontname');
            set(get(hax(2),'Ylabel'),'String','Phase (radians)')
            ypiticks(hax(2),-250:1:250);
            set(hleg1,'color','w'); set(hleg2,'color','w','fontname',origfont);
            if length(num2str(zobj,'%0.2f'))>=9, 
                figtext(0.5,1.5,['z = ',num2str(zobj,'%0.2e')]);
            else
                figtext(0.5,1.5,['z = ',num2str(zobj,'%0.2f')]);
            end
            hold off;
        end
        
       % PARSE CONTROL KEY for program actions to be performed after beam propagation calculation.
       % Store (push) onto the stack the z-location in the beam.
       if sum(strcmp(controlkey,'push'))
           bd=bd+1;
           qstore(:,bd)=q;
           zstore(:,bd)=zobj;
       end
       
    end
    cline=fgetl(cf);
end

% Clean up
try fclose(cf); end
figure(trainfig);  hold off;