At Ohmienvelope

Synchrotron radiation affects the beam distribution through the diffusion and damping effects. For the diffusion effect, we get a diffusion matrix:

(1)
\begin{align} D=\pmatrix{0&0&0&0&0&0\cr 0&0&0&0&0&0\cr 0&0&0&0&0&0\cr 0&0&0&0&0&0\cr 0&0&0&0&0&0\cr 0&0&0&0&0&d} \end{align}

with

(2)
\begin{align} d=\frac{55}{48\sqrt{3}}\alpha_0\frac{\gamma^5}{\rho^3}\left(\frac{\hbar}{mc}\right)^2 \end{align}

In terms of betatron coordinates, we transform with the dispersion matrix and get

(3)
\begin{align} D_\beta = {\cal B} D {\cal B}^{T} = d\pmatrix{\eta_x^2 & \eta_x\eta_x' & \eta_x\eta_y & \eta_x\eta_y' & 0 & -\eta_x \cr \eta_x\eta_x' & \eta_x'^2 & \eta_x'\eta_y & \eta_x\eta_y' & 0 & -\eta_x'\cr \eta_x\eta_y & \eta_x'\eta_y & \eta_y^2 & \eta_y\eta_y' & 0 & -\eta_y \cr \eta_x\eta_y' & \eta_x'\eta_y' & \eta_y\eta_y' & \eta_y'^2 & 0 & -\eta_y' \cr 0&0&0&0&0&0\cr -\eta_x & -\eta_x' & -\eta_y & -\eta_y' & 0 & 1} \end{align}

Now, to get the local diffusion coefficient, we compute

(4)
$$d_a = Tr(G_a D_a)$$

We do not need to transform the invariant or the diffusion matrix into betatron coordinates in order to compute this.

In AT, we need to compute the diffusion matrix. First, define an element. e.g.
b1s=esrf{15}

b1s =

FamName: 'B1S'
Length: 0.2927
PolynomA: [0 0 0]
PolynomB: [0 0 0]
MaxOrder: 1
NumIntSteps: 10
BendingAngle: 0.0059
EntranceAngle: -0.0432
ExitAngle: 0.0491
K: 0
PassMethod: 'BndMPoleSymplectic4E2Pass'
Class: 'Bend'
BetaCode: 'DI'
FullGap: 0
FringeInt: 0
Energy: 6.0400e+09

Next, call findmpoleraddiffmatrix on this element with an orbit (can be zero offset for a dipole). Giving a non-zero offset in a quadrupole will result in some diffusion.

findmpoleraddiffmatrix(b1s,[0 0 0 0 0 0]')

ans =

1.0e-12 *

-0.0000 -0.0000 0 0 0.0000 -0.0000
-0.0000 0.0000 0 0 0.0008 -0.0000
0 0 0 0 0 0
0 0 0 0 0 0
0.0000 0.0008 0 0 0.7762 -0.0000
-0.0000 -0.0000 0 0 -0.0000 -0.0000

There's also this paper, we wrote recently. There is the apparent emittance, the projected emittance, and the eigenemittance.
The original paper describing the AT code is here. This paper describes the coupling formalism used.
The ohmienvelope command is very useful in AT. It computes the equilibrium beam sizes for the lattice due to radiation damping and diffusion. It is based on this paper by Ohmi et. al. In order to use it properly, one needs a lattice that contains pass methods that
include radiation, and also an RF cavity must be included. Here is the code:

function [ENVELOPE, RMSDP, RMSBL] = ohmienvelope(RING,RADELEMINDEX,varargin)
%OHMIENVELOPE calculates equilibrium beam envelope in a
% circular accelerator using Ohmi's beam envelope formalism [1]
% [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994)
%
% [ENVELOPE, RMSDP, RMSBL] = ONMIENVELOPE(RING,RADELEMINDEX,REFPTS)
%
% ENVELOPE is a structure with fields
% Sigma   - [SIGMA(1); SIGMA(2)] - RMS size [m] along
%           the principal axis of a tilted ellips
%           Assuming normal distribution exp(-(Z^2)/(2*SIGMA))
% Tilt    - Tilt angle of the XY ellips [rad]
%           Positive Tilt corresponds to Corkscrew (right)
%           rotatiom of XY plane around s-axis
% R       - 6-by-6 equilibrium envelope matrix R
%
% RMSDP   - RMS momentum spread
% RMSBL   - RMS bunch length[m]


One can use ohmienvelope by invoking the function atx:

function [lindt,pm]=atx(ring,varargin)
%ATX                computes and displays global information
%
%LINDATA=ATX(RING,DPP)
%
%RING:        AT structure
%DPP:        relative energy deviation (default: 0)
%
%LINDATA is a MATLAB structure array with fields
%
% From atlinopt:
%
%   ElemIndex   - ordinal position in the RING
%   SPos        - longitudinal position [m]
%   ClosedOrbit - closed orbit column vector with
%                 components x, px, y, py (momentums, NOT angles)
%   Dispersion  - dispersion orbit position vector with
%                 components eta_x, eta_prime_x, eta_y, eta_prime_y
%                 calculated with respect to the closed orbit with
%                 momentum deviation DP
%   M44         - 4x4 transfer matrix M from the beginning of RING
%                 to the entrance of the element for specified DP [2]
%   A           - 2x2 matrix A in [3]
%   B           - 2x2 matrix B in [3]
%   C           - 2x2 matrix C in [3]
%   gamma       - gamma parameter of the transformation to eigenmodes
%   mu          - [ mux, muy] horizontal and vertical betatron phase
%   beta        - [betax, betay] vector
%   alpha       - [alphax, alphay] vector
%
% From ohmienvelope:
%
%   beam66      - 6x6 equilibrium beam matrix
%   emit66    - 6x6 emittance projections on x and y + energy spread
%   beam44      - intersection of beam66 for dpp=0
%   emit44    - emittances of the projections of beam44 on x and y
%   modemit    - [emitA emitB] emittance of modes A and B (should be constant)
%
% See also: ATREADBETA ATLINOPT ATMODUL


atx uses the atradon function to turn the radiation on in the ring. This is given by

function [ring2,radindex,cavindex]=atradon(ring,cavipass,bendpass,quadpass)
%
%
%RING:        initial AT structure
%CAVIPASS:    pass method for cavities (default ThinCavityPass)
%BENDPASS:    pass method for cavities (default BndMPoleSymplectic4RadPass)
%QUADPASS:    pass method for cavities (default: nochange)
%
%                     and cavities

global GLOBVAL
if nargin < 4, quadpass=[]; end
if nargin < 3, bendpass=[]; end
if nargin < 2
cavipass='CavityPass';
end

if ~isfield(GLOBVAL,'E0')
error('Energy not defined in GLOBVAL');
end

if ~isempty(cavipass)
cavindex=findcells(ring,'Frequency');
if isempty(cavindex)
warning('No cavity found in the structure')
else
disp(['Cavities located at position ' num2str(cavindex)]);
end
ring2=setcellstruct(ring,'PassMethod',cavindex,cavipass);
else
ring2=ring;
end

if ~isempty(bendpass)
dipindex=findcells(ring2,'BetaCode','DI');
if isempty(dipindex), warning('No dipole found in the structure'); end
ring2=setcellstruct(ring2,'PassMethod',dipindex,bendpass);
else
dipindex=[];
end

if isempty(dipindex), warning('No quadrupole found in the structure'); end
else
end

disp([num2str(length(radindex)) ' elements switched to include radiation']);


We call this using, e.g.

[esrfRad,radindex,cavindex]=atradon(esrf);


The cavity also needs to be turned on:

esrf2=setcellstruct(esrf2,'HarmNumber',cavindex,getcellstruct(esrf2,'HarmNumber',cavindex)/periods);


In OhmiEnvelope, the code needs to compute the diffusion matrix. This uses

findmpoleraddiffmatrix(esrf{15},orbit(:,15))

void FindElemB(double *orbit_in, double le, double irho, double *A, double *B,
double *pt1, double* pt2,double *PR1, double *PR2,
double entrance_angle,     double exit_angle,
double fringe_int1, double fringe_int2, double full_gap,
int max_order, int num_int_steps,
double E0, double *BDIFF)


First we need to create a lattice structure with some coupling.
esrf=build_simple(0,'s13s20thick.str');
creates an uncoupled ideal machine from a beta file. Next, we tilt a quadrupole:
esrftilt=atsettilt(esrf,6,.01);

lindata=atx(esrftilt(:),0,1:length(esrftilt(:)));
tilt=cat(1,lindata.tilt)*360/(2*pi); //Convert to degrees.
spos=findspos(esrf,1:length(esrftilt));
tiltdat=[spos tilt];
save('tiltdat.txt','tiltdat','-ASCII')
plot(spos,tilt)

Next, we extract the beam moment matrix.
sig=cat(3,lindata.beam66)
Now sig is a 6x6x1648 variable. It is the collection of all the sigma matrices at each of the lattice elements. We write it to a file as a Matlab variable
save('beamsizes.mat','sig')

We could then read it into Mathematica via
sig = Import["~/matlab/beamsizes.mat"]
In Mathematica, we could get the sigma matrix at the first element with
sig[[1, 1, 1 ;; 6, 1 ;; 6]]

projection of coupled beam onto x-y axis at center of straight section. IBS has not yet been included.