Description of calccoupling

0001 function [Tilt, Eta, Ratio, ENV, DP, DL] = calccoupling
0002 %CALCCOUPLING - Calculates the coupling and tilt of the AT model
0003 %  [Eta, Tilt, EmittanceRatio, ENV, DP, DL] = calccoupling
0004 %
0005 %  OUTPUTS
0006 %  1. Tilt - Tilts of the emittance ellipse [radian]
0007 %  2. Eta - Emittance
0008 %  3. EmittanceRatio - median(EpsX) / median(EpsX)
0009 %  4-6. ENV, DP, DL - Output of ohmienvelope
0010 %
0011 %  NOTES
0012 %  1. If there are no outputs, the coupling information will be plotted.
0013 %  2. It can be helpful the draw the lattice on top of a graph (hold on; drawlattice(Offset, Height);)
0014 
0015 %
0016 %  Written by James Safranek
0017 
0018 
0019 global THERING
0020     
0021 
0022 C0 = 299792458;      % Speed of light [m/s]
0023 
0024 ati = atindex(THERING);
0025 L0 = findspos(THERING, length(THERING)+1);
0026 
0027 % HarmNumber = 936;
0028 % THERING{ati.RF}.Frequency = HarmNumber*C0/L0;
0029 % THERING{ati.RF}.PassMethod = 'CavityPass';
0030 % for i = ati.RF
0031 %     THERING{i}.Frequency = THERING{i}.HarmNumber*C0/L0;
0032 %     THERING{i}.PassMethod = 'CavityPass';
0033 % end
0034 
0035 
0036 [PassMethod, ATIndex, FamName, PassMethodOld, ATIndexOld, FamNameOld] = setradiation('On');
0037 
0038 CavityState = getcavity;
0039 setcavity('On');
0040 
0041 
0042 % Get all the AT elements that add radiation
0043 BendCell = findmemberof('BEND');
0044 iBend = family2atindex(BendCell);
0045 for ii = 1:length(iBend)
0046     if size(iBend{ii},2) > 1
0047         iBend{ii} = sort(iBend{ii}(:));
0048         nanindx = find(isnan(iBend{ii}));
0049         iBend{ii}(nanindx) = [];
0050     end
0051 end
0052 iBend = cell2mat(iBend(:));
0053 
0054 
0055 QuadCell = findmemberof('QUAD');
0056 iQuad = family2atindex(QuadCell);
0057 for ii = 1:length(iQuad)
0058     if size(iQuad{ii},2) > 1
0059         iQuad{ii} = sort(iQuad{ii}(:));
0060         nanindx = find(isnan(iQuad{ii}));
0061         iQuad{ii}(nanindx) = [];
0062     end
0063 end
0064 iQuad = cell2mat(iQuad(:));
0065 
0066 
0067 SextCell = findmemberof('SEXT');
0068 iSext = family2atindex(SextCell);
0069 for ii = 1:length(iSext)
0070     if size(iSext{ii},2) > 1
0071         iSext{ii} = sort(iSext{ii}(:));
0072         nanindx = find(isnan(iSext{ii}));
0073         iSext{ii}(nanindx) = [];
0074     end
0075 end
0076 iSext = cell2mat(iSext(:));
0077 
0078 
0079 RadiationElemIndex = sort([iBend(:); iQuad(:); iSext(:)]');
0080 %RadiationElemIndex(find(isnan(RadiationElemIndex))) = [];
0081 
0082 [ENV, DP, DL] = ohmienvelope(THERING, RadiationElemIndex, 1:length(THERING)+1);
0083 
0084 sigmas = cat(2, ENV.Sigma);
0085 Tilt = cat(2, ENV.Tilt);
0086 spos = findspos(THERING, 1:length(THERING)+1);
0087 
0088 [TwissData, tune, chrom]  = twissring(THERING, 0, 1:length(THERING)+1, 'chrom');
0089 
0090 
0091 % The the passmethods back
0092 setpassmethod(ATIndexOld, PassMethodOld);
0093 setcavity(CavityState);
0094 
0095 
0096 % Calculate tilts
0097 Beta = cat(1,TwissData.beta);
0098 spos = cat(1,TwissData.SPos);
0099 
0100 Eta = cat(2,TwissData.Dispersion);
0101 EpsX = (sigmas(1,:).^2-Eta(1,:).^2*DP^2)./Beta(:,1)';
0102 EpsY = (sigmas(2,:).^2-Eta(3,:).^2*DP^2)./Beta(:,2)';
0103 
0104 % RMS tilt
0105 TiltRMS = norm(Tilt)/sqrt(length(Tilt));
0106 EtaY = Eta(3,:);
0107 
0108 EpsX0 = mean(EpsX);
0109 EpsY0 = mean(EpsY);
0110 %EpsX0 = median(EpsX);
0111 %EpsY0 = median(EpsY);
0112 Ratio = EpsY0 / EpsX0;
0113 
0114 
0115 % Fix imaginary data
0116 % ohmienvelope seems to return complex when very close to zero
0117 if ~isreal(sigmas(1,:))
0118     % Sigma is positive so this should be ok
0119     sigmas(1,:) = abs(sigmas(1,:));
0120 end
0121 if ~isreal(sigmas(2,:))
0122     % Sigma is positive so this should be ok
0123     sigmas(2,:) = abs(sigmas(2,:));
0124 end
0125     
0126 
0127 if nargout == 0
0128     fprintf('   RMS Tilt = %f [degrees]\n', (180/pi) * TiltRMS);
0129     fprintf('   RMS Vertical Dispersion = %f [m]\n', norm(EtaY)/sqrt(length(EtaY)));
0130     fprintf('   Mean Horizontal Emittance = %f [nm rad]\n', 1e9*EpsX0);
0131     fprintf('   Mean Vertical   Emittance = %f [nm rad]\n', 1e9*EpsY0);
0132     fprintf('   Emittance Ratio = %f%% \n', 100*Ratio);
0133     
0134     %figure(1);
0135     gcf;
0136     clf reset;
0137     subplot(2,2,1);
0138     plot(spos, Tilt*180/pi, '-')
0139     set(gca,'XLim',[0 spos(end)])
0140     ylabel('Tilt [degrees]');
0141     title(sprintf('Rotation Angle of the Beam Ellipse  (RMS = %f)', (180/pi) * TiltRMS));
0142     xlabel('s - Position [m]');
0143 
0144     %figure(2);
0145     subplot(2,2,3);
0146     [AX, H1, H2] = plotyy(spos, 1e6*sigmas(1,:), spos, 1e6*sigmas(2,:));
0147 
0148     %set(H1,'Marker','.')
0149     set(AX(1),'XLim',[0 spos(end)]);
0150     %set(H2,'Marker','.')
0151     set(AX(2),'XLim',[0 spos(end)]);
0152     title('Principal Axis of the Beam Ellipse');
0153     set(get(AX(1),'Ylabel'), 'String', 'Horizontal [\mum]'); 
0154     set(get(AX(2),'Ylabel'), 'String', 'Vertical [\mum]'); 
0155     xlabel('s - Position [m]');
0156     
0157     FontSize = get(get(AX(1),'Ylabel'),'FontSize');
0158 
0159     %figure(3);
0160     subplot(2,2,2);    
0161     plot(spos, 1e9 * EpsX);
0162     title('Horizontal Emittance');
0163     ylabel(sprintf('\\fontsize{16}\\epsilon_x  \\fontsize{%d}[nm rad]', FontSize));
0164     xlabel('s - Position [m]');
0165     xaxis([0 L0]);
0166     
0167     %figure(4);
0168     subplot(2,2,4);
0169     plot(spos, 1e9 * EpsY);
0170     title('Vertical Emittance');
0171     ylabel(sprintf('\\fontsize{16}\\epsilon_y  \\fontsize{%d}[nm rad]', FontSize));
0172     xlabel('s - Position [m]');
0173     xaxis([0 L0]);
0174    
0175     h = addlabel(.75, 0, sprintf('     Emittance Ratio = %f%% ', 100*Ratio), 10);
0176     set(h,'HorizontalAlignment','Center');
0177     
0178     orient landscape;
0179 end
0180
calccoupling

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
