diff --git a/atmat/atphysics/ParameterSummaryFunctions/atsummary.m b/atmat/atphysics/ParameterSummaryFunctions/atsummary.m index 004fc222c..1784d31ad 100644 --- a/atmat/atphysics/ParameterSummaryFunctions/atsummary.m +++ b/atmat/atphysics/ParameterSummaryFunctions/atsummary.m @@ -73,19 +73,18 @@ dp=TD(1).ClosedOrbit(5); % For calculating the synchrotron integrals - [I1d,I2d,I3d,I4d,I5d,I6,~] = DipoleRadiation(ring,TD); + [I1d,I2d,I3d,I4d,I5d,I6,~] = ElementRadiation(ring,TD); [I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,TD); smm.integrals=[I1d+I1w,I2d+I2w,I3d+I3w,I4d+I4w,I5d+I5w,I6]; - + if is6d alphac=mcf(atdisable_6d(ring),dp); eloss=atgetU0(ring,'method','tracking'); % eV - offmom = (abs(rev_frequency*harm_number - freq) > 2); else alphac=mcf(ring,dp); eloss=1.0e9*Cgamma/2/pi*smm.e0.^4*smm.integrals(2); % eV - offmom = (abs(dp) > 1.e-4); end + etac=1/gamma/gamma- alphac; smm.compactionFactor = alphac; smm.etac = etac; @@ -100,19 +99,9 @@ smm.damping(3) = 2 + smm.integrals(4)/smm.integrals(2); smm.radiation = 1.0e-9*eloss; % GeV - if offmom - warning('AT:NoEmit','\n%s\n%s\n%s',... - 'For off-momentum lattices, the equilibrium emittance',... - 'cannot be computed because the contribution of quadrupoles is missing.',... - 'To avoid this warning, use ">> warning(''off'',''AT:NoEmit'')"'); - smm.naturalEnergySpread = NaN; - smm.naturalEmittance = NaN; - alphac2=NaN; - else - smm.naturalEnergySpread = gamma*sqrt(Cq*smm.integrals(3)/(2*smm.integrals(2) + smm.integrals(4))); - smm.naturalEmittance = Cq*gamma.^2*smm.integrals(5)/(smm.integrals(2)-smm.integrals(4)); - alphac2 = smm.integrals(1)/smm.circumference; - end + smm.naturalEnergySpread = gamma*sqrt(Cq*smm.integrals(3)/(2*smm.integrals(2) + smm.integrals(4))); + smm.naturalEmittance = Cq*gamma.^2*smm.integrals(5)/(smm.integrals(2)-smm.integrals(4)); + alphac2 = smm.integrals(1)/smm.circumference; % Damping times cf=smm.radiation/2/smm.revTime/smm.e0; diff --git a/atmat/atphysics/ParameterSummaryFunctions/ringpara.m b/atmat/atphysics/ParameterSummaryFunctions/ringpara.m index 0148e9177..c4904fffc 100644 --- a/atmat/atphysics/ParameterSummaryFunctions/ringpara.m +++ b/atmat/atphysics/ParameterSummaryFunctions/ringpara.m @@ -60,12 +60,7 @@ beta_c = betar*cspeed; [ringdata,lindata]=atlinopt6(ring,1:length(ring)+1,varargin{:},'get_chrom'); - if is6d - offmom = (abs(frev*harm - freq_rf) > 2); - else - offmom = (abs(lindata(1).ClosedOrbit(5)) > 1.e-4); - end - [I1d,I2d,I3d,I4d,I5d,~,Iv] = DipoleRadiation(ring,lindata); + [I1d,I2d,I3d,I4d,I5d,~,Iv] = ElementRadiation(ring,lindata); [I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,lindata); I1=I1d+I1w; I2=I2d+I2w; @@ -74,7 +69,7 @@ I5=I5d+I5w; if ~isempty(Ux) U0 = Ux*1e6; %convert MeV to eV - fprintf('dipole radiation loss: %4.5f keV\n', U0/1000.); + fprintf('Radiation loss: %4.5f keV\n', U0/1000.); elseif is6d U0 = atgetU0(ring,'method','tracking'); else @@ -83,20 +78,9 @@ Jx = 1-I4/I2; Jy = 1.00; Je = 2+I4/I2; - if offmom - warning('AT:NoEmit','\n%s\n%s\n%s',... - 'For off-momentum lattices, the equilibrium emittance',... - 'cannot be computed because the contribution of quadrupoles is missing.',... - 'To avoid this warning, use ">> warning(''off'',''AT:NoEmit'')"'); - emittx = NaN; - sigma_E = NaN; - alphac = NaN; - else - emittx = Cq*gamma^2*I5/(I2-I4); - sigma_E = gamma*sqrt(Cq*I3/(2*I2+I4)); - alphac = I1/circ; - end - + emittx = Cq*gamma^2*I5/(I2-I4); + sigma_E = gamma*sqrt(Cq*I3/(2*I2+I4)); + alphac = I1/circ; % minimum emittance due to radiation 1/gamma cone (Handbook, Chao, eq23, pag 211) emitty_lim = 13/55*Cq/Jy/I2*Iv; diff --git a/atmat/atphysics/Radiation/DipoleRadiation.m b/atmat/atphysics/Radiation/DipoleRadiation.m index c7f07c4f9..2c2accf5c 100644 --- a/atmat/atphysics/Radiation/DipoleRadiation.m +++ b/atmat/atphysics/Radiation/DipoleRadiation.m @@ -7,105 +7,5 @@ % R.H. Helm, M.J. Lee, P.L. Morton and M. Sands % SLAC-PUB-1193, March 1973 - -angle=atgetfieldvalues(ring,'BendingAngle'); -isdipole=isfinite(angle) & (angle~=0); -vini=lindata([isdipole;false])'; -vend=lindata([false;isdipole])'; -[di1,di2,di3,di4,di5,di6,div]=cellfun(@diprad,ring(isdipole),num2cell(vini),num2cell(vend)); -I1=sum(di1); -I2=sum(di2); -I3=sum(di3); -I4=sum(di4); -I5=sum(di5); -I6=sum(di6); -Iv=sum(div); - - function [di1,di2,di3,di4,di5,di6,div]=diprad(elem,dini,dend) - betax0 = dini.beta(1); - betaz0 = dini.beta(2); - alphax0 = dini.alpha(1); - eta0 = dini.Dispersion(1); - etap0 = dini.Dispersion(2); - - ll = elem.Length; - theta = elem.BendingAngle; - rho = ll / theta; - rho2 = rho * rho; - K = getfoc(elem); - kx2 = K + 1.0/rho2; - kz2 = -K; - eps1 = tan(elem.EntranceAngle) / rho; - eps2 = tan(elem.ExitAngle) / rho; - - eta3 = dend.Dispersion(1); - alphax1 = alphax0 - betax0*eps1; - gammax1 = (1.0 + alphax1 * alphax1) / betax0; - alphaz1 = dini.alpha(2) + betaz0*eps1; - alphaz2 = dend.alpha(2) - dend.beta(2)*eps2; - gammaz1 = (1.0 + alphaz1*alphaz1) / betaz0; - etap1 = etap0 + eta0*eps1; - etap2 = dend.Dispersion(2) - eta3*eps2; - - h0 = gammax1*eta0*eta0 + 2.0*alphax1*eta0*etap1 + betax0*etap1*etap1; - - if kx2 ~= 0.0 - if kx2 > 0.0 % Focusing - kl = ll * sqrt(kx2); - ss = sin(kl) / kl; - cc = cos(kl); - else % Defocusing - kl = ll * sqrt(-kx2); - ss = sinh(kl) / kl; - cc = cosh(kl); - end - eta_ave = (theta - (etap2 - etap1)) / kx2 / ll; - bb = 2.0 * (alphax1*eta0 + betax0*etap1) * rho; - aa = -2.0 * (alphax1*etap1 + gammax1*eta0) * rho; - h_ave = h0 + (aa * (1.0-ss) + bb * (1.0-cc) / ll ... - + gammax1 * (3.0-4.0*ss+ss*cc) / 2.0 / kx2 ... - - alphax1 * (1.0-cc)^2 / kx2 / ll ... - + betax0 * (1.0-ss*cc) / 2.0 ... - ) / kx2 / rho2; - else - eta_ave = 0.5 * (eta0 + eta3) - ll*ll / 12.0 / rho; - hp0 = 2.0 * (alphax1 * eta0 + betax0 * etap1) / rho; - h2p0 = 2.0 * (-alphax1*etap1 + betax0/rho - gammax1*eta0) / rho; - h_ave = h0 + hp0*ll/2.0 + h2p0*ll*ll/6.0 ... - - alphax1*ll^3/4.0/rho2 ... - + gammax1*ll^4/20.0/rho2; - end - if kz2 ~= 0.0 - bz_ave=(gammaz1 + kz2*betaz0 + (alphaz2-alphaz1)/ll)/2/kz2; - else - bz_ave=betaz0-alphaz1*ll + gammaz1*ll*ll/3; - end - - di1 = eta_ave * ll / rho; - di2 = ll / rho2; - di3 = ll / abs(rho) / rho2; - di4 = eta_ave * ll * (2.0*K+1.0/rho2) / rho ... - - (eta0*eps1 + eta3*eps2)/rho; - di5 = h_ave * ll / abs(rho) / rho2; - di6 = kz2^2*eta_ave^2*ll; - div = abs(theta)/rho2*bz_ave; - end - - function K=getfoc(elem) - if isfield(elem,'PolynomB') && length(elem.PolynomB) >= 2 - K = elem.PolynomB(2); - if isfield(elem,'K') && elem.K ~= K - warning('AT:InconsistentK',... - 'Values in K and PolynomB(2) are different. Using PolynomB(2)'); - end - elseif isfield(elem,'K') - K = elem.K; - if K ~= 0 - warning('AT:InconsistentK','PolynomB(2) is missing. Using K'); - end - else - K = 0; - end - end -end +[I1,I2,I3,I4,I5,I6,Iv] = ElementRadiation(ring,lindata,'UseQuadrupole', false); diff --git a/atmat/atphysics/Radiation/ElementRadiation.m b/atmat/atphysics/Radiation/ElementRadiation.m new file mode 100644 index 000000000..5e2449a52 --- /dev/null +++ b/atmat/atphysics/Radiation/ElementRadiation.m @@ -0,0 +1,147 @@ +function [I1,I2,I3,I4,I5,I6,Iv] = ElementRadiation(ring,lindata,varargin) +%ELEMENTRADIATION Compute the radiation integrals in dipoles +%and quadrupoles + +% Analytical integration from: +% +% EVALUATION OF SYNCHROTRON RADIATION INTEGRALS +% R.H. Helm, M.J. Lee, P.L. Morton and M. Sands +% SLAC-PUB-1193, March 1973 + +[quadon, ~]=getoption(varargin,'UseQuadrupoles', true); +angle=atgetfieldvalues(ring,'BendingAngle'); +isdipole=isfinite(angle) & (angle~=0); + +if quadon + class=atgetcells(ring,'Class','Quadrupole'); + pb = atgetfieldvalues(ring,'PolynomB',{2}); + isquadrupole=class & (pb~=0); + iselement = isdipole | isquadrupole; +else + iselement = isdipole; +end +vini=lindata([iselement;false])'; +vend=lindata([false;iselement])'; + +[di1,di2,di3,di4,di5,di6,div]=cellfun(@elrad,ring(iselement),num2cell(vini),num2cell(vend)); +I1=sum(di1); +I2=sum(di2); +I3=sum(di3); +I4=sum(di4); +I5=sum(di5); +I6=sum(di6); +Iv=sum(div); + + function [di1,di2,di3,di4,di5,di6,div]=elrad(elem,dini,dend) + betax0 = dini.beta(1); + betaz0 = dini.beta(2); + alphax0 = dini.alpha(1); + eta0 = dini.Dispersion(1); + etap0 = dini.Dispersion(2); + try + theta = elem.BendingAngle; + catch + xpi = dini.ClosedOrbit(2)/(1+dini.ClosedOrbit(5)); + xpo = dend.ClosedOrbit(2)/(1+dini.ClosedOrbit(5)); + theta = xpi-xpo; + end + if abs(theta)<1.0e-7 + di1 = 0.0; + di2 = 0.0; + di3 = 0.0; + di4 = 0.0; + di5 = 0.0; + di6 = 0.0; + div = 0.0; + return + end + try + ange = elem.EntranceAngle; + catch + ange = 0.0; + end + try + ango = elem.ExitAngle; + catch + ango = 0.0; + end + ll = elem.Length; + rho = ll / theta; + rho2 = rho * rho; + K = getfoc(elem); + kx2 = K + 1.0/rho2; + kz2 = -K; + eps1 = tan(ange) / rho; + eps2 = tan(ango) / rho; + + eta3 = dend.Dispersion(1); + alphax1 = alphax0 - betax0*eps1; + gammax1 = (1.0 + alphax1 * alphax1) / betax0; + alphaz1 = dini.alpha(2) + betaz0*eps1; + alphaz2 = dend.alpha(2) - dend.beta(2)*eps2; + gammaz1 = (1.0 + alphaz1*alphaz1) / betaz0; + etap1 = etap0 + eta0*eps1; + etap2 = dend.Dispersion(2) - eta3*eps2; + + h0 = gammax1*eta0*eta0 + 2.0*alphax1*eta0*etap1 + betax0*etap1*etap1; + + if kx2 ~= 0.0 + if kx2 > 0.0 % Focusing + kl = ll * sqrt(kx2); + ss = sin(kl) / kl; + cc = cos(kl); + else % Defocusing + kl = ll * sqrt(-kx2); + ss = sinh(kl) / kl; + cc = cosh(kl); + end + eta_ave = (theta - (etap2 - etap1)) / kx2 / ll; + bb = 2.0 * (alphax1*eta0 + betax0*etap1) * rho; + aa = -2.0 * (alphax1*etap1 + gammax1*eta0) * rho; + h_ave = h0 + (aa * (1.0-ss) + bb * (1.0-cc) / ll ... + + gammax1 * (3.0-4.0*ss+ss*cc) / 2.0 / kx2 ... + - alphax1 * (1.0-cc)^2 / kx2 / ll ... + + betax0 * (1.0-ss*cc) / 2.0 ... + ) / kx2 / rho2; + else + eta_ave = 0.5 * (eta0 + eta3) - ll*ll / 12.0 / rho; + hp0 = 2.0 * (alphax1 * eta0 + betax0 * etap1) / rho; + h2p0 = 2.0 * (-alphax1*etap1 + betax0/rho - gammax1*eta0) / rho; + h_ave = h0 + hp0*ll/2.0 + h2p0*ll*ll/6.0 ... + - alphax1*ll^3/4.0/rho2 ... + + gammax1*ll^4/20.0/rho2; + end + if kz2 ~= 0.0 + bz_ave=(gammaz1 + kz2*betaz0 + (alphaz2-alphaz1)/ll)/2/kz2; + else + bz_ave=betaz0-alphaz1*ll + gammaz1*ll*ll/3; + end + + di1 = eta_ave * ll / rho; + di2 = ll / rho2; + di3 = ll / abs(rho) / rho2; + di4 = eta_ave * ll * (2.0*K+1.0/rho2) / rho ... + - (eta0*eps1 + eta3*eps2)/rho; + di5 = h_ave * ll / abs(rho) / rho2; + di6 = kz2^2*eta_ave^2*ll; + div = abs(theta)/rho2*bz_ave; + end + + function K=getfoc(elem) + if isfield(elem,'PolynomB') && length(elem.PolynomB) >= 2 + K = elem.PolynomB(2); + if isfield(elem,'K') && elem.K ~= K + warning('AT:InconsistentK',... + 'Values in K and PolynomB(2) are different. Using PolynomB(2)'); + end + elseif isfield(elem,'K') + K = elem.K; + if K ~= 0 + warning('AT:InconsistentK','PolynomB(2) is missing. Using K'); + end + else + K = 0; + end + end +end + diff --git a/pyat/at/lattice/lattice_object.py b/pyat/at/lattice/lattice_object.py index 9940790d3..7fad31d3b 100644 --- a/pyat/at/lattice/lattice_object.py +++ b/pyat/at/lattice/lattice_object.py @@ -15,7 +15,7 @@ import copy import numpy import math -from typing import Optional, Union, Tuple +from typing import Optional, Union if sys.version_info.minor < 9: from typing import Callable, Iterable, Generator SupportsIndex = int @@ -529,6 +529,31 @@ def copy(self) -> "Lattice": def deepcopy(self) -> "Lattice": """Returns a deep copy of the lattice""" return copy.deepcopy(self) + + def slice_elements(self, refpts, slices: Optional[int] = 1) \ + -> "Lattice": + """Create a new lattice by slicing the elements at refpts + + Keyword arguments: + size=None: Length of a slice. Default: computed from the + range and number of points: ``size = (s_max-s_min)/slices``. + slices=1: Number of slices in the specified range. Ignored if + size is specified. Default: no slicing + + Returns: + newring: New Lattice object + """ + def slice_generator(_): + check = get_bool_index(self, refpts) + for el, ok in zip(self, check): + if ok and (slices>1): + frac = numpy.ones(slices) / slices + for elem in el.divide(frac): + yield elem + else: + yield el + + return Lattice(slice_generator, iterator=self.attrs_filter) def slice(self, size: Optional[float] = None, slices: Optional[int] = 1) \ -> "Lattice": diff --git a/pyat/at/physics/radiation.py b/pyat/at/physics/radiation.py index 24d48b230..a2ed451fb 100644 --- a/pyat/at/physics/radiation.py +++ b/pyat/at/physics/radiation.py @@ -1,13 +1,15 @@ """ Radiation and equilibrium emittances """ +from __future__ import annotations from math import sin, cos, tan, sqrt, sinh, cosh, pi import numpy -from typing import Tuple +from typing import Union from scipy.linalg import inv, det, solve_sylvester from at.lattice import Lattice, check_radiation, Refpts -from at.lattice import Dipole, Wiggler, DConstant, Multipole, QuantumDiffusion -from at.lattice import get_refpts, get_value_refpts +from at.lattice import Dipole, Wiggler, DConstant +from at.lattice import Quadrupole, Multipole, QuantumDiffusion +from at.lattice import get_refpts, get_value_refpts, frequency_control from at.lattice import uint32_refpts, set_value_refpts from at.tracking import lattice_pass from at.physics import find_orbit6, find_m66, find_elem_m66, Orbit @@ -197,19 +199,20 @@ def propag(m, cumb, orbit6): return data0, r66data, data +@frequency_control def get_radiation_integrals(ring, dp: float = None, twiss=None, **kwargs)\ - -> Tuple[float, float, float, float, float]: + -> tuple[float, float, float, float, float]: r"""Computes the 5 radiation integrals for uncoupled lattices. Parameters: ring: Lattice description - dp: Momentum deviation. Ignored if radiation is ON twiss: Linear optics at all points (from :py:func:`.linopt6`). If ``None``, it will be computed. Keyword Args: - dct (float): Path lengthening. If specified, ``dp`` is - ignored and the off-momentum is deduced from the path lengthening. + dp (float): Momentum deviation. Defaults to :py:obj:`None` + dct (float): Path lengthening. Defaults to :py:obj:`None` + df (float): Deviation of RF frequency. Defaults to method (Callable): Method for linear optics: :py:obj:`~.linear.linopt2`: no longitudinal motion, no H/V coupling, @@ -225,21 +228,28 @@ def get_radiation_integrals(ring, dp: float = None, twiss=None, **kwargs)\ i5 (float): :math:`I_5 \quad [m^{-1}]` """ - def dipole_radiation(elem: Dipole, vini, vend): + def element_radiation(elem: Union[Dipole, Quadrupole], vini, vend): """Analytically compute the radiation integrals in dipoles""" beta0 = vini.beta[0] alpha0 = vini.alpha[0] eta0 = vini.dispersion[0] etap0 = vini.dispersion[1] + theta = getattr(elem, 'BendingAngle', None) + if theta is None: + xpi = vini.closed_orbit[1]/(1+vini.closed_orbit[4]) + xpo = vend.closed_orbit[1]/(1+vend.closed_orbit[4]) + theta = xpi-xpo + if abs(theta)<1.0e-7: + return numpy.zeros(5) + angin = getattr(elem, 'EntranceAngle', 0.0) + angout = getattr(elem, 'ExitAngle', 0.0) ll = elem.Length - theta = elem.BendingAngle rho = ll / theta rho2 = rho * rho k2 = elem.K + 1.0 / rho2 - eps1 = tan(elem.EntranceAngle) / rho - eps2 = tan(elem.ExitAngle) / rho - + eps1 = tan(angin) / rho + eps2 = tan(angout) / rho eta3 = vend.dispersion[0] alpha1 = alpha0 - beta0 * eps1 gamma1 = (1.0 + alpha1 * alpha1) / beta0 @@ -349,8 +359,8 @@ def harm(coef, h, phi): raise ValueError('length of Twiss data should be {0}' .format(len(ring) + 1)) for (el, vini, vend) in zip(ring, twiss[:-1], twiss[1:]): - if isinstance(el, Dipole) and el.BendingAngle != 0.0: - integrals += dipole_radiation(el, vini, vend) + if isinstance(el, (Dipole, Quadrupole)): + integrals += element_radiation(el, vini, vend) elif isinstance(el, Wiggler) and el.PassMethod != 'DriftPass': integrals += wiggler_radiation(el, vini) return tuple(integrals) diff --git a/pyat/at/physics/ring_parameters.py b/pyat/at/physics/ring_parameters.py index 33c7e1062..85b31734e 100644 --- a/pyat/at/physics/ring_parameters.py +++ b/pyat/at/physics/ring_parameters.py @@ -132,7 +132,6 @@ def radiation_parameters(ring: Lattice, dp: Optional[float] = None, rp.E0 = E0 rp.U0 = U0 emitx = Cq * gamma2 * rp.i5 / Jx / rp.i2 - rp.emittances = numpy.array([emitx, nan, nan]) alphac = rp.i1 / circumference etac = 1.0/gamma2 - alphac rp.phi_s = (pi - asin(U0 / voltage)) if U0 <= voltage else nan @@ -144,6 +143,7 @@ def radiation_parameters(ring: Lattice, dp: Optional[float] = None, rp.J = damping_partition_numbers rp.sigma_e = sqrt(Cq * gamma2 * rp.i3 / Je / rp.i2) rp.sigma_l = beta * abs(etac) * circumference / 2.0 / pi / nus * rp.sigma_e + rp.emittances = numpy.array([emitx, nan, rp.sigma_e*rp.sigma_l]) ringtunes, _ = numpy.modf(ring.periodicity * ringdata.tune) if len(ringtunes) < 3: rp.tunes = numpy.concatenate((ringtunes, (nus,)))