initial

The optics_toolkit code taken from
http://mercury.pr.erau.edu/~greta9a1/downloads/index.html

the older version is also available at mathwork web site
http://www.mathworks.com/matlabcentral/fileexchange/15459-basic-paraxial-optics-toolkit

5 files changed, 377 insertions, 0 deletions

diff --git a/opticstrain/demo_optics.txt b/opticstrain/demo_optics.txt
new file mode 100644
index 0000000..28ed4f8
--- /dev/null
+++ b/opticstrain/demo_optics.txt
@@ -0,0 +1,28 @@
+% Scatterometer Optics Chain Model                                                                                                                                                          
+% Author: A. Gretarsson, Jan. 2010.                                                                                                                                                         
+% Length Unit: millimeters                                                                                                                                                                  
+%                                                                                                                                                                                           
+% Position                                     Optic                                          Values                                         Program Control                                
+%                                                                                                                                                                                           
+0                                              index                                          1                                                                                             
+0                                              lambda                                         0.001064                                       pref                                           
+0                                              q                                              0  80.8243                                                                                    
+0                                              opticdiam                                      12.5                                                                                          
+0                                              mode                                           0 0                                                                                           
+740.8                                          lens                                           200                                                                                           
+924.8                                          lens                                           -100                                           profile                                        
+1349.8                                         cavityR                                        800 -1.8748e-4 Inf   1000   0.9487   0.9899    push profile                                   
+1349.8                                         cavityT                                                                                       profile                                        
+1349.8                                         cavityI                                                                                       profile                                        
+2500                                           lens                                           -250                                                                                          
+2700                                           window                                          1.46 30                                       push                                           
+3000                                           beamstop                                                                                                                                     
+%                                                                                                                                                                                           
+2770                                           cavityT                                        200 0 Inf 1000 0.9 0.99                        down pop                                       
+3300                                           lens                                           50                                                                                            
+3500                                           beamstop                                                                                                                                     
+%                                                                                                                                                                                           
+2500                                           lens                                           100                                            up up pop                                      
+2600                                           lens                                           300                                                                                           
+2800                                           lens                                           100                                                                                           
+3500                                           beamstop                                                                                                                                     
diff --git a/opticstrain/demo_optics.xlsx b/opticstrain/demo_optics.xlsx
new file mode 100644
index 0000000..0c02768
--- /dev/null
+++ b/opticstrain/demo_optics.xlsx
diff --git a/opticstrain/opticstrain help.docx b/opticstrain/opticstrain help.docx
new file mode 100644
index 0000000..c75e348
--- /dev/null
+++ b/opticstrain/opticstrain help.docx
diff --git a/opticstrain/opticstrain help.pdf b/opticstrain/opticstrain help.pdf
new file mode 100644
index 0000000..de1d0ff
--- /dev/null
+++ b/opticstrain/opticstrain help.pdf
diff --git a/opticstrain/opticstrain.m b/opticstrain/opticstrain.m
new file mode 100644
index 0000000..6d9ef0e
--- /dev/null
+++ b/opticstrain/opticstrain.m
@@ -0,0 +1,349 @@
+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;

+

+