atreadbeta

ATREADBETA                reads a BETA file

ring=ATREADBETA(fname,cavipass,bendpass,quadpass,multipass)

FILENAME:    BETA file
CAVIPASS:    pass method for cavities (default IdentityPass)
BENDPASS:    pass method for dipoles (default BndMPoleSymplectic4E2Pass)
QUADPASS:    pass method for quadrupoles (default QuadMPoleFringePass)
MULTIPASS:    pass method for sextupoles (default StrMPoleSymplectic4Pass)

[superp,periods]=ATREADBETA(fname,cavipass,bendpass,quadpass)
        returns only one superperiod and the number of superperiods

See also: ATX, ATLINOPT, ATMODUL

0001 function [superp,periods]=atreadbeta(filename,cavipass,bendpass,quadpass)
0002 %ATREADBETA                reads a BETA file
0003 %
0004 %ring=ATREADBETA(fname,cavipass,bendpass,quadpass,multipass)
0005 %
0006 %FILENAME:    BETA file
0007 %CAVIPASS:    pass method for cavities (default IdentityPass)
0008 %BENDPASS:    pass method for dipoles (default BndMPoleSymplectic4E2Pass)
0009 %QUADPASS:    pass method for quadrupoles (default QuadMPoleFringePass)
0010 %MULTIPASS:    pass method for sextupoles (default StrMPoleSymplectic4Pass)
0011 %
0012 %[superp,periods]=ATREADBETA(fname,cavipass,bendpass,quadpass)
0013 %        returns only one superperiod and the number of superperiods
0014 %
0015 %See also: ATX, ATLINOPT, ATMODUL
0016 
0017 global GLOBVAL
0018 persistent fpath
0019 
0020 if isempty(fpath), fpath=getenv('DBETA'); end
0021 
0022 if nargin < 5, multipass='StrMPoleSymplectic4Pass'; end
0023 if nargin < 4, quadpass='QuadMPoleFringePass'; end
0024 if nargin < 3, bendpass='BndMPoleSymplectic4Pass'; end
0025 if nargin < 2, cavipass='IdentityPass'; end
0026 if nargin < 1, filename=''; end
0027 
0028 if isempty(filename)
0029     [fname,fpath]=uigetfile('*.str','BETA structure',[fpath filesep]);
0030     if ~ischar(fname), error('ReadBeta:NoFile','No file selected'); end
0031     filename=fullfile(fpath,fname);
0032 end
0033 
0034 fid=fopen(filename,'rt');
0035 if fid < 0
0036     error('ReadBeta:couldNotReadFile','Unable to read file ''%s''. No such file',filename);
0037 end
0038 
0039 betadelim(fid,'LIST OF ELEMENTS');
0040 GLOBVAL.E0=1E9;
0041 line=fgetl(fid);
0042 nb_elems=sscanf(line,'%d');
0043 elemtable=struct();
0044 cavilist={};
0045 for el=1:nb_elems
0046     nextelem=readelem(fid,cavipass,bendpass,quadpass,multipass);
0047     try
0048         elemtable.(nextelem.FamName)=nextelem;
0049     catch %#ok<CTCH>
0050         nextelem.FamName=['x' nextelem.FamName];
0051         elemtable.(nextelem.FamName)=nextelem;
0052     end
0053     if isfield(nextelem,'Class') && strcmp(nextelem.Class,'RFCavity')
0054         cavilist{end+1}=nextelem; %#ok<AGROW>
0055     end
0056 end
0057 if isempty(cavilist)% Prepare a default cavity
0058     cavilist{1}=atrfcavity('RFCAV',0,0,0,1,GLOBVAL.E0,cavipass);
0059 end
0060 eledict=fieldnames(elemtable);
0061 disp(['Elements processed (' num2str(nb_elems) ' elements)']);
0062 
0063 betadelim(fid,'STRUCTURE');
0064 line=fgetl(fid);
0065 nb_stru=sscanf(line,'%d');
0066 superp=cell(nb_stru,1);
0067 
0068 dipelem=[];         % Pending dipole element (waiting for exit face)
0069 anglein=0;          % Current dipole face angles
0070 angleout=0;
0071 lff=0;              % Current dipole fringe field extension
0072 displ=zeros(1,3);   % Current misalignment vector
0073 srot=0;             % Current element rotation
0074 
0075 id_stru=0;
0076 for el=1:nb_stru
0077     elcode=fscanf(fid,'%s',1);              % Select element in the table
0078     try
0079         elnum=str2double(elcode);
0080         if isfinite(elnum)
0081             elcode=eledict{elnum};
0082         end
0083         nextelem=elemtable.(elcode);
0084     catch %#ok<CTCH>
0085         error('ReadBeta:BadElem',['Cannot identify element ' elcode]);
0086     end
0087     switch nextelem.BetaCode                % Process the element
0088         case 'CO'
0089             if isempty(dipelem)             % Entrance face
0090                 anglein=nextelem.Angle;
0091                 lff=nextelem.Lff;
0092             else
0093                 angleout=nextelem.Angle;    % Exit face
0094                 id_stru=id_stru+1;          % create immediately in case of 2 adjacent CO elems
0095                 superp{id_stru}=atelem(dipelem,'EntranceAngle',anglein,...
0096                     'ExitAngle',angleout,'FullGap',0,'FringeInt',lff);
0097                 anglein=0;
0098                 angleout=0;
0099                 lff=0;
0100                 dipelem=[];
0101             end
0102         case 'RO'
0103             srot=srot+nextelem.Srot;
0104         case 'DE'
0105             displ=displ+nextelem.Displacement;
0106         otherwise
0107             if ~isempty(dipelem)
0108                 id_stru=id_stru+1;
0109                 superp{id_stru}=atelem(dipelem,'EntranceAngle',anglein,...
0110                     'ExitAngle',angleout,'FullGap',0,'FringeInt',lff);
0111                 anglein=0;
0112                 angleout=0;
0113                 lff=0;
0114                 dipelem=[];
0115             end
0116             if srot ~= 0
0117                 srollmat=mkSRotationMatrix(srot);
0118                 nextelem.R1=srollmat;
0119                 nextelem.R2=srollmat';
0120             end
0121             if max(abs(displ)) > 0
0122                 nextelem.T1([1 3 5])=displ;
0123                 nextelem.T2([1 3 5])=-displ;
0124             end
0125             if strcmp(nextelem.BetaCode,'DI')
0126                 dipelem=nextelem;
0127             else
0128                 id_stru=id_stru+1;
0129                 superp{id_stru}=nextelem;
0130             end
0131     end
0132 end
0133 superp(id_stru+1:end)=[];
0134 disp(['Structure processed (' num2str(nb_stru) ' elements)']);
0135 
0136 nper=fscanf(fid,'%d',1);
0137 fclose(fid);
0138 
0139 cavities=find(atgetcells(superp,'Class','RFCavity'));
0140 if isempty(cavities)        % add implicit cavity if necessary
0141     superp{end+1}=cavilist{1};
0142     cavities=length(superp);
0143 end
0144 superp=atsetfieldvalues(superp,'Energy',GLOBVAL.E0);
0145 
0146 if nargout >= 2
0147     periods=nper;
0148     for i=cavities'            % set cavity frequency
0149         superp{i}=tunecavity(superp{i},...
0150             findspos(superp,id_stru+1),1,nper);
0151     end
0152 else
0153     for i=cavities'            % set cavity frequency
0154         superp{i}=tunecavity(superp{i},...
0155             findspos(superp,id_stru+1),nper,nper);
0156     end
0157     if nper > 1
0158         superp=repmat(superp,1,nper);
0159     end
0160 end
0161 
0162 
0163 evalin('base','global GLOBVAL');
0164 
0165 function cav=tunecavity(cav,clength,ncell,ntot)
0166 frev=2.99792458E8/clength;
0167 if cav.HarmNumber > 1
0168     harm=ceil(cav.HarmNumber/ntot);
0169 else
0170     harm=round(1.66678*clength);
0171 end
0172 cav.Frequency=frev*harm;
0173 cav.HarmNumber=ncell*harm;
0174 
0175 function newelem=readelem(fid,cavipass,bendpass,quadpass,multipass)
0176 
0177 global GLOBVAL
0178 
0179 line=fgetl(fid);
0180 next=1;
0181 [elname,count,errmess,nl]=sscanf(line(next:end),'%s',1); %#ok<ASGLU>
0182 next=next+nl;
0183 [code,count,errmess,nl]=sscanf(line(next:end),'%s',1); %#ok<ASGLU>
0184 next=next+nl;
0185 params=sscanf(line(next:end),'%f')';
0186 params((length(params)+1):3)=0;
0187 switch (code)
0188     case 'SD'
0189         newelem=atdrift(elname,params(1));
0190     case 'QP'
0191         newelem=atquadrupole(elname,params(1),params(2),quadpass);
0192     case 'DE'
0193         newelem=atelem(atmarker(elname),'Displacement',params(1:3));
0194     case 'RO'
0195         newelem=atelem(atmarker(elname),'Srot',params(1));
0196     case 'CO'
0197         newelem=atelem(atmarker(elname),'Angle',params(1),'Lff',params(3));
0198     case 'DI'
0199         strength=-params(3)/params(2)/params(2);
0200         newelem=atsbend(elname,params(1)*params(2),params(1),strength,bendpass);
0201     case 'SX'
0202         if params(1) < 0.001
0203             code='LD3';
0204             newelem=atthinmultipole(elname,[],[0 0 params(1)*params(2)]);
0205         else
0206             newelem=atsextupole(elname,params(1),params(2),multipass);
0207         end
0208     case 'LD'
0209         order=params(2)/2;
0210         polb=[];
0211         polb(1,order)=params(1);
0212         code=[code int2str(order)];
0213         newelem=atthinmultipole(elname,[],polb);
0214     case 'LT'
0215         order=params(2)/2;
0216         pola=[];
0217         pola(1,order)=params(1);
0218         code=[code int2str(order)];
0219         newelem=atthinmultipole(elname,pola,[]);
0220     case 'PU'
0221         newelem=atmonitor(elname);
0222     case 'KI'
0223         if params(3) > 0, code='CHV'; end
0224         newelem=atcorrector(elname,0,[params(1) params(2)],'IdentityPass');
0225     case 'CA'
0226         GLOBVAL.E0=params(3);
0227         newelem=atrfcavity(elname,0,abs(params(1)),0,params(2),params(3),cavipass);
0228     otherwise
0229         newelem=atmarker(elname);
0230 end
0231 newelem.BetaCode=code;
0232 
0233 function betadelim(fid,code)
0234 while true
0235     while true
0236         line=fgetl(fid);
0237         if ~ischar(line), error('ReadBeta:EndOfFile','Encountered unexpected end of file.'); end
0238         if ~isempty(strfind(line,'***')), break; end
0239     end
0240     if ~isempty(strfind(line,code)), break; end
0241 end