diff options
| author | Eugeniy Mikhailov <evgmik@gmail.com> | 2014-05-13 15:35:26 -0400 |
|---|---|---|
| committer | Eugeniy Mikhailov <evgmik@gmail.com> | 2014-05-13 15:58:56 -0400 |
| commit | d1b62466d81b164656588ac1c55acff15ee9ea43 (patch) | |
| tree | cce78f90b0768361c4a268b946848842e1e47456 /opticstrain/opticstrain.m | |
| download | optics_toolkit-d1b62466d81b164656588ac1c55acff15ee9ea43.tar.gz optics_toolkit-d1b62466d81b164656588ac1c55acff15ee9ea43.zip | |
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
Diffstat (limited to 'opticstrain/opticstrain.m')
| -rw-r--r-- | opticstrain/opticstrain.m | 349 |
1 files changed, 349 insertions, 0 deletions
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;
+
+
|