From d1b62466d81b164656588ac1c55acff15ee9ea43 Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Tue, 13 May 2014 15:35:26 -0400 Subject: 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 --- opticstrain/opticstrain.m | 349 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 349 insertions(+) create mode 100644 opticstrain/opticstrain.m (limited to 'opticstrain/opticstrain.m') 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 %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; + + -- cgit v1.2.3